Skip to content

Lagrangian Interpolation

Definition

定义： $$ L_{n,i} (x) = \prod_{j=0, j \neq i}^n \frac {x-x_j} {x_i-x_j} $$ 也就是说：\(L_{n,i}\) 就是一个 \(n\) 阶多项式。它在 \(x_0, x_1, \dots, x_{i-1}, x_{i+1}, \dots, x_n\) 处为 0，在 \(x_i\) 处为 1。

定义： $$ P_n(x) = \sum_{i=0}^n f(x_i) L_{n,i}(x) $$ 也就是说：\(P_n\) 就是一个 \(n\) 阶多项式。它在 \(x_0, x_1, \dots x_{n-1}, x_n\) 处分别为 \(f(x_0), f(x_1), \dots, f(x_{n-1}), f(x_n)\)。

Properties

  1. 如果 n+1 个点的横坐标不相同，那么，有且只有一个 n 阶多项式，能够经过这 n+1 个点
    • 假如有两个 n 阶多项式都经过这 n+1 个点，那么，两式之差必然以这 n+1 个点为零点。因此，两式之差的 n+1 个点要么是 1，

Remainder Analysis

前提：

  1. 我们只讨论 \([x_0, x_n]\) 里面的误差
  2. \(f(x)\) 足够光滑

定义： $$ R_n(x) = f(x) - P_n(x) $$

分析：

由于 \(\forall i \in \set{0, \dots, n}: R_n(x_i) = 0\)

从而：\(R_n(x) = K(x) \prod_{i=0}^n (x - x_i)\)

任取 \(x' \in [x_0, x_n]\) 且 \(\forall i \in \set{0, \dots, n}: x' \neq x_i\)

令 \(g_n(t) = R_n(t) - K(x') \prod_{i=0}^n (t-x_i) = (K(t) - K(x')) \prod_{i=0}^n (t-x_i)\)

