5. Conjugate Gradient Methods

개요

Large linear systems을 풀기에 가장 효율적인 방법이다 (이 목적으로 개발된 것이기도 하다). 그리고 이를 nonlinear optimization problem을 풀기위해 활용할 수 있다.

Linear system을 풀기 위해 사용될 때는 (Gaussian elimination을 대신하는 개념, linear conjugate gradient) 성능이 coefficient matrix의 eigenvalue의 분포에 따라 달라지는데, transform하거나, precondition 하는 방법으로 convergence를 향상시킬 수 있다.

Nonlinear conjugate gradient method로 large scale nonlinear optimization problem을 풀 수 있고, 이 때 장점은 steepest descent 보다 빠르고, matrix storage를 필요로 하지 않는다는 점이다.

The linear conjugate gradient method

Linear system of equation을 interative하게 푼다. 이 때, $A$는 $n\times n$ symmetric positive definite matrix이다. 이는 optimization 문제로 바꿔서 표현할 수 있다.

처음의 linear equation의 해와 위의 optimization 문제의 해는 같다. $\phi$의 gradient은 linear system의 residual과 같다.

Nonzero vectors $\{p_{0},p_{1},\dots,p_{l}\}$이 symmetric positive definite matrix $A$에 대해서 conjugate하다는 것은 다음을 만족하는 것이다.

이 vector들은 linearly independent하다는 것을 볼 수 있다. Conjugate gradient method의 장점은 conjugacy를 가지는 vector들의 set을 경제적으로 만들 수 있다는 것이다. Conjugacy를 가지고 있으면 $\phi(\cdot)$을 $n$ step 동안 각 direction에 따라 minimize 할 수 있다. 이를 보이기 위해 다음의 conjugate direction method를 살펴보자. 시작점 $x_{0}\in \mathbb{R}^{n}$과 conjugate direction $\{ p_{0},p_{1},\dots,p_{l-1}\}$에 대해 sequence $\{x_{k}\}$를 다음과 같이 얻을 수 있다.

이 때, $\alpha_{k}$는 $\phi(\cdot)$의 $x_{k}+\alpha p_{k}$ 를 따라 얻은 one-dimensional minimizer이고

이다.

Theorem 5.1

아무 시작점 $x_{0}\in \mathbb{R}^{n}$에 대해, conjugate direction method를 통해 얻어진 sequence $\{x_{k}\}$는 linear system의 solution $x^{\ast}$에 최대 $n$ step 안에 converge한다.

• 증명

  $\{p_{i}\}$는 서로 linearly independent 하기 때문에 $\mathbb{R}^{n}$을 span한다. 따라서

  $x^{\ast}-x_{0}=\sum^{n-1}_{i=0}\sigma_{i}p_{i}$

  로 표현할 수 있다. 양변의 앞에 $p^{T}_{k}A$를 곱하고 정리하면

  $\sigma_{k}=\frac{p^{T}_{k}A(x^{\ast}-x_{0})}{p^{T}_{k}Ap_{k}}$

  을 얻을 수 있다. $x_{k}$가 conjugate direction method를 통해 얻어지기 때문에

  $x_{k}=x_{0}+\sum^{k-1}_{i=0}\alpha_{i}p_{i}$

  이고, 양변의 앞에 $p^{T}_{k}A$를 곱하고 정리하면

  $p^{T}_{k}A(x_{k}-x_{0})=0$

  이고,

  $p^{T}_{k}A(x^{\ast}-x_{0})=p^{T}_{k}A(x^{\ast}-x_{k})=p^{T}_{k}(b-Ax_{k})=-p^{T}_{k}r_{k}$

  따라서 $\alpha_{k}=\sigma_{k}$이다.
  

Conjugate direction method의 $k$ iteration 후 얻어진 $x_{k}$는 $\text{span}\{p_{0},p_{1},\dots,p_{k-1}\}$에서 $\phi(\cdot)$이 minimize되는 점이다.

Theorem 5.2

아무 시작점 $x_{0}\in\mathbb{R}^{n}$에 대해, sequence $\{x_{k}\}$가 conjugate direction method를 통해 얻어진다고 가정하자. 그러면,

$r^{T}_{k}p_{i}=0 \quad \text{for }i=0,1,\dots,k-1$

이고, $x_{k}$는 $\phi(x)=\frac{1}{2}x^{T}Ax-b^{T}x$의 set

