Leon's Blog - 逐次线性插值

Lagrange算法是在n维线性空间中找了一组形式简单对称的基函数，使得插值问题容易计算，一般情况下这个算法已经足够好用了。但存在一个缺点是，如果由于精度不满足要求而临时增加一个插值节点时，所有基函数都得重新计算。我们想，那能不能找到一种改进的算法来避免这种资源浪费呢？实际上这就需要逐次线性插值法来解决，也就是说将高次插值问题逐步转化为线性插值，或者说将两组由k个节点算得的Lagrange插值作为形式节点，再进行一次线性插值来得到k+1个节点时的值。

先做如下符号约定：令$I_{i_1 i_2 \cdots i_n}(x)$表示函数$f(x)$关于节点$i_1, i_2, \cdots ,i_n$的n-1次插值多项式，$I_{i_k}(x)$是零次多项式，$I_{i_k}(x)=f(x_{i_k})=y_{i_k}$，$i_1,\cdots, i_n$均为非负整数。下面通过一个例子来说明“形式节点”的线性插值是如何做到的，演示最简单的情形：由2个节点增加到3个节点的插值问题，也就是由线性插值到抛物线插值。已知线性插值

$$I_{0,1}=\frac{x-x_1}{x_0-x_1}y_0+\frac{x-x_0}{x_1-x_0}y_1=L_1(x),$$

$$I_{1,2}=\frac{x-x_1}{x_2-x_1}y_2+\frac{x-x_2}{x_1-x_2}y_1,$$

求节点为$x_0, x_1,x_2$的$I_{0,1,2}=L_2(x)$.

根据n阶Lagrange插值多项式为n次齐次多项式，可以写出如下组合关系：

$$I_{0,1,2}=a(x-x_2)I_{0,1}(x)+b(x-x_0)I_{1,2}(x),$$

满足$I_{0,1,2}(x_i)=y_i$，则由$y_0和y_1$可以得到$a=\frac{1}{x_0-x_2}$， $b=\frac{1}{x_2-x_0}$，可以验证对于$y_2$自洽成立。从而，

$$I_{0,1,2}(x)=\frac{x-x_2}{x_0-x_2}I_{0,1}(x)+\frac{x-x_0}{x_2-x_0}I_{1,2}(x)$$

此式即可以看成由两点$(x_0,I_{0,1})，(x_2,I_{1,2})$的线性插值。

对于三个点的插值可以写成$I_{0,1,2}(x)=\frac{(x-x_0)I_{1,2}(x)-(x-x_2)I_{0,1}(x)}{x_2-x_0}$；一般地，有

$$I_{0,1,\cdots,n}(x)=\frac{(x-x_0)I_{1,2,\cdots,n}(x)-(x-x_2)I_{0,1,\cdots,n-1}(x)}{x_2-x_0},n=0,1,\cdots$$

更一般地，可以任意选取插值节点的顺序，公式表达上用$x_{i_k}$来代替$x_k$，即$I_{i_0,i_1,\cdots,i_n}(x)$。上式可以直观地接受，这里我不打算在数学上详细求证（其实证明起来也并不复杂，已然是次数不超过k的多项式，只需证明其穿过相应节点即可）。把递推关系式写成上面这种形式，是为了给出在一般化的递推关系，具体操作可以是多种多样的，但在实际运用中一般只选用Aitken方法或者Neville方法。也就是说，下面介绍的这两种插值方法实属同一原理，只是在具体的递推操作中选用不同的流程。

Aitken方法

如果已有$0,1,\cdots,k$节点，添加第$l$节点，把上述改写成

$$I_{0,1,\cdots,k,l}(x)=I_{0,1,\cdots,k}(x)+\frac{I_{0,1,\cdots,k-1,l}(x)-I_{0,1,\cdots,k}(x)}{x_l-x_k}(x-x_k)$$

