区块链中的数学（八十四）-- 重心插值（barycentric evaluation）

写在前面

很久前，写了篇文章解释拉格朗日插值，参见拉格朗日插值法  本文介绍一种优化的插值法：重心拉格朗日插值法。

拉格朗日插值

1.朴素的拉格朗日插值

这里再简单介绍下拉格朗日插值法： 假设有$k+1$个数据点$(x_0, y_0), (x_1, y_1), \dots, (x_k, y_k)$，其中所有$x_j$互异。目标是构造次数不超过$k$的多项式$P(x)$满足$P(x_j) = y_j$。

拉格朗日基多项式$l_j(x)$定义为： $lj(x) := \prod{i=0}^{k} \frac{x - x_i}{x_j - x_i}$ 满足： $\bullet \,\, l_j(x_j) = 1$ $\bullet \,\, l_j(xi) = 0 \quad (i \neq j)$ 插值多项式为基多项式的加权和： $P(x) = \sum{j=0}^k y_j \cdot l_j(x)$ 计算效率：

1. 每个$l_j(x)$需$O(k)$次乘法
2. 整体$P(x)$需$O(k^2)$次运算（包含分子分母运算）

3. 增删数据点时需完全重新计算

2. 重心插值法第一型

定义零化多项式（vanishing  polynomial）：

其导数在$xj$处的值（注：如果你对此处导数不理解，但最终能明白推导过程中的等式替换也OK，导数出处来自于第三个引用资料）： $l(x) := \prod{i=0}^{k} (x - x_i)$ 基多项式可表示为： $l'(xj) = \prod{\begin{subarray}{c} i=0 \ i\ne j \end{subarray}}^{k}(x_j - x_i)$

重心权重