$\{x|x=x_{0}+\text{span}\{p_{0},p_{1},\dots,p_{k-1}\}\}$

에서의 minimizer이다.

• 증명

  먼저 점 $\tilde{x}$가 set $\{x|x=x_{0}+\text{span}\{p_{0},p_{1},\dots,p_{k-1}\}\}$에서 $\phi$를 minimize한다와 $i=0,1,\dots,k-1$에 대해 $r(\tilde{x})^{T}p_{i}=0$ 이다가 동치임을 보이자. $h(\sigma)=\phi(x_{0}+\sigma_{0}p_{0}+\cdots+\sigma_{k-1}p_{k-1})$ (이 때, $\sigma=(\sigma_{0},\dots,\sigma_{k-1})^{T}$을 정의하자. $h(\sigma)$는 strictly convex quadratic이기 때문에,

  $\frac{\partial h(\sigma^{\ast})}{\partial \sigma_{i}}=0 \quad i=0,1,\dots,k-1$

  을 만족하는 unique minimizer $\sigma^{\ast}$가 존재한다. Chain rule에 따라

  $\nabla \phi(x_{0}+\sigma^{\ast}_{0}p_{0}+\cdots + \sigma^{\ast}_{k-1}p_{k-1})^{T}p_{i}=0\quad i=0,1,\dots,k-1$

  이다. 이 때, $\nabla \phi(x)=Ax-b=r(x)$로 정의했었으니, minimizer $\tilde{x}=x_{0}+\sigma^{\ast}_{0}p_{0}+\cdots + \sigma^{\ast}_{k-1}p_{k-1}$에 대해서, $r(\tilde{x})^{T}p_{i}=0$임을 보인 것이다.

  이제 $x_{k}$가 $r^{T}_{k}p_{i}=0$ ($i=0,1,\dots,k-1$)을 만족함을 보이자. $k=1$인 경우, $x_{1}=x_{0}+\alpha_{0}p_{0}$이 $p_{0}$을 따라 $\phi$를 minimize한다는 점에서 $r^{T}_{1}p_{0}=0$이다. 이제 $i=0,1,\dots ,k-2$에 대해서 $r^{T}_{k-1}p_{i}=0$이라 가정하자.

  $r_{k}=r_{k-1}+\alpha_{k-1}Ap_{k-1}$

  의 앞에 $p^{T}_{k-1}$을 곱하면

  $p^{T}_{k-1}r_{k}=p^{T}_{k-1}r_{k-1}+\alpha_{k-1}p^{T}_{k-1}Ap_{k-1}=0$

  $\alpha_{k-1}$의 정의에 의해 $r^{T}_{k}p_{k-1}=0$을 얻는다. $i=0,1,\dots, k-2$인 $p_{i}$에 대해서는

  $p^{T}_{i}r_{k}=p^{T}_{i}r_{k-1}+\alpha_{k-1}p^{T}_{i}Ap_{k-1}=0$

  가정에 의해 $p^{T}_{i}r_{k-1}=0$ 이고, $p_{i}$의 conjugacy에 의해 $p^{T}_{i}Ap_{k-1}=0$이다. 따라서 $i=0,1,\dots,k-1$에 대해 $r^{T}_{k}p_{i}=0$이다.

$\{p_{i}\}$는 conjugacy를 가지기만 하면 어떤 것이든 상관이 없다. 가장 쉽게 생각해볼 수 있는 건 Eigenvector일텐데 계산 복잡도가 커서 현실적이지 못하다. Gram-Schmidt orthogonalization을 이용하는 것도 마찬가지이다.

Conjugate gradient method는 conjugate direction method의 한 종류이다. Conjugate vector $p_{k}$를 만드는 과정에서 $p_{k-1}$만을 이용하여 만드는데, $p_{k}$는 $p_{0},\dots,p_{k-2}$와 conjugate하다. 이 때문에, conjugate gradient method는 memory 사용과 computational cost 측면에서 우수하다. 각 direction $p_{k}$는

을 통해서 구해지는데, 이 때, $\beta_{k}$는 $p_{k-1}$와 $p_{k}$가 $A$에 대해서 conjugate 하도록 정한다.

$p_{0}$은 $x_{0}$에서의 steepest direction으로 잡는다. Conjugate direction method 처럼 각 direction에 대해 one-dimensional minimization을 하게 된다.

이제, $p_{0},p_{1},\dots ,p_{n-1}$이 서로 conjugate하고, residual $r_{i}$들이 서로 orthogonal하고, search direction $p_{k}$와 residual $r_{k}$가 다음과 같이 정의되는 $r_{0}$의 Krylov subspace of degree $k$에 포함됨을 보이자.

Theorem 5.3

Conjugate gradient method에 의해 $k$번째 생성된 점이 solution $x^{\ast}$가 아니라 가정하자. 다음의 네 성질이 성립한다.

$r^{T}_{k}r_{i}=0 \quad \text{for }i=0,1,\dots,k-1 \\ \text{span}\{r_{0},r_{1},\dots,r_{k}\}=\text{span}\{r_{0}, Ar_{0},\dots,A^{k}r_{0}\} \\ \text{span}\{ p_{0},p_{1},\dots, p_{k}\}=\text{span}\{r_{0},Ar_{0},\dots,A^{k}r_{0}\} \\ p^{T}_{k}Ap_{i}=0 \quad \text{for }i=0,1,\dots,k-1$

따라서, sequence $\{x_{k}\}$는 $x^{\ast}$로 최대 $n$ step 안에 수렴한다.

• 증명

  귀납적으로 증명하자. 2, 3번째 성질은 $k=0$에서 성립함을 쉽게 알 수 있다. 그리고 4번째 성질은 $k=1$일 때 성립한다. 이제 $k$에서 이 세 성질이 성립하면 $k+1$에서도 성립함을 보이자.

  2번째 성질이 성립함을 보이기 위해 먼저 좌변의 set이 우변의 set에 포함됨을 보이자. 먼저 $k$에서 성립함을 가정했으므로

  $r_{k}\in\text{span}\{r_{0}, Ar_{0},\dots A^{k}r_{0}\} \\ p_{k}\in\text{span}\{r_{0}, Ar_{0},\dots A^{k}r_{0}\}$

  이다. 앞에 $A$를 곱하면

  $Ap_{k}\in\text{span}\{Ar_{0}, A^{2}r_{0},\dots A^{k+1}r_{0}\}$

  이 때, $r_{k+1}=r_{k}+\alpha_{k}Ap_{k}$이므로,

  $r_{k+1}\in\text{span}\{r_{0}, Ar_{0},\dots A^{k+1}r_{0}\}$

  따라서

  $\text{span}\{r_{0},r_{1},\dots,r_{k+1}\}\sub\text{span}\{r_{0}, Ar_{0},\dots,A^{k+1}r_{0}\}$

  역이 성립함을 보이자. 먼저 $k$에서 성립함을 가정했으므로

  $A^{k}r_{0}\in\text{span}\{ p_{0},p_{1},\dots, p_{k}\}$

  앞에 $A$를 곱하면

  $A^{k+1}r_{0}\in\text{span}\{ Ap_{0},Ap_{1},\dots, Ap_{k}\}$

  업데이트 식에 의해 $Ap_{i}=(r_{i+1}-r_{i})/\alpha_{i}$ ($i=0,1,\dots, k$)이니까 결국

  $A^{k+1}r_{0}\in \text{span}\{r_{0},r_{1},\dots,r_{k+1}\}$

  이다. 따라서

  $\text{span}\{r_{0}, Ar_{0},\dots,A^{k+1}r_{0}\} \sub \text{span}\{r_{0},r_{1},\dots,r_{k+1}\}$

  2번째 성질에 대해서 증명하였다.

  이제 3번째 성질에 대해서 증명하자.

  $\text{span}\{p_{0},p_{1},\dots,p_{k},p_{k+1}\}=\text{span}\{p_{0},p_{1},\dots,p_{k},r_{k+1}\}\\=\text{span}\{r_{0}, Ar_{0},\dots,A^{k}r_{0},r_{k+1}\} \\=\text{span}\{r_{0},r_{1},\dots,r_{k},r_{k+1}\}\\=\text{span}\{r_{0}, Ar_{0},\dots,A^{k+1}r_{0}\}$

  이제 4번째 성질에 대해서 증명하자. Direction $p_{k+1}$를 구하는 식에 $Ap_{i}$를 곱하면

  $p^{T}_{k+1}Ap_{i}=-r^{T}_{k+1}Ap_{i}+\beta_{k+1}p^{T}_{k}Ap_{i} \tag{5.21}$

  이다. 여기서 $\beta_{k+1}$의 정의에 따라 $i=k$일 때, 오른쪽 term이 0이 된다. $i\leq k-1$에 대해서 따져보자. 먼저 $k$에 대해 네 성질들이 성립함을 가정했으므로 $p_{0},p_{1},\dots,p_{k}$는 conjugate 하고 Theorem 5.2에 의해서

  $r^{T}_{k+1}p_{i}=0 \quad \text{for }i=0,1,\dots,k$

  또,

  $Ap_{i}\in \text{ span}\{Ar_{0}, A^{2}r_{0},\dots,A^{i+1}r_{0}\}\\\sub \text{span}\{p_{0},p_{1},\dots,p_{i+1}\}$

  이므로 이 두 사실을 결합하면

  $r^{T}_{k+1}Ap_{i}=0\quad \text{for }i=0,1,\dots,k-1$

  따라서 (5.21)의 우변의 두 term이 모두 0가 되어 $k+1$에 대해 성립함을 보였다.

  마지막으로 1번 성질은 비 귀납적인 방법으로 증명한다. Direction set이 conjugate하므로 모든 $i=0,1,\dots,k-1$과 아무 $k=1,2,\dots,n-1$에 대해서 $r^{T}_{k}p_{i}=0$이다.

  $p_{i}=-r_{i}+\beta_{i}p_{i-1}$

  이므로, 모든 $i=0,1,\dots,k-1$에 대해서 $r_{i}\in\text{span}\{p_{i},p_{i-1}\}$이다. 따라서 모든 $i=1,\dots,k-1$에 대해서 $r^{T}_{k}r_{i}=0$이다. $p_{0}$의 정의에 의해 $r^{T}_{k}r_{0}=-r^{T}_{k}p_{0}=0$이다. 따라서 증명이 성립한다.

Theorem 5.2와 5.3의 결과를 이용하여 더 경제적인 conjugate gradient method를 도출할 수 있다. 먼저 Theorem 5.2의 $r^{T}_{k}p_{i}=0 \quad \text{for }i=0,1,\dots,k-1$을 이용해서 $\alpha_{k}$의 업데이트 식을 다음으로 대체할 수 있다.

두번째로, $\alpha_{k}Ap_{k}=r_{k+1}-r_{k}$임을 이용해서 $\beta_{k+1}$식을 단순화 할 수 있다.

요약하자면, 일반적으로 사용하는 conjugate gradient method는 다음의 알고리즘과 같다.

계산복잡도가 낮고, 메모리 효율적이나 문제의 크기가 작은 경우 round error가 적은 Gaussain elimination이나 SVD 같은 방식을 고려하는게 낫다.

기본적으로 정확히 계산했을 때 (수치해석적 오차를 고려하지 않았을 때) $n$ step 안에 solution을 찾는데, 만약에 $A$의 eigenvalue가 특정한 성질을 가지고 있다면 이보다 더 빠를 수 있다. 이에 대해 얘기하기 위해 Theorem 5.2를 살짝 다르게 바라보자.

Algorithm 5.2의 $x_{k+1}$ 업데이트 식과 theorem 5.3의 세번째 성질로부터

이 때, $\gamma_{i}$는 임의의 상수이다. $\gamma_{0},\gamma_{1},\dots,\gamma_{k}$를 coefficient로 가지는 degree $k$의 polynomial $P^{\ast}_{k}(\cdot)$을 정의하면

따라서 $x_{k+1}$의 식을 $P^{\ast}_{k}(\cdot)$을 이용해서 표현할 수 있다.

이제, 첫 $k$ step이 Krylov subspace $\mathcal{K}(r_{0};k)$으로 제한된 모든 가능한 method 중 Algorithm 5.2가 weighted norm measure $\|\cdot\|_{A}$ ($\|z\|_{A}=z^{T}Az$)로 측정된 거리 관점에서 $k$ step 후 최적임을 보이자.

Theorem 5.2에 의해 $x_{k+1}$은 set $x_{0}+\text{span}\{p_{0},p_{1},\dots,p_{k}\}$ ($x_{0}+\text{span}\{r_{0},Ar_{0},\dots,A^{k}r_{0}\}$과 같다.)에서 $\phi$를 minimize한다. (따라서 $\|x-x^{\ast}\|^{2}_{A}$도 minimize한다.). Polynomial $P^{\ast}_{k}$는 다음의 문제를 가능한 모든 degree $k$의 polynomial 공간에서 최소화하는 것이다.

이 섹션동안 이 optimality property를 반복적으로 이용할것이다.

이므로

이다. $0<\lambda_{1}\leq \lambda_{2}\leq\cdots\leq\lambda_{n}$을 $A$의 eigenvalue라 하고, $v_{1},v_{2},\dots,v_{n}$을 각 eigenvalue에 해당하는 eigenvector라고 하자. 그러면

Eigenvector가 $\mathbb{R}^{n}$ 전체를 span하므로 임의의 상수 $\xi_i$에 대해

이다. $A$의 eigenvector는 polynomial $P_{k}(A)$의 eigenvector이기도 하다. 그러므로,

이를 이용하면,

$\|z\|^{2}_{A}=z^{T}Az=\sum^{n}_{i=1}\lambda_{i}(v^{T}_{i}z)^{2}$이므로

CG로 얻은 polynomial $P^{\ast}_{k}$는 $\|\cdot\|_{A}$ 관점에서 최적이므로

$\lambda_{i}[1+\lambda_{i}P_{k}(\lambda_{i})]^{2}$ 중 가장 큰 term을 밖으로 빼내면

따라서 CG method의 convergence rate를 다음의 값을 가지고 나타낼 수 있다.

바꿔말하면, 이 식을 가능한 최소화하는 polynomial $P_{k}$를 탐색한다. 몇몇 경우 explicit하게 찾을 수 있다.

Theorem 5.4

만약 $A$가 $r$개에 서로 구분되는 eigenvalue를 가지고 있다면, CG는 최대 $r$ iteration안에 solution을 찾는다.

비슷한 추론 방식으로 얻을 수 있는 CG method의 다른 유용한 성질로 다음이 있다.

Theorem 5.5

만약 $A$가 eigenvalue $\lambda_1\leq \lambda_2\leq\cdots\leq\lambda_n$을 가지면,

$\|x_{k+1}-x^\ast\|^2_A\leq \left( \frac{\lambda_{n-k}-\lambda_1}{\lambda_{n-k}+\lambda_1} \right)^2\|x_0 -x^\ast\|^2_A$

이다. (증명 생략. (5.33)에서 유도됨.)

Theorem 5.5가 어떻게 CG의 동작들 예측하는지 보자.

그림처럼 Eigenvalue들이 cluster 되어있는 경우(특히 $n-m$개의 eigen value가 1주변에 모여있음에 주목)를 가정하자. Theorem 5.5에 따르면 $m+1$ step후 conjuagte gradient 알고리즘에 의해

따라서, $m+1$ step만으로도 충분히 좋은 solution을 얻을 수 있다.

CG의 convergence를 $A$의 Euclidean condition number

을 이용하여 표현하면

인데, 일반적으로 over estimation이 심하지만 $\lambda_1$과 $\lambda_n$ 밖에 모를 때는 유용하게 사용할 수 있다.

Preconditioning

Eigenvalue의 분포에 따라 CG의 수렴이 빨라질 수 있음을 알았다. 그렇다면 A를 적절히 transform해서 이런 성질을 유도할 수 없을까? 이를 위해 nonsingular matrix $C$를 이용해 $x$를 $\hat{x}$로 바꾼다.

그러면 처음 우리가 linear system을 optimization으로 바꾸기 위해 사용했던 식 $\phi$가

로 바뀐다. 그러면 CG의 convergence rate은 $A$대신 $C^{-T}AC^{-1}$의 eigenvalue에 따를 것이다. Eigenvalue가 잘 cluster 되거나, condition number가 작도록 $C$를 고르면 된다. 이런 transform을 explicit하게 하지는 않고 \hat{x}에 대해 CG를 풀고 transform을 되돌릴 수 있다. matrix $M=C^TC$를 이용하여

$M=I$일 때 standard CG이다. 참고로 orthogonality property는 이렇게 표현된다.

Preconditioned CG에서는 (5.39d)의 계산이 추가됐다.

Practical Preconditioners

통용되는 최고의 preconditioning 방법이 있지는 않다. $M$의 효과, $M$의 저장과 계산복잡도, (5.39d)의 계산 복잡도를 고려해야한다.

PDE 풀기 위한 방법은 연구 많이 됐다는데 잘 이해 못해서 생략.

General-purpose preconditioner 중에 제일 대중적인 건 incomplete Cholesky 가 있다. 정확한 Cholesky decomposition $A=LL^T$를 하는 대신 $L$보다 sparse한 approximation $\tilde{L}$을 이용하여 $A\approx \tilde{L}\tilde{L}^T$을 만드는 방식이다. $C=\tilde{L}^T$로 잡으면 $M=\tilde{L}\tilde{L}^T$이고

가 돼서 $C^{-T}AC^{-1}$의 eigen value들이 1 주위로 모이게 된다. M을 직접적으로 계산하지 않고 $\tilde{L}$을 저장하고 $My=r$을 두번의 triangular substitution으로 푼다. $\tilde{L}$의 sparsity가 A랑 비슷하므로 $My=r$을 푸는 비용은 Ap를 계산하는 것과 비슷하다.

단, matrix $A$가 충분히 positive definite하지 않아 $\tilde{L}$을 구하기 위해서는 diaognal element의 값을 증가시켜야 할 수도 있고, sparsity condition 때문에 numerically instable 할 수 있다.

Nonlinear conjugate gradient method

CG를 앞서 언급한 $\phi$ 말고 일반적인 (nonlinear까지) function을 minimize하는 데 사용할 수는 없을까?

Fletcher-Reeves Method

CG 알고리즘에서 1) step length $\alpha_k$를 line search를 통해 찾고, 2) residual $r$을 $f$의 gradient로 수정한 버전이다.

