Classe « AAA »

Informations générales

Héritage

builtins.object
    _BarycentricRational
        AAA

Définition

class AAA(_BarycentricRational):

help(AAA)

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.

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 :math:`|f(z) - r(z)|_\infty` over `points` in the successive iterations
    of AAA.

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 :math:`m` (at which point there are :math:`m` terms in the both the
numerator and denominator of the approximation), the
rational approximation in the AAA algorithm takes the barycentric form

.. math::

    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 :math:`z_1,\dots,z_m` are real or complex support points selected from
`x`, :math:`f_1,\dots,f_m` are the corresponding real or complex data values
from `y`, and :math:`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 :math:`z_{m+1}` from the remaining
unselected `x`, 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, except at the
support points where approximation exactly interpolates the supplied data.

In the second part of each iteration, the weights :math:`w_j` are selected to solve
the least-squares problem

.. math::

    \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 :math:`z_j`,
which is the closest support point to the pole :math:`a` with residue
:math:`\alpha`, is removed if the following is satisfied

.. math::

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

where :math:`\tilde{f}` is the geometric mean of `support_values`.


References
----------
.. [1] 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
:math:`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.