定义权重常数： $w_j := \frac{1}{l'(xj)} = \left[ \prod{\substack{i=0 \ i \neq j}}^{k} (x_j - xi) \right]^{-1}$ 插值公式转化为： $P(x) = l(x) \sum{j=0}^{k} \frac{w_j y_j}{x - x_j}$ 优势：预计算权重$w_j$（$O(k^2)$时间）后，单点求值仅需$O(k)$ 时间， 关键是只需计算一次，后面可以缓存起来使用。

3. 重心插值法第二型（标准形式）

当$yj=1$时（常数函数， 注：这里借助了常数函数特例，辅助后面公式推导），有恒等式： $1 = l(x) \sum{j=0}^{k} \frac{w_j}{x - xj}$ 将第一型公式P(x)除以上式得标准形式： $P(x) = \frac{\sum{j=0}^k \frac{w_j y_j}{x - xj}}{\sum{j=0}^k \frac{w_j}{x - x_j}}$ 优势：

1. 稳定性：避免计算$l(x）$的可能带来大幅值波动
2. 高效性：单点求值保持$O(k)$复杂度

3. 可扩展性：增删点时权重局部更新

4. 计算示例

给定数据点： $\bullet \,\, (x_0, y_0) = (1, 3)$ $\bullet \,\, (x_1, y_1) = (2, 8)$ $\bullet \,\, (x_2, y_2) = (4, 6)$ 求二次多项式$P(x)$在$x = 3$处的值。

步骤1：计算重心权重

权重公式： $wj = \prod{\substack{i=0 \ i \neq j}}^{k} (x_j - x_i)^{-1}$

$\bullet \,\, w_0 = \frac{1}{(1-2)(1-4)} = \frac{1}{(-1)(-3)} = \frac{1}{3}$ $\bullet \,\, w_1 = \frac{1}{(2-1)(2-4)} = \frac{1}{(1)(-2)} = -\frac{1}{2}$ $\bullet \,\, w_2 = \frac{1}{(4-1)(4-2)} = \frac{1}{(3)(2)} = \frac{1}{6}$ 得权重向量： $w = \left[ \frac{1}{3}, \, -\frac{1}{2}, \, \frac{1}{6} \right]$

步骤2：计算 

标准形式： $P(3) = \frac{\sum_{j=0}^2 \frac{w_j y_j}{3 - xj}}{\sum{j=0}^2 \frac{w_j}{3 - x_j}}$ 分子计算： $分子 = \frac{w_0 y_0}{3 - 1} + \frac{w_1 y_1}{3 - 2} + \frac{w_2 y_2}{3 - 4}$ $=\frac{(1/3) \cdot 3}{2} + \frac{(-1/2) \cdot 8}{1} + \frac{(1/6) \cdot 6}{-1}$ $=\frac{1}{2} - 4 - 1 = -4.5$

分母计算：

$分母=\frac{w_0}{3-1} + \frac{w_1}{3-2} + \frac{w_2}{3-4}$ $=\frac{1/3}{2} + \frac{-1/2}{1} + \frac{1/6}{-1}$ $=\frac{1}{6} - \frac{1}{2} - \frac{1}{6} = -0.5$ 最终结果： $P(3) = \frac{-4.5}{-0.5} = 9$

验证

通过方程组解得多项式： $P(x) = -2x^2 + 11x - 6$ 代入$x = 3$：

结果一致。 $P(3) = -2(9) + 11(3) - 6 = -18 +33 -6 =9$

有限域中的插值

基础设定

• 有限域：$\mathbb{F}$（特征为素数$p)$
• 单位根：$\omega \in \mathbb{F} 满足 \omega^n = 1其中 (n = 2^k）$
• 插值点集：$x_j = \omega^j (j = 0, 1, \ldots, n - 1）$
• 零化多项式：$l(x) = \prod_{j=0}^{n-1} (x - \omega^j) = x^n - 1$

1. 导数计算

$l'(x) = \frac{d}{dx}(x^n - 1) = nx^{n-1}$

2. 在插值点处的导数值

$l'(x_j) = l'(\omega^j) = n(\omega^j)^{n-1} = n\omega^{j(n-1)}$

3. 单位根性质利用

由$\omega^n = 1$可得： $\omega^{j(n-1)} = \omega^{jn - j} = (\omega^n)^j \cdot \omega^{-j} = 1^j \cdot \omega^{-j} = \omega^{-j}$

4. 权重闭式解

$w_j = \frac{1}{l'(x_j)} = \frac{1}{n\omega^{-j}} = \frac{\omega^j}{n}$

重心插值公式简化

1. 标准第二型公式

$P(z) = \frac{\sum_{j=0}^{n-1} \frac{w_j yj}{z - \omega^j}}{\sum{j=0}^{n-1} \frac{w_j}{z - \omega^j}}$

2. 代入权重闭式解

$P(z) = \frac{\sum_{j=0}^{n-1} \frac{\left( \frac{\omega^j}{n} \right) yj}{z - \omega^j}}{\sum{j=0}^{n-1} \frac{\omega^j / n}{z - \omega^j}} = \frac{\sum_{j=0}^{n-1} \frac{\omega^j yj}{z - \omega^j}}{\sum{j=0}^{n-1} \frac{\omega^j}{z - \omega^j}}$

3. 分母的闭式简化

考虑分母求和： $D(z) = \sum{j=0}^{n-1} \frac{\omega^j}{z - \omega^j}$ 通过代数变换： $D(z) = \sum{j=0}^{n-1} \frac{\omega^j}{z - \omega^j}$ $= \sum{j=0}^{n-1} \left( -1 + \frac{z}{z - \omega^j} \right)$ $= -n + z \sum{j=0}^{n-1} \frac{1}{z - \omega^j}$

利用单位根性质可得： $\sum_{j=0}^{n-1} \frac{1}{z - \omega^j} = \frac{n z^{n-1}}{z^n - 1}$ 因此： $D(z) = -n + z \cdot \frac{n z^{n-1}}{z^n - 1} = -n + \frac{n z^n}{z^n - 1} = \frac{n}{z^n - 1}$

4. 最终简化公式

$P(z) = \frac{z^n - 1}{n} \sum_{j=0}^{n-1} \frac{\omega^j y_j}{z - \omega^j}$

边界情况

1. 当$z = \omega^k$时

$\lim_{z \to \omega^k} P(z) = y_k$

2. 当$z^n = 1$但$z \ne \omega^k$ 

此时$z$是$d$次单位根（$d|n$），需要：

1. 计算$z$ 的阶$m$
2. 生成$m$阶循环子群$ ${z^k \mid k = 0, 1, \dots, m - 1}$
3. 在子群上重新插值并求值

小结

结合重心权重的拉格朗日插值法往往是实际工程应用中更常用的算法，文中说的闭式是指closed form 有限域中的乘法子群。碰巧看到V神最近也在学习这个东西（见参考二）， what a coincidence!

References

https://github.com/Plonky3/Plonky3/blob/main/interpolation/src/lib.rs

https://hackmd.io/@vbuterin/barycentric_evaluation

https://people.maths.ox.ac.uk/trefethen/barycentric.pdf