$f$가 strongly convex quadratic이고 $\alpha_k$가 exact minimizer면 linear CG가 된다. FR은 $f$의 값과 gradient 만 계산하면 되고 matrix operation도 필요 없고 vector 몇개만 저장하면 돼서 좋다.

다만 $\alpha_k$를 어떻게 계산하냐에 따라서 $p_{k+1}$이 descent direction이 아닐 수 있다. (5.41b)에 $\nabla f_{k}$를 곱하면

이다. 만약 line search가 정확했다면 $\alpha_{k-1}$이 $p_{k-1}$ 방향으로 minimizer라 $\nabla f^T_k p_{k-1}=0$이고 따라서 $\nabla f^T_kp_k<0$이고 $p_k$는 descent direction이다. Line search가 정확하지 않은 경우에는 $p_k$가 descent direction이라고 확신할 수 없는데, 이를 $\alpha_k$가 strong Wolfe condition을 만족하도록 요구하면 해결된다.

Polak-Ribiere Method와 변형

Fletcher-Reeves에서 $\beta_k$를 바꾼 변형들이 여럿 있는데 그 중 하나가 Polak-Ribiere이다.

$f$가 strongly convex quadratic이고 line search가 exact할 때 FR과 동일하다. Line search가 정확하지 않을 때 PR이 더 robust하고 효율적인 편이라 한다.

