前回
筑波大学の機械学習講座：課題のPythonスクリプト部分を作りながらsklearnを勉強する (7) 最急降下法を自作
 https://github.com/legacyworld/sklearn-basic

課題 4.3 最急降下法と確率的最急降下法

解説は第5回(1) 24分30秒あたり
前回は再急降下法しか出来なかったので、今回は確率的再急降下法を実装してみる。
プログラム自体はそこまで大きくは変わらない。
数学的には以下のような感じ。

\lambda = 正則化パラメータ,
\beta = \begin{pmatrix} \beta_0 \\ \beta_1\\ \vdots \\ \beta_m \end{pmatrix},
y = \begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_N \end{pmatrix},
X = \begin{pmatrix}
1&x_{11}&x_{12}&\cdots&x_{1m}\\
1&x_{21}&x_{22}&\cdots&x_{2m}\\
\vdots\\
1&x_{N1}&x_{N2}&\cdots&x_{Nm}
\end{pmatrix}\\ \\
\beta^{t+1} = \beta^{t}(1-2\lambda\eta) - \eta\frac{1}{N}x_i^T(x_i\beta^t-y_i)

これまで勾配の計算に全てのデータを利用していたのをランダムに選んだ$x_i,y_i$のみで計算を行う。
ちょっとひっかかったのが、NumpyのTransposeの仕様。
1次元ベクトルの場合transpose(.T)を使うとそのまま返ってきてしまうため、.reshape(-1,1)を使う必要がある。
こちらを参照。
https://note.nkmk.me/python-numpy-transpose/

a_1d = np.arange(3)
print(a_1d)
# [0 1 2]

print(a_1d.T)
# [0 1 2]

a_col = a_1d.reshape(-1, 1)
print(a_col)
# [[0]
#  [1]
#  [2]]

ソースコードはこちら。

Homework_4.3SGD.py

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn import preprocessing
from sklearn.base import BaseEstimator
from sklearn.model_selection import cross_validate
import statsmodels.api as sm

class MyEstimator(BaseEstimator):
    def __init__(self,ep,eta,l):
        self.ep = ep
        self.eta = eta
        self.l = l
        self.loss = []
    # fit()を実装
    def fit(self, X, y):
        self.coef_ = self.stochastic_grad_desc(X,y)
        # fit は self を返す
        return self

    # predict()を実装
    def predict(self, X):
        return np.dot(X, self.coef_)

    def shuffle(self,X,y):
        r = np.random.permutation(len(y))
        return X[r],y[r]

    def stochastic_grad_desc(self,X,y):
        m = len(y)
        loss = []
        # 特徴量の種類
        dim = X.shape[1]
        # betaの初期値
        beta = np.ones(dim).reshape(-1,1)
        eta = self.eta
        l = self.l
        X_shuffle, y_shuffle = self.shuffle(X,y)
        # T回改善されなければ終了
        T = 100
        # 改善されない回数
        not_improve = 0
        # 目的関数最小値初期値
        min = 10 ** 9
        while True:
            for Xi,yi in zip(X_shuffle,y_shuffle):
                loss.append((1/(2*m))*np.sum(np.square(np.dot(X,beta)-y)))
                beta = beta*(1-2*l*eta) - eta*Xi.reshape(-1,1)*(np.dot(Xi,beta)-yi)
                if loss[len(loss)-1] < min:
                    min = loss[len(loss)-1]
                    min_beta = beta
                    not_improve = 0
                else:
                    # 目的関数の最小値が更新されない場合
                    not_improve += 1
                    if not_improve >= T:
                        break
            # 全サンプル終わったがT回以内に最小値が変わっている場合再度ループ
            if not_improve >= T:
                self.loss = loss
                break
        return min_beta

