Least Squares Method

最小二乗法と多項式近似：離散データから連続な真理を抽出するアルゴリズム

1. イントロダクション：多項式近似の目的

観測や実験から得られた $N$ 個の離散データセット $(x_1, y_1), (x_2, y_2), \dots, (x_N, y_N)$ があるとします。これらのデータの背後にある傾向を最もよく表す $M$ 次の多項式 $P(x)$ を導き出すことが、多項式近似の目的です。

P(x) = a_0 + a_1 x + a_2 x^2 + \dots + a_M x^M = \sum_{j=0}^{M} a_j x^j

最小二乗法（Least Squares Method）では、各データ点における「実際の値 $y_i$」と「多項式の予測値 $P(x_i)$」のズレ（残差）の二乗和 $S$ を最小化するような係数ベクトル $\boldsymbol{a} = (a_0, a_1, \dots, a_M)^T$ を探索します。

S = \sum_{i=1}^{N} (y_i - P(x_i))^2

2. 正規方程式の導出とLU分解への帰着

二乗和 $S$ が最小となる極値を見つけるため、$S$ を各係数 $a_k \ (k = 0, 1, \dots, M)$ で偏微分し、それが $0$ になる条件を求めます。

$$ \frac{\partial S}{\partial a_k} = -2 \sum_{i=1}^{N} x_i^k \left( y_i - \sum_{j=0}^{M} a_j x_i^j \right) = 0 $$

これを展開して $\boldsymbol{a}$ について整理すると、以下のような行列表現による正規方程式（Normal Equations）が得られます。

\begin{pmatrix} N & \sum x_i & \sum x_i^2 & \cdots & \sum x_i^M \\ \sum x_i & \sum x_i^2 & \sum x_i^3 & \cdots & \sum x_i^{M+1} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \sum x_i^M & \sum x_i^{M+1} & \sum x_i^{M+2} & \cdots & \sum x_i^{2M} \end{pmatrix} \begin{pmatrix} a_0 \\ a_1 \\ \vdots \\ a_M \end{pmatrix} = \begin{pmatrix} \sum y_i \\ \sum x_i y_i \\ \vdots \\ \sum x_i^M y_i \end{pmatrix}

連立一次方程式としての求解：
一見複雑な非線形の問題に見えますが、左辺の行列の各要素（$\sum x_i^k$ など）や右辺のベクトルは、与えられたデータ群から一意に計算できる単なる定数です。
つまり、この問題は最終的に $\boldsymbol{Ax} = \boldsymbol{b}$ という標準的な連立一次方程式に帰着します。プログラムで実装する際は、構築した定数行列 $\boldsymbol{A}$ に対してLU分解（LU Decomposition）を適用し、前進代入・後退代入を行うことで、未知の係数 $\boldsymbol{a}$ を高速かつ決定論的に算出することができます。

3. 具体的な計算例（4次多項式）

実際に、ノイズを含んだ以下の7つの観測データから、最小二乗法を用いて4次多項式（$M=4$）を導出するプロセスを見てみましょう。

$i$	$x_i$ (入力)	$y_i$ (観測値)
1	-3.0	47.5
2	-2.0	7.5
3	-1.0	-0.2
4	0.0	1.2
5	1.0	-0.8
6	2.0	-4.5
7	3.0	-0.1

まず、正規方程式を構成するために必要な各成分（$\sum x_i^k$ および $\sum x_i^k y_i$）を計算します。今回、入力 $x_i$ が $-3$ から $3$ まで対称に分布しているため、奇数乗の和（$\sum x_i, \sum x_i^3, \dots$）はすべて $0$ になります。

$N = 7$
$\sum x_i^2 = (-3)^2 + \dots + 3^2 = 28$
$\sum x_i^4 = (-3)^4 + \dots + 3^4 = 196$
$\sum x_i^6 = 1588, \quad \sum x_i^8 = 13636$
$\sum y_i = 47.5 + 7.5 - 0.2 + 1.2 - 0.8 - 4.5 - 0.1 = 50.6$
$\sum x_i y_i = -167.4$
$\sum x_i^2 y_i = 437.6$
$\sum x_i^3 y_i = -1381.8$
$\sum x_i^4 y_i = 3886.4$

これらの定数を正規方程式の行列 $\boldsymbol{A}\boldsymbol{a} = \boldsymbol{b}$ に当てはめると、以下の $5 \times 5$ の連立一次方程式が組み上がります。

\begin{pmatrix} 7 & 0 & 28 & 0 & 196 \\ 0 & 28 & 0 & 196 & 0 \\ 28 & 0 & 196 & 0 & 1588 \\ 0 & 196 & 0 & 1588 & 0 \\ 196 & 0 & 1588 & 0 & 13636 \end{pmatrix} \begin{pmatrix} a_0 \\ a_1 \\ a_2 \\ a_3 \\ a_4 \end{pmatrix} = \begin{pmatrix} 50.6 \\ -167.4 \\ 437.6 \\ -1381.8 \\ 3886.4 \end{pmatrix}

あとは、この係数行列をLU分解などのアルゴリズムで解くだけです。計算の結果、以下の係数ベクトルが得られます。

導出された4次多項式

$P(x) = 1.057 + 1.037x - 2.012x^2 - 1.006x^3 + 0.501x^4$

以下のグラフは、観測データ（赤い点）と、LU分解によって導き出された4次多項式の近似曲線（青い線）をプロットしたものです。ノイズの影響を平滑化しつつ、データの背後にある傾向を見事に捉えていることが視覚的に確認できます。