PR에서는 strong Wolfe condition이 $p_k$가 descent direction임을 보장해주지 않는다. 이 문제를 해결한 버전이 PR+ 이다.

Hestenes-Stiefel 도 $f$가 strongly convex quadratic이고 line search가 exact할 때 FR과 동일하다.

HS는 연속된 search direction들이 line segment $[x_{k}, x_{k+1}]$의 average Hessian에 대해 conjuagte 하도록 만드는 것에서 유도됐다. Average Hessian은

Taylor’s theorem에 의하면 $\nabla f_{k+1}=\nabla f_k +\alpha_{k}\bar{G}_kp_k$ 인데 $p_{k+1}=-\nabla f_{k+1}+\beta_{k+1}p_k$ 꼴로 표현되는 어느 direction이든, $p^T_{k+1}\bar{G}_kp_k=0$ 조건을 달면 $\beta_{k+1}$이 HS의 꼴이 되어야 한다.

나중에 보게되겠지만, 모든 $k\geq 2$에 대해

를 만족하는 모든 $\beta_{k}$는 global convergence를 보장한다. 이 성질을 이용하기 위해 PR을 수정하면

이를 FR-PR이라 한다.

최근에 (이라지만 각각 1999년A nonlinear conjugate gradient method with a strong global convergence property, 2005년A new conjugate gradient method with guaranteed descent and an efficient line search이다. 아무래도 CG가 1950년대 나와서 상대적 최근인듯) 나온 방식 두개가 있는데 PR과 견줄만하다 한다.