显然，$I_{0,1,\cdots,k,l}(x_i)=I_{0,1,\cdots,k}(x_i)=f(x_i),i=0,1,\cdots,k-1$ 成立。而$x=x_k时，I_{0,1,\cdots,k}(x_k)=I_{0,1,\cdots,k,l}(x_k)=f(x_k)$，且$x=x_l时，I_{0,1,\cdots,k,l}(x_k)=I_{0,1,\cdots,k}(x_l)+\frac{f(x_l)-I_{0,1,\cdots,k}(x_l)}{x_l-x_k}(x_l-x_k)=f(x_l)$。也就说这样的插值多项式满足插值条件。这种递推式就是Aitken逐次线性插值多项式。

我们可以由上面的递推式看到，计算时利用k=0到k=n-1逐次求得所需的插值多项式。例如k=0时，为线性插值，k=1时插值节点为$x_0,x_1,x_l$三点，公式为

$$I_{0,1,l}(x)=I_{0,1}(x)+\frac{I_{0,l}(x)-I_{0,1}(x)}{x_l-x_1}(x-x_1)$$

现在我们来直观地体会一下Aitken方法是如何操作的：

Aitken interpolation

Neville方法

或者改写成

$$I_{0,1,\cdots,k,k+1}(x)=I_{0,1,\cdots,k}(x)+\frac{I_{1,\cdots,k,k+1}(x)-I_{0,1,\cdots,k}(x)}{x_{k+1}-x_0}(x_{k+1}-x_0)$$

此即Neville方法，直观的计算过程如下：

Neville interpolation

可以看出若增加一个节点，便多计算一行，前面的结果可以重复利用，这便于快速应对精度不足的情形。另外，对于选取哪些节点来完成插值也比较自由，有利于自动选并逐步比较精度。Lagrange插值和逐次线性插值的算法复杂度都是$O(n^2)$，细致的分析给出逐次线性插值要比Lagrange插值略快几倍。

Python编程实现

对于递推公式表达的问题，可以使用递归算法来编程实现，例如我们把Neville方法中的递推公式写成对称的两项，考虑到通用性我们用$i_j$代替$j$，

$$I_{i_0,i_1,\cdots,i_k,i_{k+1}}(x)=I_{i_0,i_1,\cdots,i_k}(x)\frac{x-x_{i_{k+1}}}{x_{i_0}-x_{i_{k+1}}}+I_{i_1,\cdots,i_k,i_{k+1}}(x)\frac{x-x_{i_0}}{x_{i_{k+1}}-x_{i_0}},$$

核心是我们用相邻的两项来计算更高一次的插值。 Python代码如下：

def neville(xVal, yVal, x):
"""
Polynomial Interpolation: Neville's Method    
"""
if len(yVal) == 1:
    return yVal[0]
else:
    I_0 = aitken(xVal[1:], yVal[1:], x)*(x - xVal[0])/(xVal[-1] - xVal[0])
    I_n = aitken(xVal[:-1], yVal[:-1], x)*(x - xVal[-1])/(xVal[0] - xVal[-1])
    return I_0 + I_n

# example:
xVal = [4, 9, 16]
yVal = [2, 3, 4]
print(neville(xVal, yVal, 11))
#3.33333..

有时候我觉得递归算法是不直观的实施方式，我们可以直接对照Neville方法的计算表来写代码，对于n个节点，我们更新n次列表，每次计算形式节点数目会减少一个，直到获取目标值。这里我重画了计算示意表格，便于理解数值计算的代码实现。

Neville interpolation

我们可以用二重循环来完成Neville插值（记得Lagrange插值也是二重循环实现，所以复杂度都是$O(n^2)$），python代码如下：

def neville2(xVal, yVal, x):
    """
    Polynomial Interpolation: Neville's Method
    Another implement   
    """
    n = len(xVal)
    y = yVal.copy()
    for m in range(1,n):
        for i in range(n-m):
            y[i] = ((x - xVal[i+m])*y[i]+(xVal[i]-x)*y[i+1])/(xVal[i]-xVal[i+m])

    return y[0]

对于Aitken方法的Python实现，如果有时间再写。