分段多項式與樣條#
1D 插值常式 在前一節討論的，透過建構特定的分段多項式來運作：插值範圍透過所謂的斷點分割成多個區間，並且每個區間上都有特定的多項式。然後，這些多項式片段在斷點處以預定義的平滑度匹配：三次樣條的二階導數，單調插值器的一階導數等等。
次數為 \(k\) 的多項式可以被認為是 \(k+1\) 個單項式基底元素 \(1, x, x^2, \cdots, x^k\) 的線性組合。在某些應用中，考慮替代（如果形式上等價）基底會很有用。scipy.interpolate 中實作了兩個常用的基底，分別是 B 樣條 (BSpline) 和 Bernstein 多項式 (BPoly)。B 樣條通常用於例如非參數迴歸問題，而 Bernstein 多項式則用於建構 Bezier 曲線。
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 工廠函數從資料建構。
最後，Bernstein 多項式表示為 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)\) 的根 \(\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. ]]