多项式插值 | 学而知之

插值是一种通过已知的数据点，求新数据点的过程或方法。比如我们已经通过采样、或者实验等方式，获取到了 $n+1$ 个不同的采样点 $(x_i, y_i)$，其中 $x_i$ 为输入，$y_i$ 为输出。现在遇到了一个新的输入 $x_t$，如何预测对应的输出 $y_t$ 呢？插值就是一种能解决此问题的方法。而多项式由于形式简单，便于计算，是最常用的插值函数。使用多项式作为插值函数，就是多项式插值。

本文将介绍多项式插值的基本原理，以及常用的插值方法，其中包括了拉格朗日插值法和牛顿插值法。对于拉格朗日插值和牛顿插值，本文不仅介绍了原理，还给出了 Python 代码的实现。虽然代码基于 numpy，但是并没有使用 numpy 的高级特性。如果读者想替换为标准的 Python list，也是比较容易的。

1. 多项式求值

使用多项式插值，必然涉及到多项式求值问题：已知一个 $n$ 次多项式的表达式 $f(x) = a_0 + a_1 \cdot x + a_2 \cdot x^2 + ... + a_n \cdot x^n$，再给一个特定的点 $x$，如何计算对应的 $f(x)$ 呢？当然可以直接按多项式的定义去求值，但是代码比较繁琐。这里使用霍纳法则，转换一下多项式的形式：

$f(x) = (((a_n \cdot x + a_{n-1})x + ... + a_2)x + a_1)x + a_0$ 在使用霍纳法则迭代计算多项式时，代码会比较简洁，时间复杂度为 $O(n)$。

# 霍纳法则计算多项式，a 是多项式的系数，x 是给定的点
def Polynomial(a, x):
    r = 0
    for i in range(len(a)-1, -1, -1):
        r = x * r + a[i]
    return r

2. 多项式插值基本思路

回到本文刚开始的问题：已知 $n+1$ 个不同的采样点 $(x_i, y_i)$，其中 $x_i$ 为输入，$y_i$ 为输出。现在遇到了一个新的输入 $x_t$，如何预测对应的输出 $y_t$ 呢？

首先寻找一个多项式函数 $f(x)$，使得 $f(x)$ 恰好经过这 $n+1$ 个点。然后假定对于任意的输入 $x$，输出均为 $f(x)$。最后对于给定的输入 $x_t$，按照假定条件，输出 $y_t = f(x_t)$。

3. 多项式插值原理

定理1： 对于 $n+1$ 个不同的坐标$(x_0, y_0)$，$(x_1, y_1)$...，$(x_n, y_n)$。存在一个唯一的不超过 $n$ 次的多项式，恰好经过这 $n+1$ 个点。

证明： 使用待定系数法，设多项式为： $f(x) = a_0 + a_1 \cdot x + a_2 \cdot x^2 + ... + a_n \cdot x^n$，由于 $f(x)$ 经过这 $n+1$ 个点。将 $(x_i, y_i)$ 坐标带入 $f(x)$ 的表达式，可构建如下方程组 \[ f(x_0) = a_0 + a_1 \cdot x_0 + a_2 \cdot x_0^2 + ... + a_n \cdot x_0^n = y_0 \\ f(x_1) = a_0 + a_1 \cdot x_1 + a_2 \cdot x_1^2 + ... + a_n \cdot x_1^n = y_1 \\ ... \\ f(x_n) = a_0 + a_1 \cdot x_n + a_2 \cdot x_n^2 + ... + a_n \cdot x_n^n = y_n \] 写成矩阵形式 $\mathbf{X a} = \mathbf{y}$，其中 $\mathbf{a} = (a_0, a_1, a_2, ..., a_n)^\mathrm{T}$，$\mathbf{y} = (y_0, y_1, y_2, ..., y_n)^\mathrm{T}$，$\mathbf{X}$ 是一个 $(n+1) \times (n+1)$ 的矩阵。