从而：\(g_n(t)\) 至少有 \(x_0, \dots, x_n, x'\) 共 n+2 个零点。

从而：\(\exists \xi_x \in (x_0, x_n): g_n^{(n+1)}(\xi_x) = 0\)

从而：

\[ \begin{aligned} 0 &= R_n^{(n+1)}(\xi_x) - (n+1)!K(x') \newline &= f^{(n+1)}(\xi_x) - P^{(n+1)}(\xi_x) - (n+1)!K(x') \newline &= f^{(n+1)}(\xi_x) - (n+1)!K(x') \end{aligned} \]

从而：

\[ \begin{aligned} R_n(x') &= K(x')\prod_{i=0}^n (x' - x_i) \newline &= \frac 1 {(n+1)!} f^{(n+1)}(\xi_x) \prod_{i=0}^n (x' - x_i) \newline \end{aligned} \]

如果和泰勒展开的拉格朗日余项进行对比，会发现两者的形式其实很像，除了泰勒展开后面不是 \(\prod_{i=0}^n (x' - x_i)\)，而是 \((x'-a)^{n+1}\)。

Example: sin

也就是说：如果采用均匀插值的方式，那么后面的多项式会随着插入点的增加，以 \(\mathcal O(\frac 1 {n^n})\) 的速度收敛。我们只需要保证前面的 \(f^{(n)}(x)\) 的导数不爆炸就行。

Neville Interpolation

Lift the Degree of a Polynomial

我们可以发现，令多项式 \(p_{\text{some indices}}\) 经过的点的 index 即为其下标，且阶数为 \(\#\text{ of indices} - 1\)。那么，如果已知两个 \(n\) 阶多项式 \(p_{1,2,\dots,i-1,i+1,\dots,n}, p_{1,2,\dots,j-1,j+1,\dots,n}\)，我们就可以凑出：

\[ p_{1, 2, \dots, n} = \frac {(x-x_i) p_{1,2,\dots,i-1,i+1,\dots,n} - (x-x_j) p_{1,2,\dots,j-1,j+1,\dots,n}} {x_j - x_i} \]

Interpolation

对于高阶的插值，我们可以用 neville 插值：

img

如图，上图中的 \(p_{(m, n)} := p_{m, m+1, \dots, n-1, n}\)。因此，我们有递推关系： $$ p_{m, n} = \frac{(x-x_m) p_{m+1, n} - (x - x_n)p_{m,n-1}} {x_n - x_m} $$ 因此，一开始，我们只需要把最下面的 \(p_{m}(x) := f(x_m)\) 算出来（这里没有显示），就可以层层递推上去了。

如果 \(p_{m,n}\) 是值

计算复杂度是 \(\mathcal O(\sum_{i=1}^n (n-i+1)) = O(N^2)\)，和 naïve Lagrangian interpolation 一样，其中 \(N\) 是总共的点数。

好处是

  1. 增加插值点（如上图蓝色部分所示）的时候，我们可以利用之前的结果，从而增加一个点的计算复杂度只有 \(\mathcal O(\sum_{i=1}^n 1) = \mathcal O(N)\)。
  2. 可以递归地进行计算，代码简洁。

如果 \(p_{m,n}\) 是多项式

计算复杂度是 \(\mathcal O(\sum_{i=1}^n (n-i+1)i) = \mathcal O(N^3)\)，比 naïve Lagrangian interpolation 要大。增加插值点的时候，计算复杂度也是 \(\mathcal O(\sum_{i=1}^n i) = \mathcal O(N^2)\)，并不比 naïve 的更好。

因此，如果我们不仅要多项式在某个点 \(x^\ast\) 处的值，更需要多项式本身，那么，还需要下面的牛顿插值法。

Newton Interpolation

image-20240410164922610

如图，给定 \(n+1\) 个点和 \(n+1\)​ 个系数，我们就可以构建出 newton interpolation polynomial。

  • 不过，在构造 polynomial 的时候，我们用不到第 \(n+1\) 个点，i.e. \(x_{n }\)

然后，我们将这 \(n+1\) 个点带入 polynomial，形成 \(n+1\) 个值：\(f_0, \dots, f_n\)​。值和系数的关系如上图所示。

从而，我们只需要通过 \(f_i := f(x_i)\) 反解出 \(\vec \alpha\)​ 即可。

Properties

  1. 由于 \(\forall i \geq j: N^{(i)}(x_j) = 0\)，因此矩阵是下三角矩阵 ，本质上是可以直接用来求解线性方程组的。
  2. 由插值多项式的唯一性：Newton 插值、Lagrange 插值乃至所有的多项式插值，在最终结果上是等价的。

Computation

我们通过差商来计算，差商的定义如下：

牛顿多项式计算结果如下：

\[ N(x)=[y_{0}]+[y_{0},y_{1}](x-x_{0})+\cdots +[y_{0},\ldots ,y_{k}](x-x_{0})(x-x_{1})\cdots (x-x_{k-1}) \]

Example

image-20240410173335910

  • e.g. \([x_{-1}, x_0, x_1, x_2] = ([x_{-1}, x_0, x_1] - [x_0, x_1, x_2]) / (-1 - 2) = (3 - 2.5) / -3 = - 1/6\)

可以发现，

  • 差商表的递推顺序，和 Neville 的递推顺序完全一样。差别就在于递推函数
  • 由于差商即系数，因此建表只需要 \(\mathcal O(N^2)\)，增加一个点只需要 \(\mathcal O(N)\)。从而，Newton 多项式的复杂度，和 Neville 的值的复杂度一样。（具体对比见下）

Comparison of Neville and Newton

Newton vs Neville vs 待定系数: Matrix form

  • 待定系数法 vs Lagrangian interpolation vs Newton interpolation
    • 待定系数法的基底是 \(\set{1, x, x^2, \dots, x^n}\)，矩阵是稠密矩阵
    • Lagrangian basis 是 \(\set{L_{n,0}, L_{n,1},\dots, L_{n,n}}\)，矩阵是单位矩阵
    • Newton basis 是 \(\set{N^{(0)},N^{(1)},\dots,N^{(n)}}\)​，矩阵是下三角矩阵
  • 我们可以认为：我们通过将基底变换成 Lagrangian 和 Newton，使得这个矩阵变成了更容易求解的形式

差商表 vs Neville 表的递推复杂度

  • 差商表是 \([m\sim n] = \frac{[m+1 \sim n] - [m \sim n-1]}{n-m}\)​
    • \([\dots]\) 是系数（值），本质上也是多项式，但是多项式的各系数比例是固定的，因此可以用单一系数表示。
  • Neville 是 \(p_{m, n} = \frac{(x-x_m) p_{m, n-1} - (x - x_n)p_{m+1,n}} {x_n - x_m}\)
    • \(p_{\dots}\) 是值或者多项式。如果 \(p_{\dots}\) 是多项式，由于多项式的各系数比例是不固定的，因此必须用数组表示。从而造成了 \(\mathcal O(N^3)\) 的复杂度。

Hermite Interpolation

Hermite Interpolation 结合了 Naïve Interpolation 以及 Taylor Expansion 两个多项式近似方法。也就是：在每一个插值点上，我们不仅要考虑零阶导数（i.e. 函数值），也要考虑一阶及以上的导数。

Example

假设我们希望通过 \(f(x_0), f(x_1), f(x_2), f'(x_1)\) 进行插值。那么，4 个约束，就对应 3 阶多项式。 $$ P_3 = f(x_0) h_0(x) + f(x_1)h_1(x) + f(x_2) h_2(x) + f'(x_1) \widehat h_1(x) $$ 其中： $$ \begin{aligned} &h_i(x_j) = \delta_{ij}, h_1'(x_1) = 0 \\ &\widehat h_1(x_j) = 0, \widehat h_1'(x_1) = 1 \end{aligned} $$

  • 或者用下面的表格来形象说明（x 代表考虑这一项，i.e. 我希望函数在这里满足这样的条件）
\(x_0\) \(x_1\) \(x_2\)
\(f\) x (\(h_1\)) x (\(h_2\)) x (\(h_3\))
\(f'\) x (\(\widehat h_1\))
  • 每一个函数（如上表），在
    • 自己的 x 处为 1
    • 在其它的 x 处为 0
    • 在空白处不用管

对于 \(h_0(x)\)，由于 \(h_1'(x_1) = 0\)，因此就是 \(x=x_1\) 为双重根： $$ h_i(x) = C_0(x-x_1)^2(x-x_2) $$ 同样：\(h_2(x) = C_2(x-x_0)(x-x_1)^2\)。

对于 \(h_1(x)\)，由于 \(h_1(x_1) = 1, h_1'(x_1) = 0\)： $$ h_1(x) = (Ax+B)(x-x_0)(x-x_2) $$ 对于 \(\widehat h_1(x)\)，则显然是： $$ \widehat h_1(x) = C(x-x_0)(x-x_1)(x-x_2) $$ 最后，我们通过待定系数法解出其中的系数。注意 \(h_1(x)\) 有两个系数。

Error

误差也和 Lagrange 是类似的： $$ \begin{aligned} R_3(x) &= f(x) - P_3(x) \\ &= K(x)(x-x_0)(x-x_1)^2(x-x_2) \implies \\ K(x) &= \frac {f^{(4)}(\xi_x)}{4!} \end{aligned} $$

二重 Hermite 插值

常用的插值方法，就是使用所有插值点的函数值和一阶导数值。 $$ H_{2 n+1}(x)=\sum_{i=0}^{n} y_{i} A_{i}(x)+\sum_{i=0}^{n} y_{i}^{\prime} B_{i}(x)\ $$ 容易证明： $$ A_i(x) = (Ax+B)\prod_{j = 0 \land j \neq i}^n (x-x_j)^2 $$

  • 因为除了 \(x_i\) 以外，所有都是重根；\(x_i\)​ 处没有根，但是导数为 0。
    • 所以无法直接求出在 \(x_i\) 处的式子，必须用 \(A, B\) 两个待定量
\[ B_i(x) = C(x-x_i) \prod_{j = 0 \land j \neq i}^n (x-x_j)^2 \]
  • 因为除了 \(x_i\) 以外，所有都是重根；\(x_i\)​ 处为单根。
    • 所以可以直接求出在 \(x_i\) 处的式子，只用 \(C\) 一个待定量

Error

误差和 Lagragian 的形式类似。