1999년 방식

2005년 방식

Quadratic Termination and Restarts

Restart 라는 방식은 $n$ step 마다 $\beta_k=0$으로 설정해서 steepest descent step을 택하고 nonlinear conjugate gradient를 재시작하는 방식이다. f가 solution의 neighborhood에서 strongly convex quadratic하지만 다른 곳에선 nonquadratic하다고 가정하면, solution에 다가가다가 strongly convex quadratic 한 region에서 nonlinear CG를 재시작할 것이고, 그렇게 되면 linear CG처럼 동작하고 quadratic하게 converge 할 것이다. 다시말해,

그런데 CG는 일반적으로 $n$이 클 때 사용하는데 그런 경우 restart는 현실적으로 불가능하다. 그래서 일반적으로 restarts를 안하거나 consecutive gradient가 orthogonal에서 많이 벗어났는지 확인하고 결정한다. 일반적으로 $v=0.1$로 설정하고

Fletcher-Reeves Method 분석

아까 증명없이 언급했던 global convergence에 대해 얘기해보자.

Lemma 5.6

Strong Wolfe condition이 $0<c_2<\frac{1}{2}$로 만족하도록 하는 step length $\alpha_k$로 Fletcher-Reeves 알고리즘이 구현됐다고 가정하자. FR은 다음의 부등식을 만족하는 descent direction $p_k$를 만든다.

