7. Large-Scale Unconstrained Optimization

son·2023년 3월 6일
0

개요

크기가 너무 큰 문제는 계산, 메모리 효율적으로 풀어야한다. 현실적인 implementation을 위한 변형들에 대해 알아보자.

Inexact Newton method

Symmetric n×nn\times n linear system을 풀기 위한 Newton step pkNp^N_k를 생각해보자.

2fkpkN=fk(7.1)\nabla^2 f_kp^N_k = - \nabla f_k\tag{7.1}

이 section에서는 계산 효율적으로 pkNp^N_k를 approximation 하는 방법에 대해 알아본다. 이는 위의 식을 CG 또는 Lanczos method (Hessian의 negative curvature를 다루기 위한 변형 버전)을 통해 푸는 것을 기반으로 한다. Line search, trust-region을 모두 다룰 것이다.

Local convergence of inexact Newton methods

(7.1)을 푸는 대부분의 iterative solver는 residual에 근거해서 종료를 결정한다. (Exact Newton step pkNp^N_k의 경우 residual은 0일 것이다.)

rk=2fkpk+fk(7.2)r_k=\nabla ^2f_kp_k+\nabla f_k \tag{7.2}

여기서 pkp_k는 inexact Netwon step이다. 일반적으로, CG iteration은

rkηkfk(7.3)\|r_k\|\leq \eta_k\|\nabla f_k\|\tag{7.3}

일때 종료하며, 여기서 sequence {ηk}\{\eta_k\} (모든 k에 대해 0<ηk<10<\eta_k<1)을 forcing sequence라고 부른다. Forcing sequence의 선택이 어떻게 inexact Newton method의 convergence 속도에 영향을 미치는지 살펴보자.

Theorem 7.1

2f(x)\nabla^2 f(x^\ast)가 PD인 minimizer xx^\ast의 neighborhood에서 2f(x)\nabla^2f(x)가 존재하며 continuous 하다고 가정하자. pkp_k가 (7.3)을 만족시키고, 어떤 상수 η[0,1)\eta\in[0,1)에 대해 ηkη\eta_k\leq \eta를 가정할 때, iteration xk+1=xk+pkx_{k+1}=x_k+p_k을 생각해보자. 그러면, 만약 시작점 x0x_0xx^\ast에 충분히 가깝다면, sequence {xk}\{x_k\}xx^\ast로 수렴하며, η<η^<1\eta<\hat \eta < 1인 임의의 상수 η^\hat{\eta}에 대해 다음을 만족한다.

2f(x)(xk+1x)η^2f(x)(xkx)\|\nabla^2 f(x^\ast)(x_{k+1}-x^\ast)\|\leq \hat\eta\|\nabla^2 f(x^\ast)(x_{k}-x^\ast)\|

Hessian matrix 2f\nabla^2fxx^\ast에서 PD이고 xx^\ast의 주위에서 continuous하므로, 모든 충분히 xx^\ast에 가까운 xkx_k에 대해 (2fk)1L\|(\nabla^2 f_k)^{-1}\|\leq L을 만족시키는 양의 상수 LL이 존재한다 (nonsingular matrix AA에 대해서 I=A1AI=A^{-1}A이므로 A1(A)1\|A^{-1}\|\geq (\|A\|)^{-1}이므로). 따라서 우리는 (7.2)로부터 inexact Newton step이 다음을 만족한다는 것을 알 수 있다.

pk=(2fk)1(rkfk)L(fk+rk)2Lfk\|p_k\|=\|(\nabla^2f_k)^{-1}(r_k-\nabla f_k)\| \\\leq L(\|\nabla f_k\|+\|r_k\|)\leq 2L\|\nabla f_k\|

Taylor’s theorem을 이용하면

fk+1=fk+2fkpk+01[f(xk+tpk)f(xk)]pkdt=fk+2fkpk+o(pk)=fk(fkrk)+o(fk)=rk+o(fk)(7.5)\nabla f_{k+1}=\nabla f_k + \nabla^2 f_k p_k + \int^1_0[\nabla f(x_k+tp_k)-\nabla f(x_k)]p_kdt \\=\nabla f_k + \nabla^2 f_k p_k+o(\|p_k\|)\\=\nabla f_k-(\nabla f_k-r_k)+o(\|\nabla f_k\|)\\ = r_k + o(\|\nabla f_k\|)\tag{7.5}

Norm을 취하면

fk+1rk+o(fk)ηkfk+o(fk)(ηk+o(1))fk(7.6)\|\nabla f_{k+1}\|\leq \|r_k\|+o(\|\nabla f_k\|)\leq \eta_k\|\nabla f_k\|+o(\|\nabla f_k\|)\leq (\eta_k+o(1))\|\nabla f_k\|\tag{7.6}

xkx_kxx^\ast에 충분히 가까울 때, o(1)o(1) term은 (1η)/2(1-\eta)/2로 bounded 된다.

fk+1(ηk+(1η)/2)fk1+η2fk(7.7)\|\nabla f_{k+1}\|\leq (\eta_k + (1-\eta)/2)\|\nabla f_k\|\leq \frac{1+\eta}{2}\|\nabla f_k\| \tag{7.7}

따라서 gradient의 norm은 매 iteration마다 (1+η)/2(1+\eta)/2 만큼 줄어들게 된다.

Theorem 7.1을 증명하기 위해, smoothness assumption 하에

fk=2f(x)(xkx)+o(xkx)\nabla f_k= \nabla ^2 f(x^\ast)(x_k-x^\ast)+o(\|x_k-x^\ast\|)

임을 주목하자. 따라서 xkx_kxx^\ast에 가까우면 fk\nabla f_k는 scaled error 2f(x)(xkx)\nabla ^2 f(x^\ast)(x_k-x^\ast)와 상대적으로 작은 perturbation 만큼밖에 차이나지 않는다. 이 결과를 (7.7)에 집어넣으면 Theorem 7.1이 증명된다.

(7.6)에 의해

fk+1fkηk+o(1)\frac{\|\nabla f_{k+1}\|}{\|\nabla f_k\|}\leq \eta_k +o(1)

만약에 limkηk=0\lim_{k\rarr\infty}\eta_k=0이 되도록 forcing sequence를 선택하면,

limkfk+1fk=0\lim_{k\rarr \infty} \frac{\|\nabla f_{k+1}\|}{\|\nabla f_k\|} = 0

