三次样条插值算法

前言

该算法涉及到对数据进行二次处理，插值运算等，属于一个非常常见的算法，本篇文章就来详细介绍和梳理一番。

正文

首先就是这个名字，三次大概能猜的出来，是阶数为3的函数，那么这里的样条指的是什么呢？样条来源于早期的工程制图，是为了将一些固定点连成一条光滑的曲线，采用具有弹性的木条固定在这些点上，如图：

这里要求连续且平滑，所以每两个点之间的木条其实就可以抽象为一个函数，也就是样条函数。

样条函数

从上图中可以清晰地知道，这些表示木条的函数，在数学上把它定义为一个分段多项式函数，这里以三次多项式为例：

即每两个点之间用一个多项式来表示，这些多项式的阶数相同，但是系数不同。

样条插值的最终目的就是求出这些多项式的系数，比如如果求出了$[x_1,x_2]$之间的函数，那么在这之间任意取一个点，其值就是可以确定了。

样条插值求解

想要求解样条函数，这里必须要遵循2个原则：连续和光滑。以三次样条插值为例，从$(x_1,y_1)$到$(x_{n+1},y_{n_1})$共n+1个点，从上面的公式可知一共有n个方程，共计有4n个未知数，即$a_1、b_1、c_1、d_1、.... a_n、b_n、c_n、d_n$，从解多元一次方程的知识来说，想解开4n个系数，需要构建4n个方程。

首先就是连续性，即每个分段函数都经过其2个端点，第一分段函数经过第1和第2个点，第二段分段函数经过第2和第3个点，以此类推，n段分段函数共可以得到2n个方程：

其次，根据光滑性的原则，每2段分段函数在连接处是平滑的，即一阶倒数$f_i^1(x) = 3a_ix^2 + 2b_ix + c_i$和二阶导数$f_i^2(x) = 6a_x + 2b_i$在连接点处相等，这样除了2个端点外，把$(x_2,y_2) ... (x_n,y_n)$带入上面方程中，可以得到2(n-1)个方程：

上面一共可以得到4n-2个方程，距离解出所有参数只少了2个方程，其实少的那2个方程也明白，就是2个端点的位置，因为这里无法拟定其导数是什么。

边界条件

对于边界的情况，为了凑出2个方程，有如下几种情况：

1. 自然边界，假设第1个点与最后一个点的二阶导数为0，即

2. 假设第1段函数和第2段函数的三阶导数在端点处相等，同时第n-1段函数和第n段函数在端点处也相等，可以得到$a_1 = a_2， a_{n-1} = a_{n-2}$。

3. 假设最后一个分段函数的一阶导数与二阶导数与第一个分段的一阶导数和二阶导数相等，这种假设特别适合在周期函数中：

4. 还有一种最简单的，强制$a_1 = a_n = 0$。

直接求解法

从上面可知，我们现在就可以构建出4n个方程，来解4n个参数，比如给定5个数据点为${(0,21),(1,24),(2,24),(3,18),(4,16)}$，利用自然边界来进行计算。这里共5个数据点，包含4个样条函数，每段样条函数的表达式为:$f_i = a_ix^3 + b_ix^2 + c_ix + d_i，i = 1,2,3,4$，这里需要解出16个参数，根据连续性，中间的点都符合左右2段样条函数，可以得到8个方程：

然后根据光滑原则，每个函数的一阶导数和二阶导数在连接处相等，可以得到6个方程：

再加上自然边界的条件，可以有2个方程：

由之前的知识我们知道，这里就是解16元一次方程组，而矩阵就是方程组的简写，所以上面方程组写成矩阵形式就是：

我们可以进行求解，得到16个参数。

推导公式

上面的过程非常容易理解，尤其是求解方程组，但是问题也很明显，就是矩阵不好列出来，现在推导出一种更为方便和容易计算的过程。

现在把区间$[a,b]$分为n个区间$[(x_0,x-1),(x_1,x_2),...,(x_{n-1},x_n)]$，共有n+1个点，其中$x_0 = a，x_n = b$，现在在每一个分段$[x_i，x_{i+1}]$上，$S(x)$都是一个三次方程。

假设分段函数定义为$S_i(x) = a_i + b_i(x-x_i) + c_i(x-x+i)^2 + d_i(x-x_i)^3$，这里的i取值是0到n-1，共n个方程，为什么要定义为这种，而非普通的写法呢，肯定是有原因的。

