大学で受けた講義の復習として, 表題の内容について自分なりの理解でまとめます.

制約なし最適化問題

とは, 一般の$d$変数関数について, 以下のように定式化されるものでした.

$\mathcal{X} \subset \mathbb{R}^d$ 上の目的関数 $f:\mathcal{X}\mapsto\mathbb{R}$ の最小値を求めよ. $$ \min_{\vec{x}\in\mathcal{X}} f(\vec{x}) $$ ただし, $f$ は微分可能とする.

最大化問題は $f$ の符号を逆にすることで最小化問題になるので, 同じことです.

$f$ の定義域における制約条件（$\vec{x}$ に関する等式・不等式）が特にないので, 制約なし最適化問題です.

ここで, $f$の勾配を $\nabla f(\vec{x}) = \left(\dfrac{\partial f}{\partial x_1}, \dfrac{\partial f}{\partial x_2}, \ldots, \dfrac{\partial f}{\partial x_d} \right)^\mathrm{T}$ とすると, 以下の最適性条件

$\vec{x}^*$ が最適解ならば $\nabla f(\vec{x}^*) = \vec{0}$

が成り立ちます. 関数がなめらかな場合, 最小値をとるときは必ず微分がゼロだという話です.

もちろん微分がゼロでも最小値をとるとは限らないのですが, $f$ が凸関数である場合, 「$\vec{x}^*$ が最適解」と「 $\nabla f(\vec{x}^*) = \vec{0}$」は同値になります.（凸最適化問題）

凸関数

「$f$は凸関数である」とは,