즉, gradient norm이 0으로 Q-superlinearly converge한다. {xk}\{x_k\}xx^\ast로 superlinearly converge함도 이에 뒤따른다.

Rates of convergence
Q-superlinear convergence가 무슨 말인지 몰라서 appendix를 참조한다. Q는 quotient를 의미한다고 하고, 연속된 error의 비율을 이용하기 때문에 붙었다고 한다.

(Q-linear)
충분히 큰 모든 k에 대해, 다음을 만족시키는 상수 r(0,1)r\in(0,1)이 존재한다.

xk+1xxkxr\frac{\|x_{k+1}-x^\ast\|}{\|x_k-x^\ast\|}\leq r

(Q-superlinear)

limkxk+1xxk>x=0\lim_{k\rarr\infty}\frac{\|x_{k+1}-x^\ast\|}{\|x_k->x^\ast\|}=0

(General 하게)

Order of convergence qq와 rate of convergence μ\mu에 대해

limkxk+1xxk>xq=μ\lim_{k\rarr\infty}\frac{\|x_{k+1}-x^\ast\|}{\|x_k->x^\ast\|^q}=\mu

q=1q=1이면 linear convergence (여기서 μ(0,1)\mu\in(0,1)이면 Q-linear), q=2q=2이면 quadratic convergence이다.

Hessian 2f(x)\nabla^2 f(x)xx^\ast의 주위에서 Lipschitz continuous 하다고 가정하면 quadratic convergence도 얻을 수 있다. 이 경우, (7.5)는

fk+1=rk+O(fk2)\nabla f_{k+1}=r_k+O(\|\nabla f_k\|^2)

가 될 수 있고, ηk=O(fk)\eta_k=O(\|\nabla f_k\|)로 forcing sequence를 고르면

fk+1=O(fk2)\|\nabla f_{k+1}\| = O(\|\nabla f_k\|^2)

가 되기 때문에, gradient norm이 0으로 Q-quadratic하게 converge한다 (so does xkx_k to xx^\ast)

Theorem 7.2

Theorem 7.1의 조건들이 만족하고, inexact Newton method로 생성된 {xk}\{x_k\}xx^\ast로 converge 한다 가정하자. 그렇다면, 만약 ηk0\eta_k\rarr0일 경우 rate of convergence는 superlinear 한다. 추가로, 2f(x)\nabla^2 f(x)xx^\ast 주위에서 Lipschitz continuous하고, ηk=O(fk)\eta_k=O(\|\nabla f_k\|)이면, convergence는 quadratic하다.

Line search Newton-CG method

Truncated Newton method라고도 한다. Newton equation (7.1)에 CG를 적용해서 search direction을 계산한다. 그리고 (7.3)의 termination test를 만족하는가 본다. 그런데 CG는 positive definite system을 풀려고 나왔는데, 2fk\nabla^2 f_kxkx_k가 solution에 가깝지 않을 경우 negative eignevalue를 가질 수도 있다. 따라서 negative curvature를 가지는 direction이 생성되면 CG를 바로 종료한다. 이렇게 함으로써 pkp_k가 descent direction이 되게 된다. 또한, 가능할 경우 항상 step length αk=1\alpha_k=1로 함으로써, pure Newton method의 빠른 convergence rate이 보장된다.

즉, Line search Newton-CG는 search direction pkp_k를 계산하기위해 CG를 내부 iteration으로 사용하는 알고리즘이다. (7.1)을 다음과 같이 바꿔쓰고,

Bkp=fk(7.9)B_kp=-\nabla f_k\tag{7.9}

Search direction을 djd_j, CG가 생성하는 sequence를 zjz_j라고 하자. BkB_k가 PD면, {zj}\{z_j\}는 (7.9)의 Newton step pkNp^N_k로 converge한다. 각 iteration마다 요구되는 정확도를 나타내는 tolerance ϵk\epsilon_k를 정한다. Forcing sequence는 superlinear convergence rate을 얻기 위해 ηk=min(0.5,fk)\eta_k=\min (0.5, \sqrt{\|\nabla f_k\|})로 정할 수 있고, 다른 선택지도 있다.

Line search Newton-CG에서의 CG part와 그냥 CG의 큰 차이는 1) 여기선 시작점 z0=0z_0=0이 사용됨, 2) inexact solution에서 CG가 끝나도록 하는 tolerance ϵk\epsilon_k의 사용, 3) pkp_k가 descent direction이 되도록 보장해주는 negative curvature test djTBkdj0d^T_jB_kd_j\leq 0의 사용이다. 만약에 첫 inner iteration j=0j=0에서 negative curvature가 감지되면, 출력되는 direction pk=fkp_k=-\nabla f_k는 descent direction이면서 nonpositive curvature를 갖는다.

내부 CG iteration에 preconditiong을 도입할 수도 있다.

Line search Newton-CG는 큰 문제에 적합하지만 Hessian 2fk\nabla^2f_k가 singular에 가까울 때, direction이 너무 길고 좋지 않아, line search에 필요한 계산이 많아지고 objective function의 reduction도 작다. 해결책으로 Newton step을 normalize 하거나 negative curvature test에 threshold를 도입하는 것을 생각해볼 수 있는데 구체적인 방법을 정하기 어렵다. 저자들 의견엔 그냥 이런 경우에는 밑에서 다룰 trust-region Newton-CG를 사용하는게 낫다고 말한다.

Line search Newton-CG에서는 Hessian BkB_k를 explicit하게 알 필요 없이 Hessian-vector product인 2fkd\nabla^2 f_k d만 계산할 수 있으면 되는데, Chapter 8에서 배우겠지만 Hessian을 저장하지 않고 계산하는 방식이 있다. 이를 Hessian-free Newton method라고 부른다.

Trust-region Newton-CG method

알고리즘을 만든 사람의 이름을 따서 CG-Steihaug 라고도 부른다. Trust-region subproblem을 풀 때 CG의 변형을 사용한다.

minpRnmk(p)=fk+(fk)Tp+12pTBkp subject to pΔk\min_{p\in\mathbb{R}^n}m_k(p)=f_k+(\nabla f_k)^T p + \frac{1}{2}p^TB_kp\text{ subject to }\|p\|\leq \Delta_k


