【入門】最適化アルゴリズム(Julia)【数値計算】

MATLAB、Python、Scilab、Julia比較ページはこちら
https://www.simulationroom999.com/blog/comparison-of-matlab-python-scilab/

はじめに
Adamの更新式【再掲】
Juliaコード
処理結果
まとめ

はじめに

の、

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

を書き直したもの。

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

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}
\)

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

Juliaコード

Juliaコードは以下。

using PyPlot

function sigmoid(x)
    return 1.0 ./ (1.0 + exp.(-x))
end

function sigmoid_derivative(x)
    return sigmoid(x) .* (1.0 - sigmoid(x))
end

function meshgrid(xin,yin)
    nx=length(xin)
    ny=length(yin)
    xout=zeros(ny,nx)
    yout=zeros(ny,nx)
    for jx=1:nx
        for ix=1:ny
            xout[ix,jx]=xin[jx]
            yout[ix,jx]=yin[ix]
        end
    end
    return (x=xout, y=yout)
end

function MultilayerPerceptron()
	# データの準備
	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 in 1:epochs
	    # 順伝播
        Z1 = X * W1 .+ b1 # 隠れ層の入力
	    A1 = sigmoid.(Z1)  # 隠れ層の出力
        Z2 = A1 * W2 .+ b2 # 出力層の入力
	    A2 = sigmoid.(Z2)  # 出力層の出力

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

	    # 逆伝播
	    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 -= learning_rate * mHatW1 ./ (sqrt.(vHatW1) .+ epsilon)
	    b1 -= learning_rate * mHatb1 ./ (sqrt.(vHatb1) .+ epsilon)
	    W2 -= learning_rate * mHatW2 ./ (sqrt.(vHatW2) .+ epsilon)
	    b2 -= learning_rate * mHatb2 ./ (sqrt.(vHatb2) .+ epsilon)
	end

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

	hidden_layer_mesh = sigmoid.(X_mesh * W1 .+ b1)
	output_layer_mesh = sigmoid.(hidden_layer_mesh * W2 .+ b2)
	y_mesh = round.(output_layer_mesh)
	figure()
	decision_mesh = reshape(y_mesh, size(x1))  # 分類結果のメッシュを元のグリッドサイズに変形する
	colormap = ["#CCFFCC","#FFCCCC"]  # 各領域の色を指定する
	contourf(x1, x2, decision_mesh, levels=1, colors=colormap)
	scatter(X[y .== 1, 1], X[y .== 1, 2], 100, "r", facecolors="none",label="Class 1")  # クラス1のデータ点を赤でプロット
	scatter(X[y .== 0, 1], X[y .== 0, 2], 100, "g", facecolors="none",label="Class 0")  # クラス0のデータ点を緑でプロット
	xlabel("X1")
	ylabel("X2")
	title("XOR Classification")
	legend(loc="best")  # クラスの順序を入れ替える
	grid(true)
	show()
	
	figure()
	plot(errors)
	show()
end

MultilayerPerceptron()