ENH: interpolate: add AAA algorithm for rational approximation #21167
Conversation
[docs only] Co-authored-by: Lucas Colley <lucas.colley8@gmail.com>
Yes I could add a sentence in the notes
Done |
Sorry, I meant adding some results to this PR conversation to justify why one may want to use |
Sure, whilst both are rational approximations, they are fundamentally quite different. Pade is concerned with rational approximation of a Taylor series at a point where as AAA approximates a function on a set of points. Furthermore, Pade approximation is a quotient of polynomials which breaks down when the poles are clustered whereas AAA use a barycentric representation which turns out to be incredibly stable |
|
CI is green so this is now ready for review |
There was a problem hiding this comment.
Initial cursory review: only looked at the docs and the overall structure.
Small minor comments are mostly optional. Additionally,
@cached_property: are we sure it plays ball with no-gil? Overall I'd suggest keeping properties for basically fetching things off the object and computations can as well be methods.- Consider making the input validation stricter.
array or callable, auto raveling thezarrays are best avoided IMO. - The amount of
with np.errstate("ignore")is a bit worrying. I did not look into the math though, maybe there's no way around
[skip ci] Co-authored-by: Lucas Colley <lucas.colley8@gmail.com>
[docs only]
|
@ev-br I wanted to check whether there was anything outstanding from your initial review or any additional comments? I have removed the cached properties and removed the I would be happy to look at any upcoming interpolate prs in return! |
|
Give me a couple of days @j-bowhay to take a second look? Unless you need it in earlier, in which case I'm happy to declare it good enough :-). |
That would be great, thank you! |
There was a problem hiding this comment.
Looks good modulo a few minor points of confusion over the docs and a possible FUD about licenses. Two more bigger questions, one about the math, and one about the design.
- the AAA paper talks about some way of dealing with numerical doublets--- IIUC, if run for long enough, the algorithm will produce poles with residues of the order of machine epsilon. Do we want/need to add some step to deal with those? Or, if not, add some documentation on how to deal with these on the user side?
- the PR adds a class which does two different things: constructs a rational function with a certain recipe, and adds a callable object which represents a rational function. Have you considered splitting the latter into a separate entity (a BarycentricInterpolator? No, this name is taken by something related but different and not very perf-friendly). I personally like a pattern of having a factory function to return an instance of a class which does evaluations, knows about poles, residues and zeros etc. I'm not going to insist on exactly this, just would encourage you to give it a thought, also depending on what the future plans are. My main concern is to not paint outselves into a corner API-wise.
None of this is blocking though, and can be left for a possible follow-up. So consider this PR as approved from my side.
| from `values`, and :math:`w_1,\dots,w_m` are real or complex weights. The algorithm | ||
| then proceeds to select the next support point :math:`z_{m+1}` from the remaining | ||
| unselected `points` such that the nonlinear residual | ||
| :math:`|f(z_{m+1}) - n(z_{m+1})/d(z_{m+1})|` is maximised. The algorithm terminates |
There was a problem hiding this comment.
residual ... is minimized? Ouch, no. It selects as the m+1) -th point the one, k, from previously unselected points, which corresponds to the maximal value of the residual, is this correct?
What confused me was that I was expecting an LSQ type process under which some residual is minimized. I guess this is what you're saying, actually. Consider rephrasing?
There was a problem hiding this comment.
I will try to make the two parts of the algorithm clearer: the next support point is chosen to maximise the residual but the weights are chosen to minimise the residual
scipy/interpolate/_aaa.py
Outdated
| unselected `points` such that the nonlinear residual | ||
| :math:`|f(z_{m+1}) - n(z_{m+1})/d(z_{m+1})|` is maximised. The algorithm terminates | ||
| when this maximum is less than ``rtol * np.linalg.norm(f, ord=np.inf)``. This means | ||
| the interpolation property is only satisfied up to a tolerance. The weights are |
There was a problem hiding this comment.
But it is satisfied exactly on the subset of the data points, which were selected? Feel free to discount my confusion if you think it's just me, or rephrase.
There was a problem hiding this comment.
This is correct, will add an extra little bit
| \text{minimise}|fd - n| \quad \text{subject to} \quad | ||
| \sum_{i=1}^{m+1} w_i = 1, | ||
|
|
||
| over the unselected elements of `points`. |
There was a problem hiding this comment.
Would be great to add that the minimization is over the weights w_j
Thanks for the review!
Yes good spot, I have code written which implements this. I had planned to leave this as a follow up to keep the diff smaller but could add it here if you prefer
Yes good point I have thought about this as we may want to add other weights in the future such as Floater-Hormann. In which case my plan was to have a base class |
|
I have hopefully cleared up some points of confusion in the docs @ev-br but if not happy to polish in follow up prs! |
|
Thanks @ev-br @mdhaber @lucascolley @stefanv @rkern @izaid @nakatsukasayuji for your help! |
|
congrats Jake, great work! Happy to review adding array API standard support, even if partially blocked on And congrats on the inclusion in SciPy @nakatsukasayuji ! |
Reference issue
https://discuss.scientific-python.org/t/inclusion-of-aaa-rational-approximation-in-scipy/1282
What does this implement/fix?
This PR implements the core AAA algorithm for rational approximation. For more background please see the mailing list post linked above. There are a number of extra features that would be good to include however I would prefer to leave these as follow ups to keep the diff manageable.
Features:
Planned Features:
Additional information
cc @izaid @nakatsukasayuji