重心拉格朗日插值

上一节中，我们介绍了 Lagrange 插值多项式，证明了它的唯一性。也给出了在 $n+1$ 个采样点下，对于连续函数 $f(x)$ 的逼近误差。最后用一个示例，介绍了等距采样时的龙格现象。很多人都说 Lagrange 插值公式在数学上很优美，但是工程上很难用。主要在于它的时间复杂度很高，而且存在严重的数值不稳定问题。

实际上 Lagrange 多项式性能一点也不差，只是没有用对。2004 年，Berrut 和 Trefethen 专门发表了一篇名为 Barycentric Lagrange Interpolation 的文章重新发现了重心拉格朗日插值公式。重心拉格朗日只是在数学形式上，对传统的拉格朗日多项式做了等价变形，就获得了极佳的数值稳定性。并且计算效率也得到了提升，经过一次 $O(n^2)$ 的初始化，之后求值都只需要 $O(n)$ 的时间复杂度。这篇文章还特意讨论了龙格现象，指出这完全是因为等距采样导致的，如果换成 Chebyshev 采样就不会有问题。

本节我们来研究一下这个重心 Lagrange 插值公式及其工程实现。关于解决龙格现象的 Chebyshev 采样我们后续再专门写文章介绍。

1. 重心 Lagrange 插值公式

给定 $n + 1$ 个采样点 $\{x_0, \cdots, x_n\}$ 及其采样值 $\{y_0, \cdots, y_n\}$，要求采样点互不相同，即 $x_i \neq x_j, i \neq j$。那么我们可以写出如下 Lagrange 多项式，经过所有的采样点 $(x_i, y_i)$：

$$ \begin{equation}\label{f1} P(x) = \sum_{i = 0}^n L_i(x) = \sum_{i=0}^n y_i \prod_{i \neq k}\frac{x - x_i}{x_k - x_i} \end{equation} $$

我们对上式作如下变形

$$ \begin{equation}\label{f2} P(x) = \sum_{i = 0}^n L_i(x) = \sum_{i=0}^n \left[ \underbrace{(x - x_0)\cdots(x - x_{i-1})(x - x_{i+1}) \cdots (x - x_n)}_{\ell_i(x)} \cdot \frac{y_i}{\prod_{i \neq k} (x_k - x_i)} \right] \end{equation} $$

上式中的 $\ell_i(x)$ 中间缺少了 $(x - x_i)$ 这一项，那么我们完全可以将其写成如下的形式:

$$ \begin{equation}\label{f3} \ell_i(x) = \frac{\prod_{j=0}^n (x - x_j)}{(x - x_i)} \end{equation} $$

我们把上式中 $\prod_{j=0}^n (x - x_j)$ 记为 $\ell(x)$。那么式$(\ref{f1})$中累积求和的每一项都有一个 $\ell(x)$。可以改写成如下形式:

$$ \begin{equation}\label{f4} P(x) = \sum_{i = 0}^n L_i(x) = \ell(x)\sum_{i=0}^n \left[ \frac{1}{x - x_i} \cdot \frac{y_i}{\prod_{i \neq k} (x_k - x_i)} \right] = \ell(x)\sum_{i=0}^n \left[ \frac{1}{\prod_{i \neq k} (x_k - x_i)} \cdot \frac{y_i}{x - x_i} \right] \end{equation} $$

对于确定的采样点 $\{x_0, \cdots, x_n\}$ 而言，$\frac{1}{\prod_{i \neq k} (x_k - x_i)}$ 就是一个常数，可以将之看做是每一项的权重，用符号 $w_i$ 来表示。

$$ \begin{equation}\label{f5} P(x) = \sum_{i = 0}^n L_i(x) = \ell(x)\sum_{i=0}^n \left[ \frac{w_i}{x - x_i}\cdot y_i \right] \end{equation} $$

上式被称为第一型重心插值公式(first form of the barycentric interpolation formula)。这个式子里面权重 $w_i$ 只跟采样点$\{x_0, \cdots, x_n\}$的分布有关，跟目标函数 $f(x)$ 没有任何关系。如果对于函数 $f(x) = 1$ 进行采样，一定有 $y_i = 1$，套用上式我们可以直接写出其重心插值公式 $P_1(x)$：

$$ P_1(x) = \ell(x)\sum_{i=0}^n \left[ \frac{w_i}{x - x_i}\cdot 1 \right] $$

对常数函数插值得到的一定还是该常数，所以 $P_1(x) = 1$。我们将式$(\ref{f5})$除以$P_1(x)$就可以直接约去它们的公共因子 $\ell(x)$，如此就得到了第二型重心插值公式(second(true) form of the barycentric formula):

