7. Large-Scale Unconstrained Optimization

개요

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

Inexact Newton method

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

이 section에서는 계산 효율적으로 $p^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 $p^N_k$의 경우 residual은 0일 것이다.)

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

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

Theorem 7.1

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

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

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

Taylor’s theorem을 이용하면

Norm을 취하면

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

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

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

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

(7.6)에 의해

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

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

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

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

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

(Q-superlinear)

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

(General 하게)

Order of convergence $q$와 rate of convergence $\mu$에 대해

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

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

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

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

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

Theorem 7.2

Theorem 7.1의 조건들이 만족하고, inexact Newton method로 생성된 $\{x_k\}$가 $x^\ast$로 converge 한다 가정하자. 그렇다면, 만약 $\eta_k\rarr0$일 경우 rate of convergence는 superlinear 한다. 추가로, $\nabla^2 f(x)$가 $x^\ast$ 주위에서 Lipschitz continuous하고, $\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을 풀려고 나왔는데, $\nabla^2 f_k$는 $x_k$가 solution에 가깝지 않을 경우 negative eignevalue를 가질 수도 있다. 따라서 negative curvature를 가지는 direction이 생성되면 CG를 바로 종료한다. 이렇게 함으로써 $p_k$가 descent direction이 되게 된다. 또한, 가능할 경우 항상 step length $\alpha_k=1$로 함으로써, pure Newton method의 빠른 convergence rate이 보장된다.

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

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

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

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

Line search Newton-CG는 큰 문제에 적합하지만 Hessian $\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 $B_k$를 explicit하게 알 필요 없이 Hessian-vector product인 $\nabla^2 f_k d$만 계산할 수 있으면 되는데, Chapter 8에서 배우겠지만 Hessian을 저장하지 않고 계산하는 방식이 있다. 이를 Hessian-free Newton method라고 부른다.

Trust-region Newton-CG method

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

$d_j$는 변형된 CG iteration의 search direction이고, $z_j$는 변형된 CG가 생성하는 sequence이다.

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

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

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

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

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

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

Theorem 7.3

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

$0=\|z_0\|_2<\cdots<\|z_j\|_2<\|z_{j+1}\|_2<\cdots<\|p_k\|_2\leq \Delta_k$

• 증명

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

  $z_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이라는 내용)

  $z^T_j r_j=\sum^{j-1}_{i=0}\alpha_id^T_ir_j=0$

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

  $z^T_1 d_1=(\alpha_0d_0)^T(-r_1+\beta_1 d_0)=\alpha_0\beta_1d^T_0d_0>0$

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

  $z^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$

  $z^T_jd_j$과 $\beta_{j+1}$와 $\alpha_j$모두 양수이므로 $z^T_{j+1}d_{j+1}>0$이다.

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

  $\|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에 따라 $z_j$들은 $z_1$으로부터 $p_k$로 이어지는 path를 interpolate한다고 볼 수 있다는 걸 알았다. $B_k=\nabla^2 f_k$가 PD일 때, 이 path는 dogleg method와 비교해볼 수 있다. 두 method 모두 $-\nabla f_k$에서 시작해서 trust-region boundary condition이 위반되지 않는 범위에서 차근차근 $p^N_k$로 향한다. $B_k=\nabla^2f_k$가 PD일 때, Algorithm 7.2는 optimal 보다 적어도 절반 이상 좋은 model $m_k$의 reduction을 가져온다.

Preconditioning the trust-region Newton-CG method

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

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

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

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

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

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

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

Trust-region Newton-Lanczos method

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

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

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

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

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

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

Limited-memory quasi-Newton methods

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

Limited-memory BFGS

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

여기서 $\alpha_k$는 step length이고 $H_k$는 approximated inverse Hessian이며

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

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

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

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

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

Iteration이 $m$ 번 이상 진행되어서 계산에 필요한 vector pair들이 쌓이기 전에는 BFGS와 동일하게 작동한다 (initial Hessian이 똑같이 선택되고, L-BFGS에서 $H^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은

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

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

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

Hestenes-Stiefel과 완전히 동일하다. 여기서 모든 $k$에 대해 $\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

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

$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}$

여기서

$S_k=[s_0,\dots, s_{k-1}]\\Y_k=[y_0,\dots,y_{k-1}]$

는 $n\times k$ matrix이고

$(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\times k$ matrix이다.

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

$k>m$ 번째 iteration에서는 $\{s_i,y_i\}$의 변하는 성질을 반영하기 위해, $n\times m$ matrix $S_k$와 $Y_k$를 다음과 같이 정의하고

$B^{(k)}_0=\delta_k I$에 $m$ 번의 업데이트가 적용된다 하면,

여기서,

은 $m\times m$ matrix이다. $x_{k+1}$이 생성되고 나면 $S_k$에서 $s_{k-m}$을 지우고 $s_k$를 추가해서 $S_{k+1}$를 얻는다 ($Y_{k+1}$, $L_{k+1}$, $D_{k+1}$도 비슷하게 얻음).

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

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

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

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

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

Unrolling the update

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

여기서,

Basic matrix $B^0_k$에 대해

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

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

Sparse quasi-Newton updates

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

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

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

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

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

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

Algorithms for partially separable functions

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

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

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

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

여기서 $f_i$는 각각 몇개의 $x$에만 의존한다. 그러면,

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

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

$f_1$에 집중하자. $f_1$은 $x$의 함수지만 $x_1$과 $x_3$에만 의존한다. 이를 $f_1$의 element variable이라 부른다. Element variable을 vector로 모아놓은 것을 다음 처럼 표기한다.

이를 $x$에 대해 표현하면,

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

그러면,

우리가 가정한 $f_1$의 경우

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

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

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

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

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

이다.