$-\frac{1}{1-c_2}\leq \frac{\nabla f^T_k p_k}{\| \nabla f_k\|^2} \leq \frac{2c_2-1}{1-c_2}, \forall k=0,1,\dots$

• 증명

  $t(\xi)=(2\xi-1)/(1-\xi)$가 $[0,\frac{1}{2}]$에서 단조증가하고 $t(0)=-1$, $t(\frac{1}{2})=0$이라고 하자. $c_2\in(0, \frac{1}{2})$니까

  $-1<\frac{2c_2-1}{1-c_2}<0$

  이다. 따라서 descent condition $\nabla f^T_kp_k<0$이 성립한다.

  귀납적으로 증명하자. $k=0$에 대해서 steepest descent일 테니 $\nabla f^T_k p_k=-1$이고 부등식이 만족한다. 이제 $k$에 대해서 부등식이 만족한다 가정하면 $k+1$에 대해서

  $\frac{\nabla f^T_{k+1}p_{k+1}}{\|\nabla f_{k+1}\|^2}= -1 + \beta_{k+1}\frac{\nabla f^T_{k+1}p_{k}}{\|\nabla f_{k+1}\|^2}= -1 + \frac{\nabla f^T_{k+1}p_{k}}{\|\nabla f_{k}\|^2}$

  인데, Wolfe condition에 의해

  $|\nabla f^T_{k+1}p_k|\leq -c_2 \nabla f^T_k p_k$

  이니 잘 합치면

  $-1 + c_2\frac{\nabla f^T_{k}p_{k}}{\|\nabla f_{k}\|^2}\leq \frac{\nabla f^T_{k+1}p_{k+1}}{\|\nabla f_{k+1}\|^2}\leq -1 -c_2 \frac{\nabla f^T_{k}p_{k}}{\|\nabla f_{k}\|^2}$

  $k$에 대해 부등식이 성립한다 했으니 $\frac{\nabla f^T_{k}p_{k}}{\|\nabla f_{k}\|^2}$를 치환하면

  $-1-\frac{c_2}{1-c_2}\leq \frac{\nabla f^T_{k+1}p_{k+1}}{\|\nabla f_{k+1}\|^2} \leq -1+\frac{c_2}{1-c_2}$

  따라서 $k+1$에서도 성립한다.