\[ \begin{pmatrix} 1 & x_0 & x_0^2 & \cdot & x_0^n \\ 1 & x_1 & x_1^2 & \cdot & x_1^n \\ 1 & x_2 & x_2^2 & \cdot & x_2^n \\ \cdot &\cdot & \cdot & \cdot & \cdot \\ 1 & x_n & x_n^2 & \cdot & x_n^n \end{pmatrix} \cdot \begin{pmatrix} a_0 \\ a_1 \\ a_2 \\ \cdot \\ a_n \end{pmatrix} = \begin{pmatrix} y_0 \\ y_1 \\ y_2 \\ \cdot \\ y_n \end{pmatrix} \]

可以看出，$\mathbf{X}$ 是一个范德蒙矩阵，行列式不为 $0$，所以该方程组有唯一解：$\mathbf{a} = \mathbf{X}^{-1}\mathbf{y}$。所以存在一个唯一的多项式，经过这 $n+1$ 个点。

虽然可以使用高斯消元法对范德蒙矩阵矩阵求逆，进而求出多项式插值函数，但是该方法计算的复杂度比较高，时间复杂度为 $O(n^3)$，而且精度不高，实际很少使用。更常用的方法是拉格朗日插值和牛顿插值。由定理 1 可知，不管使用什么方法，最终的插值多项式都是一样的，只是计算的复杂程度不同而已。

4. 拉格朗日插值

拉格朗日使用构造法来构造多项式插值函数。基本思想是寻找 $n+1$ 个多项式基函数，其中第 $i$ 个基函数在 $x_i$ 处取值为 $1$，在其他 $x_j$ 处 ($j \neq i$) 取值为 $0$。即第 $i (0 \leq i \leq n)$ 个多项式基函数 $\omega_i(x)$ 满足如下特性：