$$ \begin{equation}\label{f6} P(x) = \frac{\sum_{i=0}^n \left[\frac{w_i}{x - x_i}\cdot y_i\right]}{\sum_{i=0}^n \left[\frac{w_i}{x - x_i}\right]} \end{equation} $$

根据上式，我们只需要在初始化的时候，把权重系数计算出来，然后就可以在一层循环里面，计算分子和分母的累加式。时间复杂度就是 $O(n)$ 的。

式$(\ref{f6})$的这种分式的结构具有很强的数值稳定性。因为它的分子和分母只在 $y_i$ 上有差异，计算过程中对 $\frac{w_i}{x - x_j}$ 产生的计算误差，会同时出现在分子分母上，最终的除法运算会很大程度上抵消掉的。这跟硬件上要用双绞线的差分电平来抵消共模干扰信号一个意思。此外，权重系数 $w_i$ 只跟采样点 $x_i$ 之间的相对位置有关，在当前的分式结构下，我们可以直接对所有的权重 $w_i$ 等比例缩放，而不影响最终结果。比如同除以绝对值最大的那个权重 $w_{max}$，这样可以防止某个权重特别大导致计算溢出。

参考公式$(\ref{f6})$，我们定义了类 BarycentricLagrange，其中有三个数组 mXs, mYs, mWs，分别用于保存采样点 $x_i$ 采样值 $y_i$ 和权重 $w_i$。我们在该类的构造函数中调用 InitWeights() 完成初始化，之后就可以直接调用 Evaluate() 计算插值。这两个成员函数的实现如下面代码片段所示:

        void InitWeights()
        {
            size_t n = mXs.size();
            mWs.resize(n);
            for (size_t i = 0; i < n; ++i) {
                Scalar w = 1;
                for (int k = 0; k < n; ++k) {
                    if (k == i)
                        continue;
                    w *= (mXs[i] - mXs[k]);
                }
                mWs[i] = 1 / w;
            }
            NormalizeWeights();
        }

        Scalar Evaluate(Scalar x) const
        {
            Scalar numerator = 0;
            Scalar denominator = 0;

            for (int i = 0; i < mXs.size(); ++i) {
                if (std::abs(x - mXs[i]) < SMALL_VALUE)
                    return mYs[i];

                Scalar term = mWs[i] / (x - mXs[i]);
                numerator += term * mYs[i];
                denominator += term;
            }
            return numerator / denominator;
        }

2. 新增采样点

        void AddPoint(Scalar x, Scalar y) {
            size_t n = mXs.size();
            // 1. 更新现有的权重
            for (int i = 0; i < n; ++i) {
                assert(std::abs(mXs[i] - x) > SMALL_VALUE);
                mWs[i] /= (mXs[i] - x);
            }
            // 2. 计算新采样点权重
            Scalar w = 1;
            for (int i = 0; i < n; ++i)
                w /= (x - mXs[i]);
            mXs.push_back(x); mYs.push_back(y);
            mWs.push_back(w);
            NormalizeWeights();
        }

式$(\ref{f1})$还有一点经常被人诟病：每次新增采样点的时候，都需要重新计算一次多项式系数。关于这个缺陷，我觉得还是因为传统的 Lagrange 多项式 $O(n^2)$ 的插值计算复杂度。人们觉得它太慢了，所以希望展开得到多项式的系数，这样后续再给定 $x$ 坐标插值时，就可以通过 $O(n)$ 的嵌套乘法完成了。但是式$(\ref{f1})$中采样点之间又不是互相独立的，这导致只要采样点发生了改变，就需要重新展开，增加了很多重复计算，而且过程很繁琐。所以人们更喜欢用可以增量更新的牛顿差商公式。

现在有了重心 Lagrange 多项式，我们也可以在 $O(n)$ 的时间复杂度内完成权重的更新。如右侧函数 AddPoint 所示，大体上就是两步:

更新现有的权重，给每个权重除以 $x_i - x_{new}$
为新增的点，计算权重 $1 / \prod_{i=0}^n(x_i - x_{new})$

之后我们把新增的采样点、采样值、权重记录到成员变量 mXs, mYs, mWs，并通过函数 NormalizeWeights 按照绝对值最大的权重对 mWs 缩放。

3. 完

本文对传统的拉格朗日多项式简单的做了恒等变形之后，得到了重心 Lagrange 插值公式。该公式不仅在数学上很优雅，在工程上它也是十分高效，并且极致稳定。