scipy.interpolate.

AAA#

class scipy.interpolate.AAA(x, y, *, rtol=None, max_terms=100, clean_up=True, clean_up_tol=1e-13)[source]#

AAA real or complex rational approximation.

As described in [1], the AAA algorithm is a greedy algorithm for approximation by rational functions on a real or complex set of points. The rational approximation is represented in a barycentric form from which the roots (zeros), poles, and residues can be computed.

Parameters:

x1D array_like, shape (n,): 1-D array containing values of the independent variable. Values may be real or complex but must be finite.
y1D array_like, shape (n,): Function values f(x). Infinite and NaN values of values and corresponding values of points will be discarded.
rtolfloat, optional: Relative tolerance, defaults to eps**0.75. If a small subset of the entries in values are much larger than the rest the default tolerance may be too loose. If the tolerance is too tight then the approximation may contain Froissart doublets or the algorithm may fail to converge entirely.
max_termsint, optional: Maximum number of terms in the barycentric representation, defaults to 100. Must be greater than or equal to one.
clean_upbool, optional: Automatic removal of Froissart doublets, defaults to True. See notes for more details.
clean_up_tolfloat, optional: Poles with residues less than this number times the geometric mean of values times the minimum distance to points are deemed spurious by the cleanup procedure, defaults to 1e-13. See notes for more details.

Attributes:

support_pointsarray: Support points of the approximation. These are a subset of the provided x at which the approximation strictly interpolates y. See notes for more details.
support_valuesarray: Value of the approximation at the support_points.
weightsarray: Weights of the barycentric approximation.
errorsarray: Error \(|f(z) - r(z)|_\infty\) over points in the successive iterations of AAA.

Methods

`__call__`(z)	Evaluate the rational approximation at given values.
`clean_up`([cleanup_tol])	Automatic removal of Froissart doublets.
`poles`()	Compute the poles of the rational approximation.
`residues`()	Compute the residues of the poles of the approximation.
`roots`()	Compute the zeros of the rational approximation.

Warns:

RuntimeWarning: If rtol is not achieved in max_terms iterations.

See also

FloaterHormannInterpolator: Floater-Hormann barycentric rational interpolation.
pade: Padé approximation.

Notes

At iteration \(m\) (at which point there are \(m\) terms in the both the numerator and denominator of the approximation), the rational approximation in the AAA algorithm takes the barycentric form

\[r(z) = n(z)/d(z) = \frac{\sum_{j=1}^m\ w_j f_j / (z - z_j)}{\sum_{j=1}^m w_j / (z - z_j)},\]

where \(z_1,\dots,z_m\) are real or complex support points selected from x, \(f_1,\dots,f_m\) are the corresponding real or complex data values from y, and \(w_1,\dots,w_m\) are real or complex weights.

