【例題をExcelでわかりやすく】最急降下法で単回帰の最小二乗法を解いてみる

2021年9月12日

前回の記事

【Excelでわかりやすく】勾配降下法で最小値が見つかる理由を視覚的に理解する

勾配降下法最急降下法)で最小値が見つかる様子を視覚的に示しました。

今回は、最小二乗法の問題を最急降下法で解く方法を、Excelを使って試してみたいと思います。

 

次のような10組の(x,y)のデータが得られたとして、xからyを予測する単回帰の式を求める問題を考えます。

x y
7 17
9 21
8 19
10 24
10 23
8 19
1 6
8 19
8 18
6 14

まず予測式f(x)を

f(x)=ax+b

とします。

 

10組ある(x,y)を(x1,y1)、(x2,y2)、、、(x10,y10)とすると、

(x1,y1)=(7,17)

(x2,y2)=(9,21)

・・・

(x10,y10)=(6,14)

のように書けます。

 

ですので、予測値と実際の値との誤差(残差)も次のように10個できます。

y1-(ax1+b)=17-(7a+b)

y2-(ax2+b)=21-(9a+b)

・・・

y10-(ax10+b)=14-(6a+b)

 

従って、残差の平方和は次のようになります。

{17-(7a+b)}2+{21-(9a+b)}2+・・・+{14-(6a+b)}2

この式の中のaとbを動かして、残差平方和の最小値を見つけるのが最小二乗法の問題です。

 

残差平方和はaとbの関数になっていますからQ(a,b)とすると、次のような式に一般化できます。

 

少し注意が必要なのは、この式はaとbの二次関数であるのに対して、

【Excelでわかりやすく】勾配降下法で最小値が見つかる理由を視覚的に理解する

で解説したのはxとyの二次関数であるということです。

変数をa、bにするかx、yにするかだけの違いです。

学校では2変数の二次関数というと変数をxとyにするのが一般的ですが、最小二乗法の式から来るとaとbが変数になります。

xとyは与えられたデータで、定数となるからです。

 

さて、勾配降下法(最急降下法)では初期値(x0,y0)を適当に選んで、x方向には∂f(x,y)/∂xと逆方向に、y方向には∂f(x,y)/∂yと逆方向に少しずつ動かしていけば、f(x,y)の最小値に辿り着きました。

今回は初期値(a0,b0)を適当に選んで、a方向には∂Q(a,b)/∂aと逆方向に、b方向には∂Q(a,b)/∂bと逆方向に少しずつ動かしていけば、Q(a,b)の最小値に辿り着きます。

 

それでは∂Q(a,b)/∂a∂Q(a,b)/∂bはどのように計算できるでしょうか?

nが大きいとQ(a,b)が長い式になってしまうため、この偏微分は一見難しそうに見えますが、合成関数の微分と考えれば簡単です。

{g(f(x))}’=g’(f(x))f’(x)

>> 【基本】合成関数の微分 | なかけんの数学ノート

 

これを使うと、∂Q(a,b)/∂aと∂Q(a,b)/∂bは次のようになります。

 

つまり、Q(a,b)∂Q(a,b)/∂a∂Q(a,b)/∂bも、10組のデータ(xi,yi)を代入した足し算になります。

Excelで計算してみると、次のようになります。

 

aとbは、取り合えずa=b=1としています。

 

このようにしてQ(a,b)∂Q(a,b)/∂a∂Q(a,b)/∂bは求まりますが、最急降下法ではaとbを少しずつ動かしていきます。

その時に、わざわざこのシートで別に計算するのは不便ですので、1行でこの計算を行うことにします。

 

aとbの初期値をa0=b0=10にして、これを1回目の更新とします。

次に∂Q(a,b)/∂aと∂Q(a,b)/∂bにa=b=10を代入し、

a1=a0-η∂Q(a,b)/∂a

= 10-η∂Q(10,10)/∂a

 

b1=b0-η∂Q(a,b)/∂b

= 10-η∂Q(10,10)/∂b

により、aとbを更新します。

η学習率で、0.001とします。

 

この計算を1行ずつ行うために、次のようにExcelに入力します。

1行が長いので、3つのスクリーンショットに分けています。

1枚目

2枚目

3枚目

 

このスクリーンショットでは7回目の更新までしか載っていませんが、aとbが収束するまで延々と続けるとQ(a,b)が最小値になるaとbの近似値が求まります。

その収束したかどうかの判定をJ列でやっています。

aとb共に、一つ前の値と比べて0.000001以下しか変わっていなければ収束したと判定しています。

今回の例では4,842回目の更新でa=1.95b=3.37となって収束しました。

グラフにすると次のようになります。

 

このシミュレーションでは学習率を0.001にしましたが、これを小さくすると収束するまで時間がかかります。

大きくすると早く収束しますが、やり過ぎると発散してしまいます。

 

発散とはaやbの値が更新ごとに大きく上下して、どんどん振れが大きくなってしまうことです。

今回の例では0.0017で発散してしまいました。

 

またη=0.0005まで小さくすると収束するまで8,960回もかかりました。

このことから、いい感じの学習率ηの範囲は意外と狭いことが分かります。

 

次に、初期値を変えてやってみましょう。

a=b=0を初期値とすると、4,534回目で次のように収束しました。

 

負の初期値も試してみましょう。

a=-30、b=-20とすると、5,488回目で収束しました。

 

このようにExcelで計算するとシートは大きくなってしまいますが、計算は一瞬で終わりますので、なかなか優れものだと思いました。

でも、残念ながら実務でこの方法は使われていないようです。

今回の例ではサンプルデータが10個だけでしたが、現実には何百万個ものデータを使ってパラメータを推定することも珍しくありません。

パラメータとは、今回の例ではaやbのことです。

パラメータ更新のたびに何百万個のデータを計算することは割に合わないのです。

 

そのための改良方法がいくつか知られています。

こちらで解説しています。

>> 【Excelでアルゴリズムを実装】確率的勾配降下法を使って最小二乗法を解いてみる

>> 【Excelでアルゴリズムを実装】モーメンタム法を使って最小二乗法を解いてみる

>> 【Excelでアルゴリズムを実装】RMSProp法を使って最小二乗法を解いてみる

>> 【Excelでアルゴリズムを実装】最終兵器Adamを使って最小二乗法を解いてみる