最適化問題を考える np.dot縛りで勾配降下法を実装

#最適化問題について
某先生の「ディープラーニングは最小二乗法みたい（に損失関数を最小化することでデータからモデルのパラメータチューニングを行うよう）なものですよ」という発言が話題になっていますね。（2019/4月 現在）

本稿では今流行りの(??)最小二乗法を題材にして損失関数の最適化がどういうことなのか、np計算のみで考えていきたいと思います。

注意:損失関数の議論を単純化するために値がすべて真値であることを仮定した上で議論を進めています。
注意:B4にCSの楽しさを伝えることを目的に投稿したため厳密性にかける議論が散見されます。ベイズ教徒の方ごめんなさい

#最小二乗法における損失関数の定義
最小二乗法では様々な最適化手法が存在しますが、直観的理解のため勾配降下法を用います。

さて、ここで$\mathbb{R}^{n}$を満たす$m$個の点$\{\boldsymbol{x}^{(1)},...,\boldsymbol{x}^{(m)}\}$を考えます。

この点の集合$\boldsymbol{X}$に対して

\boldsymbol{Y} = \boldsymbol{AX}\quad where\quad \boldsymbol{Y} \in \mathbb{R}^{m},\boldsymbol{A} \in \mathbb{R}^{n},\boldsymbol{X} \in \mathbb{R}^{m\times n}\dots(1)

実験等により$\boldsymbol{X}$と$\boldsymbol{Y}$が明らかになっており、上記の式が成り立つことを仮定した上で$\boldsymbol{A}$を探索するのが最小二乗法のざっくりとした理解です。

以下$\boldsymbol{A}$の探索方法について述べます。

最小二乗法を含む最適化問題の一般論としてパラメータを推定する場合には誤差を定義し、誤差を最小化するようにモデルを更新していきます。

今回はベクトル$\boldsymbol{A}$を推定するためにベクトル$\boldsymbol{w}\in \mathbb{R}^{n}$を導入し、以下の式の最小化を考えます。

Loss \quad = \quad ||\boldsymbol{Y} - \boldsymbol{wX}||^{2}_{2}\quad\dots(2)

この式では$\boldsymbol{w}$を用いて算出した$\boldsymbol{wX}$が真値である$\boldsymbol{Y}$に対してどれだけ誤差があるかを定量的に評価しています。
つまり損失関数の値が小さければ小さいほど$\boldsymbol{w}$が高精度に$\boldsymbol{A}$を予測していると言えます。

#勾配降下法の考え方
勾配降下法では定義した損失関数に対してgradientを取ることにより最小値を探索します。砕けた説明をすると勾配微分を取って勾配が大きい方向にガンガン進んで行って勾配が０になる場所を探索します。
定式化すると以下のようになります。

w^{'} = w -  \varepsilon\nabla _{\boldsymbol{w}}Loss\dots(3)

ここで$\varepsilon$は更新する幅を表します。あまりにも大きすぎると最適解まで落ちないし、細かすぎると学習が進まなかったり局所解にハマります。

損失関数の微分は以下の式になります。

\begin{align}
\nabla _{\boldsymbol{w}}Loss &= \nabla _{\boldsymbol{w}}||\boldsymbol{Y} - \boldsymbol{wX}||^{2}_{2}\dots(4)\\
&= 2(\boldsymbol{wX}-\boldsymbol{Y})\boldsymbol{X}^{T}\dots(5)
\end{align}

方向微分は式(5)を使用します。

#実装の考え方
①損失関数を計算する
②方法微分を用いてパラメーターを更新する。

この２プロセスを繰り返すことで最適なパラメータ$w^{*}$を探索します。

#実装
それではpythonを用いて勾配降下法を実装していきましょう。
まず最初に使用するサンプルデータを生成しておきます。
使用するデータは

\boldsymbol{X}\in \mathbb{R}^{30\times4}\\
\boldsymbol{A}\in \mathbb{R}^{4}\\
\boldsymbol{Y}\in \mathbb{R}^{30}\\
\boldsymbol{w}\in \mathbb{R}^{4}\\

です。また$\boldsymbol{w}$については最初の値を決めておかないといけないため乱数で出力してあります。（本当は適当ダメ）

import numpy as np
from numpy.random import randn

X = np.array(randn(4,30),dtype='float32')
A = np.array([10,30,-6,50],dtype='float32')
Y = np.dot(A,X)
w = np.random.randn(4)

次に損失関数の計算を行う関数を定義しておきます。

def loss(w,X,Y,size):
    loss = Y-np.dot(w,X)
    loss = np.sum(np.power(loss,2))
    return loss

同様にパラメーターを更新する関数を定義します。

def update_w(w,X,Y,step):
    grad = np.dot(np.dot(w,X)-Y,X.T)
    w = w - step*grad
    return w

以上の関数を用いいることにより最適化を行います。

Loss = loss(w,X,Y,30)
while Loss>0.01:
    step = 0.01
    w = update_w(w,X,Y,step)
    Loss = loss(w,X,Y,30)
print("w:{}".format(w))
print("A:{}".format(A))

実行結果は以下のようになります。

runfile('/media/lab/ADATA UFD/python/gradient_opt.py', wdir='/media/lab/ADATA UFD/python')
w:[10.01936697 29.9838687  -5.99657676 50.00292774]
A:[10,30,-6,50]

この精度に物足りなさを感じた方は最適化方法を工夫すると少し精度が上がることを体験できます。

Loss = loss(w,X,Y,30)
while Loss>0.01:
    if Loss>10:
        step = 0.01
    elif Loss > 0.1:
        step = 0.001  
    else:
        step = 0.0001
    w = update_w(w,X,Y,step)
    Loss = loss(w,X,Y,30)
print("w:{}".format(w))
print("A:{}".format(A))

runfile('/media/lab/ADATA UFD/python/gradient_opt.py', wdir='/media/lab/ADATA UFD/python')
w:[ 9.99682429 30.01050901 -6.00936462 49.98222603]
A:[10. 30. -6. 50.]

以上です。