Vandermonde linear systemとラグランジュ補間
執筆途中...
c/3e0a9138-3743-40dd-86e9-ac390cb50757
復習: Vandermonde
Vandermonde行列
列か行が等比数列を成している行列
$ \begin{pmatrix} 1 & x_1 & x_1^2 & \cdots & x_1^{n-1} \\ 1 & x_2 & x_2^2 & \cdots & x_2^{n-1} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_n & x_n^2 & \cdots & x_n^{n-1} \end{pmatrix}
Vandermonde linear system
Vandermonde行列で構成された線形方程式系
一般に$ \bm{a}は未知で$ \bm{b}は既知
$ \begin{pmatrix} 1 & x_1 & x_1^2 & \cdots & x_1^{n-1} \\ 1 & x_2 & x_2^2 & \cdots & x_2^{n-1} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_n & x_n^2 & \cdots & x_n^{n-1} \end{pmatrix} \begin{pmatrix} a_0 \\ a_1 \\ \vdots \\ a_{n-1} \end{pmatrix} = \begin{pmatrix} b_0 \\ b_1 \\ \vdots \\ b_n \end{pmatrix}
Pairwise distinctならば,行列式は非零なのでVandermonde行列は常に正則
$ \left| \begin{matrix} 1 & x_1 & x_1^2 & \cdots & x_1^{n-1} \\ 1 & x_2 & x_2^2 & \cdots & x_2^{n-1} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_n & x_n^2 & \cdots & x_n^{n-1} \end{matrix} \right| = \prod_{1 \leq i < j \leq n} (x_j - x_i) \neq 0
Vandermonde linear systemによる多項式補間
与えられた一連のデータ点 $ \{(x_1, y_1), (x_2, y_2), \ldots, (x_n, y_n)\} に対して,それらを通る多項式 $ P(x) を探すことが目的.
$ P(x) = a_0 + a_1 x + a_2 x^2 + \cdots + a_{n-1} x^{n-1}
ここで,$ \{a_0, a_1, \ldots, a_{n-1}\} は多項式の係数を指す.
このとき,各データ点 $ (x_i, y_i) について,上記の多項式 $ P(x) が $ y_i に等しくなるような方程式を立てると,以下のように Vandermonde linear systemで記述できる.
$ \begin{pmatrix} 1 & x_1 & x_1^2 & \cdots & x_1^{n-1} \\ 1 & x_2 & x_2^2 & \cdots & x_2^{n-1} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_n & x_n^2 & \cdots & x_n^{n-1} \end{pmatrix} \begin{pmatrix} a_0 \\ a_1 \\ \vdots \\ a_{n-1} \end{pmatrix} = \begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{pmatrix}
上記式を解くことで,$ \{a_0, a_1, \ldots, a_{n-1}\}が求めれば多項式での補間ができる
実験: 正弦波を多項式補間してみる
まずはvanillaな正弦波に多項式補間を適用
code: poly.py
import numpy as np
import matplotlib.pyplot as plt
x = np.linspace(0, 2*np.pi, 10) # 0から2πまでの範囲で10点
y = np.sin(x)
# Vandermonde行列
V = np.vander(x, increasing=True)
coefficients = np.linalg.solve(V, y)
def polynomial(x, coeffs):
n = len(coeffs)
return sum([coeffsi * x**i for i in range(n)]) # visualize
x_interp = np.linspace(0, 2*np.pi, 1000)
y_interp = polynomial(x_interp, coefficients)
plt.plot(x, y, 'o', label='Original data points')
plt.plot(x_interp, y_interp, label='Interpolated polynomial')
plt.legend()
plt.show()
https://gyazo.com/6edc38930541a72b343af9c474a96e65
次にホワイトノイズが加わった正弦波に多項式補間を適用
code:white_poly.py
import numpy as np
import matplotlib.pyplot as plt
x = np.linspace(0, 2*np.pi, 10)
y = np.sin(x)
# ホワイトノイズを追加
np.random.seed(0)
noise_amplitude = 0.5
y_noisy = y + noise_amplitude * np.random.randn(len(y))
V = np.vander(x, increasing=True)
coefficients_noisy = np.linalg.solve(V, y_noisy)
def polynomial(x, coeffs):
n = len(coeffs)
return sum([coeffsi * x**i for i in range(n)]) x_interp = np.linspace(0, 2*np.pi, 1000)
y_interp_noisy = polynomial(x_interp, coefficients_noisy)
plt.plot(x, y, 'o', label='Original data points')
plt.plot(x, y_noisy, 'r.', label='Noisy data points')
plt.plot(x_interp, y_interp_noisy, label='Interpolated polynomial (with noise)')
plt.legend()
plt.show()
https://gyazo.com/10e0a9c810df8211d67bbb99fba1ed3c
Vandermondeからラグランジュ補間へ
上記Vandermode linear systemをそのまま計算するのは計算量が大きすぎる・逆行列の計算による誤差が大きすぎるので,ラグランジュ基底多項式 $ L_i(x)と呼ばれる多項式を導入する.
$ L_i(x) = \prod_{j \neq i} \frac{x - x_j}{x_i - x_j} = \frac{(x - x_0)(x - x_1) \cdots (x - x_{i-1})(x - x_{i+1}) \cdots (x - x_d)}{(x_i - x_0)(x_i - x_1) \cdots (x_i - x_{i-1})(x_i - x_{i+1}) \cdots (x_i - x_d)}