Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

I followed the derivation for the Conjugate Gradient method from the documents shared below

http://en.wikipedia.org/wiki/Conjugate_gradient_method

http://www.idi.ntnu.no/~elster/tdt24/tdt24-f09/cg.pdf

I understand almost all of this, except how $\beta_k$ is derived.

From the Wiki derivation I understand we can get

$$ r_{k+1} = r_k - \alpha_k A d_k$$

reworking this

$$Ad_k = -\frac{1}{\alpha_k} (r_{k+1} - r_k)$$

And $\alpha_k$ is

$$ \alpha_k = \frac{r_k^Tr_k}{d_k^TAd_k} $$

Now, somehow all derivations I see $\beta_k$ declared to be

$$ \beta_k = \frac{r_{k+1}^Tr_{k+1}}{r_k^Tr_k} $$

and imply that it follows

$$ d_{k+1} = r_{k+1} + \beta_k d_k$$

I was not able to verify this however, my algebra must have gone wrong somewhere. Also, it seems $d_k$ and $r_k$ are assumed to be the same thing, I explain that as first direction $p_0$ is the same as the residual $r_0$ and from there we always calculate the next orthogonal direction, so $d_k$ and $r_k$ remain equal.

Any help would be appreciated,

share|improve this question
    
$p_k$ and $r_k$ are not the same for $k>0$, because you don't just use the next orthogonal direction (which would be the gradient $r_k$). You take the gradient direction $r_k$ and orthogonalize this to all the previous $p_k$ (which span the current solution space), but with respect to $A$, making them conjugate with respect to $A$ and not just orthogonal (and this way also different from $r_k$). –  Christian Rau Feb 16 '12 at 15:31
    
Ok $p_k \ne r_k$. I also should have used "A-orthogonal" instead of orthogonal. Besides those, any insights into the $\beta_k$ part of the question? –  BB_ML Feb 16 '12 at 15:39
    
This link here pg 16 www10.informatik.uni-erlangen.de/Teaching/Courses/WS2010/SiWiR/… says that "The search vectors and residuals then form one and the same subspace". Does that mean the basis of both those subspaces can be different? –  BB_ML Feb 16 '12 at 17:35
    
Of course it can. There are many bases for this sub-space, but only the $p$-basis is $A$-orthogonal, the $r$-basis isn't. –  Christian Rau Feb 16 '12 at 17:44
1  
Sort of a side note, but here is the best explanation of CG that I have found: cs.cmu.edu/~quake-papers/painless-conjugate-gradient.pdf –  Nick Alger Feb 18 '12 at 10:07
show 3 more comments

2 Answers

up vote 4 down vote accepted

There seems to be a tendency for expositions of the conjugate gradient method to get bogged down in the details of arithmetic manipulations of coefficients that hide the bigger picture. A lot of texts that explain the method in detail with formulas have already been cited in your question, your own answer and the comments; I'll try to describe on a more intuitive level how I think the method can best be understood.

The Wikipedia article, your question and the texts you link to use slightly different index offsets. In the Wikipedia article and in your question the coefficient linking $d_i$ and $d_{i+1}$ is $\beta_i$, whereas in the paper and the notes it's $\beta_{i+1}$. I'll go with the notes since I'll be referring to them every once in a while.

First, as in the notes, it's helpful to distinguish between the specific choice of the conjugate search directions and the general method of searching along conjugate directions.

So let's first assume, as in the notes, that we're given some $A$-orthogonal basis $\{d_i\}$ of search directions, that is, $d_i^\top A d_j=\delta_{ij}$, where $A$ is the positive-definite symmetric matrix of a linear system $Ax=b$ to be solved or of a quadratic function $\frac12x^\top Ax-b^\top x$ to be minimized, or the Hessian of a non-linear optimization problem.

I find it helpful to keep in mind, in addition to the space in which $x$ lives, a space related by $y=A^{1/2}x$, in which the quadratic function becomes $\frac12y^\top y-b^\top A^{-1/2}x$, the level sets become hyperspheres and $A$-orthogonality becomes orthogonality.

A key property to understand is that the first $n$ residuals and the first $n$ search directions span the same space, so it means the same thing that a vector is orthogonal to the first $n$ residuals or to the first $n$ search directions. This is because in each step the search direction is formed as the residual plus a vector from the previously spanned space.

In the $x$ space, the residuals are mutually orthogonal, and in the $y$ space, the search directions are mutually orthogonal. Because the first $n$ residuals and the first $n$ search directions span the same space, this also implies that in the $x$ space the residuals are orthogonal to all previous search directions and in the $y$ space the search directions are orthogonal to all previous residuals.

