最小二乗法

データセットを4次関数で最小二乗近似した例

最小二乗法(さいしょうにじょうほう、さいしょうじじょうほう；最小自乗法とも書く、英: least squares method)は、測定で得られた数値の組を、適当なモデルから想定される1次関数、対数曲線など特定の関数を用いて近似するときに、想定する関数が測定値に対してよい近似となるように、残差の二乗和を最小とするような係数を決定する方法、あるいはそのような方法によって近似を行うことである。

歴史

1805年にアドリアン＝マリ・ルジャンドルが出版したのが初出である。しかし、1809年にカール・フリードリヒ・ガウスが出版した際に1795年には最小二乗法を考案済みだったと主張したことで、最小二乗法の発明者が誰であるかについては不明になっている。

計算の概要

前提条件

最小二乗法では測定データy はモデル関数f (x )と誤差εの和で

$y=f(x)+\varepsilon$

と表せるとする。物理現象の測定データには、誤差が含まれ、それは系統誤差と偶然誤差を含んでいる。この内、偶然誤差は、測定における信号経路の微視的現象に由来するならば、正規分布であると期待されることが多い。また、社会調査などの誤差理由の特定が困難な場合でも誤差が正規分布になると期待する考え方もある。

誤差が正規分布に従わない場合、最小二乗法によって得られたモデル関数はもっともらしくないことに注意する必要がある。偶然誤差が正規分布していない場合、系統誤差が無視できない位大きくそれをモデル関数に含めていない場合、測定データに正規分布から大きく外れた外れ値を含む場合などが該当する。

上記を含め、最小二乗法の理論的基盤には次のような前提が設けられている。

• 測定値の誤差には偏りがない。すなわち誤差の平均値は 0 である。
• 測定値の誤差の分散は既知である。ただし測定データごとに異なる値でも良い。
• 各測定は互いに独立であり、誤差の共分散は 0 である。
• 誤差は正規分布する。
• m 個のパラメータ（フィッティングパラメータ）を含むモデル関数f が知られていて、測定量の真の値を近似誤差なく再現することのできるパラメータが存在する。

基礎的な考え方

話を簡単にするため、測定値は x, y の二次元の平面に分布するものとし、想定される分布（モデル関数）が y = f(x) の形である場合を述べる。想定している関数 f は、既知の関数 g(x) の線型結合で表されていると仮定する。すなわち、

$f(x)=\sum _{k=1}^{m}a_{k}g_{k}(x)$

例えば、g_k(x)=x^k-1 は、多項式近似であり、特に m=2 の時は $f(x)=a_{1}+a_{2}x$ という直線による近似（線形回帰）になる。

今、測定で得られた、次のような数値の組の集合があるとする。

$\{(x_{1},y_{1}),\ (x_{2},y_{2}),\ \ldots ,\ (x_{n},y_{n})\}$

これら (x, y) の分布が、y = f(x) というモデル関数に従うと仮定した時、想定される理論値は (x₁, f(x₁)), (x₂, f(x₂)), ..., (x_n, f(x_n)) ということになり、実際の測定値との残差は、各 i につき |y_i - f(x_i)| ということになる。 この残差の大きさは、xy-平面上での (x_i, y_i) と (x_i, f(x_i)) との距離でもある。

ここで、理論値からの誤差の分散の推定値は残差の平方和

$J=\sum _{i=1}^{n}(y_{i}-f(x_{i}))^{2}$

で与えられるから、J が最小になるように想定分布 f を(すなわち a_kを)、定めればよいということになる。

それには、上式は a_k を変数とする関数と見なすことができるので、J を a_k について偏微分したものをゼロと置く。こうして得られた m 個の連立方程式(正規方程式)を解き、a_k を決定すればよい。

一次方程式の場合

さらに簡単な例として、モデル関数を1次関数とし、

$f(x)=ax+b\,$

とおくと、a, b は次式で求められる。

