ArmijoのBacktracking法とは

数値最適化における反復法はある関数$f:\mathbb{R}^m \rightarrow \mathbb{R}^n$に対して、適当な探索方向$d_k$を定義して

x_k \rightarrow x_k + d_k

と更新する方法です。ただし、$d_k$は降下方向でなければならないので、

\lim_{ c \rightarrow 0 } \frac{ f(x_k+cd_k) - f(x_k) }{c} = \nabla f^T(x_k) \cdot d_k < 0

を満たす必要があります。一番簡単なのは$d_k = -\nabla f(x_k)$とする最急降下法です。

さらに降下方向にを決めたのち、各ステップの更新でステップサイズを定義してし、

x_k \rightarrow x_k + \alpha_k d_k

とすれば、更に収束性をよくすることができます。これを直線探索と言います。一番よいのは

\DeclareMathOperator*{\argmin}{arg\,min} %argmin
\alpha_k = \argmin_{\alpha\geq 0}{ f(x_k + \alpha d_k) }

とすることですが、ここでまた最適化をしなければいけないことになってしまいます。そこでこの取り方を緩めて

\exists \xi \in (0,1), \exists \alpha>0 \quad s.t. \quad f(x_k+\alpha d_k) \leq 
 f(x_k) + \xi \alpha \nabla f^T(x_k) \cdot d_k

を満たす$\alpha$でよいとします。この条件をArmijo条件と言います。

これを満たす選び方として公比$\tau \in (0,1)$の等比数列

\alpha_0 > 0 \\
\alpha_{l+1} = \tau \alpha_l \quad (l=0,1,2,\cdots,L-1)\\
f(x_k+\alpha_L d_k) \leq 
 f(x_k) + \xi \alpha_L \nabla f^T(x_k) \cdot d_k \quad (\xi \in (0,1))

です。ちなみに３行目の条件を満たす$L$は必ず存在します。仮にこのような$L$が存在しなければ、

b_l = \frac{ f(x_k+\alpha_l d_k) - f(x_k) + \xi \alpha_l \nabla f^T(x_k) \cdot d_k }{\alpha_l} \\
= \frac{ f(x_k+\alpha_l d_k) - f(x_k) }{\alpha_l} + \xi \nabla f^T(x_k) \cdot d_k \geq 0

となりますが、$\alpha_l \rightarrow 0$なので、$d_k$は降下方向であることから$l\rightarrow \infty$で

b_l \rightarrow (1-\xi) \nabla f^T(x_k) \cdot d_k < 0

となり、矛盾します。
よって、このような$\alpha_{L_k}$を使って各ステップで

x_k \rightarrow x_k + \alpha_{L_k} d_k

と更新することにします。このやり方をArmijoのBacktracking法と言います。

実装

次のような最急降下法とArjimo'sBacktrackingを合わせた最適化を行うクラスを作成しました。特にステップサイズ更新のパラメーターをいじる必要がなければ、関数$f$、$df= \nabla f$と最適化停止の閾値$\epsilon$を与えればよいです。停止条件は

|\nabla f(x_k)| < \epsilon

です。

arjimo.py

import numpy as np
from numpy import linalg as LA

class ArmijoBacktrack:
    
    def __init__(self,f,df,threshold,a0=1,xi=1e-3,r=0.5,armijo_steps=30):
        self.f = f
        self.df = df
        self.threshold = threshold
        self.a0 = a0 #ステップサイズ更新の初項
        self.xi = xi
        self.armijo_ratio = r #ステップサイズ更新の公比
        self.arjimo_steps = armijo_steps　#ステップサイズ更新の上限(>L)

    def armijo_search(self,wk,dk):
        fk = self.f(wk)
        dfk = self.df(wk)
        tau = self.a0
        for i in range(self.arjimo_steps):
            cost = self.f(wk+self.a0*tau*dk) - fk - self.xi*self.a0*tau*dfk.dot(dk)
            if cost <= 0.: return tau
            tau *= self.armijo_ratio
        print("armijo parameter not found")
        return tau

    def optim(self,w_ini,steps):
        #initial guess
        w = w_ini
        for i in range(steps):
            dw = -self.df(w)
            tau = self.armijo_search(w,dw)
            err = tau*LA.norm(dw)
            wn = w + tau*dw
            w=wn.copy()
            if err < self.threshold:
                return wn
        print("not converged")
        return wn

例えば、適当な

f(x) = |b-Ax|^2

に対して

armijo.py

from functools import partial

def f(w,A,b):
    v1 = b-A@w
    return v1.dot(v1)

def df(w,A,b):
    AA = (A.T).dot(A)
    return -2*A.T@b + 2*AA@w

steps=100000
threshold = 1e-4

f2 = partial(f,A=A,
              b=b)

df2 = partial(df,A=A,
              b=b)

w_ini = np.ones(n) + 0.1*np.random.random(n)

AB = ArmijoBacktrack(f2,df2,threshold)
AB.optim(w_ini,steps)

のように実行します。

参考文献

https://projecteuclid.org/journals/pacific-journal-of-mathematics/volume-16/issue-1/Minimization-of-functions-having-Lipschitz-continuous-first-partial-derivatives/pjm/1102995080.full