CQF预备知识：Python相关库 -- 分段多项式和样条函数 scipy.interpolate-CSDN博客

文中内容仅限技术学习与代码实践参考，市场存在不确定性，技术分析需谨慎验证，不构成任何投资建议。

分段多项式和样条函数

在上一节讨论的一维插值例程中，通过构造某些分段多项式来工作：插值范围通过所谓的断点被分割成区间，每个区间上都有一个特定的多项式。然后这些多项式片段在断点处以预定义的平滑度匹配：对于三次样条是二阶导数，对于单调插值函数是一阶导数，等等。

一个 $k$ 次多项式可以被看作是 $k + 1$ 个单项式基元素的线性组合，即 $x^2, \cdots, x^k$ 。在某些应用中，考虑替代的（尽管形式上等价的）基是有用的。在 scipy.interpolate 中实现了两种流行的基，分别是 B 样条（BSpline）和贝塞尔多项式（BPoly）。B 样条通常用于例如非参数回归问题，而贝塞尔多项式用于构造贝塞尔曲线。

PPoly 对象表示在“通常的”幂基中的分段多项式。这是 CubicSpline 实例和单调插值函数的情况。一般来说，PPoly 对象可以表示任意阶的多项式，而不仅仅是三次多项式。对于数据数组 x，断点位于数据点处，而系数数组 c 定义了 $k$ 次多项式，使得 c[i, j] 是在 x[j] 和 x[j+1] 之间的线段上 (x - x[j])**(k-i) 的系数。

BSpline 对象表示 B 样条函数 —— B 样条基元素的线性组合。这些对象可以直接实例化，或者通过 make_interp_spline 工厂函数从数据中构造。

最后，贝塞尔多项式以 BPoly 类的实例形式表示。

所有这些类都实现了一个（大部分）类似的接口，PPoly 是功能最完备的。接下来我们将考虑这个接口的主要特性，并讨论分段多项式的替代基的一些细节。

操作 `PPoly` 对象

PPoly 对象有方便的方法用于构造导数和原函数、计算积分以及寻找根。例如，我们对正弦函数进行制表，并找到其导数的根。

>>> import numpy as np
>>> from scipy.interpolate import CubicSpline
>>> x = np.linspace(0, 10, 71)
>>> y = np.sin(x)
>>> spl = CubicSpline(x, y)

现在，对样条函数求导：

>>> dspl = spl.derivative()

这里 dspl 是一个 PPoly 实例，它表示原始对象 spl 的导数的多项式近似。在固定参数处评估 dspl 等同于使用 nu=1 参数评估原始样条函数：

>>> dspl(1.1), spl(1.1, nu=1)
(0.45361436, 0.45361436)

注意，上述第二种形式是在原地评估导数，而使用 dspl 对象，我们可以找到 spl 的导数的零点：

>>> dspl.roots() / np.pi
array([-0.45480801,  0.50000034,  1.50000099,  2.5000016 ,  3.46249993])

这与 $cos⁡(x)=sin⁡′(x)\cos(x) = \sin'(x)$ 的根 $n\pi/2 + \pi\,n$ 非常吻合。注意，默认情况下它计算了在插值区间 $\leqslant x \leqslant 10$ 之外外推的根，而且外推的结果（第一个和最后一个值）要不准确得多。我们可以关闭外推，并将根查找限制在插值区间内：

>>> dspl.roots(extrapolate=False) / np.pi
array([0.50000034,  1.50000099,  2.5000016])

实际上，root 方法是更一般的 solve 方法的一个特例，后者用于找到给定常数 $y$ 的方程 $f (x) = y$ 的解，其中 $f (x)$ 是分段多项式：

>>> dspl.solve(0.5, extrapolate=False) / np.pi
array([0.33332755, 1.66667195, 2.3333271])

这与预期值 $n\pm\arccos(1/2) + 2\pi\,n$ 非常吻合。

分段多项式的积分可以使用 .integrate 方法计算，该方法接受积分的上下限。例如，我们计算完全椭圆积分 $\int_0^{\pi/2} [1 - m\sin^2 x]^{-1/2} dx$ 的近似值：

>>> from scipy.special import ellipk
>>> m = 0.5
>>> ellipk(m)
1.8540746773013719

为此，我们对被积函数进行制表，并使用单调 PCHIP 插值函数（我们也可以使用 CubicSpline）进行插值：

>>> from scipy.interpolate import PchipInterpolator
>>> x = np.linspace(0, np.pi/2, 70)
>>> y = (1 - m*np.sin(x)**2)**(-1/2)
>>> spl = PchipInterpolator(x, y)

然后进行积分：

>>> spl.integrate(0, np.pi/2)
1.854074674965991

这确实接近由 scipy.special.ellipk 计算的值。

所有分段多项式都可以用 N 维的 y 值构造。如果 y.ndim > 1，它被理解为沿插值轴（默认值为 0）排列的 1D y 值的堆叠。后者是通过 axis 参数指定的，不变量是 len(x) == y.shape[axis]。例如，我们将上面的椭圆积分示例扩展到为一系列 m 值计算近似值，使用 NumPy 的广播：

>>> from scipy.interpolate import PchipInterpolator
>>> m = np.linspace(0, 0.9, 11)
>>> x = np.linspace(0, np.pi/2, 70)
>>> y = 1 / np.sqrt(1 - m[:, None]*np.sin(x)**2)

现在 y 数组的形状是 (11, 70)，因此对于固定的 m 值，y 的值沿 y 数组的第二轴排列。

