最小二乗法の計算（理論）

POINT

• 最小二乗法の計算を解説．

最小二乗法の計算について紹介します．微分法による極値問題の一例としても良い題材です．

この記事では細かい部分を詰めることはせず，ざっくりとした計算の流れを整理することを目的とします．

【関連記事】

• 線形回帰 - Notes_JP：より詳細な話をまとめています．
• Excelで最小二乗法（近似曲線の求め方） - Notes_JP：「使う」ことをメインにした話です．Excel関数を使う方法です．

1次関数（直線）

$N$組のデータ点$(x_1,y_1),...,(x_N,y_N)$を最小二乗法で直線フィッティングすることを考えます．つまり，

\begin{aligned}
S(a,b)=\sum_{i=1}^N\bigl[y_i - (\underset{\text{直線の式}}{\underline{ax_i+b}} )\bigr]^2
\end{aligned}

を最小にするような$a,b$の値を求める問題を考えます．

---

【方針】
求める$a,b$は，$\displaystyle \frac{\partial S}{\partial a}(a,b)=0$, $\displaystyle \frac{\partial S}{\partial b}(a,b)=0$を満たします．したがって，この2式を連立して解くことで，$a,b$を計算できます．

【解】
偏微分は，合成関数の微分法を使うことで

\begin{aligned}
\frac{\partial S}{\partial a}(a,b)
&=\sum_{i=1}^N 2\left[y_i - (ax_i+b)\right]\cdot (-x_i) \\
\frac{\partial S}{\partial b}(a,b)
&=\sum_{i=1}^N 2\left[y_i - (ax_i+b)\right]
\end{aligned}

と計算できます（補足：*1）．したがって，$\displaystyle \frac{\partial S}{\partial a}(a,b)=0$, $\displaystyle \frac{\partial S}{\partial b}(a,b)=0$はそれぞれ

\begin{aligned}
\frac{\partial S}{\partial a}(a,b)&=0:\\
&\sum_{i=1}^N \left(x_i \right)^2 \cdot a +\sum_{i=1}^N x_i \cdot b = \sum_{i=1}^N x_i y_i\\
\frac{\partial S}{\partial b}(a,b)&=0: \\
&\sum_{i=1}^N x_i\cdot a+Nb=\sum_{i=1}^N y_i
\end{aligned}

となります．

ここで，それぞれ$\{x_i\}_i$, $\{y_i\}_i$を当確率で値にとる確率変数$X$, $Y$を導入しましょう（*2）．すると

• $X$の平均値：$\displaystyle\langle X\rangle=\frac{1}{N}\sum_{i=1}^N x_i$
• $Y$の平均値：$\displaystyle\langle Y\rangle=\frac{1}{N}\sum_{i=1}^N y_i$
• $X$の分散：
  
  \begin{aligned}
  \sigma_X^2
  &=\biggl\langle \bigl(X-\langle X\rangle\bigr)^2
  \biggr\rangle \\
  &=\langle X^2\rangle
  -\langle X\rangle^2
  \end{aligned}
• $X,Y$の共分散：
  
  \begin{aligned}
  \mathrm{Cov}\left(\mathrm{X,Y}\right)
  &=\biggl\langle \bigl(X-\langle X\rangle\bigr)
  \bigl(Y-\langle Y\rangle\bigr)
  \biggr\rangle \\
  &=\langle XY\rangle
  -\langle X\rangle\langle Y\rangle
  \end{aligned}

という量を考えられます．これらを使うと，得られた方程式は

\begin{aligned}
\frac{\partial S}{\partial a}(a,b)&=0:\\
& \left(\sigma_X^2 +\langle X\rangle^2 \right) a +\langle X\rangle b \\
&\qquad =\mathrm{Cov}\left(\mathrm{X,Y}\right) + \langle X\rangle\langle Y\rangle \\
\frac{\partial S}{\partial b}(a,b)&=0: \\
& \langle X\rangle a+b=\langle Y\rangle
\end{aligned}

と書き直すことができます．

これら2式を連立させて解くことで

\begin{aligned}
& y = ax + b \\
&
\begin{cases}
\, \displaystyle a =\frac{\mathrm{Cov}\left(\mathrm{X,Y}\right)}{\sigma_X^2} \\
\, \displaystyle b = -a\langle X\rangle+\langle Y\rangle
= \frac{\langle X^2\rangle \langle Y\rangle - \langle XY\rangle \langle X\rangle}{\sigma_X^2}
\end{cases}
\end{aligned}

が得られます．//

【注意】
求めた直線は

\begin{aligned}
y = a(x-\langle X\rangle) + \langle Y\rangle
\end{aligned}

という形をしています．したがって，この直線は平均値からなる点$(\langle X\rangle,\langle Y\rangle)$を通ります．

多変数の線形回帰

（上の場合を含むような）一般化した問題として，複数のパラメータについて線形な関数によるフィッティングを考えます．

フィッティング関数として

\begin{aligned}
f(x^1,…,x^m)
=\sum_{i=1}^m \theta_m x^m
\end{aligned}

を考えます（ここで，$x^i$は冪乗ではなく「単なる添字」を意味します）．この式は，ベクトル$\boldsymbol{\theta}={}^t(\theta_1,...,\theta_m)$, $\boldsymbol{x}={}^t(x^1,...,x^m)$を用いて

