Skip to main content
`Ctrl`+`K`
scipy.interpolate.

# AAA

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

  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:
    x 1D array_like, shape (n,)

      1-D array containing values of the independent variable. Values may be real or complex but must be finite.

    y 1D array_like, shape (n,)

      Function values `f(x)`. Infinite and NaN values of values and corresponding values of points will be discarded.

    rtol float, 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_terms int, optional

      Maximum number of terms in the barycentric representation, defaults to `100`. Must be greater than or equal to one.

    clean_up bool, optional

      Automatic removal of Froissart doublets, defaults to `True`. See notes for more details.

    clean_up_tol float, 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.

  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()

  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()

  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()

  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.

  Attributes:
    support_points array

      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_values array

      Value of the approximation at the `support_points`.

    weights array

      Weights of the barycentric approximation.

    errors array

      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.