$f$の定義域上の任意の$\vec{x}, \vec{x}'$に対して $$f((1-\theta)\vec{x}+\theta\vec{x}') \leq f((1-\theta)\vec{x})+f(\theta\vec{x}')$$

となることです. 「関数上のどの二点を結ぶ線分より関数の方が常に小さい」とも言えます.

また, $f$が一階微分可能な場合は, 上の定義は

$$ f(\vec{x}') \geq f(\vec{x}) + \nabla f(\vec{x})^\mathrm{T}(\vec{x}'-\vec{x}) $$

とも同値となります. 「任意の点における接線は, 常に関数自体より小さい」ということですね.

さらに, $f$ が二階微分可能な場合,

ヘッセ行列 $\nabla^{2} f(\vec{x})$の固有値が全て非負（ヘッセ行列が半正定値となる）

となることとも同値です. 一変数関数の場合の「二階微分が非負」であることに対応します.

前述の通り, $f$が凸関数である場合, 「 $\nabla f(\vec{x}^*) = \vec{0}$ ならば $\vec{x}^*$ が最適解」が成り立ちます.

(制約なし)最適化アルゴリズムは一般に「勾配がゼロになる点に収束するような点列を構成する」手続きといえるため, 凸関数であれば大域最適解が求まり嬉しいのですが, 凸でない関数だと 局所最適解しか得られません. その場合, 複数の初期値からアルゴリズムを実行して一番よいものを選ぶなどします.

勾配法（最急降下法）

代表的な最適化アルゴリズムが 勾配法（最急降下法） です. 勾配法は一般に以下の式で表されます.

適当な初期値から始め, 収束するまで以下を繰り返す : $$ \vec{x}_{k+1} \leftarrow \vec{x}_k - \varepsilon_k\nabla f(\vec{x}_k)$$

ここで $\varepsilon_k > 0$ はステップ幅と呼ばれる係数です.

要するに, $f$の勾配 $\nabla f(\vec{x}) = \left(\dfrac{\partial f}{\partial x_1}, \dfrac{\partial f}{\partial x_2}, \ldots, \dfrac{\partial f}{\partial x_d} \right)^\mathrm{T}$ は, ベクトルとしては「その点の近傍において, 最も関数値が増加する方向」を意味すると言えるので, 「各ターン, その位置から勾配の逆(減少)方向にいくらか進める」ことを繰り返すという戦略ですね.

ステップ幅 $\varepsilon_k > 0$ は大きすぎると収束しませんが, 小さすぎると収束が遅いです. これをどのように決めるかによって, 勾配法にはいくつか種類があります.

そのうちのひとつが以下で実装を示すバックトラック直線探索です.

アルミホ規準とバックトラック直線探索

バックトラック直線探索の定義をまず以下に示します.

勾配法 $$ \vec{x}_{k+1} \leftarrow \vec{x}_k - \varepsilon_k\nabla f(\vec{x}_k)$$ の各 $k$ において $\varepsilon_k $ を計算し直す.

$\varepsilon_k = 1$ と初期化したのち, アルミホ規準 : $$ f( \vec{x}_k - \varepsilon_k \nabla f(\vec{x}_k) ) - f(\vec{x}_k) \leq - \alpha \varepsilon_k || \nabla f(\vec{x}_k) ||^2 $$ が成り立つまで, $\varepsilon_k$を, $$\varepsilon_k \leftarrow \beta\varepsilon_k$$ を繰り返すことで減衰させる.

ここで, $0<\alpha, \beta<1$ は事前に定める定数である.

バックトラック直線探索とは, 勾配法の各ループでアルミホ規準というルールに従うように$\varepsilon_k$を選びなおすことで, 関数$f$の減少量に関してある保証を行うというアイデアに基づく勾配法です.

以下, アルミホ規準について見ていきます.

$$ f( \vec{x}_k - \varepsilon_k \nabla f(\vec{x}_k) ) - f(\vec{x}_k) \leq - \alpha \varepsilon_k || \nabla f(\vec{x}_k) ||^2 $$

左辺は, ある$ \varepsilon_k$が定まったとして, そのとき$ \vec{x}_{k+1} \leftarrow \vec{x}_k - \varepsilon_k\nabla f(\vec{x}_k)$と更新した場合の$f$の減少量であることがわかると思います.

問題は右辺なのですが, 結論から言うと, これは「$f$の減少量の線形予測値(の$\alpha$倍)」です.

前述の通り, 勾配法の$k$ステップにおける$f$の減少量は$\varepsilon_k$が定まれば一意に定まるので, $f$の減少量は, $\varepsilon_k$の関数

$$g(\varepsilon_k) = f( \vec{x}_k - \varepsilon_k \nabla f(\vec{x}_k) ) - f(\vec{x}_k)$$

として考えることができます. その上で, $\varepsilon_k=0$における微係数$g^{\prime}(0)$を, \begin{align} \left. \frac{\mathrm{d} g(\varepsilon_k) }{\mathrm{d} \varepsilon_k} \right|_{\varepsilon_k=0} &= \left. \left\{ - \nabla f(\vec{x})^{\mathrm{T}} \nabla f( \vec{x}_k - \varepsilon_k \nabla f(\vec{x}_k) ) \right\} \right|_{\varepsilon_k=0} \\ &= - \nabla f(\vec{x})^{\mathrm{T}} \nabla f(\vec{x}) \\ &= - || \nabla f(\vec{x}_k) ||^2 \end{align}

と計算すると, $f$の減少量$g(\varepsilon_k)$の線形予測値として, $$ \varepsilon_k g^{\prime}(0) = - \varepsilon_k || \nabla f(\vec{x}_k) ||^2 $$ を得ます.

以上のことから, アルミホ規準は,

少なくとも各ステップで, $f$の減少量の線形予測値の$\alpha < 1$倍は減少するようにする

ようなルールであることがわかります.

バックトラック直線探索の julia 実装

前置きが長くなりました. 実装しましょう. 使用するのは Julia 1.0.2 です.

今回も可視化に Plots + GR を用いることにします.

using LinearAlgebra
using Plots

gr()
ENV["PLOTS_TEST"] = "true"

また, 今回は簡単のため, 最適化を行う目的関数として,

$$ f(\vec{x}) = 10 x_1² + x_2² $$

を用いることにして, それを定義しておきます. 明らかに最適解は$(x_1, x_2)=(0,0)$でその値は$f(\vec{x}^*)=0$ですが, アルゴリズムの機能の仕方を可視化しながら観察するには単純であるほうがよいでしょう.

ついでに勾配を返す関数も定義します（都合から二変数に限ってしまいましたが）

function f(x)
    return 10x[1]^2 + x[2]^2
end

function ∇(f,x)
    h = 0.0001
    return [( f(x+[h,0]) - f(x-[h,0]) ) / 2h, 
            ( f(x+[0,h]) - f(x-[0,h]) ) / 2h]
end

勾配関数∇はあるものとして, 目的関数$f$と$\alpha, \beta$を入力として受け取って, バックトラック直線探索を行って得られた$\vec{x}$および$f(\vec{x})$の列を出力します. ここではデバッグのため各ループでも値をプリントさせるようにしました.

function BTLSGradientDescent(f, α=0.5, β=0.8)
    x = (rand(2).-0.5).*10
    i = 0
    ε = 1
    m = 0.1
        
    x_seq = zeros(0)
    f_seq = zeros(0)
    append!( x_seq, x )
    append!( f_seq, f(x) )
    print(i,":",x," f:",f(x),"\n")
    
    while abs(norm( ∇(f,x) - ∇(f,x - ε.*∇(f,x)) )) > m
        if f(x - ε.*∇(f,x)) - f(x) <= -α*ε.*norm(∇(f,x))  # Armijo
            x = x - ε.*∇(f,x)
            ε = 1
            i += 1
            print(i,":",x," f:",f(x),"\n")
            append!( x_seq, x )
            append!( f_seq, f(x) )
        else
            ε = β*ε
        end
    end
    return reshape(x_seq, 2,:), f_seq
end

実行結果は以下です. $f$が最小値$0$に収束するのがわかります.

x_seq, f_seq = BTLSGradientDescent(f)

0:[2.67815, 2.54762] f:78.21508380277415
1:[-1.92287, 2.10994] f:41.426324107210576
2:[1.3806, 1.74746] f:22.114138557399237
3:[-0.991252, 1.44725] f:11.920326667200023
4:[0.711706, 1.19861] f:6.501923402334815
5:[-0.510995, 0.992693] f:3.5966018700862796
6:[0.366888, 0.822149] f:2.0219977430959752
7:[-0.263421, 0.680905] f:1.1575369871905883
8:[0.189133, 0.563926] f:0.6757248514793868
9:[-0.135795, 0.467045] f:0.4025329889285466
10:[0.0974989, 0.386807] f:0.24467998560906534
11:[-0.0700029, 0.320354] f:0.1516308200203651
12:[0.0262084, 0.276325] f:0.08322424811647708
13:[-0.00260804, 0.245943] f:0.06055581977627555

([2.67815 -1.92287 … 0.0262084 -0.00260804; 2.54762 2.10994 … 0.276325 0.245943], [78.2151, 41.4263, 22.1141, 11.9203, 6.50192, 3.5966, 2.022, 1.15754, 0.675725, 0.402533, 0.24468, 0.151631, 0.0832242, 0.0605558])

この結果を可視化すると, 以下のようになります.

g(x, y) = begin
    10x^2 + y^2
end
a = b = -5:0.1:5
X = repeat(reshape(a, 1, :), length(b), 1)
Y = repeat(a, 1, length(b))
Z = map(g, X, Y)
p = contour(a, b, Z)
px = plot(p)

plot!(px,x_seq[1,:],x_seq[2,:],m=3,label="x")
annotate!(px,x_seq[1,:],x_seq[2,:].-0.1, [text(i,:top,8) for i in 1:length(x_seq[1,:])])

pf = plot(f_seq, m=2, label="f(x)", xticks=1:length(x_seq[1,:]), xlabel="iteration")

plot(pf,px,layout=(1,2),wsize=(1200,400))