MATLAB,Python,Scilab,Julia比較 第4章 その109【最適化アルゴリズム⑧】

バックナンバーはこちら。
https://www.simulationroom999.com/blog/compare-matlabpythonscilabjulia4-backnumber/

はじめに
登場人物
Adamの更新式【再掲】
MATLABコード
処理結果
まとめ

はじめに

Adamに至るまでの最適化アルゴリズムの系譜とそれらの依存関係とプログラムで実現する準備まで完了。
これを今回MATLABで実現する。

登場人物

博識フクロウのフクさん

イラストACにて公開の「kino_k」さんのイラストを使用しています。
https://www.ac-illust.com/main/profile.php?id=iKciwKA9&area=1

エンジニア歴8年の太郎くん

イラストACにて公開の「しのみ」さんのイラストを使用しています。
https://www.ac-illust.com/main/profile.php?id=uCKphAW2&area=1

Adamの更新式【再掲】

太郎くん

まずはAdamの更新式を再掲。

\(
\begin{eqnarray}
m_{t+1}&=&\beta_1 m_{t-1}+(1-\beta_1)\nabla J(\theta_t)\\
v_{t+1}&=&\beta_2 v_{t-1}+(1-\beta_2)(\nabla J(\theta_t))^2\\
\displaystyle\hat{m}_{t+1}&=&\frac{m_{t+1}}{1-\beta_1}\\
\displaystyle\hat{v}_{t+1}&=&\frac{v_{t+1}}{1-\beta_2}\\
\displaystyle\theta_{t+1}&=&\theta_t-\frac{\alpha}{\sqrt{\hat{v}_{t+1}}+\epsilon}\\
m_t&:&1次のモーメント\\
v_t&:&2次のモーメント\\
\hat{m}_t,\hat{v}_t&:&バイアス補正項\\
\beta_1,\beta_2&:&指数移動平均係数(\beta_1=0.9,\beta_2=0.999)
\end{eqnarray}
\)

フクさん

これを今回はMATLABで実現する。

MATLABコード

フクさん

MATLABコードは以下。

% シグモイド関数の定義
sigmoid = @(x) 1./(1 + exp(-x));

% シグモイド関数の導関数の定義
sigmoid_derivative = @(x) sigmoid(x).*(1 - sigmoid(x));

% データの準備
X = [0 0; 0 1; 1 0; 1 1]; % 入力データ
y = [0; 1; 1; 0]; % 出力データ

% ネットワークの構築
hidden_size = 4; % 隠れ層のユニット数
output_size = 1; % 出力層のユニット数
learning_rate = 0.001; % 学習率

input_size = size(X, 2);
W1 = randn(input_size, hidden_size); % 入力層から隠れ層への重み行列
b1 = randn(1, hidden_size); % 隠れ層のバイアス項
W2 = randn(hidden_size, output_size); % 隠れ層から出力層への重み行列
b2 = randn(1, output_size); % 出力層のバイアス項

% Adamのハイパーパラメータの設定
beta1 = 0.9; % モーメンタムの指数減衰率
beta2 = 0.999; % 2次モーメントの指数減衰率
epsilon = 1e-8; % 数値安定性のための小さな値

% Adam用の変数の初期化
mW1 = zeros(size(W1));
mb1 = zeros(size(b1));
mW2 = zeros(size(W2));
mb2 = zeros(size(b2));
vW1 = zeros(size(W1));
vb1 = zeros(size(b1));
vW2 = zeros(size(W2));
vb2 = zeros(size(b2));

% 学習
epochs = 20000; % エポック数

errors = zeros(epochs, 1);  % エポックごとの誤差を保存する配列

for epoch = 1:epochs
    % 順伝播
    Z1 = X * W1 + ones(size(X, 1),1)*b1; % 隠れ層の入力
    A1 = sigmoid(Z1); % 隠れ層の出力
    Z2 = A1 * W2 + b2; % 出力層の入力
    A2 = sigmoid(Z2); % 出力層の出力

    % 誤差計算（平均二乗誤差）
    error = (1/size(X, 1)) * sum((A2 - y).^2);
    errors(epoch) = error;

    % 逆伝播
    delta2 = (A2 - y) .* sigmoid_derivative(Z2);
    delta1 = (delta2 * W2') .* sigmoid_derivative(Z1);

    grad_W2 = A1' * delta2;
    grad_b2 = sum(delta2);
    grad_W1 = X' * delta1;
    grad_b1 = sum(delta1);
    
    % パラメータの更新
    gt_W1 = grad_W1;
    gt_b1 = grad_b1;
    gt_W2 = grad_W2;
    gt_b2 = grad_b2;
    
    mW1 = beta1 * mW1 + (1 - beta1) * gt_W1;
    mb1 = beta1 * mb1 + (1 - beta1) * gt_b1;
    mW2 = beta1 * mW2 + (1 - beta1) * gt_W2;
    mb2 = beta1 * mb2 + (1 - beta1) * gt_b2;
    
    vW1 = beta2 * vW1 + (1 - beta2) * (gt_W1.^2);
    vb1 = beta2 * vb1 + (1 - beta2) * (gt_b1.^2);
    vW2 = beta2 * vW2 + (1 - beta2) * (gt_W2.^2);
    vb2 = beta2 * vb2 + (1 - beta2) * (gt_b2.^2);
    
    mHatW1 = mW1 / (1 - beta1);
    mHatb1 = mb1 / (1 - beta1);
    mHatW2 = mW2 / (1 - beta1);
    mHatb2 = mb2 / (1 - beta1);
    
    vHatW1 = vW1 / (1 - beta2);
    vHatb1 = vb1 / (1 - beta2);
    vHatW2 = vW2 / (1 - beta2);
    vHatb2 = vb2 / (1 - beta2);
    
    W1 = W1 - learning_rate * mHatW1 ./ (sqrt(vHatW1) + epsilon);
    b1 = b1 - learning_rate * mHatb1 ./ (sqrt(vHatb1) + epsilon);
    W2 = W2 - learning_rate * mHatW2 ./ (sqrt(vHatW2) + epsilon);
    b2 = b2 - learning_rate * mHatb2 ./ (sqrt(vHatb2) + epsilon);
end

% 決定境界線の表示
h = 0.01; % メッシュの間隔
[x1, x2] = meshgrid(min(X(:,1))-0.5:h:max(X(:,1))+0.5, min(X(:,2))-0.5:h:max(X(:,2))+0.5);
X_mesh = [x1(:) x2(:)];

hidden_layer_mesh = sigmoid(X_mesh * W1 + ones(size(X_mesh, 1),1)*b1);
output_layer_mesh = sigmoid(hidden_layer_mesh * W2 + b2);
y_mesh = round(output_layer_mesh); % 出力を0または1に丸める

figure;
decision_mesh = reshape(y_mesh, size(x1)); % 分類結果のメッシュを元のグリッドサイズに変形する
colormap([0.8 1 0.8; 1 0.8 0.8]); % カラーマップの順序を入れ替える
contourf(x1, x2, decision_mesh);
hold on;
scatter(X(y==1, 1), X(y==1, 2), 100, 'r', 'filled'); % クラス1のデータ点を赤でプロット
scatter(X(y==0, 1), X(y==0, 2), 100, 'g', 'filled'); % クラス0のデータ点を緑でプロット
xlabel('X1');
ylabel('X2');
title('XOR Classification');
legend('', 'Class 1', 'Class 0'); % クラスの順序を入れ替える
grid

figure;
plot(errors);
hold off;