The Concept of Conjugate Gradient Descent in Python

Introduction

In this example we follow An Introduction to the Conjugate Gradient Method Without the Agonizing Pain and demonstrate few concepts in Python. I shamelessly quote the original document in few places. References to equations and figures are given in terms of the original document.

A 2D example we will work with:

The quadratic form: \[f(x) = \frac{1}{2}x^\text{T}Ax-b^\text{T}x+c\]
\[f'(x) = Ax – b\] \[Ax – b = 0\]

Evaluate the function on a small grid and draw it:

Contour plot of the same quadratic form:

The Method of Steepest Descent

The algorithm is implemented following (with slight changes) the psudocode from the Appendix B1 of the tutorial.

Here we can follow its steps:

Let us have a look at one iteration

We visualize various thingies what play a role in there. It will help to gain the intuition of how it goes. First we initialize the starting point and few parameters:

The first object to look at is residual denoted by r in the code, let’s see what exactly it stands for:

Observe now that residual actually means direction of steepest descent. The vector also has a magnitude, as you can see, but since it does not directly provide any useful information (other that some notion of how far away are we from the desired solution) we will not dig into it.

We now have the direction, the next question is how far along that direction we need to go?. The procedure to answer that is called line search and its task to take the direction and find the minimal point of the parabola formed by this aline and the quadratic bowl. We need to estimate alpha \(\in[0,1]\) which will tell how far along the line we should move. The equation to get alpha in our code is:
alpha = float(delta / (r.T * (A * r))),
which corresponds to
\(\alpha = \frac{r^\text{T}_{(0)}r_0}{r^\text{T}_{(0)}Ar_0}\) in the tutorial. Let me repeat the derivation from the tutorial to see why this will give the alpha we want (please refer to Figure 6 of the tutorial to agree with the initial premise \(r^\text{T}_{(1)}r_0 = 0\) ):
\[r^\text{T}_{(1)}r_0 = 0\]
\[(b – Ax_{(1)})^\text{T}r_0 = 0\]
\[(b – A(x_{(0)} + \alpha r_{(0)}))^\text{T}r_0 = 0\]
\[(b – Ax_{(0)})^\text{T}r_{(0)} – \alpha(Ar_{(0)})^\text{T}r_0 = 0\]
\[(b – Ax_{(0)})^\text{T}r_{(0)} = \alpha(Ar_{(0)})^\text{T}r_0\]
\[r_{0}^\text{T}r_{(0)} = \alpha r_0^\text{T}(Ar_{(0)})\]
\[\alpha = \frac{r^\text{T}_{(0)}r_0}{r^\text{T}_{(0)}Ar_0}\]

Now we know how far to got along r and we update our solution x:

And now the loop repeats: we are in some location, we compute residual, get the direction:

Update delta, and the iteration counter and we are ready to perform next line search and make the next step

Gradient descent

It might we interesting for you to compare steepest descect to what we usually refer to as gradient descent where instead of performing line search we fix alpha (learning rate) to some fixed parameter:

Change the alpha parameter to see what happens when your learning rate it too low or too high.

The Method of Conjudate Directions

As you have seen in the steepest descent method the algorihtm makes steps in same direction several times. Here’s an idea: let’s pick a set of orthoginal search directions so that if we take exactly one step of right length in each of those directions we will reach the minimum.

Conjugacy

First we introduce the notion of A-orthogonal vectors, where A is the matrix with coefficients from our orignal equation \(Ax = b\).
Take round quadratic bowl

Now take tow orhogonal vectors and draw them:

Now instread of identity matrix take a matrix that has Y axis scaled up 2 times.

To counter the scale of the matrix we scale X components of the vectors 2 times and draw them:

As you can see those vectors are not othogonal:

but they are \(A\)-othogonal under the matrix Ascaled:

Vectors \(a\) and \(b\) are \(A\)-orthogonal or conjugate iff \[a^\text{T}Ab = 0\]
The intuition is to imagine that if you stretch the contour plot so that the contours are circles and two vectors become orthogonal then they are \(A\)-orthogonal.
Ok, well, why would we care about those? Imagine that you need to find solution in non-distorted circular 2D bowl:

then you could achieve your solution in two steps using usual steepest descent:

That would be cool, wouldn’t it? Unfortunately optimization problems do not always have circular bowl error surface. So if we use conjugate directions instead of orthogonal ones then it is as-if we are solving our problem on circular bowl surface. All that is needed now is a set of conjugate (\(A\)-othogonal) directions, fortunately we can find those using conjugate Gram-Schmidt process.

Gram-Schmidt Conjugation

We will need \(n\) conjugate directions. First we take \(n\) linearly independent vectors. For example residuals we used in the steepest descent will do nicely, because we selected each next residual to be orthogonal to the previous search directions (Equation 39). Right now I will only mention that Gram-Schmidt process is the method we will use and later you will see how it works and what it does.

The Method of Conjugate Gradients

Let me generate the original error surface again.

And follow the pseudocode in B2 to implement CG method.

It works! And did reach the solution in 2 steps!

If you compare it to the steepest descent you will see that not much has changed, the only conceptually new lines are

beta = float(deltanew / deltaold)
d = r + beta * d

let’s see what they do.

The initialization is very similar, we now have d instead of r to show the direction of the descent.

Here we step into the interation of the main loop and it also looks quite familiar. I will not draw the first step since it looks exactly as the first step of the steepest descent. Only starting from the second step we will have a restricting condition: direction of our second step will have to be conjugate to the direction of the first step. But we will worry about that later, right now we boldly make our first step in the direction of the steepest descend.

The first step is done and now it is time to figure out the next direction.

This would be the direction steepest descent would take us, but this time we are bound by the rule that our next step should be A-orthogonal (conjugate) to the first step. And currently

that is not the case. So here is what we are going to do:

this \(\beta\) is a new variable, but before going into the details of what it represents let’s see what it does:

and our new direction now looks like this

and it also satisfies the condition of conjugacy

So let us go into math and see how \(\beta\) steered us to the conjugate direction.
The final equation we use to obtain \(\beta_{(i)}\) is \[\beta_{(i)} = \frac{r^\text{T}_{(i)}r_{(i)}}{r^\text{T}_{(i-1)}r_{(i-1)}}\] now let us see how we get that from Gram-Schmidt process and initial idea of building conjugate directins from the residuals.

So we have linearly independent vectors \(r_{(0)}, \ldots, r_{(n-1)}\) we call residuals. And we want to get \(d_{(0)}, \ldots, d_{(n-1)}\) conjugate vectors we call directions. To construct \(d_{(i)}\) we take \(r_{(i)}\) and subtract from it non-conjugate components of all previous \(d_{(0)}, \ldots, d_{(i-1)}\) vectors (see Figure 24).

For \(d_0\) it is easy: \(d_0 = r_0\), and for all the next ones it can be computed as \[d_{(i)} = r_{(i)} + \displaystyle\sum_{k=0}^{i-1}\beta_{ik}d_{(k)}\] where \(\beta_{ik}\) would indicate how much of each previous direction \(d_{(k)}\) should be subtracted (\(\beta\) will be a negative coefficient) from \(r_{(i)}\) For the particular direction \(d_{(i-1)}\) we compute as follows. First we multiply the whole equation by \(Ad_{(i-1)}\): \[d^\text{T}_{(i)}Ad_{(i-1)} = r^\text{T}_{(i)}Ad_{(i-1)} + \displaystyle\sum_{k=0}^{i-1}\beta_{ik}d^\text{T}_{(k)}Ad_{(i-1)}\] Next we observe that \(\forall p \in \{0, \ldots, k-2\}\) we know that \(d^\text{T}_{(p)}Ad_{(i-1)} = 0\) because of \(A\)-orthogonality. So all addends, except for the \(d^\text{T}_{(i-1)}Ad_{(i-1)}\) are 0. And the left hand side of the equation we set to 0, because this is what we want to achieve: \[0 = r^\text{T}_{(i)}Ad_{(i-1)} + \beta_{(i)(i-1)}d^\text{T}_{(i-1)}Ad_{(i-1)}\] and from that \[\beta_{(i)(i-1)} = – \frac{r^\text{T}_{(i)}Ad_{(i-1)}}{d^\text{T}_{(i-1)}Ad_{(i-1)}}\] That is already computable and in our code would look like

beta = -float((r.T * A * d) / float(d.T * A * d))

but we can simplify it further and avoid doing matrix multiplications.
Note that (page 6) \(r_{(i)} = -Ae_{(i)}\) and from that we can get a recurrence for the residual \[r_{(i+1)} = -Ae_{(i+1)} = -A(e_{(i)} + \alpha_{(i)}d_{(i)}) = r_{(i)} – \alpha_{(i)}Ad_{(i)}\] we multiply it by \(r_{(i)}\) \[r^\text{T}_{(i)}r_{(i)} = r^\text{T}_{(i)}r_{(i-1)} – \alpha_{(i-1)}r^\text{T}_{(i)}Ad_{(i-1)}\]
\[\alpha_{(i-1)}r^\text{T}_{(i)}Ad_{(i-1)} = r^\text{T}_{(i)}r_{(i-1)} – r^\text{T}_{(i)}r_{(i)}\] Since the residuals are othogonal we know \(r^\text{T}_{(i)}r_{(i-1)} = 0\) as we can leave it out: \[r^\text{T}_{(i)}Ad_{(i-1)} = – \frac{1}{\alpha_{(i-1)}}r^\text{T}_{(i)}r_{(i)}\] Now we take an equation for \(\beta_{(i)(i-1)}\) and put it in: \[\beta_{(i)(i-1)} = \frac{1}{\alpha_{(i-1)}}\frac{r^\text{T}_{(i)}r_{(i)}}{d^\text{T}_{(i-1)}Ad_{(i-1)}}\] We also had an equation for \(\alpha\): \[\alpha_{(i)} = \frac{d^\text{T}_{(i)}r_{(i)}}{d^\text{T}_{(i)}Ad_{(i)}}\] We apply it on \(i-1\), replace \(\alpha_{(i-1)}\) and get \[\frac{d^\text{T}_{(i-1)}Ad_{(i-1)}}{d^\text{T}_{(i-1)}r_{(i-1)}}\frac{r^\text{T}_{(i)}r_{(i)}}{d^\text{T}_{(i-1)}Ad_{(i-1)}} = \frac{r^\text{T}_{(i)}r_{(i)}}{d^\text{T}_{(i-1)}r_{(i-1)}}\] The final touch is to use \(d^\text{T}_{(i)}r_{(i)} = r^\text{T}_{(i)}r_{(i)}\) (see equatins 40 and 42) we get the desired \[\beta_{(i)} = \frac{r^\text{T}_{(i)}r_{(i)}}{r^\text{T}_{(i-1)}r_{(i-1)}}\]

To remind you what it was: \(\beta_{(i)}\) shows how much of the direction \(d_{(i-1)}\) we should take to keep \(d_{(i)}\) \(A\)-orthogonal (conjugate) to \(d_{(i-1)}\).

Now we increase the iteration counter and continue to the next (in that partiuclar 2D case, final) iteration.

where we will decide on how far we should go in this direction and happilty reach the minimum point:

All that is fine you say, but looking for a minumum point of a quadratic form isn’t what you usually need in practice. Next piece of this story will cast light on that.