djd_j는 변형된 CG iteration의 search direction이고, zjz_j는 변형된 CG가 생성하는 sequence이다.

Loop안의 첫 if문은 현재의 search direction에 대한 negative curvature test이다. 두번째 if문은 trust-region bound를 위반했을 때 loop을 멈추는 역할이다. 두 경우 모두 현재의 search direction과 trust-region boundary의 intersection을 통해 step pkp_k를 얻는다.

Tolerance ϵk\epsilon_k의 선택이 trust-region Newton-CG의 비용을 적게 유지하는데 중요하다. Well-behaved solution xx^\ast의 주위에서는 trust-region bound가 무시될 수 있을 것이고, trust-region Newton-CG는 Theorem 7.1과 7.2에서 다룬 inexact Newton method로 이해될 수 있을 것이다. 이런 경우 Line search Newton-CG와 비슷하게 ϵk\epsilon_k를 고르면 빠른 convergence를 얻을 수 있다.

Pure CG와 CG-Steihaug에 사용된 CG의 주요 차이점은 후자는 1) trust region bound를 위반할 때, 2) direction이 2fk\nabla^2 f_k에 대해 negative curvature를 갖는 direction일 때, 3) ϵk\epsilon_k로 정의되는 convergence tolerance를 만족할 때 iteration을 종료한다는 것이다.

z0=0z_0=0으로 initialize 하는게 중요한 부분이다. fk2ϵk\|\nabla f_k\|_2\geq \epsilon_k라면, Algorithm 7.2는 mk(pk)mk(pkC)m_k(p_k)\leq m_k(p^C_k)인 점 pkp_k에서 종료한다. 즉, model function에 대한 reduction이 Cauchy point로 얻어지는 것보다 같거나 크다는 것이다. 이를 보이기 위해 여러 경우들을 생각해보자. 먼저, 만약 d0TBkd0=(fk)TBkfk0d^T_0B_kd_0=(\nabla f_k)^TB_k\nabla f_k\leq 0이면, negative curvature test가 만족되어 Cauchy point p=Δk(fk)/fkp=-\Delta_k(\nabla f_k)/\|\nabla f_k\|을 반환한다. 이외의 경우, Algorithm 7.2는 z1z_1을 다음과 같이 정의한다.

z1=α0d0=r0Tr0d0TBkd0d0=(fk)Tfk(fk)TBkfkfkz_1=\alpha_0 d_0=\frac{r^T_0r_0}{d^T_0B_kd_0}d_0=-\frac{(\nabla f_k)^T\nabla f_k}{(\nabla f_k)^TB_k\nabla f_k}\nabla f_k

만약 z1<Δk\|z_1\| < \Delta_k라면 z1z_1은 정확히 Cauchy point이다. Algorithm 7.2의 이어지는 step들은 최종 pkp_kmk(pk)mk(z1)m_k(p_k)\leq m_k(z_1)을 만족하도록 만든다. 만약 z1Δk\|z_1\|\geq \Delta_k면 Cauchy point를 반환한다. 모든 경우에 최소 Cauchy point 만큼은 model mkm_k를 감소시키기 때문에 Algorithm 7.2는 globally convergent하다.

또 중요한 성질 중 하나는 zjz_j의 norm이 매 iteration마다 zj1z_{j-1}의 norm보다 커진다는 것이다. 이는 z0=0z_0=0으로 initialize한 결과이다. mkm_k를 넘어서 탐색해봤자 의미가 없기 때문에 boundary에 도착하면 iteration을 종료한다는 것을 함의한다.

Theorem 7.3

Algorithm 7.2에 의해 생성되는 vector의 sequence {zj}\{z_j\}는 다음을 만족한다.

0=z02<<zj2<zj+12<<pk2Δk0=\|z_0\|_2<\cdots<\|z_j\|_2<\|z_{j+1}\|_2<\cdots<\|p_k\|_2\leq \Delta_k
  • 증명

    먼저 Algorithm 7.2에 의해 생성되는 vector들의 sequence가 j0j\geq 0에 대해 zjTrj=0z^T_jr_j=0이고 j1j\geq 1에 대해 zjTdj>0z^T_j d_j>0임을 만족함을 보이자. zj+1z_{j+1}zjz_j에 대해 recursive 하게 계산되지만, 다음과 같이 풀어 쓸 수 있다.

    zj=z0+i=0j1αidi=i=0j1αidiz_j=z_0+\sum^{j-1}_{i=0}\alpha_i d_i=\sum^{j-1}_{i=0}\alpha_i d_i

    양변에 r_j를 곱하고 Theorem 5.2의 내용을 적용하면 (간략히 설명하면 CG에서는 각 iteration마다 이전 direction vector들과는 orthogonal한 새 direction에 대해 minimize 하기 때문에 residual과 지난 direction vector의 내적은 0이라는 내용)

    zjTrj=i=0j1αidiTrj=0z^T_j r_j=\sum^{j-1}_{i=0}\alpha_id^T_ir_j=0

    귀납적으로 zjTdj>0z^T_jd_j>0임을 보이자. 먼저

    z1Td1=(α0d0)T(r1+β1d0)=α0β1d0Td0>0z^T_1 d_1=(\alpha_0d_0)^T(-r_1+\beta_1 d_0)=\alpha_0\beta_1d^T_0d_0>0

    임을 알 수 있다. 이제 zjTdj>0z^T_jd_j >0이 참이면 zj+1Tdj+1>0z^T_{j+1}d_{j+1}>0도 참임을 보이면 된다.

    zj+1Tdj+1=zj+1T(rj+1+βj+1dj)=βj+1zj+1Tdj=βj+1(zj+αjdj)Tdj=βj+1zjTdj+αjβj+1djTdjz^T_{j+1}d_{j+1}=z^T_{j+1}(-r_{j+1}+\beta_{j+1}d_j)\\=\beta_{j+1}z^T_{j+1}d_j\\=\beta_{j+1}(z_j+\alpha_jd_j)^Td_j\\=\beta_{j+1}z^T_jd_j+\alpha_j\beta_{j+1}d^T_j d_j

    zjTdjz^T_jd_jβj+1\beta_{j+1}αj\alpha_j모두 양수이므로 zj+1Tdj+1>0z^T_{j+1}d_{j+1}>0이다.

    이제 Theorem을 증명하기 위한 준비가 끝났다. 만약 Algorithm 7.2가 djTBkdj0d^T_jB_kd_j\leq 0 또는 zj+12Δk\|z_{j+1}\|_2\geq \Delta_k에 의해 종료한다면, pkp_k는 가능한 최대 길이인 pk2=Δk\|p_k\|_2 = \Delta_k를 만족하도록 골라질 것이다. 다른 경우들에 대해서는 zj2<zj+12\|z_j\|_2 <\|z_{j+1}\|_2임을 보여야한다.

    zj+122=(zj+αjdj)T(zj+αjdj)=zj22+2αjzjTdj+αj2dj22\|z_{j+1}\|^2_2=(z_j+\alpha_jd_j)^T(z_j+\alpha_jd_j)=\|z_j\|^2_2+2\alpha_jz^T_jd_j+\alpha^2_j\|d_j\|^2_2

    따라서 증명이 성립한다.

