データセットを4次関数で最小二乗近似した例 最小二乗法 （さいしょうにじょうほう、さいしょうじじょうほう；最小自乗法 とも書く、英 : least squares method ）は、誤差を伴う測定値の処理において、その誤差の二乗の和を最小にするようにし、最も確からしい関係式を求める方法である。測定で得られた数値の組を、適当なモデルから想定される1次関数 、対数 曲線など特定の関数 を用いて近似 するときに、想定する関数が測定値に対してよい近似となるように、残差平方和 を最小とするような係数を決定する方法[ 1] [ 2] [ 3] [ 1] [ 2] [ 3]
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}}
m
{\displaystyle m}
a
k
{\displaystyle a_{k}}
さらに簡単な例として、モデル関数を1次関数とし、
f
(
x
)
=
a
x
+
b
{\displaystyle f(x)=ax+b\,}
とおくと、
a
,
b
{\displaystyle a,b}
a
=
n
∑
k
=
1
n
x
k
y
k
−
∑
k
=
1
n
x
k
∑
k
=
1
n
y
k
n
∑
k
=
1
n
x
k
2
−
(
∑
k
=
1
n
x
k
)
2
{\displaystyle a={\frac {n\textstyle \sum \limits _{k=1}^{n}x_{k}y_{k}-\sum \limits _{k=1}^{n}x_{k}\sum \limits _{k=1}^{n}y_{k}}{n\textstyle \sum \limits _{k=1}^{n}{x_{k}}^{2}-\left(\sum \limits _{k=1}^{n}x_{k}\right)^{2}}}}
b
=
∑
k
=
1
n
x
k
2
∑
k
=
1
n
y
k
−
∑
k
=
1
n
x
k
y
k
∑
k
=
1
n
x
k
n
∑
k
=
1
n
x
k
2
−
(
∑
k
=
1
n
x
k
)
2
{\displaystyle b={\frac {\textstyle \sum \limits _{k=1}^{n}{x_{k}}^{2}\sum \limits _{k=1}^{n}y_{k}-\sum \limits _{k=1}^{n}x_{k}y_{k}\sum \limits _{k=1}^{n}x_{k}}{n\textstyle \sum \limits _{k=1}^{n}{x_{k}}^{2}-\left(\sum \limits _{k=1}^{n}x_{k}\right)^{2}}}}
当てはめたい関数
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}}}
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)}
n
{\displaystyle n}
誤差 があらかじめ分かっている場合を考える。異なる測定方法で測定した複数のデータ列を結合する場合などでは、測定ごとに誤差が異なることはしばしばある。誤差が正規分布していると考え、その標準偏差
σ
i
(
i
=
1
,
2
,
…
,
n
)
{\displaystyle \sigma _{i}(i=1,2,\ldots ,n)}
J
′
=
∑
i
=
1
n
(
y
i
−
f
(
x
i
)
)
2
σ
i
2
{\displaystyle J'=\sum _{i=1}^{n}{\frac {(y_{i}-f(x_{i}))^{2}}{\sigma _{i}^{2}}}}
を、最小にするように f を定める方がより正確な近似を与える。
毎回の測定が独立 ならば、測定値の尤度 は
e
x
p
(
−
J
′
)
{\displaystyle exp(-J')}
J
′
{\displaystyle J'}
f
{\displaystyle f}
最尤推定値 であるとも解釈できる。また、
J
′
{\displaystyle J'}
自由度
n
−
m
{\displaystyle n-m}
カイ二乗分布 [ 7]
f
{\displaystyle f}
毎回の測定誤差が同じ場合、
J
′
{\displaystyle J'}
J
{\displaystyle J}
もし、
f
{\displaystyle f}
a
k
{\displaystyle a_{k}}
反復解法 を用いて数値的に
a
k
{\displaystyle a_{k}}
ガウス・ニュートン法 [ 8] レーベンバーグ・マーカート法 [ 9] [ 10] [ 11] [ 12] [ 13] 英 : Levenberg-Marquardt Method ）は多くの多次元非線形関数でパラメータを発散させずに効率よく収束させる（探索する）方法として知られている[ 10] [ 11] [ 12]
前提条件の節で述べたように、測定データを最小二乗法によって近似する場合、外れ値または異常値が含まれていると極端に近似の尤もらしさが低下することがある。また、様々な要因によって異常値を含む測定はしばしば得られるものである。
誤差が正規分布から極端に外れた異常値を取り除くための方法として修正トンプソン-τ法が用いられる[ 14]
ウィキメディア・コモンズには、
最小二乗法 に関連するカテゴリがあります。
ウィキブックスに
最小二乗法 関連の解説書・教科書があります。
^
m
{\displaystyle m}
^
α
2
=
min
J
(
a
)
{\displaystyle \alpha ^{2}=\min J({\boldsymbol {a}})}
G
T
G
{\displaystyle G^{\textrm {T}}G}
^ a b c 中川徹; 小柳義夫『最小二乗法による実験データ解析』東京大学出版会、1982年、30頁。ISBN 4-13-064067-4 。
^ a b c Lawson, Charles L.; Hanson, Richard J. (1995). Solving Least Squares Problems doi :10.1137/1.9781611971217 . https://epubs.siam.org/doi/10.1137/1.9781611971217 (
^ a b c Björck, Åke (1996). Numerical Methods for Least Squares Problems doi :10.1137/1.9781611971484 . https://epubs.siam.org/doi/abs/10.1137/1.9781611971484
^ a b 山本哲朗『数値解析入門』（増訂版）サイエンス社 〈サイエンスライブラリ 現代数学への入門 14〉、2003年6月。ISBN 4-7819-1038-6 。
^ Hansen, Per Christian (1987). “The truncated SVD as a method for regularization” . BIT Numerical Mathematics (Springer) 27 : 534-553. doi :10.1007/BF01937276 . https://doi.org/10.1007/BF01937276 . ( ^ 安川章. (2017). 科学実験/画像変換の近似計算に便利な 「疑似逆行列」 入門 できる人が使っている最小二乗法の一発フィット. インターフェース= Interface, 43(8), 142-146.
^ Weisstein, Eric W. "Chi-Squared Distribution." From MathWorld--A Wolfram Web Resource. mathworld
^ MAGREÑÁN, Alberto; Argyros, Ioannis (2018). A contemporary study of iterative methods: convergence, dynamics and applications . https://books.google.co.jp/books?hl=ja&lr=lang_ja ^
Weisstein, Eric W. "Levenberg-Marquardt Method." From MathWorld--A Wolfram Web Resource. mathworld
^ a b Moré Jorge J. (1978). “The Levenberg-Marquardt algorithm : Implementation and theory” . Lecture Notes in Mathematics 630 : 105-116. doi :10.1007/BFb0067700 . https://link.springer.com/chapter/10.1007/BFb0067700 . (
^ a b Yu, H., & Wilamowski, B. M. (2011). Levenberg-marquardt training. Industrial electronics handbook, 5(12), 1.
^ a b Ranganathan, Ananth (2004). “The levenberg-marquardt algorithm” (PDF). Tutoral on LM algorithm 11 (1): 101-110. https://sites.cs.ucsb.edu/~yfwang/courses/cs290i_mvg/pdf/LMA.pdf .
^ 山下信雄, 福島雅夫「Levenberg-Marquardt法の局所収束性について (最適化の数理科学) 」『数理解析研究所講究録』第1174巻、京都大学数理解析研究所、2000年10月、161-168頁、CRID 1050001201691367552 、hdl :2433/64462 ISSN 1880-2818 。 ^ Michele Rienzner (2020). Find Outliers with Thompson Tau (www MATLAB Central File Exchange. Retrieved May 17, 2020.
^ Van Huffel, S., & Vandewalle, J. (1991). The total least squares problem: computational aspects and analysis (Vol. 9). SIAM.
^ Golub, Gene H; Van Loan, Charles F (1980). “An analysis of the total least squares problem” . SIAM journal on numerical analysis (SIAM) 17 (6): 883-893. doi :10.1137/0717073 . https://doi.org/10.1137/0717073 . ^ Drygas, H. (2012). The coordinate-free approach to Gauss-Markov estimation (Vol. 40). Springer Science & Business Media.
^ Motulsky, H., & Christopoulos, A. (2004). Fitting models to biological data using linear and nonlinear regression: a practical guide to curve fitting. Oxford University Press.
要購読契約)