#scikit-leanよりワインのデータをインポートする
df= pd.read_csv('winequality-red.csv',sep=';')
# 目標値であるqualityが入っているので落としたdataframeを作る
df1 = df.drop(columns='quality')
y = df['quality'].values.reshape(-1,1)
X = df1.values
scaler = preprocessing.StandardScaler()
X_fit = scaler.fit_transform(X)
X_fit = sm.add_constant(X_fit) #最初の列に1を加える
epsilon = 10 ** (-7)
eta_list = [0.03,0.01,0.003]
loss = []
coef = []
for eta in eta_list:
    l = 10**(-5)
    test_min = 10**(9)
    while l <= 1/(2*eta):
        myest = MyEstimator(epsilon,eta,l)
        myest.fit(X_fit,y)
        scores = cross_validate(myest,X_fit,y,scoring="neg_mean_squared_error",cv=10)
        if abs(scores['test_score'].mean()) < test_min:
            test_min = abs(scores['test_score'].mean())
            loss = myest.loss
            l_min = l
            coef = myest.coef_
        l = l * 10**(0.5)
    plt.plot(loss,label=f"$\eta$={eta}")
    print(f"eta = {eta} : iter = {len(loss)}, loss = {loss[-1]}, lambda = {l_min}, TestErr = {test_min}")
    # 係数の出力　一番最初に切片が入っているので2つ目から取り出して、最後に切片を出力
    i = 1
    for column in df1.columns:
        print(column,coef[i][0])
        i+=1
    print('intercept',coef[0][0])
plt.legend()
plt.savefig("sgd.png")

解説の中で「100回連続で改善が無ければ停止」という条件が書いてあったので、それにあわせてある。
$\eta$が小さいパターンでは1599個すべて使ってもこの条件に当てはまらないため2周目に入ることもある。
結果はこんな感じ。

$\eta = 0.03$では最後の方で少し誤差が増えたりしているのがよくわかる。
最後に求められた係数等。

eta = 0.03 : iter = 298, loss = 0.29072324272824085, lambda = 0.0031622776601683803, TestErr = 0.47051639691326796
fixed acidity 0.1904239451124434
volatile acidity -0.11242984344193296
citric acid -0.00703125780915424
residual sugar 0.2092352618792849
chlorides -0.044795495356479025
free sulfur dioxide -0.018863685196341816
total sulfur dioxide 0.07447982325062003
density -0.17305138620126106
pH 0.05808006453308803
sulphates 0.13876262568557934
alcohol 0.2947134691111974
intercept 5.6501294014064145
eta = 0.01 : iter = 728, loss = 0.24203354045966255, lambda = 0.00010000000000000002, TestErr = 0.45525344581852156
fixed acidity 0.25152952212309976
volatile acidity -0.03876889927769888
citric acid 0.14059421863669852
residual sugar 0.06793602828251821
chlorides -0.0607861479963043
free sulfur dioxide 0.08441853171277111
total sulfur dioxide -0.09642176480191654
density -0.2345690991118163
pH 0.1396740265674562
sulphates 0.1449843342292861
alcohol 0.19737851967044345
intercept 5.657998427200384
eta = 0.003 : iter = 1758, loss = 0.22475918775097103, lambda = 0.00010000000000000002, TestErr = 0.44693442950748147
fixed acidity 0.2953542653448508
volatile acidity -0.12934364893075953
citric acid 0.04629080083382285
residual sugar 0.013753852832452122
chlorides -0.03688613363045954
free sulfur dioxide 0.045541235818534614
total sulfur dioxide -0.049594638329345575
density -0.17427360277645224
pH 0.13897225246491407
sulphates 0.15425590075925466
alcohol 0.26518804857692096
intercept 5.597149258230254

ランダムにシャッフルされているので、実行するごとに結果が変わっている。

過去の投稿

筑波大学の機械学習講座：課題のPythonスクリプト部分を作りながらsklearnを勉強する (1)
筑波大学の機械学習講座：課題のPythonスクリプト部分を作りながらsklearnを勉強する (2)
筑波大学の機械学習講座：課題のPythonスクリプト部分を作りながらsklearnを勉強する (3)
筑波大学の機械学習講座：課題のPythonスクリプト部分を作りながらsklearnを勉強する (4)
筑波大学の機械学習講座：課題のPythonスクリプト部分を作りながらsklearnを勉強する (5)
筑波大学の機械学習講座：課題のPythonスクリプト部分を作りながらsklearnを勉強する (6)
https://github.com/legacyworld/sklearn-basic