制約なし最適化アルゴリズムの基礎

記事の内容


今回は, 代表的な制約なし最適化アルゴリズム, 最急降下法とNewton法を実装します.

例題: Rosenblock関数

次の関数を考えます.

\[ f(x) = 100(x_2-x_1^2)^2 + (1-x_1)^2 \]

これはRosenblockの関数と呼ばれます. この関数の最小点を求めましょう. 解析的に求めると, 最小点は(1,1)です. これを数値計算で求めます.

勾配とHesse行列を計算しておきます.

\begin{align} \nabla f(x) &= \begin{bmatrix} -400(x_2-x_1^2)x_1-2(1-x_1) \\ 200(x_2-x_1^2) \end{bmatrix} \\ \nabla^2 f(x) &= \begin{bmatrix} 800x_1^2-400(x_2-x_1^2)+2 & -400x_1\\ -400x_1 & 200\end{bmatrix} \end{align}

数値実験

概要

次の制約なし最適化問題を考えます.

\[ \text{最小化}\ f(x) \]

この問題を解くためのアルゴリズムはいくつか知られています. ここでは, 最急降下法とNewton法を用います. これらの解法は直線探索アルゴリズムと呼ばれます. 一般に, 直線探索アルゴリズムは次の更新式によって極小値に収束する点列を構成します.

\[ x_{k+1} = x_k + \alpha_k p_k \]

最急降下法の場合は \[ p_k = -\nabla f(x_k) \] とし, Newton法の場合は \[ p_k = -(\nabla^2 f(x_k))^{-1}\nabla f(x_k) \] とします. また, \(\alpha_k\)は次の条件を満たすように選びます.

\[ f(x_k + \alpha_k p_k) \leq f(x_k) + c \alpha_k \nabla f(x_k)^{\mathrm{T}}p_k \]

cは定数で, \(10^{-4}\)にすることが多いようです(以下の実験では0.5としています). これを実装するために, 今回はバックトラッキングを用います. バックトラッキングは, ある値から始めて, 上の不等式を満たすまで\(\alpha\)の値を減少させていく方法です(コード2参照).

最急降下法

最急降下法を実行します. 以下が実行結果です. 初期値は\((-1.7, 1.0)\)からはじめて, 目的関数の勾配のノルムが\(10^{-5}\)を下回る, もしくは10000回を超えるまで反復を続けました. 実際には, 1396反復で収束しました. また, 真の解に収束しました. 以下の画像は, 初期値(青丸)から, 数値解(赤丸)までの軌跡を示したものです.

【コード3の実行結果】
GRADIENT DESCENT
MINIMIZER: X1 = 1.0 X2 = 1.0
|∇f(x)|: 9.72e-6
MINIMUM: 0.0
TOLERANCE: 1.0e-5
ITERATIONS: 1396
CONVERGENCE: true

下図は, 関数値と勾配のノルムの値の変化を示したものです.

【コード4の実行結果】

Newton法

最急降下法を実行します. 以下が実行結果です. 条件は最急降下法と同じです. 初期値は\((-1.7, 1.0)\)からはじめて, 目的関数の勾配のノルムが\(10^{-5}\)を下回る, もしくは10000回を超えるまで反復を続けました. 実際には, 26反復で収束しました. また, 真の解に収束しました. 以下の画像は, 初期値(青丸)から, 数値解(赤丸)までの軌跡を示したものです.

【コード5の実行結果】
NEWTON METHOD
MINIMIZER: X1 = 1.0 X2 = 1.0
|∇f(x)|: 0.0
MINIMUM: 0.0
TOLERANCE: 1.0e-5
ITERATIONS: 26
CONVERGENCE: true

下図は, 関数値と勾配のノルムの値の変化を示したものです.

【コード6の実行結果】

ソースコード

共通のコード

【Juliaコード1; インポート】
using Plots
  using LinearAlgebra