Theorem에 따라 zjz_j들은 z1z_1으로부터 pkp_k로 이어지는 path를 interpolate한다고 볼 수 있다는 걸 알았다. Bk=2fkB_k=\nabla^2 f_k가 PD일 때, 이 path는 dogleg method와 비교해볼 수 있다. 두 method 모두 fk-\nabla f_k에서 시작해서 trust-region boundary condition이 위반되지 않는 범위에서 차근차근 pkNp^N_k로 향한다. Bk=2fkB_k=\nabla^2f_k가 PD일 때, Algorithm 7.2는 optimal 보다 적어도 절반 이상 좋은 model mkm_k의 reduction을 가져온다.

Preconditioning the trust-region Newton-CG method

우리는 preconditioning이 CG를 가속시켜준다는 것을 알고있다. Preconditioning 하는 상황에 대해, Theorem 7.3을 일반화시켜 zjz_j가 weighted norm D\|D\cdot\| 에 대해 monotonically 증가한다고 할 수 있다. 우선 trust-region subproblem을 weighted norm에 대해 다시 적으면,

minpRnmk(p)=fk+(fk)Tp+12pTBkp subject to DpΔk(7.13)\min_{p\in\mathbb{R}^n}m_k(p)=f_k+(\nabla f_k)^T p + \frac{1}{2}p^TB_kp\text{ subject to }\|Dp\|\leq \Delta_k \tag{7.13}

여기서 p^=Dp\hat p=Dp로 정의하면,

g^k=DTB^k=DT(2fk)D1\hat g_k=D^{-T}\\ \hat B_k=D^{-T}(\nabla^2 f_k)D^{-1}

가 되고 (7.13)을 새로 정의한 variable로 적으면

minp^Rnfk+g^kTp^+12p^kTB^kp^ subject to p^Δ\min _{\hat p\in\mathbb{R}^n} f_k +\hat g^T_k\hat p+\frac{1}{2}\hat p^T_k\hat B_k \hat p \text{ subject to }\|\hat p\| \leq \Delta

이므로 기존의 trust-region subproblem과 똑같은 형태이다. 따라서 Algorithm 7.2를 수정하는 일 없이 바로 적용 가능하다.

Preconditioner은 여러가지 선택지가 있는데, 그 중 incomplete Cholesky factorization에 대해 얘기해보자. PD인 matrix BB에 대해 다음을 만족하는 lower triangular matrix LL을 구한다.

B=LLTRB=LL^T-R

이 때 LL에 제약조건이 붙는다. 예를 들면, BB의 lower triangular part와 같은 sparsity 구조를 가져야 한다던지, BB보다 nonzero element의 수가 적어야 한다던지. RR은 approximate factorization의 inexactness이다. Hessian 2fk\nabla^2 f_k가 indefinite한 경우에는 문제가 좀 복잡해진다. 다음 알고리즘은 trust-region Newton-CG의 preconditioner를 정의하기 위해 incomplete Cholesky와 modified Cholesky를 합쳤다.

Cholesky 분해 후에는 D=LTD=L^T로 preconditioning 하면 된다.

Trust-region Newton-Lanczos method

Algorithm 7.2의 한계는 negative curvature를 갖는 어느 direction도 모두 허용해서, model에 대한 reduction이 적은 것까지 채택한다는 것이다. 예를 들어,

minpm(p)=103p1104p12p22 subject to p1\min_p m(p)=10^{-3}p_1-10^{-4}p^2_1-p^2_2\text{ subject to }\|p\|\leq 1

문제를 생각해보자. p=0p=0일 때 steepest descent direction은 (103,0)T(-10^{-3},0)^T이고 negative curvature를 가진다. Algorithm 7.2는 이 direction과 trust-region boundary와 경계를 채택하면 model mm에 대한 reduction은 약 10310^{-3}이다. 다음 step도 비슷하게 이어질 것이다.

챕터 4에서 Hessian 2fk\nabla^2 f_k가 negative eigenvaules를 가질 때, search direction이 2fk\nabla^2 f_k의 가장 작은 eigenvalue에 해당하는 eigenvector 쪽 성분을 크게 가질 것이라는 것을 보았다. (Theorem 4.1의 1번과 그에 대한 분석 Iterative solution of the subproblem 참조, 내가 이해하기로는 p(λ)=(B+λI)1gp(\lambda)=-(B+\lambda I)^{-1}g 인데, negative eigenvalue가 존재하면 B+λIB+\lambda I가 PD여야 한다는 조건에 의해 λ\lambda가 양의 값을 가져야한다. 이 때, B+λI=Q(Λ+λI)QTB+\lambda I=Q(\Lambda +\lambda I)Q^T로 분해하면 p(λ)=j=1nqjTgλj+λqjp(\lambda)=-\sum^n_{j=1}\frac{q^T_j g}{\lambda_j +\lambda}q_j 이니, 가장 작은 eigenvalue를 가지는 eigenvector 성분이 분모 값이 작아지므로 두드러질 것이다.) 이 성질을 이용하여 algorithm이 minimizer가 아닌 stationary point에서 빠르게 벗어나도록 만들 수 있다. 여기서 앞에 괄호 안에 적힌 내용을 적용할 수 있다는 것이다.

다른 방법은 Bkp=fkB_kp=-\nabla f_k를 풀 때 CG 대신 Lanczos method를 이용하는 것이다. Lanczos method는 indefinite system에 적용 가능한 CG의 일반화 된 버전이라고 생각하면 되고, negative curvature information을 모으며 CG를 계속 진행할 수 있게 해준다.

