[Algorithm in Python] 三次样条内插算法演示

当我们需要将有序资料点视觉化，甚至是估计点之间的数值时，我们会使用内插（Interpolation）。最简单的方式就是用直接做线性内插，不过得出来的结果就是很粗糙的直线。如果要得到连续且平滑的内插值，则会需要用到曲线去做内插，比较常见的方法是用二次或三次函数。

使用Scipy即可轻鬆的实现各种曲线内插：

from scipy.interpolate import interp1dimport matplotlib.pyplot as pltimport numpy as np# Pseudo datax = np.linspace(0, 10, 11)y = np.cos(-x**2/9.0)# Cubic spline interpolationx_fine = np.linspace(0, 10, 1000)f_interpolate = interp1d(x, y, kind='cubic')y_fine = f_interpolate(x_fine)# Plotplt.plot(x, y, 'b-', alpha=0.8, label='Linear')plt.plot(x_fine, y_fine, 'r-', alpha=0.8, label='Cubic')plt.plot(x, y, 'ko', alpha=0.8, label='Data')plt.legend()plt.show()

本篇笔记使用Numpy去实作一个三次样条内插（Cubic Spline Interpolation）的算法，并着重于步骤思路演示。

原理

和线性内插一样，两个资料点间会决定一条线段，所以如果今天有n个资料点，就会有n-1条线段，而在三次内插中，每一条线段是由三次多项式代表，如第i条线段表示为：

▲两个点之间的线段都是不同的三阶多项式

我们接下来要做的事，就是解出所有线段的係数才能够得到线段的函式。一个三次多项式中会有4个未知数，也就是abcd，而假设我们有n笔资料，就会有n-1条线段，总共就会有4(n-1)个未知数，也就是说我们会需要4(n-1)条方程式去解所有未知数。
为演示方便起见，这里设一个n=3的例子：x = [0, 1, 2], y = [1, 3, 2]

所以我们就会有两条线段：

条件一：端点

最简单的条件就是我们知道线段的左右两端一定要黏在点上，所以每个线段可以产生出2个条件，如此我们可以得到2(n-1)个方程式。

左端：

右端：

条件二：平滑

接下来我们还要考虑一件事情，就是两个线段连接处要连续且平滑，这个条件透过微分达成。两个线段在同一点微分相等就代表在该点时斜率会相同，所以就会有“平”的结果，而如果要“滑”的话，我们就可以进一步用二阶微分，让该点处的斜率变化量相等而更圆滑。两个线段可以提供1个微分等式，因此一阶微分可以提供n-2个条件，二阶微分也提供n-2个，共有2(n-2)个条件。

一阶微分

二阶微分

那既然可以用二阶微分，那可否使用三阶甚至四阶微分呢？可以这么想，三次多项式在二阶微分后，会变成一条直线，而两条直线在资料点上本质上就会是个折点，会跟平滑矛盾，所以在三次多项式三阶微分以上的条件都是没有几何意义的。

▲三阶微分（Third）几何上就不会连接起来

条件三：边界

那么有了前两个条件一共4n-6个方程式，就还缺2个条件，而剩下这两个条件的花样就多了，根据不同的设计可以得到不同的几何效果。就以Scipy的scipy.interpolate.CubicSplineAPI来说，就提供了not-a-knot、periodic、clamped、natural，后面的笔记会各别简介，大致上都是拿头和尾的2条线段做操作，毕竟刚刚找微分条件的时候，第一个和最后一个资料节点没有做出贡献嘛。以下会拿natural做演示。
natural的方法就是把头跟尾2线段，做二阶微分然后让它等于0。这在几何意义上是什么意思呢，简化来说就是末端线段接到节点时，会接近一条直线，不会有曲线的情形。有兴趣深入探讨的话，可以到维基百科查阅反曲点。

边界条件

求解

至此，我们终于集齐4(n-1)个龙珠方程式啦！不过在求解之前，我们要先将上面那一大坨方程式整理成矩阵形式。矩阵虽然是个很反人类的数学符号，但和电脑却是很情投意合呢。

矩阵形式

而要求解其实只需要算出矩阵A的反矩阵，此笔记会使用Numpy去处理。

程式

步骤

在编写程式码之前，先把大致要做的步骤列出来：