Each iteration of the algorithm has two parts: the greedy selection the next support point and the computation of the weights. The first part of each iteration is to select the next support point to be added \(z_{m+1}\) from the remaining unselected x, such that the nonlinear residual \(|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, except at the support points where approximation exactly interpolates the supplied data.

In the second part of each iteration, the weights \(w_j\) are selected to solve the least-squares problem

\[\text{minimise}_{w_j}|fd - n| \quad \text{subject to} \quad \sum_{j=1}^{m+1} w_j = 1,\]

over the unselected elements of x.

One of the challenges with working with rational approximations is the presence of Froissart doublets, which are either poles with vanishingly small residues or pole-zero pairs that are close enough together to nearly cancel, see [2]. The greedy nature of the AAA algorithm means Froissart doublets are rare. However, if rtol is set too tight then the approximation will stagnate and many Froissart doublets will appear. Froissart doublets can usually be removed by removing support points and then resolving the least squares problem. The support point \(z_j\), which is the closest support point to the pole \(a\) with residue \(\alpha\), is removed if the following is satisfied

\[|\alpha| / |z_j - a| < \verb|clean_up_tol| \cdot \tilde{f},\]

where \(\tilde{f}\) is the geometric mean of support_values.

References

[1] (1,2,3)

Y. Nakatsukasa, O. Sete, and L. N. Trefethen, “The AAA algorithm for rational approximation”, SIAM J. Sci. Comp. 40 (2018), A1494-A1522. DOI:10.1137/16M1106122

[2]

J. Gilewicz and M. Pindor, Pade approximants and noise: rational functions, J. Comp. Appl. Math. 105 (1999), pp. 285-297. DOI:10.1016/S0377-0427(02)00674-X

Examples

Here we reproduce a number of the numerical examples from [1] as a demonstration of the functionality offered by this method.

>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from scipy.interpolate import AAA
>>> import warnings

For the first example we approximate the gamma function on [-3.5, 4.5] by extrapolating from 100 samples in [-1.5, 1.5].

>>> from scipy.special import gamma
>>> sample_points = np.linspace(-1.5, 1.5, num=100)
>>> r = AAA(sample_points, gamma(sample_points))
>>> z = np.linspace(-3.5, 4.5, num=1000)
>>> fig, ax = plt.subplots()
>>> ax.plot(z, gamma(z), label="Gamma")
>>> ax.plot(sample_points, gamma(sample_points), label="Sample points")
>>> ax.plot(z, r(z).real, '--', label="AAA approximation")
>>> ax.set(xlabel="z", ylabel="r(z)", ylim=[-8, 8], xlim=[-3.5, 4.5])
>>> ax.legend()
>>> plt.show()

../../_images/scipy-interpolate-AAA-1_00_00.png

We can also view the poles of the rational approximation and their residues:

>>> order = np.argsort(r.poles())
>>> r.poles()[order]
array([-3.81591039e+00+0.j        , -3.00269049e+00+0.j        ,
       -1.99999988e+00+0.j        , -1.00000000e+00+0.j        ,
        5.85842812e-17+0.j        ,  4.77485458e+00-3.06919376j,
        4.77485458e+00+3.06919376j,  5.29095868e+00-0.97373072j,
        5.29095868e+00+0.97373072j])
>>> r.residues()[order]
array([ 0.03658074 +0.j        , -0.16915426 -0.j        ,
        0.49999915 +0.j        , -1.         +0.j        ,
        1.         +0.j        , -0.81132013 -2.30193429j,
       -0.81132013 +2.30193429j,  0.87326839+10.70148546j,
        0.87326839-10.70148546j])

For the second example, we call AAA with a spiral of 1000 points that wind 7.5 times around the origin in the complex plane.

>>> z = np.exp(np.linspace(-0.5, 0.5 + 15j*np.pi, 1000))
>>> r = AAA(z, np.tan(np.pi*z/2), rtol=1e-13)

We see that AAA takes 12 steps to converge with the following errors:

>>> r.errors.size
12
>>> r.errors
array([2.49261500e+01, 4.28045609e+01, 1.71346935e+01, 8.65055336e-02,
       1.27106444e-02, 9.90889874e-04, 5.86910543e-05, 1.28735561e-06,
       3.57007424e-08, 6.37007837e-10, 1.67103357e-11, 1.17112299e-13])

We can also plot the computed poles:

>>> fig, ax = plt.subplots()
>>> ax.plot(z.real, z.imag, '.', markersize=2, label="Sample points")
>>> ax.plot(r.poles().real, r.poles().imag, '.', markersize=5,
...         label="Computed poles")
>>> ax.set(xlim=[-3.5, 3.5], ylim=[-3.5, 3.5], aspect="equal")
>>> ax.legend()
>>> plt.show()

../../_images/scipy-interpolate-AAA-1_01_00.png

We now demonstrate the removal of Froissart doublets using the clean_up method using an example from [1]. Here we approximate the function \(f(z)=\log(2 + z^4)/(1 + 16z^4)\) by sampling it at 1000 roots of unity. The algorithm is run with rtol=0 and clean_up=False to deliberately cause Froissart doublets to appear.

>>> z = np.exp(1j*2*np.pi*np.linspace(0,1, num=1000))
>>> def f(z):
...     return np.log(2 + z**4)/(1 - 16*z**4)
>>> with warnings.catch_warnings():  # filter convergence warning due to rtol=0
...     warnings.simplefilter('ignore', RuntimeWarning)
...     r = AAA(z, f(z), rtol=0, max_terms=50, clean_up=False)
>>> mask = np.abs(r.residues()) < 1e-13
>>> fig, axs = plt.subplots(ncols=2)
>>> axs[0].plot(r.poles().real[~mask], r.poles().imag[~mask], '.')
>>> axs[0].plot(r.poles().real[mask], r.poles().imag[mask], 'r.')

Now we call the clean_up method to remove Froissart doublets.

>>> with warnings.catch_warnings():
...     warnings.simplefilter('ignore', RuntimeWarning)
...     r.clean_up()
4
>>> mask = np.abs(r.residues()) < 1e-13
>>> axs[1].plot(r.poles().real[~mask], r.poles().imag[~mask], '.')
>>> axs[1].plot(r.poles().real[mask], r.poles().imag[mask], 'r.')
>>> plt.show()

../../_images/scipy-interpolate-AAA-1_02_00.png

The left image shows the poles prior of the approximation clean_up=False with poles with residue less than 10^-13 in absolute value shown in red. The right image then shows the poles after the clean_up method has been called.