증명 과정에서 strong Wolfe condition의 한가지 성질만 사용했는데 다른 성질을 이용하면 global convergence도 증명할 수 있다. 우선 여기선 $\nabla f^T_k p_k$의 bound가 step의 norm $\|p_k\|$가 커지는 속도에 제약을 건다는 점을 주목하자.

Lemma 5.6은 FR의 단점도 설명하는데, 한번 direction과 step length가 좋지 못했으면, 다음 step에서도 그럴 가능성이 크다.

$p_k$가 좋지 못하다는 것은 descent direction으로부터 멀어져 $\cos \theta_k \approx 0$이라는 것이다. Lemma 5.6의 부등식에 $\|\nabla f_k \| / \|p_k\|$를 곱하면

이 부등식에 따라 $\cos \theta_k \approx 0$은

와 동치이다. $p_k$가 gradient와 orthogonal에 가깝다면 step도 작을것이고 $x_{k+1}\approx x_k$라 한다면 $\nabla f_{k+1}\approx \nabla f_{k}$일테니

이 될 것이고, 그렇다면 $\|\nabla f_{k+1}\| \approx \|\nabla f_{k}\| \ll \|p_k\|$라서

가 된다. 새로운 search direction도 비슷할 것이다. PR의 경우는 같은 상황에서 $\beta^{PR}_{k+1}\approx 0$이 되서 사실상 restart를 하게된다. 앞서 언급한 다른 modification들의 경우도 마찬가지라 이런 문제가 생기지 않는다.

