重心插值(Barycentric interpolation )是拉格朗日插值的变换。
有时候需要通过一组多项式点值直接计算另一个不同点处的值,例如,$p(x)$ 是一个度为 2 的多项式,可以在 O(N) 时间内用 $p(0), p(1), p(2)$ 直接计算 $p(100)$ 处的值,而无须重构多项式。
回顾下拉格朗日插值
$$ P(x) = \sum^{n}_{i=1}y_iL_i(x) $$
其中 $L_i(x)$ 是拉格朗日基
$$ l_i(x)=\prod^{n-1}_{j=1,j \ne i}\frac{x-x_j}{x_i-x_j}=\frac{x-x_0}{x_i-x_0}...\frac{x-x_{i-1}}{x_i-x_{i-1}}\frac{x-x_{i+1}}{x_i-x_{i+1}}...\frac{x-x_{n-1}}{x_i-x_{n-1}} $$
分子在除了 $x_i$ 外的每个评估点都等于 0,分母在 $x_i$ 处等于分子。所以 $L_i(x)$ 在 $x_i$ 处等于 1,非 $x_i$ 处等于 0
令 $M(x)=\prod_{j=0}^{n-1}(x-x_j)=(x-x_0)(x-x_1)...(x-x_{n-1})$,$l_i(x)$ 的分母用 $d_i$(重心权重的倒数) 表示,得到
$$ L_i(x)=\frac{M(x)}{d_i(x-x_i)} $$
得到
$$ P(x) = M(x) * \sum_i\frac{y_i}{d_i(x-x_i)} $$
其中,$d_i$ 可以预计算,$M(x)$ 的时间复杂度为 O(N),所以 $L_i(x)$ 的时间复杂度为 O(N),所以 $P(x)$ 的时间复杂度为 O(N)
我们可以选择求值点以避免预计算 $d_i$,比如 $\left\{1, w, w^2, ...,w^{n-1}\right\}$,$n$ 是 2 的幂,即 $n=2^k$。有
$$ d_i=\prod_{0 \le j \ne i < n}(w^i-w^j) =(w^i-w^{i+\frac{n}{2}})\prod_{0 \le j \ne i < \frac{n}{2}}(w^i-w^j)(w^i-w^{j+\frac{n}{2}}) $$
因为 $w^{\frac{n}{2}}=-1$,所以有
$$ d_i=(w^i+w^i)\prod_{0 \le j \ne i < \frac{n}{2}}(w^i-w^j)(w^i+w^{j})=2w^i\prod_{0 \le j \ne i < \frac{n}{2}}(w^{2i}-w^{2j}) $$
重复以上步骤,得到 $\text{log}_2n$ 项乘积
$$ d_i=2w^i * 2w^{2i} * 2w^{4i} * ... * 2w^{n/2i}=n * w^{(n-1)i} $$
因为 $w^n=1$,所以有
$$ d_i=\frac{n}{w^i} $$
另外,对于 2-adic 子群,有
$$ z_H(x) = \prod_{i=0}^{n-1}(x-w^i) =x^n-1 $$
证明:因为 $w^{\frac{n}{2}}=-1$,当 $n=4$ 时,有
$$ \begin{align} & (x-w^0)(x-w^1)(x-w^2)(x-w^3) \\ &=(x-1)(x-w)(x+1)(x-w^3)\\ &=(x^2-1)(x-w)(x+w)\\ &=(x^2-1)(x^2-w^2)\\ &=(x^2-1)(x^2+1)\\ &=x^4-1 \end{align} $$
所以有
$$ M(x)=\prod_{j=0}^n(x-x_j)=(x-x_0)(x-x_1)...(x-x_{n})=x^n-1\\ $$
最终得到
$$ P(x)=\frac{x^n-1}{n} * \sum_i\frac{y_i * w^i}{x-w^i} $$
没有评论