>>> spl = PchipInterpolator(x, y, axis=1)  # 默认值是 axis=0
>>> import matplotlib.pyplot as plt
>>> plt.plot(m, spl.integrate(0, np.pi/2), '--')

>>> from scipy.special import ellipk
>>> plt.plot(m, ellipk(m), 'o')
>>> plt.legend(['`ellipk`', 'integrated piecewise polynomial'])
>>> plt.show()

B 样条：结点和系数

例如通过 make_interp_spline 调用从数据构造的 B 样条函数，是由所谓的结点和系数定义的。

以再次构造正弦函数的插值为例。B 样条实例的结点可通过其 t 属性获得：

>>> x = np.linspace(0, 3/2, 7)
>>> y = np.sin(np.pi*x)
>>> from scipy.interpolate import make_interp_spline
>>> bspl = make_interp_spline(x, y, k=3)
>>> print(bspl.t)
[0.  0.  0.  0.        0.5  0.75  1.        1.5  1.5  1.5  1.5 ]
>>> print(x)
[            0.  0.25  0.5  0.75  1.  1.25  1.5 ]

我们看到，默认情况下，结点向量是从输入数组 x 构造的：首先，它被构造为 $(k + 1)$ -正则（它在开头和结尾各附加了 $k$ 个重复的结点）；然后，输入数组的第二个和倒数第二个点被移除——这就是所谓的非结点边界条件。

一般来说，一个 $k$ 次插值样条需要 len(t) - len(x) - k - 1 个边界条件。对于具有 $(k + 1)$ -正则结点数组的三次样条，这意味着两个边界条件——或者从 x 数组中移除两个值。可以使用 make_interp_spline 的可选参数 bc_type 请求各种边界条件。

B 样条的系数可以通过 BSpline 对象的 c 属性访问：

约定是对于 len(t) 个结点，有 len(t) - k - 1 个系数。某些例程（见“平滑样条”部分）会对 c 数组进行零填充，使得 len(c) == len(t)。这些额外的系数在评估中会被忽略。

我们强调，系数是用 B 样条基给出的，而不是幂基 $\cdots, x^k$ 。

B 样条基元素

B 样条是分段多项式，表示为B 样条基元素的线性组合——而这些基元素本身是某些幂基元素 $x^m$ （ $\dots, k$ ）的线性组合。

B 样条基通常比幂基更计算稳定，并且适用于包括插值、回归和曲线表示在内的各种应用。主要特点是这些基元素是局部化的，并且在由结点数组定义的区间之外等于零。

具体来说，一个 $k$ 次的 B 样条基元素（例如三次样条的 $k = 3$ ）由 $k + 2$ 个结点定义，并且在这两个结点之外为零。为了说明这一点，我们在某个区间上绘制一系列非零基元素：

>>> k = 3      # 三次样条
>>> t = [0., 1.4, 2., 3.1, 5.]   # 内部结点
>>> t = np.r_[[0]*k, t, [5]*k]   # 添加边界结点

>>> from scipy.interpolate import BSpline
>>> import matplotlib.pyplot as plt
>>> for j in [-2, -1, 0, 1, 2]:
...     a, b = t[k+j], t[-k+j-1]
...     xx = np.linspace(a, b, 101)
...     bspl = BSpline.basis_element(t[k+j:-k+j])
...     plt.plot(xx, bspl(xx), label=f'j = {j}')
>>> plt.legend(loc='best')
>>> plt.show()

这里 BSpline.basis_element 实际上是构造一个仅有一个非零系数的样条的简写。例如，上述示例中的 j=2 元素等同于

>>> c = np.zeros(t.size - k - 1)
>>> c[-2] = 1
>>> b = BSpline(t, c, k)
>>> np.allclose(b(xx), bspl(xx))
True

如果需要，可以使用 PPoly.from_spline 方法将 B 样条转换为 PPoly 对象，该方法接受一个 BSpline 实例并返回一个 PPoly 实例。反向转换由 BSpline.from_power_basis 方法执行。然而，最好避免在基之间进行转换，因为这会累积舍入误差。

B 样条基中的设计矩阵

B 样条的一个常见应用是非参数回归。原因是 B 样条基元素的局部化特性使得线性代数变得稀疏。这是因为在一个给定的评估点上，最多有 $k + 1$ 个基元素是非零的，因此基于 B 样条构建的设计矩阵最多有 $k + 1$ 条对角线。

为了说明这一点，我们考虑一个玩具示例。假设我们的数据是一维的，并且限制在区间 $[0, 6]$ 内。我们构造一个 4-正则的结点向量，对应于 7 个数据点和三次样条（ $k = 3$ ）：

>>> t = [0., 0., 0., 0., 2., 3., 4., 6., 6., 6., 6.]

然后，取“观测值”为

并构造以稀疏 CSR 格式的设计矩阵

>>> from scipy.interpolate import BSpline
>>> mat = BSpline.design_matrix(xnew, t, k=3)
>>> mat
<Compressed Sparse Row sparse array of dtype 'float64'
        with 12 stored elements and shape (3, 7)>

这里设计矩阵的每一行对应于 xnew 数组中的一个值，每一行最多有 k+1 = 4 个非零元素；第 j 行包含在 xnew[j] 处评估的基元素：

>>> with np.printoptions(precision=3):
...     print(mat.toarray())
[[0.125 0.514 0.319 0.042 0.    0.    0.   ]
 [0.    0.111 0.556 0.333 0.    0.    0.   ]
 [0.    0.    0.125 0.75  0.125 0.    0.   ]]