收集端点相等的方程式左端右端收集二阶导数的方程式收集三阶导数的方程式收集边界条件的方程式矩阵求解计算曲线内插

Numpy演示

import matplotlib.pyplot as pltimport numpy as np# Cook pseudo datax = np.linspace(0, 10, 11)y = np.cos(-x**2/9.0)# Sure the sequence is sorted by xsort_i = np.argsort(x)x = x[sort_i]y = y[sort_i]# 4*(n-1) unknown variablesn = len(x)# Ax = bmatrix = []  # Atarget = []  # b# (1) Endpoint equality: 2*(n-1) equationsfor i in range(n-1):  # n-1 curves    for j in range(2):  # left and right        row = np.zeros(4*(n-1))        row[4*i:4*i+4] = np.array([x[i+j]**3, x[i+j]**2, x[i+j], 1])        matrix.append(row)        target.append(y[i+j])# (2) First derivative equalityfor i in range(n-2):  # n-2 equations    row = np.zeros(4*(n-1))    row[4*i:4*i+8] = np.array([3*x[i+1]**2, 2*x[i+1], 1, 0] + [-3*x[i+1]**2, -2*x[i+1], -1, 0])    matrix.append(row)    target.append(0)# (3) Second derivative equalityfor i in range(n-2):  # n-2 equations    row = np.zeros(4*(n-1))    row[4*i:4*i+8] = np.array([6*x[i+1], 2, 0, 0] + [-6*x[i+1], -2, 0, 0])    matrix.append(row)    target.append(0)# (4) Boundary condition# the leftmost noderow = np.zeros(4*(n-1))row[:4] = np.array([6*x[0], 2, 0, 0])matrix.append(row)target.append(0)# the rightmost noderow = np.zeros(4*(n-1))row[-4:] = np.array([6*x[-1], 2, 0, 0])matrix.append(row)target.append(0)# Cook martix and targetmatrix = np.stack(matrix)target = np.array(target)# (5) Solve the equationsinverse_matrix = np.linalg.inv(matrix)abcd = np.dot(inverse_matrix, target)# Define spline functiondef spline(x, a, b, c, d):    return a*x**3 + b*x**2 + c*x + d# (6) Interplationx_fine = np.linspace(0, 10, 1000)  # fine gridsy_fine = []ci = 0  # curve ifor x_i in x_fine:    if x_i > x[ci+1] and ci < (n-1):        ci += 1    a, b, c, d = abcd[4*ci:4*ci+4]    y_fine.append(spline(x_i, a, b, c, d))# Plotplt.figure(dpi=100)plt.plot(x, y, 'ko')plt.plot(x_fine, y_fine, 'r-')plt.show()

最后插分的部分提供了一个比较naive但直觉的写法。因为效能的关係，最好尽量使用Numpy做操作。

Scipy演示

如前面所介绍的，Scipy的scipy.interpolate.CubicSpline提供了三次内插，而其中bc_type引数可以指定使用什么样的边界条件。

from scipy.interpolate import CubicSplineimport matplotlib.pyplot as pltimport numpy as npx = np.linspace(0, 10, 11)y = np.cos(-x**2/9.0)y[-1] = y[0]  # make the sequence periodicx_fine = np.linspace(0, 10, 1000)bc_types = ['not-a-knot', 'periodic', 'clamped', 'natural']for s in bc_types:    f = CubicSpline(x, y, bc_type=s)    plt.plot(x_fine, f(x_fine), '-', alpha=0.8, label=s)plt.plot(x, y, 'ko')plt.grid(linestyle='--', alpha=0.3)plt.legend()plt.show()

not-a-knot: 最后与倒数第二条线段为相同线段，如此一来就少了4个方程式条件要满足。虽然API文件推荐这个模式，但笔者觉得容易造成尾端有过拟合现象。periodic: 在两末端相同的时使用，如此第一与最后一点也可以加入一阶与二阶导数条件。正如其名，推荐在已知资料点是有週期性的时候再使用。clamped: 两末端点的一阶微分为0，可以理解为让曲线的末端点与水平线相切。natural: 两末端点的二阶微分为0，可以理解为让曲线接近末端点的部分趋近直线。

参考

Lagrange Polynomial Interpolation
scipy.interpolate.CubicSpline
Legacy interface for 1-D interpolation (interp1d)
插值-样条插值