最小二乗法 (さいしょうにじょうほう、さいしょうじじょうほう;最小自乗法 とも書く、英 : least squares method )は、誤差を伴う測定値の処理において、その誤差の二乗の和を最小にするようにし、最も確からしい関係式を求める方法である。測定で得られた数値の組を、適当なモデルから想定される1次関数 、対数 曲線など特定の関数 を用いて近似 するときに、想定する関数が測定値に対してよい近似となるように、残差平方和 を最小とするような係数を決定する方法[1] [2] [3] 、あるいはそのような方法によって近似を行うことである[1] [2] [3] 。
データセットを4次関数で最小二乗近似した例
1805年 にアドリアン=マリ・ルジャンドル が出版したのが初出である。しかし、1809年 にカール・フリードリヒ・ガウス が出版した際に1795年 には最小二乗法を考案済みだったと主張したことで、最小二乗法の発明者が誰であるかについては不明になっている。
前提条件
最小二乗法では測定データ
y
{\displaystyle y}
はモデル関数
f
(
x
)
{\displaystyle f(x)}
と誤差
ε
{\displaystyle \varepsilon }
の和で
y
=
f
(
x
)
+
ε
{\displaystyle y=f(x)+\varepsilon }
と表せるとする。物理現象の測定データには、誤差 が含まれ、それは系統誤差 と偶然誤差 を含んでいる。この内、偶然誤差は、測定における信号経路の微視的 現象に由来するならば、正規分布 であると期待されることが多い。また、社会調査などの誤差理由の特定が困難な場合でも誤差が正規分布になると期待する考え方 もある。
誤差が正規分布に従わない場合、最小二乗法によって得られたモデル関数はもっともらしくないことに注意する必要がある。偶然誤差が正規分布していない場合、系統誤差が無視できない位大きくそれをモデル関数に含めていない場合、測定データに正規分布から大きく外れた外れ値 を含む場合などが該当する。
上記を含め、最小二乗法の理論的基盤には次のような前提が設けられている[1] [2] [3] 。
測定値の誤差には偏りがない。すなわち誤差の平均値 は 0 である。
測定値の誤差の分散 は既知である。ただし測定データごとに異なる値でも良い。
各測定は互いに独立 であり、誤差の共分散 は 0 である。
誤差は正規分布する。
m
{\displaystyle m}
個のパラメータ[注釈 1] (フィッティングパラメータ)を含むモデル関数
f
{\displaystyle f}
が知られていて、測定量の真の値を近似誤差なく再現することのできるパラメータが存在する。
基礎的な考え方
話を簡単にするため、測定値は
x
,
y
{\displaystyle x,y}
の二次元の平面に分布するものとし、想定される分布(モデル関数)が
y
=
f
(
x
)
{\displaystyle y=f(x)}
の形である場合を述べる。想定している関数
f
{\displaystyle f}
は、既知の関数
g
(
x
)
{\displaystyle g(x)}
の線型結合 で表されていると仮定する。すなわち、
f
(
x
)
=
∑
k
=
1
m
a
k
g
k
(
x
)
{\displaystyle f(x)=\sum _{k=1}^{m}a_{k}g_{k}(x)}
例えば、
g
k
(
x
)
=
x
k
−
1
{\displaystyle g_{k}(x)=x^{k-1}}
は、多項式近似であり、特に
m
=
2
{\displaystyle m=2}
の時は
f
(
x
)
=
a
1
+
a
2
x
{\displaystyle f(x)=a_{1}+a_{2}x}
という直線による近似(線形回帰 )になる。
今、測定で得られた、次のような数値の組の集合があるとする。
{
(
x
1
,
y
1
)
,
(
x
2
,
y
2
)
,
…
,
(
x
n
,
y
n
)
}
{\displaystyle \{(x_{1},y_{1}),\ (x_{2},y_{2}),\ \ldots ,\ (x_{n},y_{n})\}}
これら
(
x
,
y
)
{\displaystyle (x,y)}
の分布が、
y
=
f
(
x
)
{\displaystyle y=f(x)}
というモデル関数に従うと仮定した時、想定される理論値は
(
x
1
,
f
(
x
1
)
)
,
(
x
2
,
f
(
x
2
)
)
,
.
.
.
,
(
x
n
,
f
(
x
n
)
)
{\displaystyle (x_{1},f(x_{1})),(x_{2},f(x_{2})),...,(x_{n},f(x_{n}))}
ということになり、実際の測定値との残差は、各
i
{\displaystyle i}
につき
|
y
i
−
f
(
x
i
)
|
{\displaystyle |y_{i}-f(x_{i})|}
ということになる。
この残差の大きさは、
x
y
{\displaystyle xy}
平面上での
(
x
i
,
y
i
)
{\displaystyle (x_{i},y_{i})}
と
(
x
i
,
f
(
x
i
)
)
{\displaystyle (x_{i},f(x_{i}))}
との距離でもある。
ここで、理論値からの誤差の分散 の推定値は残差の平方和
J
=
∑
i
=
1
n
(
y
i
−
f
(
x
i
)
)
2
{\displaystyle J=\sum _{i=1}^{n}(y_{i}-f(x_{i}))^{2}}
で与えられるから、
J
{\displaystyle J}
が最小になるように想定分布
f
{\displaystyle f}
を(すなわち
a
k
{\displaystyle a_{k}}
を)、定めればよいということになる。
それには、上式は
a
k
{\displaystyle a_{k}}
を変数とする関数と見なすことができるので、
J
{\displaystyle J}
を
a
k
{\displaystyle a_{k}}
について偏微分したものを 0 と置く。こうして得られた
m
{\displaystyle m}
個の連立方程式(正規方程式)を解き、
a
k
{\displaystyle a_{k}}
を決定すればよい。
当てはめたい関数
f
{\displaystyle f}
は、
f
(
x
)
=
(
g
1
(
x
)
,
g
2
(
x
)
,
…
,
g
m
(
x
)
)
(
a
1
,
a
2
,
…
,
a
m
)
T
{\displaystyle f(x)=(g_{1}(x),g_{2}(x),\ldots ,g_{m}(x))(a_{1},a_{2},\ldots ,a_{m})^{\textrm {T}}}
と表すことができる。上付き添字 T は転置行列 を表す。最小にすべき関数
J
{\displaystyle J}
は
J
(
a
)
=
(
G
a
−
y
)
T
(
G
a
−
y
)
=
(
[
G
y
]
[
a
−
1
]
)
T
(
[
G
y
]
[
a
−
1
]
)
{\displaystyle {\begin{aligned}J({\boldsymbol {a}})&=(G{\boldsymbol {a}}-{\boldsymbol {y}})^{\textrm {T}}\,(G{\boldsymbol {a}}-{\boldsymbol {y}})\\&=\left({\begin{bmatrix}G&{\boldsymbol {y}}\end{bmatrix}}{\begin{bmatrix}{\boldsymbol {a}}\\-1\end{bmatrix}}\right)^{\textrm {T}}\,\left({\begin{bmatrix}G&{\boldsymbol {y}}\end{bmatrix}}{\begin{bmatrix}{\boldsymbol {a}}\\-1\end{bmatrix}}\right)\end{aligned}}}
と表される。ここに
G
{\displaystyle G}
は、
G
i
j
=
g
j
(
x
i
)
{\displaystyle G_{ij}=g_{j}(x_{i})}
なる成分を持つ
n
×
m
{\displaystyle n\times m}
行列、
y
=
(
y
1
,
y
2
,
…
,
y
n
)
T
{\displaystyle {\boldsymbol {y}}=(y_{1},y_{2},\ldots ,y_{n})^{\textrm {T}}}
、係数
a
=
(
a
1
,
a
2
,
…
,
a
m
)
T
{\displaystyle {\boldsymbol {a}}=(a_{1},a_{2},\ldots ,a_{m})^{\textrm {T}}}
である。
これの最小解
a
{\displaystyle {\boldsymbol {a}}}
は、
[
G
y
]
T
[
G
y
]
=
R
~
T
R
~
{\displaystyle {\begin{bmatrix}G&{\boldsymbol {y}}\end{bmatrix}}^{\textrm {T}}{\begin{bmatrix}G&{\boldsymbol {y}}\end{bmatrix}}={\tilde {R}}^{\textrm {T}}{\tilde {R}}}
を満たす上三角行列
R
~
{\displaystyle {\tilde {R}}}
の計算を経て、解
a
{\displaystyle {\boldsymbol {a}}}
を得ることができ、全体の計算量に無駄が少ない。下記の表式を用いると
R
~
=
[
R
Q
T
y
0
T
α
]
{\displaystyle {\tilde {R}}={\begin{bmatrix}R&Q^{\textrm {T}}{\boldsymbol {y}}\\{\boldsymbol {0}}^{\textrm {T}}&\alpha \end{bmatrix}}}
が得られ、
R
a
=
Q
T
y
{\displaystyle R{\boldsymbol {a}}=Q^{\textrm {T}}{\boldsymbol {y}}}
から係数解
a
{\displaystyle {\boldsymbol {a}}}
を求める[注釈 2] 。
また前節で述べたように J を
a
{\displaystyle {\boldsymbol {a}}}
のそれぞれの成分で偏微分して 0 と置いた
m
{\displaystyle m}
個の式(正規方程式)は行列を用いて、
G
T
G
a
=
G
T
y
{\displaystyle G^{\textrm {T}}G{\boldsymbol {a}}=G^{\textrm {T}}{\boldsymbol {y}}}
と表される。これを正規方程式 (英 : normal equation )と呼ぶ。この正規方程式を解けば係数解
a
{\displaystyle {\boldsymbol {a}}}
が求まる。
係数解
a
{\displaystyle {\boldsymbol {a}}}
の解法には以下のようないくつかの方法がある。
逆行列で正規方程式を解く
行列
G
T
G
{\displaystyle G^{T}G}
が正則行列 (つまりフルランク)である場合は、解
a
=
(
G
T
G
)
−
1
G
T
y
{\displaystyle {\boldsymbol {a}}=(G^{\textrm {T}}G)^{-1}G^{\textrm {T}}{\boldsymbol {y}}}
は一意に求まる。ただし
G
T
G
{\displaystyle G^{T}G}
の逆行列を明示的に求めることは通常は良い方法ではない。
計算量が小さい方法としてコレスキー分解 [4] (
G
T
G
=
R
T
R
{\displaystyle G^{\textrm {T}}G=R^{\textrm {T}}R}
、
R
{\displaystyle R}
は
m
×
m
{\displaystyle m\times m}
上三角行列 )による三角行列分解を経て、最終的に
R
a
=
R
−
T
G
T
y
{\displaystyle R{\boldsymbol {a}}=R^{-{\textrm {T}}}G^{\textrm {T}}{\boldsymbol {y}}}
を解けばよい。
数値的安定性 確保のためには、積
G
T
G
{\displaystyle G^{T}G}
を経ない三角行列分解が良い。すなわち以下と同じくQR分解 (直交分解)による
G
=
Q
R
{\displaystyle G=QR}
から、
R
a
=
Q
T
y
{\displaystyle R{\boldsymbol {a}}=Q^{\textrm {T}}{\boldsymbol {y}}}
を解く。
直交分解で正規方程式を解く
コレスキー分解の方法よりも計算量が大きいが、数値的に安定かつ汎用な方法として、QR分解 や特異値分解 (SVD)を用いる方法がある[4] 。これらの方法では計算の過程で積
G
T
G
{\displaystyle G^{T}G}
を必要としないため数値的安定性 が高い。また
G
T
G
{\displaystyle G^{T}G}
が正則行列でない(ランク落ちしている)場合は正規方程式の解が不定となるが、その場合でも、これらの手法では解 a のうちノルム が最も小さいものを求めることができる。特異値分解を用いる場合は、特異値のうち極めて小さい値を 0 とみなして計算することで数値計算上の大きな誤差の発生を防ぐことができる(英 : truncated SVD )[5] 。
擬似逆行列 を使う方法もあるが[6] 、計算効率が悪いため、特殊な場合(解析的な数式が必要な場合など)を除いてあまり用いられない。
多次元
想定される分布が媒介変数 t を用いて
(
x
,
y
)
=
(
f
(
t
)
,
g
(
t
)
)
{\displaystyle (x,y)=(f(t),g(t))}
の形(あるいは
f
,
g
{\displaystyle f,g}
は複数の媒介変数によって決まるとしても同様)であっても考察される。
すなわち、測定値
(
x
i
,
y
i
)
{\displaystyle (x_{i},y_{i})}
がパラメータ
t
i
{\displaystyle t_{i}}
に対する
(
f
(
t
i
)
,
g
(
t
i
)
)
{\displaystyle (f(t_{i}),g(t_{i}))}
を理論値として近似されているものと考えるのである。
この場合、各点の理論値
(
f
(
t
i
)
,
g
(
t
i
)
)
{\displaystyle (f(t_{i}),g(t_{i}))}
と測定値
(
x
i
,
y
i
)
{\displaystyle (x_{i},y_{i})}
の間に生じる残差は
(
x
i
−
f
(
t
i
)
)
2
+
(
y
i
−
g
(
t
i
)
)
2
{\displaystyle {\sqrt {(x_{i}-f(t_{i}))^{2}+(y_{i}-g(t_{i}))^{2}}}}
である。故に、残差平方和 は
∑
i
=
1
n
{
(
x
i
−
f
(
t
i
)
)
2
+
(
y
i
−
g
(
t
i
)
)
2
}
{\displaystyle \sum _{i=1}^{n}\left\{(x_{i}-f(t_{i}))^{2}+(y_{i}-g(t_{i}))^{2}\right\}}
となるから、この値が最小であるように、
f
,
g
{\displaystyle f,g}
を決定するのである。
このように、
n
{\displaystyle n}
組の
(
x
,
y
)
{\displaystyle (x,y)}
の測定値
(
x
i
,
y
i
)
(
i
=
1
,
2
,
.
.
.
,
n
)
{\displaystyle (x_{i},y_{i})(i=1,2,...,n)}
を
n
{\displaystyle n}
組の
(
x
1
,
x
2
,
.
.
.
,
x
m
)
{\displaystyle (x_{1},x_{2},...,x_{m})}
の測定値
(
x
1
i
,
x
2
i
,
.
.
.
,
x
m
i
)
(
i
=
1
,
2
,
.
.
.
,
n
)
{\displaystyle (x_{1i},x_{2i},...,x_{mi})(i=1,2,...,n)}
に拡張したものも考察することができる。
異常値の除去
前提条件の節で述べたように、測定データを最小二乗法によって近似する場合、外れ値または異常値が含まれていると極端に近似の尤もらしさが低下することがある。また、様々な要因によって異常値を含む測定はしばしば得られるものである。
誤差が正規分布から極端に外れた異常値を取り除くための方法として修正トンプソン-τ法が用いられる[14] 。