分段多项式和样条#
一维插值例程在上一节中讨论的，通过构建某些分段多项式工作：插值范围被所谓的断点划分成区间，每个区间上都有一个特定的多项式。这些多项式片段在断点处以预定义的平滑度匹配：对于三次样条是二阶导数，对于单调插值器是一阶导数等等。
阶数为 \(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[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])
这与 \(\pi/2 + \pi\,n\) 的根 \(\cos(x) = \sin'(x)\) 非常吻合。请注意，默认情况下，它计算的根是外推到插值区间 \(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，则它被理解为一堆一维 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样条是分段多项式，表示为b样条基元素的线性组合——这些基元素本身是普通单项式 \(x^m\)（其中 \(m=0, 1, \dots, k\)）的某些线性组合。
B样条的性质在文献中已有详细描述（例如，参见 BSpline 文档字符串中列出的参考文献）。出于我们的目的，只需知道B样条函数由系数数组和所谓的节点数组唯一确定即可，这些节点可能与数据点 x 重合，也可能不重合。
具体来说，阶数为 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
如果需要，B样条可以使用 PPoly.from_spline 方法转换为 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. ]]
伯恩斯坦多项式，BPoly#
对于 \(t \in [0, 1]\)，阶数为 \(k\) 的伯恩斯坦基多项式定义为
其中 \(C_k^a\) 是二项式系数，且 \(a=0, 1, \dots, k\)，因此有 \(k+1\) 个阶数为 \(k\) 的基多项式。
一个 BPoly 对象表示一个分段伯恩斯坦多项式，通过断点 x 和系数 c 来定义：c[a, j] 是在 x[j] 和 x[j+1] 之间的区间上，对应于 \(b(t; k, a)\) 的系数。
BPoly 对象的接口与 PPoly 对象的接口非常相似：它们都可以进行评估、微分和积分。
BPoly 对象的另一个附加功能是备用构造函数 BPoly.from_derivatives，它从数据值和导数构造一个 BPoly 对象。具体来说，b = BPoly.from_derivatives(x, y) 返回一个可调用对象，该对象插值所提供的值（b(x[i]) == y[i])），并具有所提供的导数（b(x[i], nu=j) == y[i][j]）。
此操作类似于 CubicHermiteSpline，但它更灵活，因为它可以在不同数据点处理不同数量的导数；即，y 参数可以是不同长度数组的列表。有关进一步讨论和示例，请参见 BPoly.from_derivatives。
基之间的转换#
原则上，分段多项式所有三种基（幂基、伯恩斯坦基和B样条）都是等价的，并且一种基中的多项式可以转换为不同的基。在不同基之间进行转换的一个原因是并非所有基都实现了所有操作。例如，寻根操作仅针对 PPoly 实现，因此要找到 BSpline 对象的根，您需要首先转换为 PPoly。有关转换的详细信息，请参见方法 PPoly.from_bernstein_basis、PPoly.from_spline、BPoly.from_power_basis 和 BSpline.from_power_basis。
然而，在浮点运算中，转换总是会带来一些精度损失。这种损失是否显著取决于具体问题，因此建议在基之间进行转换时谨慎操作。