NOTES | HOME
$$ \newcommand{\RR}{\mathbb{R}} \newcommand{\GG}{\mathbb{G}} \newcommand{\PP}{\mathbb{P}} \newcommand{\PS}{\mathcal{P}} \newcommand{\SS}{\mathbb{S}} \newcommand{\NN}{\mathbb{N}} \newcommand{\ZZ}{\mathbb{Z}} \newcommand{\CC}{\mathbb{C}} \newcommand{\HH}{\mathbb{H}} \newcommand{\ones}{\mathbb{1\hspace{-0.4em}1}} \newcommand{\alg}[1]{\mathfrak{#1}} \newcommand{\mat}[1]{ \begin{pmatrix} #1 \end{pmatrix} } \renewcommand{\bar}{\overline} \renewcommand{\hat}{\widehat} \renewcommand{\tilde}{\widetilde} \newcommand{\inv}[1]{ {#1}^{-1} } \newcommand{\eqdef}{\overset{\text{def}}=} \newcommand{\block}[1]{\left(#1\right)} \newcommand{\set}[1]{\left\{#1\right\}} \newcommand{\abs}[1]{\left|#1\right|} \newcommand{\trace}[1]{\mathrm{tr}\block{#1}} \newcommand{\norm}[1]{ \left\| #1 \right\| } \newcommand{\argmin}[1]{ \underset{#1}{\mathrm{argmin}} } \newcommand{\argmax}[1]{ \underset{#1}{\mathrm{argmax}} } \newcommand{\st}{\ \mathrm{s.t.}\ } \newcommand{\sign}[1]{\mathrm{sign}\block{#1}} \newcommand{\half}{\frac{1}{2}} \newcommand{\inner}[1]{\langle #1 \rangle} \newcommand{\dd}{\mathrm{d}} \newcommand{\ddd}[2]{\frac{\partial #1}{\partial #2} } \newcommand{\db}{\dd^b} \newcommand{\ds}{\dd^s} \newcommand{\dL}{\dd_L} \newcommand{\dR}{\dd_R} \newcommand{\Ad}{\mathrm{Ad}} \newcommand{\ad}{\mathrm{ad}} \newcommand{\LL}{\mathcal{L}} \newcommand{\Krylov}{\mathcal{K}} \newcommand{\Span}[1]{\mathrm{Span}\block{#1}} \newcommand{\diag}{\mathrm{diag}} \newcommand{\tr}{\mathrm{tr}} \newcommand{\sinc}{\mathrm{sinc}} \newcommand{\cat}[1]{\mathcal{#1}} \newcommand{\Ob}[1]{\mathrm{Ob}\block{\cat{#1}}} \newcommand{\Hom}[1]{\mathrm{Hom}\block{\cat{#1}}} \newcommand{\op}[1]{\cat{#1}^{op}} \newcommand{\hom}[2]{\cat{#1}\block{#2}} \newcommand{\id}{\mathrm{id}} \newcommand{\Set}{\mathbb{Set}} \newcommand{\Cat}{\mathbb{Cat}} \newcommand{\Hask}{\mathbb{Hask}} \newcommand{\lim}{\mathrm{lim}\ } \newcommand{\funcat}[1]{\left[\cat{#1}\right]} \newcommand{\natsq}[6]{ \begin{matrix} & #2\block{#4} & \overset{#2\block{#6}}\longrightarrow & #2\block{#5} & \\ {#1}_{#4} \hspace{-1.5em} &\downarrow & & \downarrow & \hspace{-1.5em} {#1}_{#5}\\ & #3\block{#4} & \underset{#3\block{#6}}\longrightarrow & #3\block{#5} & \\ \end{matrix} } \newcommand{\comtri}[6]{ \begin{matrix} #1 & \overset{#4}\longrightarrow & #2 & \\ #6 \hspace{-1em} & \searrow & \downarrow & \hspace{-1em} #5 \\ & & #3 & \end{matrix} } \newcommand{\natism}[6]{ \begin{matrix} & #2\block{#4} & \overset{#2\block{#6}}\longrightarrow & #2\block{#5} & \\ {#1}_{#4} \hspace{-1.5em} &\downarrow \uparrow & & \downarrow \uparrow & \hspace{-1.5em} {#1}_{#5}\\ & #3\block{#4} & \underset{#3\block{#6}}\longrightarrow & #3\block{#5} & \\ \end{matrix} } \newcommand{\cone}[1]{\mathcal{#1}} $$

Conjugate Gradient

Over time, I found myself checking the Wikipedia page for Conjugate Gradient over and over again, and finally got tired of it. So here is a complete, self-contained and hopefully correct derivation of the method, including non-standard inner products and preconditioning, up to the Conjugate Residuals variations. The alternate Lanczos formulation can be found in the notes for Krylov methods.

  1. Gradient Descent
  2. Conjugate Directions
  3. Standard Algorithm
  4. (TODO) Convergence Analysis
  5. Non-Standard Inner Product
    1. Breakdown
  6. Examples
    1. Preconditioned Conjugate Gradient
    2. Conjugate Residuals
    3. Preconditioned Conjugate Residuals
    4. Preconditioned Conjugate Residuals (alternate)
    5. Mixed Conjugate Gradient/Residuals

Gradient Descent

When a matrix \(A\) is positive definite, the solution to \(Ax = b\) is given by the unique minimizer of the following convex quadratic form:

\[E(x) = \half x^T A x - b^T x\]

The gradient descent algorithm proceeds by minimizing \(E\) along successive descent directions \(p_k\), given in this case by the negative gradient \(r_k = b - A x_k\) of \(E\):

\[x_{k+1} = x_k + \alpha_k p_k\]

where \(p_k = r_k\) is known as the steepest descent direction. \(\alpha_k\) is found by the following line-search procedure, a one-dimensional minimization problem (convex in this case):

\[\alpha_k = \argmin{\alpha \in \RR} \quad f_k(\alpha) = E\block{x_k + \alpha p_k}\]

Stationarity conditions give:

\[\begin{align} \dd f_k(\alpha) = 0 &\iff \block{x_k + \alpha p_k}^T A p_k - b^Tp_k = 0 \\ & \iff \quad p_k^T A p_k = \block{b - A x_k}^T p_k \\ & \iff \quad \alpha = \frac{\block{b - A x_k}^T p_k}{p_k^T A p_k} \end{align}\]

Geometrically, the line-search can also be interpreted in terms of the \(A\)-projection of \(x_k\) onto \(p_k\):

\[x_{k+1} = x_k - p_k \frac{p_k^TAx_k}{p_k^TAp_k} \quad + \quad p_k \frac{p_k^T b}{p_k^TAp_k}\]

At each iteration, the \(A\)-projection of \(x_k\) along \(p_k\) is replaced with that of the solution \(\inv{A}b = x^\star\) to obtain \(x_{k+1}\). Consequently, an interesting choice for descent directions \(\block{p_k}_k\) is to use mutually \(A\)-orthogonal (“\(A\)-conjugate”, or simply “conjugate”) directions, so that the solution is found in maximum \(n\) steps in exact arithmetic. This is precisely the idea behind the Conjugate Gradient (CG) method.

Conjugate Directions

The CG algorithm chooses \(p_{k+1}\) by \(A\)-orthogonalizing the gradient \(r_{k+1}\) against the previous \(p_k\), so that the family \(\block{p_k}_k\) is \(A\)-orthogonal:

\[p_{k+1} = r_{k+1} - \sum_{i=1}^k p_i \frac{p_i^TAr_{k+1}}{p_i^T A p_i}\]

We let \(\beta_{ik} = p_i^TAr_k\) if \(i < k\), \(\beta_{ii} = 1\) and \(\beta_{ik} = 0\) if \(i > k\), so that the recurrence relation becomes \(r_{k+1} = \sum_{i=1}^{i = k+1} p_i \beta_{ik}\) or, in matrix form:

\[R_k = P_k B_k\]

where the matrix \(B_k = \block{\beta_{ij}}_{i, j \leq k}\) is upper-triangular with unit diagonal (hence invertible). Now, since the \(p_i\) are now mutually \(A\)-orthogonal, and from the geometric interpretation above (or simply by induction over the update formula for \(x_k\)), the \(A\)-projection of \(x_{k+1}\) over the previous \(\block{p_i}_{i\leq k} = P_k\) gives:

\[P_k^T A x_{k+1} = P_k^T b\]

Grouping both sides, we obtain:

\[\underbrace{P_k^T}_{B_k^{-T} R_k^T} r_{k+1} = 0\]

Which implies the following important property:

\[R_k^T r_{k+1} = 0\]

That is, the \(r_k\) are mutually orthogonal. Before moving on, we also note that \(p_k^T r_k = \norm{r_k}^2\). At this point, let us rewrite the \(r_k\) update equations in block form, to get a clearer view of the situation. From each \(r_k - r_{k+1} = \alpha_k A p_k\), we obtain:

\[R_{k+1} \underbrace{\mat{ 1& \\ -1 & 1\\ & -1 & 1 \\ & & \ddots & \ddots \\ & & & -1 & 1 \\ & & & & -1 \\ }}_{\underline{D}_k} = A P_k \block{\alpha_k}_k\]

Projecting onto \(R_k\), and introducing \(\phi_k = \norm{r_k}^2\), gives:

\[\begin{align} R_k^T R_{k+1} \underline{D}_k &= R_k^T A P_k \block{\alpha_k}_k \\ \mat{ \diag\block{\phi_k} & 0_k} \underline{D}_k &= R_k^T A P_k \block{\alpha_k}_k \\ \diag\block{\phi_k} D_k &= R_k^T A P_k \block{\alpha_k}_k \\ \end{align}\]

where \(D_k\) is the upper \(k\times k\) (upper Hessenberg) block from \(\underline{D}_k\). More precisely:

\[\diag\block{\phi_k} D_k = \mat{ \phi_1 & \\ -\phi_2 & \phi_2\\ & -\phi_3 & \phi_3 \\ & & \ddots & \ddots \\ & & & -\phi_k & \phi_k \\ }\]

From this, we finally get:

\[r_i^T A p_j \alpha_j = \cases{ \phi_i & if $i = j$ \\ -\phi_i & if $j = i - 1$ \\ 0 & otherwise }\]

Since \(\alpha_k \neq 0\) (unless we are finished), \(p_i^T A r_{k+1} = 0\) for \(i < k\) and the update for \(p_{k+1}\) reduces to:

\[p_{k+1} = r_{k+1} - p_k \frac{p_k^TAr_{k+1}}{p_k^T A p_k}\]

This is called a short-recurrence relation: each \(p_k\) only needs to be conjugated with the previous one. Some alternate formula before finishing:

\[\begin{align} \beta_k &= \frac{p_k^TAr_{k+1}}{p_k^T A p_k} \\ &= -\frac{\phi_{k+1}}{\alpha_k p_k^T A p_k } = -\frac{\phi_{k+1}}{r_k^T p_k} \\ &= -\frac{\phi_{k+1}}{\phi_k} \\ \\ \alpha_k &= \frac{r_k^T p_k}{p_k^T A p_k} \\ &= \frac{\phi_k}{p_k^T A p_k} \\ \end{align}\]

Standard Algorithm

\[p_0 = r_0 = b - Ax_0\] \[\begin{align} x_{k+1} &= x_k + \alpha_k p_k &\alpha_k = \frac{r_k^Tp_k}{p_k^T A p_k} &= \frac{r_k^T r_k}{p_k^T A p_k} \\ r_{k+1} &= r_k - \alpha_k Ap_k & & \\ p_{k+1} &= r_{k+1} - \beta_k p_k &\beta_k = \frac{r_{k+1}^T Ap_k}{p_k^T A p_k} &= -\frac{ r_{k+1}^T r_{k+1} }{r_k^T r_k} \\ \end{align}\]

While the initial search direction \(p_0\) is set to the initial residual \(r_0\), any other initial direction could be chosen. Some references mention, however, that the speed of convergence can be affected ( i.e. linear instead of superlinear) when the initial residual is not chosen for \(p_0\).

The normalized \(r_k\) provide an orthogonal basis for the Krylov subspace \(\Krylov_k(A, b)\) and can be related to the Lanczos vectors. If \(p_0 = r_0\), then the \(p_k\) vectors also provide an \(A\)-orthogonal basis for \(\Krylov_k(A, b)\).

(TODO) Convergence Analysis

Non-Standard Inner Product

Let \(A\) be self-adjoint and positive definite with respect to an inner product \(M\), that is:

(or simply: let \(MA\) be symmetric positive definite). The following quadratic form is positive-definite (thus convex):

\[\begin{align} E_M(x) &= \half \inner{x, A x - b}_M \\ &= \half x^TMAx - x^TMb \\ \end{align}\]

So again, we look for its minimum. We apply the CG algorithm to the modified system \(MA x = Mb\), letting \(z_k = M r_k\) with \(r_k = b - Ax_k\) as before, and denoting the \(MA\)-conjugate directions by \(q_k\). We again obtain the following conditions for step length and conjugation:

Now we are facing a choice regarding whom should span the \(q_k\). If the \(q_k\) are spanned by the residuals \(r_k\) (i.e. the gradient of \(E_M\) with respect to \(M\)), we get:

In contrast, if the \(q_k\) are spanned by the preconditioned residuals \(z_k\) (i.e. the gradient of \(E_M\) with respect to the standard inner product), we get:

which is vanilla CG on the preconditioned system. Using the gradient with respect to \(M\) provides an extra degree of freedom in the algorithm: intuitively, the gradient descent follows directions that are orthogonal to level-sets in order to minimize the function, and orthogonality depends on the metric. Different metrics produce different trajectories across level-sets, hopefully to converge faster towards the solution. Since the second case was already described above, let us now consider the first and obtain the following algorithm:

\[p_0 = r_0 = b - Ax_0\] \[\begin{align} x_{k+1} &= x_k + \alpha_k p_k & \alpha_k &= \frac{p_k^TMr_k}{p_k^TMAp_k} &= \frac{r_k^TMr_k}{p_k^TMAp_k}\\ r_{k+1} &= r_k - \alpha_k Ap_k & & & \\ p_{k+1} &= r_{k+1} - \beta_k p_k & \beta_k &= \frac{r_{k+1}^TMAp_k}{p_k^TMAp_k} &= -\frac{ r_{k+1}^T M r_{k+1} }{r_k^T M r_k} \\ \end{align}\]

which is exactly the standard CG algorithm using \(M\) for all inner-products, as one could expect. At each iteration of the algorithm:

Breakdown

The above requires that \(MA\) be symmetric, which could happen even though \(M\) is indefinite, or not even symmetric. However, \(r_k\) is not guaranteed to be a descent direction of the quadratic form when \(M\) is not positive definite: it might happen that the quadratic form simply stagnates along \(r_{k+1} \neq 0\):

\[r_{k+1}^T z_{k+1} = r_{k+1}^T M r_{k+1} = 0\]

Should this happen, the subsequent \(p_{k+1}\) will also be a direction of stagnation since \(P_k^T M r_{k+1} = 0\) always holds, resulting in:

\[p_{k+1}^T M r_{k+1} = 0\]

At this point, the algorithm will stop making any progress since \(\alpha_{k+1}\) will be zero and \(r_{k+2} = r_{k+1}\). It seems that it is necessary and sufficient that the symmetric part of the metric is positive definite in order that the above breakdown scenario can never happen.

Examples

It is now straightforward to derive most classical symmetric Krylov methods as instances of the conjugate gradient algorithm for a suitable choice of metric \(M\).

Preconditioned Conjugate Gradient

The matrix \(\inv{B}A\) is trivially positive definite for an inner-product \(B\), and we immediately obtain the Preconditioned Conjugate Gradient algorithm. The \(p_k\) will be \(A\)-conjugate, and the energy norm will be minimized along the way, but the residuals \(r_k\) will now be \(B\)-conjugate, hopefully to follow a faster path towards the solution.

\[p_0 = r_0 = \inv{B}\block{b - Ax_0}\] \[\begin{align} x_{k+1} &= x_k + \alpha_k p_k & \alpha_k &= \frac{p_k^T B r_k}{p_k^TAp_k} \\ r_{k+1} &= r_k - \alpha_k \inv{B}Ap_k & & \\ p_{k+1} &= r_{k+1} - \beta_k p_k & \beta_k &= \frac{r_k^TAp_k}{p_k^TAp_k} \\ \end{align}\]

It is generally more convenient to use the residual for the original system instead, and introduce an extra variable \(z = \inv{B} r\) as follows:

\[r_0 = b - A x_0, \quad p_0 = z_0 = \inv{B} r_0\] \[\begin{align} x_{k+1} &= x_k + \alpha_k p_k & \alpha_k &= \frac{p_k^Tr_k}{p_k^TAp_k} &= \frac{z_k^Tr_k}{p_k^TAp_k} \\ r_{k+1} &= r_k - \alpha_k Ap_k & & \\ z_{k+1} &= \inv{B} r_{k+1} \\ p_{k+1} &= z_{k+1} - \beta_k p_k & \beta_k &= \frac{z_{k+1}^TAp_k}{p_k^TAp_k} &= -\frac{z_{k+1}^Tr_{k+1}}{z_k^T r_k} \\ \end{align}\]

Intuitively, the metric \(B\) should ease displacements where the curvature is low, and damp it where the curvature is high. By assigning a high cost to highly-curved directions, the gradient in these directions will become small (since you don’t have to move a lot to get energy decrease, due to the high associated energy cost), and conversely. In this respect, the diagonal metric \(B = \diag(A)\) approximates the curvature along each axis (similar to the Levernberg-Marquardt algorithm), and damps the corresponding displacements accordingly.

Conjugate Residuals

This is exactly CG on \(A\) with non-standard metric \(M = A^T\). The residual norm is minimized, and the residuals are \(A^T\)-conjugate.

\[p_0 = r_0 = b - Ax_0\] \[\begin{align} x_{k+1} &= x_k + \alpha_k p_k & \alpha_k &= \frac{p_k^T A^T r_k}{p_k^TA^TAp_k} &= \frac{r_k^TA^Tr_k}{p_k^TA^TAp_k}\\ r_{k+1} &= r_k - \alpha_k Ap_k & & \\ p_{k+1} &= r_{k+1} - \beta_k p_k & \beta_k &= \frac{r_{k+1}^TA^TAp_k}{p_k^TA^TAp_k} &= -\frac{ r_{k+1}^T A^T r_{k+1} }{r_k^T A^T r_k} \\ \end{align}\]

This algorithm yields the same iterates as MINRES (and is cheaper to compute). \(Ap\) can be obtained from \(Ar\) at each iteration to keep only one multiplication by \(A\) per iteration.

Preconditioned Conjugate Residuals

This is the one found on wikipedia. Again, this is CG on \(\inv{B}A\) with metric \(A^T\). The \(\inv{B}\)-norm of the residual is minimized, and the residuals are again \(A^T\)-conjugate.

\[p_0 = r_0 = \inv{B}\block{b - Ax_0}\] \[\begin{align} x_{k+1} &= x_k + \alpha_k p_k & \alpha_k &= \frac{p_k^T A^T r_k}{p_k^TA^T\inv{B} Ap_k} &= \frac{r_k^T A r_k}{p_k^TA^T\inv{B} Ap_k}\\ r_{k+1} &= r_k - \alpha_k \inv{B} Ap_k & & \\ p_{k+1} &= r_{k+1} - \beta_k p_k & \beta_k &= \frac{r_{k+1}^TA^T \inv{B} A p_k}{p_k^TA^T \inv{B} A p_k} &= -\frac{ r_{k+1}^T A r_{k+1} }{r_k^T A r_k} \\ \end{align}\]

\(Ap\) can be obtained from \(Ar\) at each iteration to keep only one multiplication by \(A\) per iteration.

Preconditioned Conjugate Residuals (alternate)

This one is CG on \(\inv{B}A\) using metric \(A^TB\). In praticular, \(A\) does not need to be symmetric, only \(A^TB\) has to be positive definite. The norm of the residual is minimized, and the residuals are \(A^TB\)-conjugate:

\[p_0 = r_0 = b - \inv{B}Ax_0\] \[\begin{align} x_{k+1} &= x_k + \alpha_k p_k & \alpha_k &= \frac{p_k^T A^T Br_k}{p_k^TA^TAp_k} &= \frac{r_k^TA^TBr_k}{p_k^TA^TAp_k}\\ r_{k+1} &= r_k - \alpha_k \inv{B}Ap_k & & & \\ p_{k+1} &= r_{k+1} - \beta_k p_k & \beta_k &= \frac{r_{k+1}^TA^TAp_k}{p_k^TA^TAp_k} &= -\frac{ r_{k+1}^T A^TB r_{k+1} }{r_k^T A^TB r_k} \\ \end{align}\]

Using an extra variable \(z = \inv{B} r\):

\[r_0 = b - Ax_0, z_0 = p_0 = \inv{B} r_0\] \[\begin{align} x_{k+1} &= x_k + \alpha_k p_k & \alpha_k &= \frac{p_k^T A^T r_k}{p_k^TA^TAp_k} &= \frac{z_k^TA^Tr_k}{p_k^TA^TAp_k}\\ r_{k+1} &= r_k - \alpha_k Ap_k & & & \\ z_{k+1} &= \inv{B}r_{k+1} & & & \\ p_{k+1} &= z_{k+1} - \beta_k p_k & \beta_k &= \frac{z_{k+1}^TA^TAp_k}{p_k^TA^TAp_k} &= -\frac{ z_{k+1}^T A^T r_{k+1} }{z_k^T A^T r_k} \\ \end{align}\]

Mixed Conjugate Gradient/Residuals

Any convex combination of \(I, A^T\) as a metric will minimize the corresponding weighted sum of energy and residual norm.