\[ \omega_i(x_k) = \left\{ \begin{aligned} &1 , i = k \\ &0 , i \neq k \end{aligned} \right. \]

构造多项式函数 $f(x) = \sum\limits_{i=1}^{n}y_i \cdot {\omega_i(x)}$，那么显然函数 $f(x)$ 满足插值点 $f(x_i) = y_i$。由于每个基函数 $\omega_i(x)$ 是一个关于 $x$ 的 $n$ 次多项式，所以 $f(x)$ 也一个关于 $x$ 的 $n$ 次多项式。关键的问题是如何寻找基函数。

对于 $i (0 \leq i \leq n)$，拉格朗日选取的第 $i$ 个基函数： $\displaystyle \omega_i(x) = \prod\limits_{j=0, j \neq i}^{j=n} \frac{(x-x_j)}{(x_i-x_j)}$。

可以看到，基函数 $\omega_i(x)$ 是一个关于 $x$ 的 $n$ 次多项式。下面验证一下是否满足基函数特性。将 $x_i$ 带入第 $i$ 个基函数可得：

$\displaystyle \omega_i(x_i) = \prod\limits_{j=0,j \neq i}^{j=n} \frac{(x_i-x_j)}{(x_i-x_j)}$，约分后，可得 $\omega_i(x_i) = 1$

对于$k \neq i$ 且 $0 \leq k \leq n$，将 $x_k$ 带入第 $i$ 个基函数可得： $\displaystyle \omega_i(x_k) = \prod\limits_{j=0,j \neq i}^{j=n} \frac{(x_k-x_j)}{(x_i-x_j)}$。

在该式中，分子中必然存在一个 $j$ 使得 $k = j$，即$x_k - x_j = 0$，那么分子的乘积也为 $0$，即 $\omega_i(x_k) = 0$，所以满足基函数特性。

拉格朗日多项式插值函数表达式为：$f(x) = \sum\limits_{i=1}^{n}y_i \cdot {\omega_i(x)}$

举例

以上的介绍有些抽象，看一个例子就会清晰一些。假设有 $3$ 个点 $(1,2)$，$(3,12)$，$(4,23)$，求一个二次函数 $y = a_0 + a_1 \cdot x + a_2 \cdot x^2$ 经过这 $3$ 个点。首先找到 3 个拉格朗日基：

$\displaystyle \omega_0(x) = \frac{x-x_1}{x_0-x_1} \cdot \frac{x-x_2}{x_0-x_2} = \frac{x-3}{1-3} \cdot \frac{x-4}{1-4} = \frac{(x-3)(x-4)}{6}$

$\displaystyle \omega_1(x) = \frac{x-x_0}{x_1-x_0} \cdot \frac{x-x_2}{x_1-x_2} = \frac{x-1}{3-1} \cdot \frac{x-4}{3-4} = -\frac{(x-1)(x-4)}{2}$

$\displaystyle \omega_2(x) = \frac{x-x_0}{x_2-x_0} \cdot \frac{x-x_1}{x_2-x_1} = \frac{x-1}{4-1} \cdot \frac{x-3}{4-3} = \frac{(x-1)(x-3)}{3}$

所以最终的表达式：

\[ \begin{aligned} f(x) &= y_0 \cdot \omega_0(x) + y_1 \cdot \omega_1(x) + y_2 \cdot \omega_2(x) \\ &= 2 \cdot \frac{(x-3)(x-4)}{6} + 12 \cdot -\frac{(x-1)(x-4)}{2} + 23 \cdot \frac{(x-1)(x-3)}{3} \end{aligned} \]

容易验证，对于以上表达式：$f(1) = 2$，$f(3) = 12$，$f(4) = 23$

Python 代码

直接使用拉格朗日基函数的定义，来求插值函数

import numpy as np

def Lagrange(X, Y):
    def f(x):
        r = 0
        for i in range(len(Y)):
            w = Y[i]
            for j in range(len(X)):
                if i == j:
                    continue
                w *= (x - X[j]) / (X[i] - X[j])
            r += w
        return r
    return f # 返回插值函数

测试：给定一个 $5$ 次多项式，以及 $6$ 个插值点，可以得到拉格朗日插值函数。再给定一个新的点，测试通过原始的多项式和插值函数计算的结果是否一样。

def test_lagrange():
    # f(x) 是一个 5 次多项式，表达式如下：
    # f(x) = 1 + 5 * x + 2 * x^2 + 4 * x^3 + 6 * x^4 + 3 * x^5
    a = np.array([1, 5, 2, 4, 6, 3])
    y = Polynomial(a, 3.5)  # 直接使用多项式计算 f(3.5)
    print(y)

    # 给定 6 个插值点
    X = np.array([1, 2, 5, 7, 9, 10])
    Y = np.zeros_like(X)
    for i in range(len(X)):
        Y[i] = Polynomial(a, X[i])

    f = Lagrange(X, Y) # 返回拉格朗日插值函数
    y = f(3.5) # 使用插值函数计算
    print(y)

test_lagrange()

拉格朗日插值虽然形式简单，便于理解，但是不具备增量特特性，即在已经计算出 $n$ 个点的插值函数后，如果再增加一个点，所有基函数都要重建。

如果不对插值多项式做化简，每次计算一个新的点时，时间复杂度为 $O(n^2)$。

5. 牛顿插值

牛顿插值法可以解决拉格朗日不具备增量计算的问题，解决方法是换了一组基函数。对于 $n+1$ 个不同的点 $(x_0, y_0)$，$(x_1, y_1)$...，$(x_n, y_n)$ ，牛顿法也是构造 $n+1$ 个基函数。对于上述每个点 $(x_i, y_i)$ ，基函数需要满足如下特性：

\[ \begin{aligned} \sum\limits_{j=0}^{i}{a_j\phi_j(x_i)} &= y_i \\ \sum\limits_{j=i+1}^{n}{a_j\phi_j(x_i)} &= 0 \end{aligned} \]

如果找到了这样一组基函数，那么最终插值多项式的表达式：$\displaystyle f(x) = \sum\limits_{i=0}^{n}a_i\phi_i(x)$，根据基函数特性可得：$f(x_i) = y_i$。

如何寻找基函数呢，仍然是构造法。牛顿构造的第 $i$ 个基函数表达式如下：

\[ \phi_i(x) = \left\{ \begin{aligned} &1 , i = 0 \\ &\prod\limits_{j=0}^{i-1}(x-x_j) , 1 \leq i \leq n \end{aligned} \right. \]

牛顿插值多项式中的 $a_i$ 是待定系数，需要满足牛顿基函数特性，并不能通过观测直接得到。为了求解 $a_i$，首先将牛顿插值多项式展开： \[ \begin{aligned} f(x) &= \phi_0(x) + \phi_1(x) + \phi_2(x) + ... + \phi_n(x) \\ &= a_0 + a_1(x-x_0) + a_2(x-x_0)(x-x_1) + ... + a_n(x-x_0)(x-x_1)...(x-x_{n-1}) \end{aligned} \]

将 $n+1$ 个不同的点的坐标 $(x_0, y_0)$，$(x_1, y_1)$...，$(x_n, y_n)$ 带入牛顿插值表达式 $f(x)$ 可得到以下方程组：

\[ \begin{array}{llll} f(x_0) &= a_0 &= y_0 \\ f(x_1) &= a_0 + a_1(x_1-x_0) &= y_1 \\ f(x_2) &= a_0 + a_1(x_2-x_0) + a_2(x_2-x_0)(x_2-x_1) &= y_2 \\ ... \\ f(x_n) &= a_0 + a_1(x_n-x_0) + a_2(x_n-x_0)(x_n-x_1) + ... + a_n(x_n-x_0)...(x_n-x_{n-1}) &= y_n \end{array} \]

上述方程组有 $n+1$ 个方程，$n+1$ 个未知量，由于这 $n+1$ 个点的坐标不同，而且各个方程之间线性无关，该方程组有唯一解。虽然可以使用高斯消元法来求解，但是使用带入法更加快捷。可以先通过第一个方程求解出 $a_0$，然后带入到第二个方程求解 $a_1$，等等以此类推直到求解出 $a_n$。

举例

和拉格朗日插值一样，假设有 $3$ 个点 $(1,2)$，$(3,12)$，$(4,23)$，求一个经过这 $3$ 个点的二次函数。

牛顿插值的二次函数的形式为：$f(x) = a_0 + a_1(x-1) + a_2(x-1)(x-3)$。带入坐标之后为：

$f(1) = a_0 = 2$

$f(3) = a_0 + a_1(3-1) = 2 + a_1(3-1) = 12 => a_1 = 5$

$f(4) = a_0 + a_1(4-1) + a_2(4-1)(4-3) = 2 + 5 \times 3 + a_2 \times 3 = 23 => a_2 = 2$

这样就把全部系数都求出来了，所以插值函数：$f(x) = 2 + 5(x-1) + 2(x-1)(x-3)$

编程实现

考察一下关于 $a_i$ 的方程： \[ \begin{aligned} f(x_i) &= a_0 \\ &+ a_1(x_i-x_0) \\ &+ a_2(x_i-x_0)(x_i-x_1) \\ &+ ... \\ &+ a_{i-1}(x_i-x_0)(x_i-x_1)...(x_i-x_{i-2}) \\ &+ a_i(x_i-x_0)(x_i-x_1)...(x_i-x_{i-1}) \end{aligned} \]

变形可得： \[ \begin{aligned} a_i &= \frac{f(x_i)-a_0}{(x_i-x_0)(x_i-x_1)...(x_i-x_{i-1})} \\ &- \frac{a_1}{(x_i-x_1)...(x_i-x_{i-1})} \\ &- ... \\ &- \frac{a_2}{(x_i-x_2)...(x_i-x_{i-1})} \\ &- \frac{a_{i-1}}{x_i-x_{i-1}} \end{aligned} \]

仿照霍纳法则，再对上述公式变形：

$a_i = \displaystyle(((\frac{f(x_i)-a_0}{x_i-x_0} - a_1) \frac{1}{x_i-x_1} - a_2) \frac{1}{x_i-x_2} - ... - a_{i-1}) \frac{1}{x_i-x_{i-1}}$

这样就把系数计算出来了。对于一个未知的点 $x$，如何计算 $f(x)$ 呢？仍然是仿霍纳法则

\[ \begin{aligned} f(x) &= a_0 + a_1(x-x_0) + a_2(x-x_0)(x-x_1) + ... + a_n(x-x_0)(x-x_1)...(x-x_{n-1}) \\ & = ((a_n(x-x_{n-1}) + a_{n-1})(x-x_{n-2}) + ... + a_1)(x-x_0) + a_0 \end{aligned} \]

Python 代码

import numpy as np

def Newton(X, Y):
    a = np.zeros_like(X)
    a = Y.copy() # 初始化 a[i] = Y[i]
    # 计算系数 a[i]
    for i in range(1, len(a)):
        for j in range(0, i):
            a[i] = (a[i] - a[j]) / (X[i] - X[j])

    # 已知系数 a[i]， 返回插值多项式 f(x)
    def f(x):
        r = a[len(a) - 1]
        for i in range(len(a) - 2, -1, -1):
            r = r * (x - X[i]) + a[i]
        return r
    return f

测试：给定一个 $5$ 次多项式，以及 $6$ 个插值点，可以得到牛顿插值函数。再给定一个新的点，测试通过原始的多项式和插值函数计算的结果是否一样。

def test_newton():
    # f(x) 是一个 5 次多项式，表达式如下：
    # f(x) = 1 + 5 * x + 2 * x^2 + 4 * x^3 + 6 * x^4 + 3 * x^5
    a = np.array([1, 5, 2, 4, 6, 3])
    y = Polynomial(a, 3.5)  # 直接使用多项式计算 f(3.5)
    print(y)

    # 给定 6 个插值点
    X = np.array([1, 2, 5, 7, 9, 10])
    Y = np.zeros_like(X)
    for i in range(len(X)):
        Y[i] = Polynomial(a, X[i])

    f = Newton(X, Y) # 返回牛顿插值函数
    y = f(3.5) # 使用插值函数计算
    print(y)

test_newton()

牛顿插值法有一个初始化阶段，时间复杂度为 $O(n^2)$，初始化完成之后，在计算新的点对应的值时，时间复杂度为 $O(n)$。而且牛顿基函数具备增量计算特性：如果已知前 $n$ 个基函数，再增加 1 个插值点时，前 $n$ 个基函数都不变，只需增加一个新的基函数即可。

差商

虽然通过以上算法可以计算出牛顿插值基的系数，但是对于系数的物理意义并没有解释，实际上第 $i$ 个插值基函数的系数，是第 $i$ 阶差商。差商的定义是递归的，其定义如下：

零阶差商：$F[x_i] = f(x_i)$

一阶差商：$\displaystyle F[x_i,x_j] = \frac{F[x_i] - F[x_j]}{x_i-x_j}$

二阶差商：$\displaystyle F[x_i,x_k,x_j] = \frac{F[x_i,x_k] - F[x_k,x_j]}{x_i-x_j}$

...

$k$ 阶差商：$\displaystyle F[x_i,x_{i+1},...,x_{i+k}] = \frac{F[x_i,x_{i+1}, ..., x_{k-1}] - F[x_{i+1},x_{i+2}, ..., x_{i+k}]}{x_i-x_{i+k}}$

为了得到插值函数 $f(x)$ 和差商的关系，可以按阶从低到高考察差商函数。

考察一阶差商函数：$\displaystyle F[x, x_0] = \frac{F[x] - F[x_0]}{x-x_0}$，变形可得： \[ F[x] = F[x_0] + (F[x] - F[x_0])(x-x_0) \]

由于零阶差商： $F[x] = f(x)$，替换之后可得： \[ f(x) = F[x_0] + (F[x] - F[x_0])(x-x_0) \]

再考察二阶差商函数：

\[ \begin{aligned} F[x, x_0, x_1] &= \frac{F[x,x_0] - F[x_0,x_1]}{x-x_1} &= \frac{\displaystyle\frac{F[x]-F[x_0]}{x-x_0} - F[x_0,x_1]}{x - x_1} \end{aligned} \]

等式左右两边同时乘以 $(x-x_0)(x-x_1)$ 之后可得：

\[ F[x, x_0, x_1](x-x_0)(x-x_1) = F[x]-F[x_0] - F[x_0,x_1](x-x_0) \] 移项并用 $f(x)$ 替换 $F[x]$ 得：

\[ f(x) = F[x_0] + F[x_0,x_1](x-x_0) + F[x, x_0, x_1](x-x_0)(x-x_1) \]

...

一步一步替换之后，最终可得： \[ \begin{aligned} f(x) &= F[x_0] \\ &+ F[x_0,x_1](x-x_0) \\ &+ F[x_0, x_1, x_2](x-x_0)(x-x_1) \\ &+ ... \\ &+ F[x_0, x_1, ..., x_{n}](x-x_0)(x-x_1)...(x-x_{n-1}) \\ &+ F[x, x_1, ..., x_{n-1}, x_{n}](x-x_1)...(x-x_{n-1})(x-x_n) \end{aligned} \] 以上是推导过程，更严格的证明可以使用数学归纳法。

把 $x_0$ 带入以上表达式可得

\[f(x_0) = F[x_0] + F[x, x_1, ..., x_{n-1}, x_{n}](x_0-x_1)...(x_0-x_{n-1})(x_0-x_n)\]

由于零阶差商：$F[x_0] = f(x_0)$，所以 $F[x, x_1, ..., x_{n-1}, x_{n}] = 0$。最终 $f(x)$ 的表达式如下：

\[ \begin{aligned} f(x) &= F[x_0] \\ &+ F[x_0,x_1](x-x_0) \\ &+ F[x_0, x_1, x_2](x-x_0)(x-x_1) \\ &+ ... \\ &+ F[x_0, x_1, ..., x_{n}](x-x_0)(x-x_1)...(x-x_{n-1}) \end{aligned} \]

通过以上推导可以发现，牛顿插值法的系数 $a_k$，其实就是 $k$ 阶差商 $F[x_0, x_1, ..., x_k]$。

Python 代码

首先计算 $0$ 阶差商，然后通过递推公式，依次计算 $1$ 阶，$2$ 阶，$3$ 阶， ...， $n$ 阶差商。差商计算出来之后，仍然使用霍纳法则计算多项式的值。

def NewtonDivDiff(X, Y):
    d = np.zeros([len(Y), len(Y)])
    # 计算 0 阶差商，这里表示为 d[i][i]
    for i in range(0, len(Y)):
        d[i][i] = Y[i]

    # k 为当前差商的阶数，先计算 1 阶，然后 2 阶，依次类推
    for k in range(1, len(Y)):
        for i in range(0, len(Y) - k):
            d[i][i+k] = (d[i][i+k-1] - d[i+1][i+k]) / (X[i] - X[i+k])

    # 已知系数 a[i]， 计算多项式 f(x)
    def f(x):
        r = d[0][len(Y) - 1]
        for i in range(len(Y) - 2, -1, -1):
            r = r * (x - X[i]) + d[0][i]
        return r
    return f

可以看到，直接使用差商的定义和前面求解的方式代码是相似的，二者的时间复杂度是一样的，都是 $O(n)$。直接求解时，代码更加简洁，空间复杂度为 $O(n)$，而使用差商定义时，需要用一个二维的数组，空间复杂度为 $O(n^2)$，要高一些。但是使用差商的定义，物理意义更加明确一些。

总结

本文介绍了一般的多项式插值、拉格朗日插值以及牛顿插值，不仅介绍了原理，还给出了 Python 代码的实现，可以供读者研究使用。虽然通过定理 1 可以知道，多项式插值最终的表达式是一样的，但是使用不同的插值方法其算法复杂程度是不一样的。一般多项式插值的算法复杂度最高，拉格朗日插值算法复杂度稍低，牛顿插值算法复杂度最低。$$