\begin{aligned}
f(\boldsymbol{x})
=\boldsymbol{x} \cdot \boldsymbol{\theta}
\end{aligned}

と表すことができます．

それでは$N$組のデータ点$(\boldsymbol{x}_1,y_1),…,(\boldsymbol{x}_N,y_N)$のフィッティングを考えましょう．$N$組の式

\begin{aligned}
y_i - \boldsymbol{x}_i \cdot \boldsymbol{\theta}
\end{aligned}

は行列$A=(a_{ij})=x^j_{\:i}$（$x^j_{\:i}$は$\boldsymbol{x}_i$の$j$成分）を用いて

\begin{aligned}
\boldsymbol{y} - A \boldsymbol{\theta}
\end{aligned}

と書き直せます．

したがって，二乗和は$\left\| \boldsymbol{y} - A \boldsymbol{\theta}\right\|^2$$={}^t(\boldsymbol{y} - A \boldsymbol{\theta})(\boldsymbol{y} - A \boldsymbol{\theta})$となります．ここでは，（データに相関がある場合にも扱えるよう）重み付きの二乗和

\begin{aligned}
S(\boldsymbol{\theta})=
{}^t(\boldsymbol{y} - A \boldsymbol{\theta}) W (\boldsymbol{y} - A \boldsymbol{\theta})
\end{aligned}

を最小化することを考えます（*3）．

1変数の場合と同じように$\displaystyle\frac{\partial S}{\partial \theta_i}(\boldsymbol{\theta})=0$$\,(i=1,...,m)$をみたす$\boldsymbol{\theta}$が求めるものです．これは簡単に計算することができ（*4）

\begin{aligned}
\frac{\partial S}{\partial \theta_i}(\boldsymbol{\theta})&=0: \\
&-2\:{}^tA W (\boldsymbol{y} - A \boldsymbol{\theta})=0
\end{aligned}

となります．したがって

\begin{aligned}
\boldsymbol{\theta}=\left({}^tA WA \right)^{-1} \left( {}^tA W\boldsymbol{y} \right)
\end{aligned}

が求める値です．

参考文献

データ・誤差解析の基礎

*1:$g(x)=x^2$, $h_i(a,b)=y_i - (ax_i+b)$により，$S(a,b)=\sum_i g\left(h_i(a,b)\right)$と表されるので，合成関数の微分法より

\begin{aligned} &\frac{\partial S}{\partial a}(a,b) \\ &=\sum_{i=1}^N \frac{\mathrm{d} g}{\mathrm{d} x}\left(h_i(a,b)\right) \frac{\partial h_i}{\partial a}(a,b) \\ &=\sum_{i=1}^N 2\left[y_i - (ax_i+b)\right]\cdot (-x_i) \end{aligned}

と計算できます．$\partial S/\partial b$についても同様です．

*2:あるいは，

\begin{aligned} & \langle X\rangle =\frac{1}{N}\sum_{i=1}^N x_i \\ & \Delta x_i =x_i - \langle X\rangle \\ & \langle f(\Delta X,\Delta Y) \rangle =\frac{1}{N}\sum_{i=1}^N f(\Delta x_i, \Delta y_i) \end{aligned}

という量を「定義」したと考えても良いです（$f$は関数）．この場合，最終的な結果は

\begin{aligned} \sigma_X^2 &\rightarrow \langle (\Delta X)^2\rangle \\ \mathrm{Cov}\left(\mathrm{X,Y}\right) &\rightarrow \langle \Delta X \Delta Y\rangle \\ \end{aligned}

と置き換えたものになります．

*3:$W=(w_{ij})$は重みを表す対称行列です．ベクトル$\boldsymbol{a}$に対して

\begin{aligned} {}^t\boldsymbol{a} W\boldsymbol{a} =\sum_{i,j} a_i w_{ij} a_j \end{aligned}

となります．

*4:関数$f(\boldsymbol{\theta})=[{}^t\boldsymbol{a}(\boldsymbol{\theta})] W [\boldsymbol{b}(\boldsymbol{\theta})]$を考えましょう．ベクトル・行列の成分で表せば

\begin{aligned} f(\boldsymbol{\theta})=\sum_{jk}a_j w_{jk}b_k \end{aligned}

なので，積の微分法より

\begin{aligned} \frac{\partial f}{\partial \theta_i} &=\sum_{jk}w_{jk}\left( \frac{\partial a_j}{\partial \theta_i}\ b_k +a_j \frac{\partial b_k}{\partial \theta_i} \right) \\ &={}^t \left( \partial \boldsymbol{a} / \partial \theta_i \right) W \boldsymbol{b} +{}^t \boldsymbol{a} W \left( \partial \boldsymbol{b} / \partial \theta_i \right) \end{aligned}

となります．特に，$\boldsymbol{a}=\boldsymbol{b}$の場合には

\begin{aligned} \frac{\partial f}{\partial \theta_i} &={}^t \left( \partial \boldsymbol{a} / \partial \theta_i \right) (W + {}^t W)\boldsymbol{a} \end{aligned}

と変形することができます（スカラー量は転置をとっても値を変えないので，${}^t\boldsymbol{a} W\boldsymbol{b}$$={}^t\left({}^t\boldsymbol{a} W\boldsymbol{b}\right)$$={}^t\boldsymbol{b} {}^tW\boldsymbol{a}$.）．