scipy.interpolate.

AAA#

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

AAA 实数或复数有理逼近。

如 [1] 中所述，AAA 算法是一种贪婪算法，用于在实数或复数点集上通过有理函数进行逼近。有理逼近以重心形式表示，可以从中计算根（零点）、极点和残差。

参数:

x一维类数组，形状 (n,): 包含自变量值的 1 维数组。值可以是实数或复数，但必须是有限的。
y一维类数组，形状 (n,): 函数值 f(x)。将丢弃 values 的无穷大和 NaN 值以及 points 的对应值。
rtol浮点数，可选: 相对容差，默认为 eps**0.75。如果 values 中的一小部分条目远大于其余部分，则默认容差可能太宽松。如果容差太紧，则逼近可能包含 Froissart 双峰，或者算法可能完全无法收敛。
max_terms整数，可选: 重心表示中的最大项数，默认为 100。必须大于或等于 1。
clean_up布尔值，可选: 自动删除 Froissart 双峰，默认为 True。有关更多详细信息，请参见注释。
clean_up_tol浮点数，可选: 残差小于此数乘以 values 的几何平均值乘以到 points 的最小距离的极点被清理程序视为伪极点，默认为 1e-13。有关更多详细信息，请参见注释。

警告:

RuntimeWarning: 如果在 max_terms 次迭代中未达到 rtol。

另请参见

FloaterHormannInterpolator: Floater-Hormann 重心有理插值。
pade: Padé 逼近。

注释

在迭代 \(m\)（此时逼近的分子和分母中都有 \(m\) 项）时，AAA 算法中的有理逼近采用重心形式

\[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)},\]

其中 \(z_1,\dots,z_m\) 是从 x 中选择的实数或复数支持点，\(f_1,\dots,f_m\) 是来自 y 的对应的实数或复数数据值，\(w_1,\dots,w_m\) 是实数或复数权重。

算法的每次迭代都有两个部分：贪婪选择下一个支持点和计算权重。每次迭代的第一部分是从剩余的未选择的 x 中选择要添加的下一个支持点 \(z_{m+1}\)，使得非线性残差 \(|f(z_{m+1}) - n(z_{m+1})/d(z_{m+1})|\) 最大化。当此最大值小于 rtol * np.linalg.norm(f, ord=np.inf) 时，算法终止。这意味着插值属性仅在容差范围内满足，除了在逼近完全插值所提供数据的支持点之外。

在每次迭代的第二部分中，选择权重 \(w_j\) 来解决最小二乘问题

\[\text{最小化}_{w_j}|fd - n| \quad \text{受限于} \quad \sum_{j=1}^{m+1} w_j = 1,\]

在 x 的未选择元素上。

使用有理逼近的一个挑战是存在 Froissart 双峰，这些双峰要么是残差极小的极点，要么是靠得很近足以几乎相互抵消的极点-零点对，请参阅 [2]。AAA 算法的贪婪性质意味着 Froissart 双峰很少见。但是，如果 rtol 设置得太紧，则逼近将停滞，并且会出现许多 Froissart 双峰。Froissart 双峰通常可以通过删除支持点然后解决最小二乘问题来删除。如果满足以下条件，则删除最接近残差为 \(\alpha\) 的极点 \(a\) 的支持点 \(z_j\)

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

其中 \(\tilde{f}\) 是 support_values 的几何平均值。

参考文献

[1] (1,2,3)

Y. Nakatsukasa、O. Sete 和 L. N. Trefethen，“用于有理逼近的 AAA 算法”，SIAM J. Sci. Comp. 40 (2018)，A1494-A1522。DOI:10.1137/16M1106122

[2]

J. Gilewicz 和 M. Pindor，Pade 逼近和噪声：有理函数，J. Comp. Appl. Math. 105 (1999)，第 285-297 页。DOI:10.1016/S0377-0427(02)00674-X

示例

在这里，我们重现了 [1] 中的许多数值示例，以演示此方法提供的功能。

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

对于第一个示例，我们通过从 [-1.5, 1.5] 中的 100 个样本进行外推来逼近 [-3.5, 4.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

我们还可以查看有理逼近的极点及其残差

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

对于第二个示例，我们使用复平面中围绕原点缠绕 7.5 次的 1000 个点的螺旋来调用 AAA。

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

我们看到 AAA 需要 12 个步骤才能收敛，并且具有以下误差

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

我们还可以绘制计算出的极点

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

现在，我们使用 [1] 中的示例演示如何使用 clean_up 方法删除 Froissart 双峰。这里，我们通过在 1000 个单位根处对其进行采样来逼近函数 \(f(z)=\log(2 + z^4)/(1 + 16z^4)\)。该算法以 rtol=0 和 clean_up=False 运行，以故意导致出现 Froissart 双峰。

>>> 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.')

现在，我们调用 clean_up 方法删除 Froissart 双峰。

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

左图显示了逼近 clean_up=False 之前的极点，其中以红色显示绝对值小于 10^-13 的残差极点。右图显示了调用 clean_up 方法后的极点。

属性:

support_points数组: 逼近的支持点。这些是提供的 x 的子集，在该子集中，逼近严格插值 y。有关更多详细信息，请参见注释。
support_values数组: 在 support_points 处的逼近值。
weights数组: 重心逼近的权重。
errors数组: 在 AAA 迭代过程中，\(|f(z) - r(z)|_\infty\) 在点上的误差。

方法

`__call__`(z)	在给定值处评估有理逼近。
`clean_up`([cleanup_tol])	自动移除 Froissart 重复极点。
`poles`()	计算有理逼近的极点。
`residues`()	计算逼近极点的残差。
`roots`()	计算有理逼近的零点。