那么其一阶倒数就是$S^`_i(x) = b_i +2c_i(x-x_i) + 3d_i(x-x_i)^2$ ，二阶导数就是$S^{``}_i(x) = 2c_i + 6d_i(x-x_i)$。

上面的构造方法中，每一段的样条函数都和其开始点有关，这么构造是有一定道理的，接着分析。

每个样条函数肯定都经过它的起始点，所以$S_i(x) = a_i + b_i(x_i-x_i) + c_i(x_i-x+i)^2 + d_i(x_i-x_i)^3 = y_i$，可得$a_i = y_i$【方程A】。

然后是每个样条函数都经过它的终点，这里令$h_i = x_{i+1} - x_i$表示步长，由$S_i(x_{i+1}) = y_{i+1}$，推导出：$a_i + h_ib_i + h_i^2c_i + h_i^3d_i = y_{i+1}$。

每段样条函数的一阶导数在其连接处相等，所以第i段样条函数的导数，在其$x_{i+1}$上即$S_i^`(x_{i+1}) = b_i + 2c_i(x_{i+1} - x_i) + 3d_i(x_{i+1} - x_i)^2 = b_i + 2c_ih + 3d_ih^2$等于第i+1段上即$S^`_{i+1}(x_{i+1}) = b_{i+1} + 2c_i(x_{i+1} - x_{i+1}) + 3d_i(x_{i+1} - x_{i+1})^2 = b_{i+1}$，可以得到$b_i + 2c_ih + 3d_ih^2 = b_{i+1}$。

最后是二阶导数在连接处相等，可以推导出$2c_i + 6h_id_i = 2c_{i+1}$。

这里的推导到目前还是蛮有技巧性的，主要是引入了步长h，现在我们观察一下二阶导数，我们可以令$m_i = S^{``}_i(x_i) = 2c_i$【方程B】，所以上面二阶导推出的公式可以写成$m_i + 6h_id_i = m_{i+1}$，可以得到：

【方程C】

现在$a_i、c_i、d_i$都是y与m的表达方式，将其带入到$a_i + h_ib_i + h_i^2c_i + h_i^3d_i = y_{i+1}$，可以得到:

【方程D】

然后将a、b、c、d再带入$b_i + 2c_ih + 3d_ih^2 = b_{i+1}$，可以得到：

【最终方程】

这样我们就可以构建出以m为未知数的线性方程组，不得不说这里的做法非常巧妙。

具体代码

有了方程组就可以更加方便的列出矩阵了，直接按步骤计算即可：

//待修正的曲线，zeroX是需要插值修正后的x坐标
void correctRealSpectrum(QList<double> *realX,
                         QList<double> *realY,
                         QList<double> *zeroX) {
    int s = realX->size();
    const int size = s;
    double* x = new double[size];
    double* y = new double[size];
    std::vector<double> xx, yy;
    
    for (int i = 0; i < size; i++) {
        x[i] = realX->at(i);
        y[i] = realY->at(i);
        xx.push_back(zeroX->at(i));
    }
    int size_xx = xx.size();
    
    std::vector<double> dx;
    std::vector<double> dy;
        for (int i = 0; i < size - 1; i++) {
            //x的步长
            double temp_dx = x[i + 1] - x[i];
            dx.push_back(temp_dx);
            //y的步长
            double temp_dy = y[i + 1] - y[i];
            dy.push_back(temp_dy);
        }

        //构建一个sizexsize的矩阵，全部赋予0，是参数矩阵
        Eigen::MatrixXd H = Eigen::MatrixXd::Random(size, size);
        for (int i = 0; i < size; i++) {
            for (int j = 0; j < size; j++) {
                H(i, j) = 0;
            }
        }
        
        //构建一个列向量，是结果向量
        Eigen::VectorXd Y(size);
        for (int i = 0; i < size; i++) {
            Y(i) = 0;
        }
        
        //构建一个列向量，是待求参数向量
        Eigen::VectorXd M(size);
        for (int i = 0; i < size; i++) {
            M(i) = 0;
        }

        //2个端点为1
        H(0, 0) = 1;
        H(size - 1, size - 1) = 1;
        //按上面方程依次赋值
        for (int i = 1; i < size - 1; i++) {
            H(i, i - 1) = dx[i - 1];
            H(i, i) = 2 * (dx[i - 1] + dx[i]);
            H(i, i + 1) = dx[i];
            //注意这里的倍数是3倍
            Y(i) = 3 * (dy[i] / dx[i] - dy[i - 1] / dx[i - 1]);
        }
        try {
        //几种求解方式
//            M = H.colPivHouseholderQr().solve(Y);
            M = H.llt().solve(Y);
//              M = H.bdcSvd(Eigen::ComputeThinU|Eigen::ComputeThinV).solve(Y);
//            M = H.inverse() * Y;
//            M = H.bdcSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(Y);
        } catch (std::exception e) {
            logError() <<  "calculate error" << QString(e.what());
        }

        //每一段样条的4个参数
        std::vector<float> ai, bi, ci, di;
        for (int i = 0; i < size - 1; i++) {
            ai.push_back(y[i]);
            di.push_back((M(i + 1) - M(i)) / (3 * dx[i]));
            bi.push_back(dy[i] / dx[i] - dx[i] * (2 * M(i) + M(i + 1)) / 3);
            ci.push_back(M(i));
        }

        //进行插值
        for (int i = 0; i < size_xx; i++) {
            int k = 0;
            for (int j = 0; j < size - 1; j++) {
                if (xx[i] >= x[j] && xx[i] < x[j + 1]) {
                    k = j;
                    break;
                } else if (xx[i] == x[size - 1]) {
                    k = size - 1;
                }
            }
            double temp = ai[k] + bi[k] * (xx[i] - x[k]) + M(k) * pow((xx[i] - x[k]), 2) + di[k] * pow((xx[i] - x[k]), 3);
            yy.push_back(temp);
        }

    delete[] x;
    delete[] y;
    for (int i = 0; i < zeroX->length(); i++) {
        realX->replace(i, zeroX->at(i));
        realY->replace(i, yy.at(i));
    }
}

总结

样条插值算法总体来说理解不算困难，就是这个推导过程比较复杂，利用巧妙的构建函数来达到最后可以构建方便的矩阵进行求解。

更多算法相关总结移步专栏： juejin.cn/column/7415…