【python】最小二乗法の理論と実装を解説！

1次関数で近似
2次関数で近似
おわりに

最小二乗法をご存じでしょうか？
測定データの近似等でよく用いられる手法です。
大学の時に勉強していたことがほとんどふっとんでるので簡単な理論から復習がてら書きます…

1次関数で近似

簡単のためまずはよくある直線近似からやっていきます。
データを以下の形で近似します。

$y = ax + b$

サンプルデータ

今回は以下のデータの近似を行います。
洗浄装置の昇温データです。

経過時間 x [分]	0	5	10	15	20	25	30
液体温度 y [度]	20.0	24.2	25.5	30.5	32.4	36.5	39.3

電気ヒータで水を温めているため、液体温度と経過時間は以下のような関係になるはずなので一次関数で近似できると想定できます。

$T - T_{0} = \dfrac{Pt}{mc}$

ここで $T$ は液体温度、 $T_{0}$ は初期温度、 $P$ は電力、 $t$ は経過時間、 $m$ は水の質量、 $c$ は比熱です。

実際にプロットしてみましょう。

import numpy as np
import matplotlib.pyplot as plt

# 測定データ
x = np.array([0, 5, 10, 15, 20, 25, 30])
y = np.array([20, 24.2, 25.5, 30.5, 32.4, 36.5, 39.3])

# 測定データのプロット
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.scatter(x, y, label='measure data')
ax.legend(loc='best')

グラフで見ても直線に近似できそうです。

最小二乗法での近似

各時刻と温度のデータをそれぞれ $x_{0}, x_{1}, ..., x_{n}, y_{0}, y_{1}, ..., y_{n}$ とすると近似値と測定値の差(残差)の2乗和 $S$ は以下の式で計算できる。

$\displaystyle{ S = \sum_{i=1}^n (y_{i} - ax_{i}-b)^2 }$

上の $S$ が最小になるような $a, b$ を求めるのが最小二乗法です。
線形問題であれば、 $S$ が最小値となるとき、 $S$ を各 $a, b$ で偏微分したものが0になる性質があるので式を整理すると以下のようにあらわすことができます。

$A^T Au = A^T v$

ここで $A, u, v$ はそれぞれ以下のような行列です。

$A = \begin{pmatrix} x_{1} & 1 \\ x_{2} & 1 \\ \vdots & \vdots \\ x_{n} & 1 \end{pmatrix}, u = \begin{pmatrix} a \\ b \end{pmatrix}, v = \begin{pmatrix} y_{0} \\ y_{1} \\ \vdots \\ y_{n} \end{pmatrix}$

$u$ をもとめたいので、上の式を変換して以下を解けばOKです。

$u = (A^T A)^{-1} A^T v$

Pythonで試してみる

A = np.vstack((x, np.ones(len(x)))).T
v = y.T

# 最小二乗近似
u = np.linalg.inv(A.T @ A) @ A.T @ v

# 近似式のプロット
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.scatter(x, y, label='measure data')
ax.plot(x, x * u[0] + u[1], 'r', label='least square')
ax.legend(loc='best')

グラフで見るとうまく近似できてそうです。

プログラム説明

A = np.vstack((x, np.ones(len(x)))).T

np.vstackやnp.hstackで行列の結合をしています。
np.arrayは.Tをつけることで転置行列に変換できます。

u = np.linalg.inv(A.T @ A) @ A.T @ v

np.linalg.inv()で逆行列変換、@は行列の内積です

2次関数で近似

変更するところ

2次関数であっても基本的には同様です
近似式は以下のようになります

$y = ax^2 + bx + c$

この場合、上記の $A, u, v$ を以下のようにすることで1次関数と同様の方法で近似できます。

$A = \begin{pmatrix} x_{1}^2 & x_{1} & 1 \\ x_{2}^2 & x_{2} & 1 \\ \vdots & \vdots & \vdots \\ x_{n}^2 & x_{n} & 1 \end{pmatrix}, u = \begin{pmatrix} a \\ b \\ c \end{pmatrix}, v = \begin{pmatrix} y_{0} \\ y_{1} \\ \vdots \\ y_{n} \end{pmatrix}$

便利なライブラリ

numpyには便利な関数も用意されているので上記の計算式を覚えてなくても大丈夫です。下記はどちらも同じ結果が得られます。

# 最初に説明したやり方
u = np.linalg.inv(A.T @ A) @ A.T @ v
# np.linalg.lstsqの利用
u, residuals, rank, s = np.linalg.lstsq(A, v, rcond=None)

返り値は、u: 係数のリスト、residuals: 残差二乗和の合計、rank: Aのランク、s: Aの特異値

2次関数近似やってみましょう

# 測定データ
x = np.array([0, 5, 10, 15, 20, 25, 30])
y = 0.1 * x ** 2 - 1 * x + 20 + np.random.normal(0, 3, (len(x)))

# 最小二乗近似
A = np.vstack((x ** 2, x, np.ones(len(x)))).T
v = y.T
u, residuals, rank, s = np.linalg.lstsq(A, v, rcond=None)

# 近似式のプロット
x_graph = np.linspace(0, 30, 100)
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.scatter(x, y, label='measure data')
ax.plot(x_graph, x_graph ** 2 * u[0] + x_graph * u[1] + u[2], \
        'r', label='least square')
ax.legend(loc='best')