$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={\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 は、

$f(x)=(g_{1}(x),g_{2}(x),\ldots ,g_{m}(x))(a_{1},a_{2},\ldots ,a_{m})^{\textrm {T}}$

と表すことができる。上付き添字 T は転置行列を表す。最小にすべき関数 J は

${\begin{aligned}J({\boldsymbol {a}})&=(G{\boldsymbol {a}}-{\boldsymbol {y}})^{\textrm {T}}\,(G{\boldsymbol {a}}-{\boldsymbol {y}})\\&=({\begin{bmatrix}G&{\boldsymbol {y}}\end{bmatrix}}{\begin{bmatrix}{\boldsymbol {a}}\\-1\end{bmatrix}})^{\textrm {T}}\,({\begin{bmatrix}G&{\boldsymbol {y}}\end{bmatrix}}{\begin{bmatrix}{\boldsymbol {a}}\\-1\end{bmatrix}})\end{aligned}}$

と表される。ここにG は、$G_{ij}=g_{j}(x_{i})$ なる成分を持つn×m行列、${\boldsymbol {y}}=(y_{1},y_{2},\ldots ,y_{n})^{\textrm {T}}$、係数${\boldsymbol {a}}=(a_{1},a_{2},\ldots ,a_{m})^{\textrm {T}}$ である。

これの最小解${\boldsymbol {a}}$は、${\begin{bmatrix}G&{\boldsymbol {y}}\end{bmatrix}}^{\textrm {T}}{\begin{bmatrix}G&{\boldsymbol {y}}\end{bmatrix}}={\tilde {R}}^{\textrm {T}}{\tilde {R}}$を満たす上三角行列${\tilde {R}}$の計算を経て、解${\boldsymbol {a}}$を得ることができ、全体の計算量に無駄が少ない。下記の表式を用いると${\tilde {R}}={\begin{bmatrix}R&Q^{\textrm {T}}{\boldsymbol {y}}\\{\boldsymbol {0}}^{\textrm {T}}&\alpha \end{bmatrix}}$が得られ、$R{\boldsymbol {a}}=Q^{\textrm {T}}{\boldsymbol {y}}$から係数解${\boldsymbol {a}}$を求める。

また前節で述べたように J を${\boldsymbol {a}}$のそれぞれの成分で偏微分してゼロと置いた m 個の式（正規方程式）は行列を用いて、

$G^{\textrm {T}}G{\boldsymbol {a}}=G^{\textrm {T}}{\boldsymbol {y}}$

と表される。これを正規方程式 (normal equation) と呼ぶ。この正規方程式を解けば係数解${\boldsymbol {a}}$が求まる。

係数解${\boldsymbol {a}}$の解法には以下のようないくつかの方法がある。

• 逆行列で正規方程式を解く

  行列 G^T G が正則行列（つまりフルランク）である場合は、解${\boldsymbol {a}}=(G^{\textrm {T}}G)^{-1}G^{\textrm {T}}{\boldsymbol {y}}$は一意に求まる。ただしG^T G の逆行列を明示的に求めることは通常は良い方法ではない。

計算量が小さい方法としてコレスキー分解（$G^{\textrm {T}}G=R^{\textrm {T}}R$、$R$はm×m上三角行列）による三角行列分解を経て、最終的に$R{\boldsymbol {a}}=R^{-{\textrm {T}}}G^{\textrm {T}}{\boldsymbol {y}}$を解けばよい。

数値的安定性確保のためには、積 G^T G を経ない三角行列分解が良い。すなわち以下と同じくQR分解（直交分解）による$G=QR$から、$R{\boldsymbol {a}}=Q^{\textrm {T}}{\boldsymbol {y}}$を解く。

• 直交分解で正規方程式を解く

  コレスキー分解の方法よりも計算量が大きいが、数値的に安定かつ汎用な方法として、QR分解や特異値分解 (SVD) を用いる方法がある。これらの方法では計算の過程で積 G^T G を必要としないため数値的安定性が高い。また G^T G が正則行列でない（ランク落ちしている）場合は正規方程式の解が不定となるが、その場合でも、これらの手法では解 a のうちノルムが最も小さいものを求めることができる。特異値分解を用いる場合は、特異値のうち極めて小さい値を0とみなして計算することで数値計算上の大きな誤差の発生を防ぐことができる (truncated SVD)。

• 擬似逆行列を使う方法もあるが、計算効率が悪いため、特殊な場合（解析的な数式が必要な場合など）を除いてあまり用いられない。

拡張

多次元

想定される分布が媒介変数 t を用いて (x, y) = (f(t), g(t)) の形(あるいは f, g は複数の媒介変数によって決まるとしても同様)であっても考察される。

すなわち、測定値 (x_i, y_i) がパラメータ t_i に対する (f(t_i), g(t_i)) を理論値として近似されているものと考えるのである。

この場合、各点の理論値 (f(t_i), g(t_i)) と測定値 (x_i, y_i) の間に生じる残差は

${\sqrt {(x_{i}-f(t_{i}))^{2}+(y_{i}-g(t_{i}))^{2}}}$

である。故に、残差平方和は

$\sum _{i=1}^{n}\left\{(x_{i}-f(t_{i}))^{2}+(y_{i}-g(t_{i}))^{2}\right\}$

となるから、この値が最小であるように、f, g を決定するのである。

このように、n 組の (x , y ) の測定値 (x_i , y_i ) (i = 1, 2, ... , n ) をn 組の (x₁ , x₂ , ... , x_m ) の測定値 (x_1i , x_2i , ... , x_mi ) (i = 1, 2, ... , n ) に拡張したものも考察することができる。

測定の誤差が既知の場合

n 回の測定における誤差があらかじめ分かっている場合を考える。異なる測定方法で測定した複数のデータ列を結合する場合などでは、測定ごとに誤差が異なることはしばしばある。誤差が正規分布していると考え、その標準偏差 $\sigma _{i}(i=1,2,\ldots ,n)$ で、誤差の大きさを表す。すると、誤差が大きい測定より、誤差が小さい測定の結果により重みをつけて近似関数を与えるべきであるから、

$J'=\sum _{i=1}^{n}{\frac {(y_{i}-f(x_{i}))^{2}}{\sigma _{i}^{2}}}$

を、最小にするように f を定める方がより正確な近似を与える。

毎回の測定が独立ならば、測定値の尤度は exp(-J') に比例する。そこで、上記の J' を最小にする f は、最尤推定値であるとも解釈できる。また、J' は自由度 n-m のカイ二乗分布に従うので、それを用いてモデル f の妥当性を検定することもできる。

毎回の測定誤差が同じ場合、J' を最小にするのは J を最小にするのと同じ意味になる。

非線形最小二乗法

詳細は「非線形最小二乗法」を参照

もし、f が、a_k の線型結合で表されないときは、正規方程式を用いた解法は使えず、反復解法を用いて数値的に a_k の近似値を求める必要がある。例えば、ガウス・ニュートン法やレーベンバーグ・マーカート法が用いられる。とくにLevenberg-Marquardt法は多くの多次元非線形関数でパラメータを発散させずに効率よく収束させる（探索する）方法として知られている。

異常値の除去

前提条件の節で述べたように、測定データを最小二乗法によって近似する場合、外れ値または異常値が含まれていると極端に近似の尤もらしさが低下することがある。また、様々な要因によって異常値を含む測定はしばしば得られるものである。

誤差が正規分布から極端に外れた異常値を取り除くための方法として修正トンプソン－τ法が用いられる。

関連項目

• 総最小二乗法 (Total least squares)
• ガウス＝マルコフの定理
• 曲線あてはめ