jj step 후, Lanczos method는 n×jn\times j matrix QjQ_j를 만드는데, QjQ_j의 column들은 서로 orthogonal하며 Lanczos method에 의해 생성되는 Krylov subspace를 span한다. QjQ_jQjTBQj=TjQ^T_jBQ_j=T_j (여기서 TjT_j는 tridiagonal 행렬) 가 되는 성질을 가지고 있다. Tridiagonal하다는 점을 이용해서 basis QjQ_j의 range 안에서 trust-region subproblem의 approximate solution을 찾을 수 있다. e1=(1,0,0,,0)Te_1=(1,0,0,\dots,0)^T일 때, 그렇게 함으로써 다음의 문제를 풀고

minwRjfk+e1TQj(fk)e1Tw+12wTTjw subject to wΔk\min_{w\in\mathbb{R}^j}f_k+e^T_1 Q_j (\nabla f_k)e^T_1 w+\frac{1}{2}w^TT_jw \text{ subject to } \|w\|\leq \Delta_k

trust-region subproblem의 approximate solution은 pk=Qjwp_k=Q_jw로 정의한다. TjT_j가 tridiagonal이므로 위 problem은 Tj+λIT_j+\lambda I를 factoring하고 챕터 4의 Iterative solution of the subproblem의 방법처럼 풀면 된다.

Limited-memory quasi-Newton methods

큰 문제의 Hessian matrix가 계산 비용이 너무 크거나 sparse 하지 않을 때 사용한다. Hessian의 n×nn\times n approximation을 모두 저장하는 대신, nn 길이를 갖는 vector들을 몇개 저장한다. 그럼에도 불구하고 꽤 괜찮은 (linear 이지만) rate of convergence를 보여준다. 여기서는 L-BFGS를 주로 다룰것이다. L-BFGS의 메인 아이디어는 Heesian approximation에 최근 수 iteration에서 얻은 curvature 정보만 사용하는 것이다.

Limited-memory BFGS

먼저 BFGS를 표현하는 수식들로부터 시작하자. 각 BFGS의 step은

xk+1=xkαkHkfkx_{k+1}=x_k-\alpha_kH_k\nabla f_k

여기서 αk\alpha_k는 step length이고 HkH_k는 approximated inverse Hessian이며

Hk+1=VkTHkVk+ρkskskTρk=1ykTskVk=IρkykskTsk=xk+1xkyk=fk+1fkH_{k+1}=V^T_k H_k V_k + \rho_k s_k s^T_k\\ \rho_k = \frac{1}{y^T_ks_k}\\V_k = I-\rho_k y_k s^T_k \\ s_k=x_{k+1}-x_k\\ y_k= \nabla f_{k+1}-\nabla f_k

HkH_k를 모두 저장하는 대신 mm개(일반적으로 3~20개)의 {si,yi}\{s_i, y_i\} vector pair 만을 저장하자 (이제 정확히는 modified 버전 HkH_k를 저장하는 것이다). 이를 이용해서 HkH_k대신 HkfkH_k\nabla f_k를 바로 계산한다. 각 step마다 새 vector pair를 저장하고 제일 오래된 것은 버린다.

이제 update 과정을 조금 더 자세히 살펴보자. kk번째 iteration에서 먼저 initial Hessian approximation Hk0H^0_k를 고른다 (BFGS와 다르게 각 iteration마다 initial approximation이 달라도 된다). BFGS의 Hk+1H_{k+1} 구하는 식을 반복한다고 보면,

Hk=(Vk1TVkmT)Hk0(VkmVk1)+ρkm(Vk1TVkm+1T)skmskmT(Vkm+1Vk1)+ρkm+1(Vk1TVkm+2T)skm+1skm+1T(Vkm+2Vk1)++ρk1sk1sk1TH_k=(V^T_{k-1}\cdots V^T_{k-m})H^0_k(V_{k-m}\cdots V_{k-1})\\+\rho_{k-m}(V^T_{k-1}\cdots V^T_{k-m+1})s_{k-m}s^T_{k-m}(V_{k-m+1}\cdots V_{k-1})\\+\rho_{k-m+1}(V^T_{k-1}\cdots V^T_{k-m+2})s_{k-m+1}s^T_{k-m+1}(V_{k-m+2}\cdots V_{k-1}) \\ + \cdots \\ +\rho_{k-1} s_{k-1}s^T_{k-1}

이다. 이를 이용해서 HkfkH_k\nabla f_k를 구하는 recursive한 방법을 표현하면,

Hk0qH^0_kq 계산을 빼면 총 4mn4mn 번의 곱셈 연산이 필요하다. 만약 Hk0H^0_k가 diagonal이면 nn 번의 추가적인 곱셈 연산이 필요하다. 계산 복잡도가 낮을 뿐만 아니라 initial matrix Hk0H^0_k가 나머지 계산과 분리되어 있어 각 iteration마다 자유롭게 고를 수 있다는 장점도 있다. 효과적인 방법 중 하나는 Hk0=γkIH^0_k=\gamma_k I로 잡는 것이라고 한다. 이 때,

γk=sk1Tyk1yk1Tyk1\gamma_k=\frac{s^T_{k-1}y_{k-1}}{y^T_{k-1}y_{k-1}}

γk\gamma_k는 가장 최근 search direction에 대한 true Hessian matrix의 size를 추정하는 scaling factor이다. 이 방법을 사용하면 search direction pkp_k가 잘 scaled 되게 보장해주고 대부분의 iteration에서 step length αk=1\alpha_k=1이 선택되도록 만든다.

Iteration이 mm 번 이상 진행되어서 계산에 필요한 vector pair들이 쌓이기 전에는 BFGS와 동일하게 작동한다 (initial Hessian이 똑같이 선택되고, L-BFGS에서 Hk0=H0H^0_k=H_0을 선택하는 경우).

L-BFGS의 가장 큰 단점은 ill-conditioned problem (특히 eigenvalue가 넓게 퍼져 있는 문제)에서 수렴이 느리다는 것이다.

Relationship with conjugate gradient methods

Limited-memory method는 nonlinear conjugate gradient method를 개선하고자 하는 시도해서 진화했고, 초기 구현은 quasi-Newton method 보다는 CG에 닮았다. 이 두 class 간의 관계가 memoryless BFGS iteration의 근간이 된다.

