三次样条回归

样条函数

给定 $n$ 个参考点，拟合出一条分段平滑的三次样条曲线，每个参考点之间通过一个三次样条函数表示，其中第 $i$ 个样条函数以及其各阶导数的表达式如下，这样的参考点间样条函数共有 $n-1$ 个。
$\begin{aligned} f_i(x) &= a_i(x-x_i)^3 + b_i(x-x_i)^2 + c_i(x-x_i) + y_i + d_i \\ f^{'}_i(x) &= 3a_i(x-x_i)^2 + 2b_i(x-x_i) + c_i \\ f^{''}_i(x) &= 6a_i(x-x_i) + 2b_i \end{aligned}$

样条函数的约束

与样条插值不同的是，样条拟合不要求所有的拟合出的样条曲线一定通过参考点，假定第i个参考点坐标可以表示为: $\begin{bmatrix}x_i & y_i\end{bmatrix}^{T}$ ，因此 $y_i=f_i(x_i)$ 未必成立。由于连续性和平滑性的要求，相邻两个样条函数之间存在如下约束：
$\forall i \in \{2 \ldots n-1\}$ :
$\begin{aligned} f_{i}(x_{i}) - f_{i-1}(x_{i}) &= 0 \\ f^{'}_{i}(x_{i}) - f^{'}_{i-1}(x_{i}) &= 0 \\ f^{''}_{i}(x_{i}) - f^{''}_{i-1}(x_{i}) &= 0 \end{aligned}$

于此同时，满足如下边界条件：
$\begin{aligned} f^{''}_{1}(x_{1}) &= 0 \\ f^{''}_{n-1}(x_{n}) &= 0 \end{aligned}$

代价函数

样条拟合的主要目标是在满足上述约束的前体现，找到一个连续平滑且尽可能接近参考点的样条集合，下面构造误差向量：
$\begin{aligned} F_{i} &= \begin{bmatrix}f_i & f_i^{'} & f_i^{''}\end{bmatrix}^{T} \\ Y_{i} &= \begin{bmatrix}y_i & 0 & 0\end{bmatrix}^{T} \\ E_i^i &= F_i(x_i) - Y_i \\ E_i^{i+1} &= F_i(x_{i+1}) - Y_{i+1} \\ \end{aligned}$
其中 $E_i^i$ 表示第 $i$ 样条在 $i$ 参考点出的误差向量，其中 $E_i^{i+1}$ 表示第 $i$ 样条在 $i+1$ 参考点出的误差向量，通过误差向量，可以看到，样条拟合的优化目标是使得 $f_i(x)$ 靠近 $y$ ，同时，一阶导和二阶导的误差项是用来使两者尽可能较小，从而达到减小曲线变化的剧烈程度（曲率）的目的，使得回归过程对于采样点的噪声更加鲁棒。

接着考虑约束：
$\begin{aligned} G_i &= F_{i-1}(x_i)-F_i(x_i) \qquad \forall i \in \{2 \ldots n-1\}\\ G_1 &= \begin{bmatrix} g_1 \end{bmatrix} = \begin{bmatrix} f_1^{''}(x_1) \end{bmatrix} \\ G_n &= \begin{bmatrix} g_n \end{bmatrix}= \begin{bmatrix} f_{n-1}^{''}(x_n) \end{bmatrix} \\ \phi_{i} &= \begin{bmatrix}\phi_i^1 & \phi_i^2 & \phi_i^3\end{bmatrix}^{T} \\ \phi_1 &= [ \phi_1^1 ] \\ \phi_n &= [ \phi_n^1 ] \\ \end{aligned}$
其中 $G_i$ 表示样条间约束关系（连续、平滑）； $G_1$ ， $G_n$ 分别表示边缘约束； $\phi$ 表示拉格朗日乘子（Lagrange multiplier）。

综上，约束优化问题可以表示为：
$\begin{aligned} O(X) = &\sum_{i=1}^{n-1}({E_i^i}^TWE_i^i + {E_i^{i+1}}^TWE_i^{i+1}) +\sum_{i=1}^{n}\phi_{i}^TG_i \end{aligned}$
其中 $X$ 定义为：
$\begin{aligned} \beta_i &= \begin{bmatrix}a_i & b_i & c_i & d_i\end{bmatrix}^{T} \\ X &= \begin{bmatrix}\phi_1^T & \beta_1^T & \phi_2^T & \beta_2^T &\ldots & \phi_{n-1}^T & \beta_{n-1}^T & \phi_n^T\end{bmatrix}^{T} \end{aligned}$
$O$ 对 $X$ 求导和二阶导：
$\begin{aligned} \frac{\partial O}{\partial X} = \begin{bmatrix} \frac{\partial O}{\partial X_1} & \ldots & \frac{\partial O}{\partial X_m} \end{bmatrix}^{T} \end{aligned}$
$\begin{aligned} H = \frac{\partial^2 O}{\partial X^2} = \frac{\partial}{\partial X}(\frac{\partial O}{\partial X})=\begin{bmatrix} \frac{\partial^2 O}{\partial X_1^2} & \frac{\partial^2 O}{\partial X_1\partial X_2} & \ldots & \frac{\partial O}{\partial X_1\partial X_m} \\ \frac{\partial^2 O}{\partial X_2\partial X_1} & \frac{\partial^2 O}{\partial X_2^2} & \ldots & \frac{\partial O}{\partial X_2\partial X_m} \\ \vdots & \vdots & \ddots &\vdots \\ \frac{\partial^2 O}{\partial X_m\partial X_1} & \frac{\partial^2 O}{\partial X_m\partial X_2} & \ldots & \frac{\partial^2 O}{\partial X_m^2} \end{bmatrix} \end{aligned}$

求解

现在显然是一个二次优化问题，求解过程如下：
$\begin{aligned} B &= -\left.\frac{\partial O}{\partial X}\right| _{X=0} \\ B &= \begin{bmatrix} B_1^{\phi} & B_1^{\beta} & B_2^{\phi} & B_2^{\beta} & \ldots & B_{n-1}^{\phi} & B_{n-1}^{\beta} & B_n^{\phi} \end{bmatrix}^T \\ B_i^{\phi} &= -\left.\frac{\partial O}{\partial \phi_i}\right| _{\phi_i=0} = \begin{bmatrix} y_i-y_{i-1} && 0 && 0 \end{bmatrix}^T\\ B_i^{\beta} &= -\left.\frac{\partial O}{\partial \beta_i}\right| _{\beta_i=0} = \begin{bmatrix} 2w_1(y_{i+1}-y_i)(x_{i+1}-x_i)^3 + 6w_2y_{i+1}^{'}(x_{i+1}-x_i)^2 + 12w_3y_{i+1}^{''}(x_{i+1}-x_i) \\ 2w_1(y_{i+1}-y_i)(x_{i+1}-x_i)^2 + 4w_2y_{i+1}^{'}(x_{i+1}-x_i) + 4w_3(y_i^{''} + y_{i+1}^{''}) \\ 2w_1(y_{i+1} - y_i)(x_{i+1}-x_i) + 2w_2(y_i^{'} + y_{i+1}^{'}) \\ 2w_1(y_{i+1} - y_i) \end{bmatrix} \end{aligned}$

$\begin{aligned} HX = B \\ X = H^{-1}B \end{aligned}$
其中 $H$ 矩阵的形式如下：

$\begin{aligned} H = \begin{bmatrix} 0 & T_1^T \\ T_1 & A_1 & R_1 \\ & R_1^T & 0 & T_2^T \\ & & T_2 & A_2 & R_2 \\ & & & R_2^T & 0 & T_3^T \\ & & & & \vdots & \vdots & \vdots \\ & & & & & R_{n-2}^T & 0 & T_{n-1}^T \\ & & & & & & T_{n-1} & A_{n-1} & R_{n-1} \\ & & & & & & & R_{n-1}^T & 0 \end{bmatrix} \\ \end{aligned}$

$\begin{aligned} A_i^{11} &= 2w_1(x_{i+1}-x_i)^6+18w_2(x_{i+1}-x_i)^4+72w_3(x_{i+1}-x_i)^2 \\ A_i^{12} &= 2w_1(x_{i+1}-x_i)^5+12w_2(x_{i+1}-x_i)^3+24w_3(x_{i+1}-x_i) \\ A_i^{13} &= 2w_1(x_{i+1}-x_i)^4+6w_2(x_{i+1}-x_i)^2 \\ A_i^{14} &= 2w_1(x_{i+1}-x_i)^3 \\ \\ A_i^{21} &= 2w_1(x_{i+1}-x_i)^5+12w_2(x_{i+1}-x_i)^3+24w_3(x_{i+1}-x_i) \\ A_i^{22} &= 2w_1(x_{i+1}-x_i)^4+8w_2(x_{i+1}-x_i)^2+16w_3 \\ A_i^{23} &= 2w_1(x_{i+1}-x_i)^3+4w_2(x_{i+1}-x_i) \\ A_i^{24} &= 2w_1(x_{i+1}-x_i)^2 \\ \\ A_i^{31} &= 2w_1(x_{i+1}-x_i)^4+6w_2(x_{i+1}-x_i)^2 \\ A_i^{32} &= 2w_1(x_{i+1}-x_i)^3+4w_2(x_{i+1}-x_i) \\ A_i^{33} &= 2w_1(x_{i+1}-x_i)^2+4w_2 \\ A_i^{34} &= 2w_1(x_{i+1}-x_i) \\ \\ A_i^{41} &= 2w_1(x_{i+1}-x_i)^3 \\ A_i^{42} &= 2w_1(x_{i+1}-x_i)^2 \\ A_i^{43} &= 2w_1(x_{i+1}-x_i) \\ A_i^{44} &= 4w_1 \\ \\ T_i &= \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & -2 \\ 0 & -1 & 0 \\ -1 & 0 & 0 \end{bmatrix} \\ R_i &= \begin{bmatrix} (x_{i+1}-x_i)^3 & 3(x_{i+1}-x_i)^2 & 6(x_{i+1}-x_i) \\ (x_{i+1}-x_i)^2 & 2(x_{i+1}-x_i) & 2 \\ (x_{i+1}-x_i) & 1 & 0 \\ 1 & 0 & 0 \end{bmatrix} \\ T_1 &= \begin{bmatrix} 0 & 2 &0 & 0 \end{bmatrix}^T \\ R_{n-1} &= \begin{bmatrix} 6(x_n-x_{n-1}) & 2 &0 & 0 \end{bmatrix}^T \end{aligned}$
比较稀疏，可以用稀疏矩阵求逆的方式求解。