こんにちは、(株)日立製作所 Lumada Data Science Lab. の中川です。普段は人工知能を制御に適用する研究に従事しています。近年、機械学習が注目される中、機械学習理論および機械学習を使った技術開発環境は急速に進歩すると共に、多くの方がデータサイエンスに関わるようになってきました。すでにデータサイエンスに携わっている方や、これからデータサイエンスに関わってみようと思っている方の中で、理論の大切さをあらためて知りたいあるいは感じたいという方がいらっしゃいましたら、それをできるだけわかりやすく伝えられたら、という思いから基本的な内容で記事を書きます。今回は、線形重回帰の回帰係数値を様々な方法で推定し、各方法の特徴を調べてみます。
1. はじめに
線形重回帰の回帰係数の値を決定するとき、どのような手法で求めるのが良いのでしょうか。一般に、線形重回帰の回帰係数（$\hat{a}_1,\hat{a}_2,\cdots,\hat{a}_n$）の値は、（1）で示される評価関数$J$の値を最小化する値として決定します。
$$
\begin{array}
.&\,J=\sum_{i=1}^{m}\,\bigl(\,y(i)-\hat{y}(i)\,\bigl)^2\hspace{0.5cm}・・・(1)\\
&\hat{y}(i)=\hat{{\bf\theta}}^T\cdot{\bf x}(i)\\
&\hspace{1cm}=\hat{a}_1x_1(i)+\hat{a}_2x_2(i)+\cdots+\hat{a}_nx_n(i)+\hat{b}\hspace{0.5cm}・・・・(2)\\
&\hat{{\bf \theta}}=(\hat{a}_1,\hat{a}_2,\cdots,\hat{a}_n,\hat{b})^T\\
&{\bf x}(i)=(\,x_1(i),x_2(i),\cdots,x_n(i),1\,)^T\\
&\mspace{1cm}\\
&\hspace{1cm}x_j(i):入力\\
&\hspace{1cm}y(i):出力\\
&\hspace{1cm}\hat{y}(i):出力予測値\\
&\hspace{1cm}a_j:回帰係数\\
&\hspace{1cm}i:1,2,\cdots,m\\
\end{array}
$$
$J$の値を最小化する$\hat{{\bf\theta}}$は、最小二乗法で決めるのが一般的ですが、正規方程式と呼ばれる連立方程式を直接解いて求める方法と、${\bf x}(i)$と$y(i)$が得られる度に$\hat{{\bf \theta}}$を逐次更新する逐次最小二乗法があります。前者はオフライン最小二乗法、後者はオンライン最小二乗法と呼ばれることもあります。一方で、$J$は$\hat{\bf \theta}$の関数なので、$J$を最小化する$\hat{\bf \theta}$は、最適化法によっても求めることができます。今回は、最適化法として、最急降下法、ガウス・ニュートン法、PSO（Particle Swarm Optimization）を取り上げて、前述の２つの最小二乗法と合わせて５つの手法について考えてみます。
2. 各方法の紹介
2-1. 正規方程式を解析的に解く（オフライン最小二乗法）
前章で紹介しました５つの方法についてそれぞれ簡単に御紹介します。
（3）式は正規方程式と呼ばれ、（1）式において、$\partial{J}/\partial\hat{\bf \theta}=0$とすることで求められます。（3）式は連立一次方程式であり、これを解くことで、$\hat{\bf \theta}$が求まります。
$$
\begin{array}
\,&\hat{\bf \theta}=(X^TX)^{-1}X^T{\bf y}\hspace{0.5cm}・・・(3)\\
\mspace{1cm}\\
&X=\left(
\begin{array}{ccc}
\,{\bf x}(1)&,&1\,\\
\,{\bf x}(2)&,&1\,\\
\vdots&,&\vdots\\
\,{\bf x}(m)&,&1\,
\end{array}
\right)\\
\mspace{1cm}\\
&\,{\bf y}=(\,y(1),y(2),\cdots,y(m)\,)^T\\
\end{array}
$$
本方式は、入出力データを一括処理して（反復処理を行わず）、$\hat{\bf \theta}$を求めることから、オフライン最小二乗法あるいは一括最小二乗法とも呼ばれることがあります。
2-2. 逐次最小二乗法（オンライン最小二乗法）
（4）式は逐次最小二乗法を示しています。（4）式の反復処理を行うことで、$\hat{\bf \theta}$が求まります。
$$
\begin{array}
\,&\hat{\bf\theta}(k)=\hat{\bf \theta}(k-1)+\frac{P(k-1){\bf\phi}(k)}{1+{\bf\phi}(k)^T P(k-1){\bf\phi}(k)}{\bf e}(k)\hspace{0.5cm}・・・(4)\\
&\mspace{1cm}\\
&{\bf\phi}(k)=(\,{\bf x}(k),1\,)^T\\
&{\bf e}(k)={\bf y}(k)-{\bf\phi}(k)^T\hat{{\bf\theta}}(k-1)\\
&P(k)=P(k-1)-\frac{P(k-1){\bf\phi}(k){\bf\phi}(k)^T P(k-1)}{1+{\bf\phi}(k)^T P(k-1){\bf\phi}(k)}\\
&\mspace{1cm}\\
&\hspace{1cm}\hat{\bf\theta}(0):任意の値\\
&\hspace{1cm}k:1,2,\cdots,m\\
\end{array}
$$
（4）式は、（3）式を変形することで求まります。本方式は、入出力データが新たに得られる度に$\hat{\bf \theta}$を更新することから、オンライン最小二乗法とも呼ばれることがあります。
2-3. 最急降下法
（5）式は最急降下法を示しています。（5）式の反復処理を行うことで、$\hat{\bf \theta}$が求まります。
$$
\begin{array}
\,&\hat{\bf\theta}(k)=\hat{\bf\theta}(k-1)-\beta_{k}\nabla{J_{k}}\hspace{0.5cm}・・・(5)\\
&\mspace{1cm}\\
&\hspace{1cm}\beta_{k}:ステップ幅\\
&\hspace{1cm}k:1,2,\cdots\\
\end{array}
$$
本方式は、勾配法の一種であり、$-\nabla{J_{k-1}}$は$J$が最小となる方向を示しています。$\beta_{k-1}$は、更新量を調整する係数で、一般に、1次元の最適化により決定します。本手法は、評価関数の１次微分値を用いて最適化を行いますので、１次最適化と呼ばれる手法に属します（「１次元の最適化」と「１次最適化」は意味が異なりますので、注意しましょう）。
2-4. ガウス・ニュートン法
（6）式はガウス・ニュートン法を示しています。（6）式の反復処理を行うことで、$\hat{\bf \theta}$が求まります。
$$
\begin{array}
\,&\hat{\bf\theta}(k)=\hat{\bf\theta}(k-1)-(K^TK)^{-1}K^T{\bf e}(k)\hspace{0.5cm}・・・(6)\\
\mspace{1cm}\\
&K=\left(\begin{array}{cccc}
\,\frac{\partial e_1}{\partial \hat{a}_1},&\frac{\partial e_1}{\partial \hat{a}_2 },&\cdots,&\frac{\partial e_1}{\partial \hat{a}_n},&\frac{\partial e_1}{\partial \hat{b}}\,\\
\,\frac{\partial e_2}{\partial \hat{a}_1},&\frac{\partial e_2}{\partial \hat{a}_2 },&\cdots,&\frac{\partial e_2}{\partial \hat{a}_n},&\frac{\partial e_2}{\partial \hat{b}}\,\\
\vdots&\vdots&\ddots&\vdots&\vdots\\
\,\frac{\partial e_m}{\partial \hat{a}_1},&\frac{\partial e_m}{\partial \hat{a}_2 },&\cdots,&\frac{\partial e_m}{\partial \hat{a}_n},&\frac{\partial e_m}{\partial \hat{b}}\,\\
\end{array}\right)\\
\mspace{1cm}\\
&{\bf e}(k)=(\,e_1(k),e_2(k),\cdots,e_m(k)\,)^T\\
&e_i=y(i)-\hat{y}_k(i)\\
&\hat{y}_k(i)=\hat{{\bf\theta}}_k^T\cdot{\bf x}(i)\\
&\hat{{\bf\theta}}_k=(\,\hat{a}_1(k),\hat{a}_2(k),\cdots,\hat{a}_n(k),\hat{b}(k)\,)^T\\
&k:1,2,\cdots\\
\end{array}
$$
ガウス・ニュートン法は、ニュートン法を改良したもので、評価関数が２次形式で表されるときに適用できます。ニュートン法では、ヘッセ行列という要素が2階微分値である行列を用いますが、ガウス・ニュートン法では近似により要素が１階微分値となります。参考までに、（7）式にニュートン法を示します。
$$
\begin{array}
\,&\hat{\bf\theta}(k)=\hat{\bf\theta}(k-1)-H^{-1}\nabla{J_{k}}\hspace{0.5cm}・・・(7)\\
\mspace{1cm}\\
&H=\left( \begin{array}{cccc}\,\frac{\partial^2 J}{\partial \hat{a}_1^2},&\frac{\partial^2 J}{\partial \hat{a}_1 \partial \hat{a}_2},&\cdots,&\frac{\partial^2 J}{\partial \hat{a}_1 \partial \hat{a}_n},&\frac{\partial^2 J}{\partial \hat{a}_1 \partial \hat{b}}\,\\
\frac{\partial^2 J}{\partial \hat{a}_2 \partial \hat{a}_1},&\frac{\partial^2 J}{\partial \hat{a}_2^2},&\cdots,&\frac{\partial^2 J}{\partial \hat{a}_2 \partial a_n},&\frac{\partial^2 J}{\partial \hat{a}_2 \partial \hat{b}} \,\\
\vdots&\vdots&\ddots&\vdots&\vdots\\
\frac{\partial^2 J}{\partial \hat{a}_n \partial \hat{a}_1},&\frac{\partial^2 J}{\partial \hat{a}_n \partial \hat{a}_2},&\cdots,&\frac{\partial^2 J}{\partial \hat{a}_n^2},&\frac{\partial^2 J}{\partial \hat{a}_n \partial \hat{b}} \,\\
\end{array} \right)\\
\end{array}
$$
ニュートン法は、評価関数を探索点の回りで２次関数近似し、その極値（最小値もしくは最大値）に探索点を更新することで大域最適解に収束させるもので、２次最適化と呼ばれる手法に属します。ニュートン法の近似手法であるガウス・ニュートン法も同じく２次最適化に属します。
2-5. ＰＳＯ
（8）式はPSOのアルゴリズムを示しています。（8）式の反復処理を行うことで、$\hat{\bf \theta}$が求まります。
$$
\begin{array}
\,&\hat{\bf\theta}_p(k)=\hat{\bf\theta}_p(k-1)+{\bf v}_p(k)\hspace{0.5cm}・・・(8)\\
\mspace{1cm}\\
&\hat{{\bf \theta}}_p(k)=(\,\hat{a}_{1,p}(k),\hat{a}_{2,p}(k),\cdots,\hat{a}_{n,p}(k),\hat{b}_p(k)\,)^T\\
&{\bf v}_p(k)=(\,v_{1,p}(k),v_{2,p}(k),\cdots,v_{n,p}(k),v_{b,p}(k) \,)^T\\
&\hat{a}_{j,p}:\hat{a}_j最適化時のp番目の粒子\\
&\hat{b}_{p}:\hat{b}最適化時のp番目の粒子\\
&v_{j,p}(k)=w{\cdot}v_{j,p}(k-1)+C_1{\cdot}rand_{1,j,p}\cdot\bigl(\,pbest_{j,p}(k-1)-\hat{a}_{j,p}(k-1)\,\bigl)\hspace{0.3cm}\\
&\hspace{3cm}+\,C_2{\cdot}rand_{2,j,p}\cdot\bigl(\,gbest_{j,p}(k-1)-\hat{a}_{j,p}(k-1)\,\bigl)\hspace{0.5cm}・・・(9)\\
&v_{b,p}(k)=w{\cdot}v_{b,p}(k-1)+C_1{\cdot}rand_{1,b,p}\cdot\bigl(\,pbest_{b,p}(k-1)-\hat{b}_p(k-1)\,\bigl)\hspace{0.6cm}\\
&\hspace{3cm}+\,C_2{\cdot}rand_{2,b,p}\cdot\bigl(\,gbest_{b,p}(k-1)-\hat{b}_p(k-1)\,\bigl)\hspace{0.5cm}・・・(10)\\
&w,C_1,C_2:定数\\
&rand_{1,j,p},\,rand_{2,j,p},\,rand_{1,b,p},\,rand_{2,b,p}:0～1の一様乱数\\
&p:粒子番号\\
&k:1,2,\cdots\\
\end{array}
$$
PSOは、メタヒューリスティクスと呼ばれる非凸（多峰性）の解空間を効率的に探索するための手法の一つで、鳥や魚などの群れのような集団での探索行動に基づいた手法です。（8）式にありますように、複数の粒子が反復演算毎に調整される速度${\bf v}_p$に基づいて、解空間を動き回ることで最適解を探索します。（9）式および（10）式の第一項は「慣性」、第二項は「自己認識」、第三項は「社会認識」とも言われ、各項のベクトル和として、各粒子の探索する方向が更新されます。
3. 最適化結果
3-1. 評価条件
線形重回帰とその各パラメータ値を下記のように定義します。
$$
\begin{array}
\,&y=\theta^T\cdot{\bf x}\hspace{0.5cm}・・・(11)\\
&\;\;\,=a_1x_1+a_2x_2+a_3x_3+a_4x_4+b\\
&{\bf \theta}=(\,a_1,a_2,a_3,a_4,b\,)^T\\
&{\bf x}=(\,x_1,x_2,x_3,x_4,1\,)^T\\
&(\,a_1,a_2,a_3,a_4,b\,)=(\,2.6,-7.2,-1.6,4.8,2\,)\\
\end{array}
$$
（11）式で定義される線形重回帰式を元に、（12）式で定義される$y(i)$と$x(i)$を用いて、前章でご紹介しました５つの推定方法で$\theta$を推定します。
$$
\begin{array}
\,&y(i)={\bf\theta}^T\cdot{\bf x}(i)+d(i)\hspace{0.5cm}・・・(12)\\
&d(i):ノイズ（平均が0、分散が4の正規分布）\\
&i=1,2,・・・,m\\
\end{array}
$$
3-2. 正規方程式を解析的に解く（オフライン最小二乗法）
（3）式で示されます正規方程式を解いて得られた$\hat{\theta}$の値は、下記となりました。推定値$\hat{\theta}$は、小数点以下５桁目を四捨五入し、小数点以下４桁までの値を評価しました。
$$
\begin{array}
\,&\hat{a}_1=2.6002\\
\,&\hat{a}_2=-7.1995\\
\,&\hat{a}_3=-1.6005\\
\,&\hat{a}_4=4.8002\\
\,&\hat{b}=2.0019\\
\end{array}
$$
真値（$a_1,a_2,a_3,a_4,b$）との誤差は、ノイズ（$d$）によるものです。
3-3. 逐次最小二乗法（オンライン最小二乗法）
逐次最小二乗法を用いたときの反復演算回数に対する推定値のプロフィールを下図に示します。なお、$\hat{\theta}$の初期値$\hat{\theta}_0$は、$\hat{\theta}_0=[0, 0, 0, 0, 0]$としました。
反復回数500回で、$\hat{b}$以外は、ほぼ収束していることがわかります。反復回数500回における$\hat{\theta}$の値は、下記です。
$$
\begin{array}
\,&\hat{a}_1=2.5996\\
\,&\hat{a}_2=-7.2060\\
\,&\hat{a}_3=-1.6031\\
\,&\hat{a}_4=4.8085\\
\,&\hat{b}=2.1899\\
\end{array}
$$
正規方程式を解いて得た$\hat{\theta}$と比較しますと、誤差がありますが、さらに反復回数をm(=146410)まで増やしますと、得られた$\hat{\theta}$の値は、下記となりました。
$$
\begin{array}
\,&\hat{a}_1=2.6002\\
\,&\hat{a}_2=-7.1995\\
\,&\hat{a}_3=-1.6005\\
\,&\hat{a}_4=4.8002\\
\,&\hat{b}=2.0019\\
\end{array}
$$
反復回数がmのときは、正規方程式を解いて得た$\hat{\theta}$と同値となっています。正規方程式は、m組のデータを一括して扱って$\hat{\theta}$を与えますが、逐次最小二乗法では、データ一組ずつ与えて$\hat{\theta}$を更新していきます。逐次最小二乗法においてm組のデータを与えれば、逐次最小二乗法により求めた$\hat{\theta}$と正規方程式を解いて得た$\hat{\theta}$は、同値になっています。
3-4. 最急降下法
最急降下法を用いたときの反復演算回数に対する推定値のプロフィールを下図に示します。なお、$\hat{\theta}_0=[0, 0, 0, 0, 0]$です。（5）式における$\beta_k$は、ブレント法と呼ばれる一次元の最適化方法を適用しました。
反復回数１回で、$\hat{b}$以外は、ほぼ収束していることがわかります。反復回数１回における$\hat{\theta}$の値は、下記です。
$$
\begin{array}
\,&\hat{a}_1=2.6002\\
\,&\hat{a}_2=-7.1995\\
\,&\hat{a}_3=-1.6005\\
\,&\hat{a}_4=4.8002\\
\,&\hat{b}=1.2512\\
\end{array}
$$
反復回数２回では、$\hat{b}$も変化し、ほぼ収束していることがわかります。反復回数５回における$\hat{\theta}$の値は、下記です。
$$
\begin{array}
\,&\hat{a}_1=2.6002\\
\,&\hat{a}_2=-7.1995\\
\,&\hat{a}_3=-1.6005\\
\,&\hat{a}_4=4.8002\\
\,&\hat{b}=2.0019\\
\end{array}
$$
反復回数５回では$\hat{\theta}$は、正規方程式を解いて得た$\hat{\theta}$と同値となっています。
最急降下法では、ブレント法により、一次元の最適化を行いますので、反復処理が入れ子構造になっていると言えます。下図は、反復回数１回目におけるブレント法による$\beta_{k,j}$（ $j$：ブレント法における反復回数）の最適化の過程での$\hat{\theta}$のプロフィ－ルを示しています。（$\hat{a}_1,\hat{a}_2,\hat{a}_3,\hat{a}_4$）が$\hat{b}$に優先して最適化されていることがわかります。
下図は、反復回数２回目におけるブレント法による$\beta_{k,j}$の最適化の過程での$\hat{\theta}$のプロフィ－ルを示しています。上図で示しましたように反復回数１回目において、（$\hat{a}_1,\hat{a}_2,\hat{a}_3,\hat{a}_4$）が最適化されているので、反復回数２回目では、$\hat{b}$が主として最適化されていることがわかります。
これらの結果から、解空間は３次元空間で表せば、いわば、V字谷のような形状をしていると推測されます。（$\hat{a}_1,\hat{a}_2,\hat{a}_3,\hat{a}_4$）の現在値（ベクトル）と最適解（ベクトル）の差違は、V字谷の高さであり、$\hat{b}$の軸は、V字谷の奥行き方向の軸であると推測されます。したがって、まずは、谷底に向かって、すなわち（$\hat{a}_1,\hat{a}_2,\hat{a}_3,\hat{a}_4$）の最適化が進み、谷底に達した後、谷底の中でもより低い地点を目指して$\hat{b}$の最適化が進んだと考えられます。
下図は、$\beta_k$はブレント法を用いず一定値（$\beta_k=1$）として、最急降下法により最適化した場合の結果です。すなわち、反復回数１回あたりに進む距離は一定です。反復演算の当初は（$\hat{a}_1,\hat{a}_2,\hat{a}_3,\hat{a}_4$）の最適化が進み、その後、$\hat{b}$の最適化が進んでいることがわかります。特に、$\hat{b}$の最適化には、（$\hat{a}_1,\hat{a}_2,\hat{a}_3,\hat{a}_4$）の最適化より、反復演算を要しており、前述の谷底の傾斜が緩やかであることを示唆しています。
3-5. ガウス・ニュートン法
ガウス・ニュートン法を用いたときの反復演算回数に対する推定値のプロフィールを下図に示します。なお、$\hat{\theta}_0=[0, 0, 0, 0, 0]$です。反復回数１回で、すべての推定値が、ほぼ収束していることがわかります。
反復回数１回における$\hat{\theta}$の値は、下記です。
$$
\begin{array}
\,&\hat{a}_1=2.6002\\
\,&\hat{a}_2=-7.1995\\
\,&\hat{a}_3=-1.6005\\
\,&\hat{a}_4=4.8002\\
\,&\hat{b}=2.0019\\
\end{array}
$$
正規方程式を解いて得た$\hat{\theta}$と同値となっています。先述しましたように、ガウス・ニュートン法は、探索点のまわりで評価関数を２次関数で近似して、近似した２次関数の最小点に探索点を更新させます。（1）式で示される評価関数の場合、解空間は凸となるため、ガウス・ニュートン法は、優れた収束性を示します。
3-6. ＰＳＯ
PSOを用いたときの反復演算回数に対する推定値のプロフィールを下図に示します。なお、ハイパーパラメータ値は、下記としました。
$$
\begin{array}
\,&w=0.729\\
\,&C_1=1.4955\\
\,&C_2=1.4955\\
\,&粒子数：2000\\
\,&反復回数：500\\
\end{array}
$$
前述の手法と比較しますと、最適解に最短で向かうのではなく、ある程度ランダム性のある探索をしつつ、徐々に最適解に収束する動きになっていることがわかります。反復回数１００回で、すべての推定値が、ほぼ収束していることがわかります。図中にはありませんが、さらに反復回数を増やし、反復回数５００回（十分収束したと考えられる反復回数）のときの$\hat{\theta}$の値は、下記となりました。
$$
\begin{array}
\,&\hat{a}_1=2.6015\\
\,&\hat{a}_2=-7.1986\\
\,&\hat{a}_3=-1.6026\\
\,&\hat{a}_4=4.8001\\
\,&\hat{b}=2.0735\\
\end{array}
$$
正規方程式を解いて得た$\hat{\theta}$と比較しますと、やや差違のある結果となり、理論上の最適解には収束しませんでした。ただし、各ハイパーパラメータをチューニングすることで、最適解への収束性は高まるものと考えます。
4. まとめ
以上を簡単にまとめます。
＜ 正規方程式を解析的に解く ＞
- 反復計算を行わず、全データを一括して処理して推定値を求める。解析的に解いて
いるので、誤差二乗和を最小化する意味での推定値として精度はもっとも良い。
- 新たなデータが今後得られることはない、あるいは、得られても一括して計算し直
してもよい場合に有効。
- データ数が増えるにつれて、行列のサイズが大きくなり、計算負荷が増える。
＜ 逐次最小二乗法 ＞
- 正規方程式を反復演算により解く。データが得られる都度、反復演算により推定値
を更新するので、全てのデータに対して反復演算を行えば、正規方程式を解析的に
解く場合と同等の精度が得られる。
- 全てのデータが一括して得られず、かつデータが得られる度に逐次推定値を更新し
たい場合に有効。
- 反復演算一回あたりの計算負荷は比較的小さい。
＜ 最急降下法 ＞
- データの持つ情報は解空間として現れており、一次最適化（反復演算）により解空
間内の最適解を求める。正規方程式を解析的に解く場合とほぼ同等の精度が得られ
る。
- 新たなデータが今後得られることはない、あるいは、得られても一括して計算し直
してもよい場合に有効。
- 今回のケースでは解空間が凸となるので、初期値に依存せず最適解に収束する。同
様に今回のケースのように勾配が緩やかな次元があると収束が遅くなる。
- 計算負荷は比較的少ないが、一次元最適化のアルゴリズム（計算負荷）にも拠る。
＜ ガウス・ニュートン法 ＞
- データの持つ情報は解空間として現れており、二次最適化（反復演算）により解空
間内の最適解を求める。正規方程式を解析的に解く場合とほぼ同等の精度が得られ
る。
- 新たなデータが今後得られることはない、あるいは、得られても一括して計算し直
してもよい場合に有効。
- 今回のケースでは解空間が凸となるので、初期値に依存せず最適解に収束する。二
次収束なので収束も速い。
- 評価関数が二次形式の時だけ用いることが可能。データ数が増えるにつれて、行列
のサイズが大きくなり、計算負荷が増える。
＜ ＰＳＯ ＞
- データの持つ情報は解空間として現れており、複数の粒子が一定の法則に基づいて
解空間内を移動することで最適解を求める。
- 新たなデータが今後得られることはない、あるいは、得られても一括して計算し直
してもよい場合に有効。
- 初期値依存性は比較的低く、非凸の解空間においても最適解に収束する可能性があ
る。
- ハイパーパラメータが多く、最適化が必要。計算負荷は比較的高い。
以上から、線形重回帰の回帰係数を推定する場合、特段の理由がない限り、もっとも精度が良い「正規方程式を解析的に解く」もしくは同等の精度が得られる「逐次最小二乗法」を選択するのが適切と考えます。「最急降下法」を用いても良いですが、反復演算打ち切り条件などにより、多少の誤差を想定しておく必要があります。また、「ガウス・ニュートン法」を用いても良いですが、ニュートン法の近似手法であり、また、「最急降下法」と同じく最適化法なので、多少の誤差を想定しておく必要があります。なんらかの理由がない限りは、正規方程式を解く方が良いと思います。PSOは、解空間が非凸（多峰）の場合に有効な手法で、線形重回帰の回帰係数の推定手法としては、今回紹介したその他の手法と比較すると、敢えてPSOを選ぶ理由は少ないのではないかと考えます。
5. おわりに
今回は、線形重回帰の回帰係数値を５つの方法で推定し、各方法の特徴を調べてみました。線形重回帰に限らず機械学習では、なんらかの損失関数を最適化する意味でモデルのパラメータ値を決定することが多々あると思います。パラメータ値を推定（最適化）する最適な方法を選択し、また、得られたパラメータ値の持つ意味を合理的に検討できるようになるためには、様々な推定方法、最適化方法を知っておくことが有用です。今回は、線形重回帰の回帰係数値の推定を例にご紹介しましたが、興味を持って頂けましたら幸いです。
今回も、理論の大切さをあらためて知りたいあるいは感じたいという方がいらっしゃいましたら、それをできるだけわかりやすく伝えられたら、という思いから記事を書かせて頂きました。より詳しく知りたいという方は、参考文献などをご参考頂ければと思います。
参考文献
永田靖，棟近雅彦：多変量解析法入門
金谷健一：これなら分かる最適化数学
天谷賢治：工学のための最適化手法入門
足立修一：システム同定の基礎
穴井宏和：数理最適化の実践ガイド
石亀篤司，安田惠一郞：群れの知能 : particle Swarm Optimization, 知能と情報（日本知能情報ファジィ学会誌）, Vol. 20, No.6, pp.829-839 (2008)