ラグランジュ補間とその応用(Shamirの秘密分散法)

はじめに

最近、数値解析の講義を受講したluna_moonlightです。この記事では、ラグランジュ補間のPythonによる実装とその応用例について解説していきたいと思います。

多項式補間

ラグランジュ補間の前にまず多項式補間について解説していきます。
$i \neq j \Rightarrow x_i \neq x_j$を満たす$N + 1$個の点$(x_0, y_0), (x_1, y_1), \dots, (x_N, y_N)$から、この$N + 1$個の点を通る$N$次多項式$p_N(x) = a_0 + a_1 x + a_2 x ^ 2 + \dots + a_N x^N$を求めることを$N$次の多項式補間と呼びます。

直感的に分かるかもしれませんが、実はこのような$N$次多項式は存在して唯一つであることが知られています。

証明の概略

補間多項式$p_N(x) = a_0 + a_1 x + a_2 x ^ 2 + \dots + a_N x^N$に$i \neq j \Rightarrow x_i \neq x_j$を満たす$N + 1$個の点のデータ$(x_0, y_0), (x_1, y_1), \dots, (x_N, y_N)$を代入すると、以下のような連立方程式が得られる。

$$\left\{ \begin{alignedat}{6} y_0 &= &a_0 &+ &a_1 x_0 &+ &a_2 x_0^2 &+ &\dots &+ &a_N x_0^N\\ y_1 &= &a_0 &+ &a_1 x_1 &+ &a_2 x_1^2 &+ &\dots &+ &a_N x_1^N\\ \vdots\\ y_N &= &a_0 &+ &a_1 x_N &+ &a_2 x_N^2 &+ &\dots &+ &a_N x_N^N\\ \end{alignedat} \right.$$

以下のような行列とベクトルを定義すると、

$$V = \begin{pmatrix} 1 & x_0 & x_0^2 & \dots & x_0^N\\ 1 & x_1 & x_1^2 & \dots & x_1^N\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ 1 & x_N & x_N^2 & \dots & x_N^N\\ \end{pmatrix}, \bm{y} = \begin{pmatrix} y_0\\ y_1\\ \vdots\\ y_N \end{pmatrix}, \bm{a} = \begin{pmatrix} a_0\\ a_1\\ \vdots\\ a_N \end{pmatrix}$$

上の連立方程式は$\bm{y} = V\bm{a}$と表すことができる。

ここで、行列$V$はVandermonde行列と呼ばれ、その行列式は以下のようになる。

$$|V| = \prod_{i < j}(x_i - x_j)$$

$i \neq j \Rightarrow x_i \neq x_j$なので、$|V| \neq 0$となるため、$\bm{y} = V\bm{a}$の解は存在して唯一つである。
よって、$i \neq j \Rightarrow x_i \neq x_j$を満たす$N + 1$個の点$(x_0, y_0), (x_1, y_1), \dots, (x_N, y_N)$を通る$N$次多項式は存在して唯一つである。

3点の例

それでは、試しに以下のような3点$(0, 1), (3, 4), (4, 2)$で多項式補間をしてみましょう。

2次の多項式補間なので、$p_2(x) = a + b x + c x ^ 2$の$a, b, c$を求めていきましょう。3点のデータを代入すると以下のような連立方程式が得られます。

$$\left\{ \begin{alignedat}{3} 1 &= a\\ 4 &= a &+ 3b &+ 9c\\ 2 &= a &+ 4b &+ 16c\\ \end{alignedat} \right.$$

これを頑張って解くと、$a = 1, b = \dfrac{13}{4}, c = -\dfrac{3}{4}$となるため、$p_2(x) = 1 + \dfrac{13}{4} x -\dfrac{3}{4} x ^ 2$であることがわかりました。
実際にグラフを描画してみると以下のようになり、きちんと3点を通るグラフが求められたことがわかります。

N + 1点の例

次は$N + 1$点$(x_0, y_0), (x_1, y_1), \dots, (x_N, y_N)$で多項式補間をしてみましょう。

$N$次の多項式補間なので、$p_N(x) = a_0 + a_1 x + a_2 x ^ 2 + \dots + a_N x^N$を求めていきましょう。$N + 1$点のデータを代入すると以下のような連立方程式が得られます。

$$\left\{ \begin{alignedat}{6} y_0 &= &a_0 &+ &a_1 x_0 &+ &a_2 x_0^2 &+ &\dots &+ &a_N x_0^N\\ y_1 &= &a_0 &+ &a_1 x_1 &+ &a_2 x_1^2 &+ &\dots &+ &a_N x_1^N\\ \vdots\\ y_N &= &a_0 &+ &a_1 x_N &+ &a_2 x_N^2 &+ &\dots &+ &a_N x_N^N\\ \end{alignedat} \right.$$

このような連立方程式を解くためには掃き出し法と呼ばれる手法を使います。実際に、掃き出し法を行っている様子を見てみましょう。

掃き出し法は一般の連立一次方程式で使えるため、係数を以下のように変更して解説していきます。

$$\left\{ \begin{alignedat}{6} y_0 &= &a_0 x_{0, 0} &+ &a_1 x_{0, 1} &+ &a_2 x_{0, 2} &+ &\dots &+ &a_N x_{0, N}\\ y_1 &= &a_0 x_{1, 0} &+ &a_1 x_{1, 1} &+ &a_2 x_{1, 2} &+ &\dots &+ &a_N x_{1, N}\\ \vdots\\ y_N &= &a_0 x_{N, 0} &+ &a_1 x_{N, 1} &+ &a_2 x_{N, 2} &+ &\dots &+ &a_N x_{N, N}\\ \end{alignedat} \right.$$

まず、第一行をもとにそれ以外の行の第一列の係数を0にします。

基準の行の係数が0だと、0除算が発生し計算が続行できないため、非ゼロの係数を持つ行と入れ替えてください。この入れ替え操作をピボッティングと呼びます。

$$\left\{ \begin{alignedat}{6} y_0 &= &a_0 x_{0, 0} &+ &a_1 x_{0, 1} &+ &a_2 x_{0, 2} &+ &\dots &+ &a_N x_{0, N}\\ y_1' &= &0 &+ &a_1 x_{1, 1}' &+ &a_2 x_{1, 2}' &+ &\dots &+ &a_N x_{1, N}'\\ \vdots\\ y_N' &= &0 &+ &a_1 x_{N, 1}' &+ &a_2 x_{N, 2}' &+ &\dots &+ &a_N x_{N, N}'\\ \end{alignedat} \right.$$

次に第二行をもとにそれ以外の行の第二列の係数を0にします。

$$\left\{ \begin{alignedat}{6} y_0'' &= &a_0 x_{0, 0} &+ &0 &+ &a_2 x_{0, 2}'' &+ &\dots &+ &a_N x_{0, N}''\\ y_1' &= &0 &+ &a_1 x_{1, 1}' &+ &a_2 x_{1, 2}' &+ &\dots &+ &a_N x_{1, N}'\\ \vdots\\ y_N'' &= &0 &+ &0 &+ &a_2 x_{N, 2}'' &+ &\dots &+ &a_N x_{N, N}''\\ \end{alignedat} \right.$$

これを繰り返していって、対角成分のみが残るように変形していきます。

$$\left\{ \begin{alignedat}{6} y_0^* &= &a_0 x_{0, 0} &+ &0 &+ &0 &+ &\dots &+ &0\\ y_1^* &= &0 &+ &a_1 x_{1, 1}' &+ &0 &+ &\dots &+ &0\\ \vdots\\ y_N^* &= &0 &+ &0 &+ &0 &+ &\dots &+ &a_N x_{N, N}^*\\ \end{alignedat} \right.$$

最後に、各行を$x$で割り算することで答えが求まります。

$$\begin{split} a_0 &= \dfrac{y_0^*}{x_{0, 0}}\\ a_1 &= \dfrac{y_1^*}{x_{1, 1}'}\\ \vdots\\ a_N &= \dfrac{y_N^*}{x_{N, N}^*}\\ \end{split}$$

from fractions import Fraction

x = [0, 3, 4]
y = [1, 4, 2]
N = len(x)

x = [Fraction(i) for i in x]
y = [Fraction(i) for i in y]

A = [[x[j] ** i for i in range(N)] for j in range(N)]

for i in range(N):
    # ピボッティング
    for j in range(i, N):
        if A[j][i] != 0:
            A[i], A[j] = A[j], A[i]
            break
    
    # 消去操作
    for j in range(N):
        if i == j:
            continue
        v = -A[j][i] / A[i][i]
        for k in range(N):
            A[j][k] += A[i][k] * v
        y[j] += y[i] * v

ans = [y[i] / A[i][i] for i in range(N)]
for i, a in enumerate(ans):
    print(f"a_{i} = {a}")

この方法だと、データ点の数を$N$として$O(N^3)$で答えが求まります。

ラグランジュ補間

実は、先ほどの連立方程式を解かなくても補間多項式は以下のように求めることができます。

$$p_N(x) = \sum_{i=0}^N y_i l_i(x)$$

ただし、$l_i(x)$は以下の式で与えられます。

$$l_i(x) = \prod_{j=0, j \neq i}^N \dfrac{(x - x_j)}{(x_i - x_j)} = \dfrac{(x - x_0)(x - x_1)\dotsb \textcolor{red}{(x - x_{j - 1})(x - x_{j + 1})} \dotsb(x - x_N)}{(x_i - x_0)(x_i - x_1)\dotsb \textcolor{red}{(x_i - x_{j - 1})(x_i - x_{j + 1})} \dotsb(x_i - x_N)}$$

このような方法を$N$次のラグランジュ補間と呼びます。

正当性の証明の概略

$$l_i(x_j) = \delta_{i, j} = \left\{ \begin{split} 1 (i = j)\\ 0 (i \neq j)\\ \end{split} \right.$$

$l_i(x_j)$は上のような関係になるので、$p_N(x_i) = y_i$となり、$p_N(x)$は与えられた点を通る補間多項式になっている。

さらに補間多項式の一意性により、多項式補間の結果と一致する。

3点の例

それでは、試しに多項式補間のときと同じ例の3点$(0, 1), (3, 4), (4, 2)$でラグランジュ補間をしてみましょう。

$$\begin{aligned} l_0(x) &= \dfrac{(x - 3)(x - 4)}{(0 - 3)(0 - 4)} &&= 1 - \dfrac{7}{12} x + \dfrac{1}{12} x^2\\ l_1(x) &= \dfrac{(x - 0)(x - 4)}{(3 - 0)(3 - 4)} &&= \dfrac{4}{3} x - \dfrac{1}{3} x^2\\ l_2(x) &= \dfrac{(x - 0)(x - 3)}{(4 - 0)(4 - 3)} &&= -\dfrac{3}{4} x + \dfrac{1}{4} x^2\\ \end{aligned}$$

よって、

$$\begin{aligned} p_2(x) &= 1 \cdot (1 - \dfrac{7}{12} x + \dfrac{1}{12} x^2) + 4 \cdot (\dfrac{4}{3} x - \dfrac{1}{3} x^2) + 2 \cdot (-\dfrac{3}{4} x + \dfrac{1}{4} x^2)\\ &= 1 + \dfrac{13}{4} x -\dfrac{3}{4} x ^ 2 \end{aligned}$$

すなわち、$p_2(x) = 1 + \dfrac{13}{4} x -\dfrac{3}{4} x ^ 2$となり、多項式補間の結果と一致しました。

N + 1点の例

次は$N + 1$点$(x_0, y_0), (x_1, y_1), \dots, (x_N, y_N)$でラグランジュ補間をしてみましょう。

愚直

愚直に掛け算して$l_i(x)$を計算して$y_i$をかけた結果を足し合わせて$p_N(x)$を求めます。

from fractions import Fraction

x = [0, 3, 4]
y = [1, 4, 2]
N = len(x)

x = [Fraction(i) for i in x]
y = [Fraction(i) for i in y]

ans = [0] * N
for i in range(N):
    dp = [0] * (N + 1)
    dp[1] = 1
    v = 1
    for j in range(N):
        if i == j:
            continue
        for k in reversed(range(1, N + 1)):
            dp[k] = dp[k - 1] - dp[k] * x[j]
        v *= x[i] - x[j]
    
    for j in range(N):
        ans[j] += y[i] / v * dp[j + 1]

for i, a in enumerate(ans):
    print(f"a_{i} = {a}")

この方法だと、データ点の数を$N$として$O(N^3)$で答えが求まります。

愚直に計算するとは言っていますが、実は今回の実装にはいくつかの重要なテクニックを使っています。
$\textrm{dp}_{i, j}$を$l_k(x)$の$(x - a_i)$の掛け算まで行ったときの$x^j$の係数(分子は$j$に依らないので無視)とすると、$\textrm{dp}_{i, j}$は以下のような漸化式で求まります。

初期状態: $\textrm{dp}_{0, j} = \left\{\begin{aligned}&1 &(j = 0)\\&0 &(j \gt 0)\\\end{aligned}\right.$

更新式: $\textrm{dp}_{i, j} = \left\{\begin{aligned}&- a_i \textrm{dp}_{i - 1, j} &(j = 0)\\&\textrm{dp}_{i - 1, j - 1} - a_i \textrm{dp}_{i - 1, j} &(j \gt 0)\\\end{aligned}\right.$

このように、目的の問題を部分問題(漸化式)に分割して目的の問題を解く方法を、動的計画法やDPと呼びます。

今回の実装では、もう一つのテクニックを活用しています。
$\textrm{dp}_{i, j}$の更新式をみると、$i$の結果を求めるために$i - 1$の結果しか利用していません。さらに、$j$の結果の方も$j - 1$(前)の結果しか利用していないため、$j$を後ろから更新していけば、$\textrm{dp}_{i, j}$を二次元ではなく、$\textrm{dp}_j$のように一次元に圧縮して管理することができます。

このように、計算途中の情報すべてを保持せずに、必要な情報をその場で上書きしながら計算することで空間計算量を削減するDPをインラインDPと呼びます。

高速化

$l_i(x)$分母の式を式変形すると以下のようになることがわかります。

$$\prod_{j=0, j \neq i}^N (x - x_j) = \dfrac{\displaystyle\prod_{j=0}^N (x - x_j)}{(x - x_i)}$$

どの$i$に対しても$\displaystyle\prod_{j=0}^N (x - x_j)$は共通なのでこれを事前に計算しておき、$(x - x_i)$で割る操作を高速に行うことができれば、ラグランジュ補間も高速に行えます。

$\displaystyle\prod_{j=0}^N (x - x_j)$は、先ほどと同様にDPを使えば求めることができます。少し考えると、このDPは掛け算する順番を変えても問題ないことがわかります。そこで、$(x - x_i)$の掛け算を最後に行ったとしてDPの更新を逆に適用する、すなわち以下のような更新式を適用することで$(x - x_i)$の掛け算をなかったことに$\bigl( (x - x_i)$で割り算$\bigr)$することができます。このようにDPを戻す($\textrm{dp}_i$から$\textrm{dp}_{i - 1}$を復元する)方法を戻すDPと呼びます。

$$\textrm{dp}_{i - 1, j} = \left\{ \begin{aligned} &-\dfrac{\textrm{dp}_{i, j}}{x_i} &(j = 0)\\ &-\dfrac{\textrm{dp}_{i, j} - \textrm{dp}_{i - 1, j - 1}}{x_i} &(j \gt 0)\\ \end{aligned} \right.$$

$x_i = 0$のときは0除算が発生してしまうため、以下の更新式を使う必要があります。

$$\textrm{dp}_{i - 1, j - 1} = \textrm{dp}_{i, j} (j > 0)$$

from fractions import Fraction

x = [0, 3, 4]
y = [1, 4, 2]
N = len(x)

x = [Fraction(i) for i in x]
y = [Fraction(i) for i in y]

dp = [0] * (N + 2)
dp[1] = 1
for i in range(N):
    for j in reversed(range(1, N + 2)):
        dp[j] = dp[j - 1] - dp[j] * x[i]

ans = [0] * N
for i in range(N):
    v = 1
    for j in range(N):
        if i == j:
            continue
        v *= x[i] - x[j]
    
    dp2 = [0] * (N + 1)
    for j in range(1, N + 1):
        if x[i] == 0:            
            dp2[j] = dp[j + 1]
        else:
            dp2[j] = -(dp[j] - dp2[j - 1]) / x[i]
    
    for j in range(N):
        ans[j] += y[i] / v * dp2[j + 1]

for i, a in enumerate(ans):
    print(f"a_{i} = {a}")

この方法だと、データ点の数を$N$として$O(N^2)$で答えが求まります。

実行時間

参考までに、データ点の数による実行時間の違いを載せておきます。

ここで載せたプログラムはFractionを使って厳密解を求めており、データの点の数が多くなるほど四則演算の計算時間が増えてしまって傾向が分かりにくいため、下の表はFractionを使わない場合の実行時間を載せています。
標準のfloatでは大きな数字を扱うことができないため、Fractionの代わりにDecimalを使っています。

	$10^1$	$10^2$	$10^3$	$10^4$	$10^5$
多項式補間	$2.83 \cdot 10^{-4} \mathrm{s}$	$0.243 \mathrm{s}$	$291 \mathrm{s}$	-	-
ラグランジュ補間(愚直)	$2.51 \cdot 10^{-4} \mathrm{s}$	$0.200 \mathrm{s}$	$202 \mathrm{s}$	-	-
ラグランジュ補間(戻すDP)	$1.31 \cdot 10^{-4} \mathrm{s}$	$1.26 \cdot 10^{-2} \mathrm{s}$	$1.30 \mathrm{s}$	$128 \mathrm{s}$	$1.28 \cdot 10^4 \mathrm{s}$

厳密解

Fractionを使って厳密解を求めた場合の実行時間も載せておきます。

	$10^1$	$10^2$	$10^3$
多項式補間	$2.25 \cdot 10^{-3} \mathrm{s}$	$2.28 \mathrm{s}$	$1.72 \cdot 10^{4} \mathrm{s}$
ラグランジュ補間(愚直)	$1.76 \cdot 10^{-3} \mathrm{s}$	$2.04 \mathrm{s}$	$3.86 \cdot 10^{3} \mathrm{s}$
ラグランジュ補間(戻すDP)	$9.98 \cdot 10^{-4} \mathrm{s}$	$0.153 \mathrm{s}$	$556 \mathrm{s}$

応用例

ラグランジュ補間の補間以外の面白い応用例として、Shamirの秘密分散法を解説していきます。

秘密分散法とは、データをいくつかのシェアに分割して複数人(パーティ)の間で共有する暗号化手法です。ここで重要なことは、一定数のシェアを集めることで初めて元のデータを復元できる、すなわち一定数未満のシェアからは元のデータに関する情報は得られないということです。

Shamirの秘密分散法について解説する前に、加法的秘密分散法について解説します。
秘密情報$s$を$n$人で分割することを考えます。秘密情報提供者が$s = s_1 + s_2 + \dots + s_n$となるランダムなシェア$s_i$を$i$さんに配布します。このようなシェアを作成すると、すべてのシェアが集まれば秘密情報$s$が復元できます。また、シェアが1つでも欠けている場合は、秘密情報$s$はどんな値も取りうるため復元することができません。

Shamirの秘密分散法で秘密情報$s$を$n$人で分割することを考えます。秘密情報提供者がランダムな$k(k < n)$次多項式$f(x) = s + a_1 x + \dots + a_k x^k$を作成し、シェア$f(i)$を$i$さんを配布します。このようなシェアを作成すると、$k + 1$個以上のシェアが集まればラグランジュ補間を用いて$k$次多項式$f$を復元し、$f(0) = s$を求めることで秘密情報$s$が復元することがます。また、シェアが$k$個以下しか集まらない場合は、あらゆる$k$次多項式が考えられ、秘密情報$s$はどんな値も取りうるため復元することができません。

このように、加法的秘密分散法ではすべてのシェアを集まらないと秘密情報を復元できないのに対し、Shamirの秘密分散法では閾値以上のシェアが集まれば秘密情報を復元することができます。さらにShamirの秘密分散法では、複数のシェアを受け取ることを許容することで、任意のアクセス構造(役職者が2人以上集まらないと復元できないなど)を実現できます。

まとめ

この記事では、ラグランジュ補間のPythonによる実装とその応用例としてShamirの秘密分散法について解説しました。ラグランジュ補間は単なる関数補間にとどまらず、暗号分野でも重要な役割を果たしていることがわかりました。