Nonlinear CG의 하나인 Hestenes-Stiefeil로부터 시작하자. Search direction은

pk+1=fk+1+fk+1TykykTpkpk=(IskykTykTsk)fk+1=H^k+1fk+1p_{k+1}=-\nabla f_{k+1} + \frac{\nabla f^T_{k+1} y_k}{y^T_k p_k} p_k = -\left(I - \frac{s_ky^T_k}{y^T_ks_k}\right)\nabla f_{k+1}=-\hat H_{k+1}\nabla f_{k+1}

이 식은 quasi-Newton과 닮았지만 matrix H^k+1\hat H_{k+1}이 symmetric, PD가 아니다. H^k+1TH^k+1\hat H^T_{k+1}\hat H_{k+1}로 symmetric화 할 수 있지만 secant equation은 만족하지 못하고, singular하다. Symmetric하고 PD이고 secant equation을 만족시키는 iteration matrix는

Hk+1=(IskykTykTsk)(IykskTykTsk)+skskTykTskH_{k+1}=\left(I - \frac{s_ky^T_k}{y^T_ks_k}\right)\left(I - \frac{y_ks^T_k}{y^T_ks_k}\right) + \frac{s_ks^T_k}{y^T_ks_k}

이는 L-BFGS update 식이 Hk=IH_{k}=I일 때 적용된 것과 동일하다. 항상 이전 Hessian approximation이 identity matrix로 reset된다는 의미에서 위 식으로 업데이트 되는 알고리즘들을 “memoryless” BFGS 라고 할 수 있다.

Memoryless BFGS에 exact line search를 고려하면 CG와의 관계가 더 잘 드러난다. 식을 써보면,

pk+1=Hk+1fk+1=fk+1+fk+1TykykTpkpkp_{k+1}=-H_{k+1}\nabla f_{k+1}=-\nabla f_{k+1}+\frac{\nabla f^T_{k+1}y_k}{y^T_kp_k} p_k

Hestenes-Stiefel과 완전히 동일하다. 여기서 모든 kk에 대해 fk+1Tpk=0\nabla f^T_{k+1}p_k=0이라 가정하면 Hestenes-Stiefel에서 Polak-Ribiere가 된다.

General limited-memory updating

Quasi-Newton matrix를 compact (또는 outer product) form으로 나타냄으로써 모든 quasi-Newton update formula와 inverse에 대해 효율적인 implementation을 유도할 수 있다. 여기서는 correction pair를 연속적으로 update하는 경우만 다룬다 (반대되는 방식으로는 storage가 가득차면 모두 버리고 새로 프로세스를 시작하는 것이 있다. 실험 결과 효율이 떨어진다 한다.).

Compact representation of BFGS updating

Theorem 7.4

B0B_0이 symmetric, PD라고 하고, kk vector pair {si,yi}i=0k1\{s_i, y_i\}^{k-1}_{i=0}siTyi>0s^T_iy_i>0을 만족한다 가정하자. BkB_kB0B_0kk번 BFGS update를 적용하여 얻어진다고 하자. 그러면,

Bk=B0[B0SkYk][SkTB0SkLkLkTDk]1[SkTB0YkT](7.24)B_k=B_0-\left[\begin{matrix}B_0S_k & Y_k \end{matrix}\right]\left[ \begin{matrix}S^T_kB_0S_k & L_k \\ L^T_k & -D_k\end{matrix} \right]^{-1}\left[\begin{matrix}S^T_kB_0 \\ Y^T_k\end{matrix}\right]\tag{7.24}

여기서

Sk=[s0,,sk1]Yk=[y0,,yk1]S_k=[s_0,\dots, s_{k-1}]\\Y_k=[y_0,\dots,y_{k-1}]

n×kn\times k matrix이고