Global Convergence

우선 논의를 위해 objective function이 다음의 성질들을 가진다고 가정하자.

1. Level set $\mathcal{L}=\{x|f(x)\leq f(x_0)\}$이 bounded 되어 있다.
2. $\mathcal{L}$의 어떤 open neighborhood $\mathcal{N}$에 대해 objective function $f$가 Lipschitz continuously differentiable하다.

가정에 의하면

을 만족하는 상수 $\bar{\gamma}$가 존재한다.

Zoutendijk’s theorem에 따라 Wolfe condition을 만족하면

이다. 이 것을 이용하면 주기적으로 restart하는 알고리즘들의 global convergence를 증명할 수 있다. $k_i$를 $i$번째 재시작할 때의 index라고 하고, 재시작 할때는 steepest direction임을 (따라서 $\cos\theta_k=1$) 이용하면 $\lim_{j\rarr \infty} \|\nabla f_{k_j}\|=0$이므로

이다.

이제 restart하지 않은 경우를 분석해보자.

Theorem 5.7

앞의 두 가정이 성립하고 FR이 $0<c_1<c_2<\frac{1}{2}$로 strong Wolfe condition을 만족한다 하면,

$\liminf_{k\rarr \infty} \|\nabla f_k\|=0$

• 증명

  귀류적으로 증명하자. 먼저 충분히 큰 모든 $k$에 대해서

  $\|\nabla f_k\|\geq \gamma$

  를 만족하는 $\gamma>0$이 존재한다고 가정하자. Lemma 5.6의 부등식을 Zoutendijk’s condition에 대입하면

  $\sum^\infty_{k=0}\frac{\|\nabla f_k\|^4}{\|p_k\|^2}<\infty$

  Wolfe condition과 Lemma 5.6의 부등식에 의해

  $|\nabla f^T_{k}p_{k-1}|\leq -c_2\nabla f^T_{k-1}p_{k-1}\leq \frac{c_2}{1-c_2}\|\nabla f_{k-1}\|^2$

  FR 알고리즘의 계산에 의해 (일부 과정 생략)

  $\|p_k\|^2\leq \left(\frac{1+c_2}{1-c_2}\right)\|\nabla f_k\|^2+(\beta^{FR}_k)^2\|p_{k-1}\|^2$

  이 부등식을 반복해서 적용하고 $c_3=(1+c_2)/(1-c_2)\geq 1$로 정의하면

  $\|p_k\|^2\leq c_3\|\nabla f_k\|^4 \sum^k_{j=0}\|\nabla f_j\|^{-2}$

  이다. Assumtion에 의한 bound와 귀류법을 위해 만든 가정에 의해

  $\|p_k\|^2\leq \frac{c_3 \bar{\gamma}^4}{\gamma^2}k$

  이기 때문에, positive 상수 $\gamma_4$에 대해

  $\sum^\infty_{k=1}\frac{1}{\|p_k\|^2}\geq \gamma_4 \sum^\infty_{k=1}\frac{1}{k}$

  이다. 그런데

  $\sum^\infty_{k=1}\frac{1}{\|p_k\|^2}<\infty$

  이므로

  $\sum^\infty_{k=1}\frac{1}{k}<\infty$

  여야 한다. 이는 사실이 아니므로 모순이다.
  

일반적인 (nonconvex까지) 함수의 경우에는 Theorem 5.7같은 증명이 불가능하다. PR이 FR보다 잘 동작한다는 점을 생각해봤을 때 이는 다소 충격적이다. 다음 theorem은 ideal line search인 경우에도 PR이 solution point에 다가가지 못하고 영원히 cycle 돌수도 있음을 보여준다.

Theorem 5.8

Ideal line search를 사용하는 PR을 생각해보자. Gradient의 sequence $\{\|\nabla f_k\|\}$이 0으로부터 bounded away 되도록 만드는 twice continuously differentiable한 function $f:\mathbb{R}^3\rarr \mathbb{R}$과 시작점 $x_0\in\mathbb{R}^3$이 존재한다.

PR+를 이용하면 Theorem 5.8의 문제를 피할 수 있고 global convergence가 증명 가능하다고 한다.