分段多项式和样条#
前面章节讨论的 1D 插值例程通过构建特定的分段多项式来工作：插值范围通过所谓的断点分成间隔，并且每个间隔上都有一个特定的多项式。然后，这些多项式片段在断点处以预定义的平滑度匹配：三次样条的二阶导数，单调插值器的一阶导数，依此类推。
次数为 \(k\) 的多项式可以被认为是 \(k+1\) 个单项式基元素 \(1, x, x^2, \cdots, x^k\) 的线性组合。在某些应用中，考虑替代（如果形式上等效）的基会很有用。在 scipy.interpolate 中实现的两种流行的基是 B 样条 (BSpline) 和伯恩斯坦多项式 (BPoly)。B 样条通常用于例如非参数回归问题，而伯恩斯坦多项式用于构造贝塞尔曲线。
PPoly 对象以“通常”的幂基表示分段多项式。 CubicSpline 实例和单调插值器都是这种情况。通常，PPoly 对象可以表示任意阶的多项式，而不仅仅是三次多项式。对于数据数组 x，断点位于数据点处，并且系数数组 c 定义次数为 \(k\) 的多项式，使得 c[i, j] 是 (x - x[j])**(k-i) 在 x[j] 和 x[j+1] 之间的段上的系数。
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 的导数的近似多项式。使用 nu=1 参数评估原始样条，等效于在固定参数下评估 dspl
>>> 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)\) 的根 \(\pi/2 + \pi\,n\) 非常一致。请注意，默认情况下，它计算的根外推到插值区间 \(0 \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])
这与 \(\pm\arccos(1/2) + 2\pi\,n\) 的预期值非常一致。
可以使用 .integrate 方法计算分段多项式的积分，该方法接受积分的下限和上限。 例如，我们计算完全椭圆积分 \(K(m) = \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，则将其理解为 1D y 值的堆叠，这些值沿插值轴（默认值为 0）排列。后者通过 axis 参数指定，并且不变的是 len(x) == y.shape[axis]。 例如，我们扩展上面的椭圆积分示例，以使用 NumPy 广播来计算一系列 m 值的近似值
>>> 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) # the default is 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 样条：节点和系数#
b 样条函数（例如，通过 make_interp_spline 调用从数据构造）由所谓的节点和系数定义。
作为说明，让我们再次构造正弦函数的插值。 节点作为 BSpline 实例的 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(bspl.c)
7
约定是，对于 len(t) 个节点，有 len(t) - k - 1 个系数。某些例程（参见平滑样条部分）对 c 数组进行零填充，使得 len(c) == len(t)。这些额外的系数在求值时会被忽略。
我们强调，系数是以b-样条基给出的，而不是 \(1, x, \cdots, x^k\) 的幂基。
B 样条基元素#
B 样条是分段多项式，表示为b 样条基元素的线性组合，这些基元素本身是常见单项式的特定线性组合，\(x^m\)，其中 \(m=0, 1, \dots, k\)。
b 样条基通常比幂基在计算上更稳定，并且适用于各种应用，包括插值、回归和曲线表示。主要特点是这些基元素是局部化的，并且在由节点数组定义的区间之外等于零。
具体而言，一个 k 次的 b 样条基元素（例如，三次样条的 k=3）由 \(k+2\) 个节点定义，并且在这些节点之外为零。为了说明，绘制在某个区间上的一组非零基元素
>>> k = 3 # cubic splines
>>> t = [0., 1.4, 2., 3.1, 5.] # internal knots
>>> t = np.r_[[0]*k, t, [5]*k] # add boundary knots
>>> 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.]
接下来，将“观测值”取为
>>> xnew = [1, 2, 3]
并在稀疏 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. ]]