(Lk)i,j={si1Tyj1 if i>j0otherwiseDk=diag[s0Ty0,,sk1Tyk1](L_k)_{i,j}=\begin{cases}s^T_{i-1}y_{j-1} \text{ if } i>j \\ 0 \qquad\quad \text{otherwise}\end{cases} \\ D_k=\text{diag}[s^T_0y_0,\dots,s^T_{k-1}y_{k-1}]

k×kk\times k matrix이다.

siTyi>0,i=0,1,,k1s^T_i y_i>0, i=0,1,\dots,k-1 조건이 (7.24)의 가운데 matrix가 nonsingular하도록 보장해준다. L-BFGS처럼 최대 mm개의 최신 correction pair {si,yi}\{s_i, y_i\}를 갖고 오래된 pair는 버린다. Bk0=δkIB^0_k=\delta_k I (여기서 δk=1/γk\delta_k=1/\gamma_k)로 선택하면 첫 mm iteration 동안의 Theorem 7.4의 식을 바로 사용할 수 있다.

k>mk>m 번째 iteration에서는 {si,yi}\{s_i,y_i\}의 변하는 성질을 반영하기 위해, n×mn\times m matrix SkS_kYkY_k를 다음과 같이 정의하고

Sk=[skm,,sk1]Yk=[ykm,,yk1]S_k=[s_{k-m}, \dots , s_{k-1}]\\Y_k=[y_{k-m},\dots,y_{k-1}]

B0(k)=δkIB^{(k)}_0=\delta_k Imm 번의 업데이트가 적용된다 하면,

Bk=δkI[δkSkYk][δkSkTSkLkLkTDk]1[δkSkTYkT](7.29)B_k=\delta_k I-\left[\begin{matrix}\delta_kS_k & Y_k \end{matrix}\right]\left[ \begin{matrix}\delta_kS^T_kS_k & L_k \\ L^T_k & -D_k\end{matrix} \right]^{-1}\left[\begin{matrix}\delta_kS^T_k \\ Y^T_k\end{matrix}\right]\tag{7.29}

여기서,

(Lk)i,j={skm1+iTykm1+j if i>j0otherwiseDk=diag[skmTykm,,sk1Tyk1](L_k)_{i,j}=\begin{cases}s^T_{k-m-1+i}y_{k-m-1+j} \text{ if } i>j \\ 0 \qquad\qquad\qquad\qquad \text{otherwise}\end{cases} \\ D_k=\text{diag}[s^T_{k-m}y_{k-m},\dots,s^T_{k-1}y_{k-1}]

m×mm\times m matrix이다. xk+1x_{k+1}이 생성되고 나면 SkS_k에서 skms_{k-m}을 지우고 sks_k를 추가해서 Sk+1S_{k+1}를 얻는다 (Yk+1Y_{k+1}, Lk+1L_{k+1}, Dk+1D_{k+1}도 비슷하게 얻음).

(7.29) 식의 가운데 matrix의 dimension이 2m2m으로 작기 때문에 factorization의 계산 비용이 무시할만하다. 위 그림에서 볼 수 있듯이 correction을 길고 narrwo한 matrix들의 outer product로 볼 수 있다.

BkB_k의 limited-memory update는 2mn+O(m3)2mn+O(m^3) operation이 필요하고 matrix-vector 곱 BkvB_kv(4m+1)n+O(m2)(4m+1)n+O(m^2) 곱셈연산이 필요하다. mm이 작을 때는 BkB_k를 바로 다루는 것이 경제적이다.

BkB_k의 approximation이 trust-region method나, 특히 bound-constrained, general-constrained optimization에서 사용될 수 있다. L-BFGS-B은 bound constraint이 있는 큰 nonlinear optimization을 풀기 위해 compact limited-memory approximation을 활용한다. 이 경우, constraint gradients에 의해 정의된 subspace에 대한 BkB_k의 projection이 반복적으로 계산되어야 한다.

물론, compact representation을 inverse Hessain approximation에 적용할 수도 있다. 또, SR1에 적용될 수도 있는데, 표기를 Theorem 7.4에서 빌려오면

Bk=B0+(YkB0Sk)(Dk+Lk+LkTSkTB0Sk)1(YkB0Sk)TB_k=B_0+(Y_k-B_0S_k)(D_k+L_k+L^T_k-S^T_kB_0S_k)^{-1}(Y_k-B_0S_k)^T

Limited-memory SR1로 바꾸는 것도 limited-memory BFGS랑 같은 방식으로 할수 있다.

Unrolling the update

Limited-memory update를 더 간단하게 구현할 수는 없을까? Limited-memory BFGS update는 compact representation에 기반한 방법보다 꽤 비싸다는 것을 보자. Direct BFGS를 다시 써보면

Bk+1=BkakakT+bkbkTB_{k+1}=B_k-a_ka^T_k+b_kb^T_k

여기서,

ak=Bksk(skTBksk)1/2bk=yk(ykTsk)1/2a_k=\frac{B_ks_k}{(s^T_kB_ks_k)^{1/2}}\\b_k=\frac{y_k}{(y^T_ks_k)^{1/2}}

Basic matrix Bk0B^0_k에 대해

Bk=Bk0+i=kmk1[bibiTaiaiT]B_k=B^0_k+\sum^{k-1}_{i=k-m}[b_ib^T_i-a_ia^T_i]

Vector pair {ai,bi}\{a_i, b_i\}{si,yi}\{s_i, y_i\}로부터 복원될 수 있는데,

위 과정을 살펴보면, 각 iteration kk마다 aia_i는 vector pair {skm,ykm}\{s_{k-m}, y_{k-m}\}에 종속적이라 모두 재계산되어야 하는 반면, bib_i의 경우는 bk1b_{k-1}bjTsk1b^T_j s_{k-1}만 재계산하면 된다. 이 사실을 고려해 계산을 수정하면 여전히 대략 32m2n\frac{3}{2}m^2n 연산이 limited-memory matrix를 결정하기 위해 필요한데 compact form은 2mn2mn multiplication만을 필요로한다.

Sparse quasi-Newton updates

Quasi-Newton approximation BkB_k가 true Hessian과 같은 (또는 비슷한) sparsity pattern을 갖도록 하는 방법에 대해서 알아보자. 이렇게 함으로써 필요한 storage는 줄이고, 정확도는 오히려 향상될 수 있다.

우리가 어느 점에서 Hessian의 어느 component가 nonzero인지 안다고 가정하자. 즉, 우리는 다음의 set의 내용을 안다.

Ω={(i,j)[2f(x)]ij0 for some x in the domain of f}\Omega=\{(i,j)|[\nabla^2f(x)]_{ij}\neq 0\text{ for some }x\text{ in the domain of } f\}

또, 현재의 Hessian approximation BkB_k의 exact Hessian의 nonzero structure을 반영한다고 가정하자. 즉, (i,j)Ω(i,j)\notin \Omega에 대해 (Bk)ij=0(B_k)_{ij}=0이다. BkB_kBk+1B_{k+1}로 업데이트 하면서 1) secant condition이 만족하고, 2) 같은 sparsity pattern을 가지고, 3) BkB_k와 최대한 가까운 Bk+1B_{k+1}를 찾고 싶다. 이를 최적화 문제로 표현하면,

minBBBkF2=(i,j)Ω[Bij(Bk)ij]2subject to Bsk=yk,B=BT,Bij=0 for (i,j)Ω\min_B\|B-B_k\|^2_F=\sum_{(i,j)\in\Omega}[B_{ij}-(B_k)_{ij}]^2\\\text{subject to } Bs_k=y_k,\, B=B^T,\, B_{ij}=0\text{ for }(i,j)\notin \Omega

이 최적화 문제의 solution Bk+1B_{k+1}은 sparsity pattern이 Ω\Omegan×nn\times n linear system을 풀어서 얻을 수 있다. Bk+1B_{k+1}을 얻고나면 이를 trust-region method에 이용해서 xk+1x_{k+1}을 얻을 수 있다. 참고로 Bk+1B_{k+1}은 PD가 아닐 수 있다. 그런데 사실 이런 방식은 scale invariance하지 못하고 실제 해보면 성능이 떨어진다는 단점이 있다.

다른 방법으로 secant equation을 relax해서 마지막 step에 strictly 만족하도록 하기보다는 최근 몇 step에 approximately 만족하도록 하는 것이다. SkS_k, YkY_k를 limited-memory BFGS 때처럼 최근 mm개에 대해 저장한 것으로 정의하면,

minBBSkYkF2subject to B=BT,Bij=0 for (i,j)Ω\min_B\|BS_k-Y_k\|^2_F\\\text{subject to }B=B^T,\,B_{ij}=0\text{ for }(i,j)\notin \Omega

