ラグランジュ補間

ラグランジュの補間公式とその応用例 | 高校数学の美しい物語

できること

具体的な係数はわかんないけど、$f(x)=a_0+a_1x+a_2x^2+a_3x^3+...+a_Dx^D$ みたいに $x$ の $D$ 次多項式で表現できることはわかっている関数がある。

また、$D+1$ 個の入力と解の組合せ $(x_0,f(x_0)), (x_1,f(x_1)), ..., (x_D, f(x_D))$ がわかっている。（$x_i$ は互いに異なる）

ここである $N$ が1つ与えられた時、$f(N)$ は $O(D^2)$ で計算することができる。（値が巨大で多倍長演算が必要になるとかは別として）

また、$x_i$ の並びがある条件を満たせば、$O(D)$ で計算することもできる。

係数 $a_0,a_1,...,a_D$ は $O(D^2)$ で復元できるので、複数の $N$ に対して $f(N)$ が必要な場合は係数を復元する。
（$O(D(\log{D})^2)$ で出来るという話もある？）

※ラグランジュ補間は、上記のような厳密解を求める用途でなくとも、 $D+1$ 個の点を全て通る曲線の方程式を $D$ 次多項式と仮定してスムーズに結ぶ用途などにも使える。
ただ点を増やすと過適合によって暴れやすいので、一般的にはスプライン補間など別のを使うことが多い印象。

やりかた

$D$ 次多項式は係数が $D+1$ 個ある。
$f(x)$ に $x_i$ をそれぞれ代入すると、$D+1$ 個の連立方程式ができる。

$D+1$ 個の未知数に対して $D+1$ 個の式があるため、$f(x)$ は今ある情報で一意に決まる。

f(x) = a + b*x + c*x^2
f(x) は (1, 1)  (2, 3)  (3, 11)  を通る
  ↓
 1 = a +  b +  c
 3 = a + 2b + 4c
11 = a + 3b + 9c

補足

だが、実際に連立方程式や逆行列を解くのは $D$ が大きいと難しい。

そこで上手いこと式を作ると、（一見ややこしいものの）機械的に処理しやすい形で $f(x)$ を定義できる。

$\displaystyle f(x)= \sum_{i=0}^D{f(x_i) \prod_{j=0..D,j \neq i}\frac{x-x_j}{x_i-x_j}}$

                       (x-2) * (x-3)
f(x) = f(1) * -----------------------
                       (1-2) * (1-3)
     
               (x-1) *         (x-3)
     + f(2) * -----------------------
               (2-1) *         (2-3)
     
               (x-1) * (x-2)
     + f(3) * -----------------------
               (3-1) * (3-2)

例えば $x=1$ の時、一番上の分数は $1$ になり、他の分数は $0$ になる。
$x=2$ の時、真ん中の分数は $1$ になり、他は $0$ になる。

このように、与えられた $D+1$ 点を必ず通過する $D$ 次多項式を機械的に作ることができる。
$x$ に $N$ など好きな値を代入すると、$f(N)$ が求まる。

計算量は、分数1つあたりに $D$ 項のかけ算割り算、それを $D+1$ 回計算するので、$O(D^2)$ となる。

特殊な場合における高速化

求めたい $N$ が1点、サンプル点が等差数列

$D+1$ 個のサンプル点 $x_i$ が等差数列 $x_i=ai+b$ の場合は $O(D)$ で求められる。

簡単に、$(x_0,x_1,...,x_D)=(0,1,...,D)$ である場合を考える。

この形なら、分母の計算に $x_i$ ごとに $O(D)$ かかっていたところ、前計算によって $O(1)$ でできるようになる。

分子

$D$ や $N$ が相応に大きい場合、桁数が増えすぎてオーバーフローが発生する。
そういう場合、答えを素数で割ったあまりで複数個求めれば復元できるので、$\mod P$ で求めることを前提とする。

簡単には、全ての積 $(N-x_0)(N-x_1)...(N-x_D)$ をあらかじめ計算しておいて、 $i$ の場合の分子を求める際は $(N-x_i)$ で割ればいい（$N-x_i$ が0になるようなら、$N$ はサンプル点である）。

これなら $(N-x_i)$ の逆数計算にそれぞれ $\log{P}$ ずつ必要になり、全体の計算量は $O(D \log{P})$ となる。

もう少し頑張れば、左からの累積積、右からの累積積をそれぞれ計算しておくことで、

i        0    1           2                 3
Left     1  (N-x0)  (N-x0)(N-x1)  (N-x0)(N-x1)(N-x2)  ...

i                     D-3                      D-2          D-1   D
Right  ...  (N-x[D-2])(N-x[D-1])(N-xD)  (N-x[D-1])(N-xD)  (N-xD)  1

$x_i$ に対する分子は $Left[i] \times Right[i]$ で求まるため、前計算 $O(D)$、クエリ $O(1)$ になる。

分母

例えば $D=8$ で $i=3$ の時、分母は以下のようになる。

\begin{align} &(3-0)(3-1)(3-2)(3-4)(3-5)(3-6)(3-7)(3-8) \\ &= 3 \cdot 2 \cdot 1 \cdot -1 \cdot -2 \cdot -3 \cdot -4 \cdot -5 \\ &= 3! \times 5! \times (-1)^5 \end{align}

$i=4$ の時は

\begin{align} &(4-0)(4-1)(4-2)(4-3)(4-5)(4-6)(4-7)(4-8) \\ &= 4 \cdot 3 \cdot 2 \cdot 1 \cdot -1 \cdot -2 \cdot -3 \cdot -4 \\ &= 4! \times 4! \times (-1)^4 \end{align}

より一般的に書くと、$i!(D-i)!(-1)^{D-i}$ となる。

$i$ をずらしても常に階乗にまとめられるのが等間隔になっている利点であり、階乗の逆数を前計算しておけばよい。

これで分子・分母ともに $O(1)$、つまり分数1個あたりの計算が $O(1)$ になったことで、全体が $O(D)$ でできるようになった。

サンプル点が等比数列

詳しく調べていないが、$O(D \log{D})$ でできるらしい。

多項式かどうか

ラグランジュ補間がすごいのはわかったけど、与えられた問題の答えが多項式かどうかってどうやってわかるの？

基本的には以下の2つが出てくることが多いのかな……。
他にもあるかも。

二項係数

答えが二項係数で表せる場合。

$f(x)={}_xC_r$（$r$ は適当な定数）と表現できる場合、$f(x)$ は $x$ の $r$ 次多項式である。

まぁ、階乗を展開すればそのまんま。（分子が $r$ 個の項の積になる）

$$ {}_xC_r = \frac{x(x-1)(x-2)...(x-r+1)}{r(r-1)(r-2)...1}$$

総和

答えが何らかの多項式の総和を取った形で表せる場合。

$x$ の $D$ 次多項式 $g(x)$ があって、その総和 $\displaystyle f(x)=\sum_{i=1}^xg(i)$ は、$x$ の $D+1$ 次多項式となる。

目次