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)
のように実行します。