Each search ends at a minimum along the search direction, or equivalently at a point where the gradient is orthogonal to the search direction. Since the residual is the (negative) gradient, that immediately implies $r_{i+1}\perp d_i$, and it's rather misleading that the notes rederive this with complicated arithmetic as (23) instead of noting that $\alpha_k$ was chosen so as to make this true. It then follows by induction that $r_{i+1}\perp d_j$ for all $j\le i$, since $r_{i+1}=r_i+A(x_{i+1}-x_i)$, $r_i$ is orthogonal to all $r_j$ by the induction hypothesis, $x_{i+1}-x_i$ is a vector along the search direction $d_i$ and $Ad_i$ is orthogonal to all $d_j$ (by the $A$-orthogonality of the search directions), and thus to all $r_j$.

All this holds for any given $A$-orthogonal basis $\{d_i\}$ of search directions. Now comes the part where we choose this basis as we go along, based on the residuals. So far it's been a one-way street from the given search directions to the residuals, using the coefficients $\alpha_k$ to choose the local minimum in each step and make the residuals orthogonal; now we also build the search directions from the residuals, using the coefficients $\beta_k$ to make the search directions $A$-orthogonal.

This is the part that often seems a bit magical in the expositions, but it's rather straightforward if you keep in mind the above. That a single $A$-orthogonalization suffices is a not-quite-so-magical consequence of the orthogonality of the residuals. Just like the mutual orthogonality of the residuals was ensured by induction for earlier residuals and by the choice of $\alpha_k$ for the current pair of residuals, so the mutual $A$-orthogonality of the search directions is ensured by induction for earlier search directions and by the choice of $\beta_k$ for the current pair of search directions.

We want to use the residual $r_i$ to construct a new search direction $d_i$ that's $A$-orthogonal to all previous $d_j$. As described in the notes, we could do a full Gram-Schmidt-style $A$-orthogonalization to all previous $d_j$, but that's not necessary because $r_i$ is already $A$-orthogonal to all previous $d_j$ except $d_{i-1}$ (to which it's orthogonal). This is because the difference between successive residuals $r_{j+1}$ and $r_j$ is proportional to $Ad_j$, and $r_i$ is orthogonal to $r_{j+1}$ and $r_j$ (and thus $A$-orthogonal to $d_j$) unless $i=j+1$. So there's only a single search direction $d_{i-1}$ to $A$-orthogonalize to, and this determines the parameter $\beta_i$. We want $d_i=r_i+\beta_id_{i-1}$ to be $A$-orthogonal to $d_{i-1}$. Since $Ad_{i-1}$ is proportional to $r_i-r_{i-1}$, that means

$$(r_i-r_{i-1})^\top(r_i+\beta_id_{i-1})\stackrel{!}{=}0\;.$$

Then using the orthogonality of $r_i$ to $r_{i-1}$ and $d_{i-1}$ yields

$$\beta_i=\frac{r_i^\top r_i}{r_{i-1}^\top d_{i-1}}\;,$$

and using $d_{i-1}=r_{i-1}+\beta_{i-1}d_{i-2}$ and $r_{i-1}\perp d_{i-2}$ this can be rewritten as

$$\beta_i=\frac{r_i^\top r_i}{r_{i-1}^\top r_{i-1}}\;.$$

I hope I've contributed a little to making things a bit clearer with a minimal amount of arithmetic; I've probably repeated a lot of things that you'd already understood, but I thought it would be valuable to have this here for the future even if you'd already answered your original question in the meantime.

share|improve this answer
    
Very cool; I changed the accepted answer from mine to this one. Cheers. –  BB_ML Feb 22 '12 at 16:24
add comment

This link

http://www.bruguier.com/notes/ConjugateGradient.pdf

pg 4-6 has the best info on this. The key is using the recurrence relationship on the residuals for the numerator (to make $A$ disappear), and using the definition of $\alpha_i$ for the denominator. It's also best to think of the whole $\beta$ business as a piece of a modified Gram-Schmidt.

In this link

http://caig.cs.nctu.edu.tw/course/NM07S/slides/chap7_2.pdf

pg 23-24, it is proved that $d_k = r_k$, that is residual equals search direction, so that's how they were able to use use a pure $r$ based equation in the code.

share|improve this answer
    
When I open that second PDF file I don't get pages up to 23-24, and I can't find anything about $d_k=r_k$. Did you post the right link? –  joriki Feb 20 '12 at 11:25
    
Wrong link; corrected. –  BB_ML Feb 22 '12 at 21:58
    
There's no proof or mention of $d_k=r_k$ in that document, and that equation doesn't hold. (If it did, we'd have $\beta_k=0$, and there would be no need to introduce $\beta_k$ and to distinguish between $d_k$ and $r_k$.) Perhaps you're referring to the fact that $d_{k-1}$ can be replaced by $r_{k-1}$ in the calculation of $\beta_k$? See my answer for a derivation of that fact; this doesn't mean that $d_k=r_k$. –  joriki Feb 23 '12 at 15:38
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.