【Juliaコード2; 関数の用意】
#バックトラッキング
function back_tracking(α₀::Float64, ρ::Float64, c::Float64, xk::Array{Float64,1}, pk::Array{Float64,1}, ∇fk::Array{Float64,1}, f)
    #初期化
    α = α₀;

    #更新
    while f(xk.+α*pk) > f(xk) + c*α*(∇fk'*pk)
        α = ρ*α;
    end
    return α
end


#最急降下法
function my_steepest_descent(ϵ::Float64, max_iter::Int64, x₀::Array{Float64,1}, selectα::String,  f, ∇f, print_flag=false)
    #初期化
    xk = x₀;
    αk = 1.0;
    fk = f(x₀);
    ∇fk = ∇f(x₀);
    pk = -∇fk;

    #保存用
    xk_array = zeros(length(x₀), max_iter);
    αk_array = zeros(max_iter-1);
    norm∇fk_array = zeros(max_iter);
    normpk_array = zeros(max_iter-1);
    fk_array = zeros(max_iter);

    #初期値の保存
    xk_array[:,1] = xk;
    norm∇fk_array[1] = norm(∇fk);
    fk_array[1] = fk;

    #反復回数
    n_iter = 1;

    #収束判定
    is_conv = false;

    #結果保存用
    ret = Dict{String,Any}();
    ret["n_iter"] = n_iter;
    ret["max_iter"] = max_iter;
    ret["x_opt"] = xk;
    ret["xk_array"] = xk_array;
    ret["αk_array"] = αk_array;
    ret["pk_array"] = normpk_array;
    ret["∇fk_array"] = norm∇fk_array;
    ret["fk_array"] = fk_array;

    #反復計算
    for k in 2:max_iter
        #局所最小値ならば抜け出す
        if norm(∇fk) < ϵ
            is_conv = true;
            break
        end

        #更新方向のベクトル
        pk = -∇fk;
        normpk_array[k-1] = norm(pk);

        #探索幅の更新
        if selectα == "back_tracking"
            αk = back_tracking(1.0, 0.5, 0.5, xk, pk, ∇fk, f);
        elseif selectα == "uniform"
            αk = αk;
        end
        αk_array[k-1] = αk;

        #解の更新
        xk = xk .+ αk*pk;
        xk_array[:,k] = xk;

        #関数値の更新
        ∇fk = ∇f(xk);
        norm∇fk_array[k] = norm(∇fk);

        #関数値の保存
        fk_array[k] = f(xk);

        #カウンタの更新
        n_iter += 1;
    end
    #更新
    ret["xk_array"] = xk_array[:, 1:n_iter];
    ret["αk_array"] = αk_array[1:n_iter-1];
    ret["pk_array"] = normpk_array[1:n_iter-1];
    ret["∇fk_array"] = norm∇fk_array[1:n_iter];
    ret["fk_array"] = fk_array[1:n_iter];
    ret["x_opt"] = xk;
    ret["n_iter"] = n_iter;
    ret["convergence"] = is_conv;

    #結果を表示する
    if print_flag == true
        println("GRADIENT DESCENT");
        print("MINIMIZER: ");
        for i in 1:length(x₀)
            print("X$(i) = $(round(xk[i], digits=3)) ");
        end
        print("\n")
        println("|∇f(x)|: $(round(norm(∇fk), digits=8))");
        println("MINIMUM: $(round(f(xk), digits=3))");
        println("TOLERANCE: $(ϵ)");
        println("ITERATIONS: $(n_iter)");
        println("CONVERGENCE: $(is_conv)");
    end
    return ret
end

#Newton法
function my_newton(ϵ::Float64, max_iter::Int64, x₀::Array{Float64,1}, selectα::String,  f, ∇f, Hf, print_flag=false)
    #初期化
    xk = x₀;
    αk = 1.0;
    fk = f(x₀);
    ∇fk = ∇f(x₀);
    Hfk = Hf(x₀);
    pk = -Hfk\∇fk;

    #保存用
    xk_array = zeros(length(x₀), max_iter);
    αk_array = zeros(max_iter-1);
    norm∇fk_array = zeros(max_iter);
    normpk_array = zeros(max_iter-1);
    fk_array = zeros(max_iter);

    #初期値の保存
    xk_array[:,1] = xk;
    norm∇fk_array[1] = norm(∇fk);
    fk_array[1] = fk;

    #反復回数
    n_iter = 1;

    #収束判定
    is_conv = false;

    #結果保存用
    ret = Dict{String,Any}();
    ret["n_iter"] = n_iter;
    ret["max_iter"] = max_iter;
    ret["x_opt"] = xk;
    ret["xk_array"] = xk_array;
    ret["αk_array"] = αk_array;
    ret["pk_array"] = normpk_array;
    ret["∇fk_array"] = norm∇fk_array;
    ret["fk_array"] = fk_array;

    #反復計算
    for k in 2:max_iter
        #局所最小値ならば抜け出す
        if norm(∇fk) < ϵ
            is_conv = true;
            break
        end

        #更新方向のベクトル
        pk = -Hfk\∇fk;
        normpk_array[k-1] = norm(pk);

        #探索幅の更新
        if selectα == "back_tracking"
            αk = back_tracking(1.0, 0.5, 0.5, xk, pk, ∇fk, f);
        elseif selectα == "uniform"
            αk = αk;
        end
        αk_array[k-1] = αk;

        #解の更新
        xk = xk .+ αk*pk;
        xk_array[:,k] = xk;

        #関数値の更新
        ∇fk = ∇f(xk);
        Hfk = Hf(xk);
        norm∇fk_array[k] = norm(∇fk);

        #関数値の保存
        fk_array[k] = f(xk);

        #カウンタの更新
        n_iter += 1;
    end
    #更新
    ret["xk_array"] = xk_array[:, 1:n_iter];
    ret["αk_array"] = αk_array[1:n_iter-1];
    ret["pk_array"] = normpk_array[1:n_iter-1];
    ret["∇fk_array"] = norm∇fk_array[1:n_iter];
    ret["fk_array"] = fk_array[1:n_iter];
    ret["x_opt"] = xk;
    ret["n_iter"] = n_iter;
    ret["convergence"] = is_conv;

    #結果を表示する
    if print_flag == true
        println("NEWTON METHOD");
        print("MINIMIZER: ");
        for i in 1:length(x₀)
            print("X$(i) = $(round(xk[i], digits=3)) ");
        end
        print("\n")
        println("|∇f(x)|: $(round(norm(∇fk), digits=8))");
        println("MINIMUM: $(round(f(xk), digits=3))");
        println("TOLERANCE: $(ϵ)");
        println("ITERATIONS: $(n_iter)");
        println("CONVERGENCE: $(is_conv)");
    end
    return ret
end

#目的関数
function f(x)
    return 100*(x[2]-x[1]^2)^2 + (1-x[1])^2
end

#勾配
function ∇f(x)
    return [-400*(x[2]-x[1]^2)*x[1]-2*(1-x[1]); 200*(x[2]-x[1]^2)]
end

#Hesse行列
function Hf(x)
    return [800*x[1]^2-400*(x[2]-x[1]^2)+2 -400*x[1]; -400*x[1] 200]
end

#プロット
function result_plot(figname::String, x₀::Array{Float64,1}, x_opt::Array{Float64,1}, xk_array::Array{Float64,2}, n_iter::Int64, top::Union{Float64,Int64}, right::Union{Float64,Int64}, bottom::Union{Float64,Int64}, left::Union{Float64,Int64},  f)
    xx = range(left,right,step=0.01);
    yy = range(bottom,top,step=0.01);
    zz = [f([x,y]) for y in yy, x in xx];
    p = contour(xx, yy, zz);
    for j in 1:n_iter-1
        plot!(xk_array[1,j:j+1], xk_array[2,j:j+1], color=:gray) ;
    end
    scatter!(xk_array[1,:], xk_array[2,:], markershape=:cross, markercolor=:gray, markersize=3, legend=false);
    scatter!([x₀[1]], [x₀[2]], color=:blue, label=["initial"]);
    scatter!([x_opt[1]], [x_opt[2]], color=:red, label=["initial"]);
    #savefig(p, figname);
end

最急降下法の実行

【Juliaコード3; 最急降下法の実行とプロット】
#パラメータと結果
x₀ = [-1.7; 1.0];
ret = my_steepest_descent(1e-5, 10000, x₀, "back_tracking", f, ∇f, true);
x_opt = ret["x_opt"];
xk_array = ret["xk_array"];
n_iter = ret["n_iter"];

#プロット
result_plot("fig1.png", x₀, x_opt, xk_array, n_iter, 2, 2, -2, -2, f);
【Juliaコード4; 最急降下法の結果をプロット】
p1 = plot(1:n_iter, ret["fk_array"], ylabel="f(x)", title="Objective function", yaxis=:log, legend=false, linewidth=1);
p2 = plot(1:n_iter, ret["∇fk_array"], ylabel="∇f(x)", title="Gradient of objective function", yaxis=:log, legend=false, linewidth=0.5);
p = plot(p1, p2, layout=(1,2), xlabel="iter", xlim=[1,n_iter], size=(700,400));
#savefig(p, "fig2.png");

Newton法の実行

【Juliaコード5; Newton法の実行】
#パラメータと結果
x₀ = [-1.7; 1.0];
ret = my_newton(1e-5, 10000, x₀, "back_tracking", f, ∇f, Hf, true);
x_opt = ret["x_opt"];
xk_array = ret["xk_array"];
n_iter = ret["n_iter"];

#プロット
result_plot("fig3.png", x₀, x_opt, xk_array, n_iter, 4, 2, 0, -2, f);
【Juliaコード6; Newton法の結果をプロット】
p1 = plot(1:n_iter, ret["fk_array"], ylabel="f(x)", title="Objective function", yaxis=:log, legend=false, linewidth=1);
p2 = plot(1:n_iter, ret["∇fk_array"], ylabel="∇f(x)", title="Gradient of objective function", yaxis=:log, legend=false, linewidth=0.5);
p = plot(p1, p2, layout=(1,2), xlabel="iter", xlim=[1,n_iter], size=(700,400));
#savefig(p, "fig4.png");