이 convex optimization은 solution이 있지만 계산하기 쉽지 않다. 또 생성된 Hessian approximation이 singular하거나 poorly condition 되어있을 수 있다. 앞에 방식보단 낫지만 큰 문제에서는 성능이 썩 좋지는 않다.

Algorithms for partially separable functions

Separable하다는 것은 함수가 더 간단한 독립적으로 최적화 될 수 잇는 함수들의 합으로 분해될 수 있다는 것을 말한다. 예를들어,

f(x)=f1(x1,x3)+f2(x2,x4,x6)+f3(x5)f(x)=f_1(x_1, x_3)+f_2(x_2,x_4,x_6)+f_3(x_5)

변수가 겹치지 않음을 주목하자. 일반적으로 더 낮은 차원의 최적화를 여러번 하는게 높은 차원 최적화 한번 하는것 보다 경제적이다.

일반적으로 큰 objective function f:RnRf: \mathbb{R}^n\rarr \mathbb{R}은 separable하지 않은데, 그래도 element function이라는 것의 sum으로 적을 수 있다. 각 element function은 많은 수의 linearly independent direction으로 움직였을 때 영향받지 않는다는 성질이 있다. 이런 성질을 만족할 때 ff가 partially separble하다고 말한다. Hessian 2f\nabla^2 f가 sparse한 모든 함수는 partially separable하다. Partially separable하면 어떻게 도움되는지 살펴보자.

가장 간단한 형태의 partial separability는 objective function이 다음과 같이 쓰여질 수 있을 때다.

f(x)=i=1nefi(x)f(x)=\sum^{ne}_{i=1}f_i(x)

여기서 fif_i는 각각 몇개의 xx에만 의존한다. 그러면,

f(x)=i=1nefi(x)2f(x)=i=1ne2fi(x)\nabla f(x)=\sum^{ne}_{i=1}\nabla f_i(x)\\\nabla^2 f(x)=\sum^{ne}_{i=1}\nabla^2 f_i(x)

이런 경우에 전체 Hessian 2f(x)\nabla^2 f(x)를 approximate하는 것보다 효율적인 방법이 있을까? Quasi-Newton approximation이 이런 구조를 잘 활용한다면 가능하다. 다음의 objective function을 예를 들어 설명하자.

f(x)=(x1x32)2+(x2x42)2+(x3x22)2+(x4x12)2=f1(x)+f2(x)+f3(x)+f4(x)f(x)=(x_1-x_3^2)^2+(x_2-x_4^2)^2+(x_3-x_2^2)^2+(x_4-x_1^2)^2\\=f_1(x)+f_2(x)+f_3(x)+f_4(x)

각 element function fif_i의 Hessian 은 4×44\times 4 sparse, singular, nonzero element가 4개 있는 matrix이다.

f1f_1에 집중하자. f1f_1xx의 함수지만 x1x_1x3x_3에만 의존한다. 이를 f1f_1의 element variable이라 부른다. Element variable을 vector로 모아놓은 것을 다음 처럼 표기한다.

x[1]=[x1x3]x_{[1]}=\left[\begin{matrix} x_1 \\ x_3 \end{matrix}\right]

이를 xx에 대해 표현하면,

x[1]=U1xU1=[10000010]x_{[1]}=U_1 x\\U_1=\left[\begin{matrix} 1 & 0& 0& 0 \\ 0 &0 &1&0 \end{matrix}\right]

그리고, function ϕ1\phi_1을 다음과 같이 정의하자.

ϕ1(z1,z2)=(z1z22)2\phi_1(z_1,z_2)=(z_1-z^2_2)^2

그러면,

f1(x)=ϕ1(U1x)f1(x)=U1Tϕ1(U1x)2f1(x)=U1T2ϕ1(U1x)U1f_1(x)=\phi_1(U_1x)\\ \nabla f_1(x)=U^T_1\nabla \phi_1(U_1x)\\ \nabla^2 f_1(x)=U^T_1\nabla^2\phi_1(U_1x)U_1

우리가 가정한 f1f_1의 경우

2ϕ1(U1x)=[24x34x312x324x1]2f1(x)=[204x3000004x3012x324x100000]\nabla^2\phi_1(U_1x)=\left[\begin{matrix} 2&-4x_3 \\-4x_3 &12x^2_3-4x_1\end{matrix}\right]\\\nabla^2 f_1(x)= \left[\begin{matrix} 2 & 0 & -4x_3 & 0 \\ 0 & 0& 0& 0\\ -4x_3 & 0&12x^2_3-4x_1&0\\ 0 & 0& 0& 0\end{matrix}\right]

Matrix U1U_1을 compactifying matrix라고 부르는데, 위 두 matrix를 비교해보면 알 수 있듯이 derivative 정보를 작은 dimension에 저장할 수 있게 해준다.

그래서 핵심 아이디어는 quasi-Newton approximation 2f1\nabla^2 f_1을 저장하는 대신에 2ϕ1\nabla^2\phi_12×22\times 2 quasi-Newton approximation B[1]B_{[1]}을 저장해두자는 것이다. B[1]B_{[1]}을 업데이트 하기 위해 xx에서 x+x^+로 가는 step 후 다음 정보를 저장한다.

s[1]=x[1]+x[1]y[1]=ϕ1(x[1]+)ϕ1(x[1])s_{[1]}=x^+_{[1]}-x_{[1]}\\y_{[1]}=\nabla \phi_1(x^+_{[1]})-\nabla \phi_1(x_{[1]})

그리고 B[1]+B^+_{[1]}를 얻기 위해 BFGS나 SR1을 사용한다. 그러면

B[1]2ϕ1(U1x)=2ϕ1(x[1])B_{[1]}\approx \nabla^2\phi_1(U_1x)=\nabla^2\phi_1(x_{[1]})

이고, 2f1\nabla^2f_1을 얻기 위해서는

2f1(x)U1TB[1]U1\nabla ^2f_1(x)\approx U^T_1B_{[1]}U_1

로 transformation해서 얻을 수 있다. 이를 각 element function에 대해 적용하면

f(x)=i=1neϕi(Uix)B=i=1neUiTB[i]Uif(x)=\sum^{ne}_{i=1}\phi_i(U_ix) \\ B=\sum^{ne}_{i=1}U^T_iB_{[i]}U_i

이다.

0개의 댓글