1.1. Setup¶

First, we import the necessary libraries. We’ll be using numpy for numerical operations and matplotlib for plotting.

1.2. Generate an Example Problem¶

As mentioned before, the CG method requires the matrix $A$ to be symmetric and positive-definite. A simple way to generate such a matrix is to create a random matrix $M$ and compute $A = M^T M$. This construction guarantees that $A$ is symmetric and positive semi-definite. If $M$ is non-singular, $A$ will be positive-definite.

We will set up a linear system $Ax = b$ where we know the true solution x_true beforehand. This allows us to easily verify the correctness of our implementations.

1.3. The “Naive” Approach: Direct Solver¶

Before diving into iterative methods, let’s solve the system using a standard direct solver from NumPy, np.linalg.solve. This function typically uses an LU decomposition, which is very accurate and robust for small to moderately sized dense matrices.

This gives us a baseline for both speed and accuracy. For a 1000x1000 matrix, the direct solver should be quite fast, but as we scale up, its computational cost ($O(n^3)$) and memory requirements ($O(n^2)$) will become significant.

1.4. The Idea behind the Conjugate Gradient Method¶

The CG method combines the concept of steepest descent with Krylow subspaces:

First, we start at some initial point $x_0$, and take a step in the direction of steepest descent.

• This step takes place in the first, one-dimensional Krylow-subspace

• Ideally, this step takes us the optimal amount of length along that direction

Then, for each further step, we find the new direction of maximum descent, but we only consider directions that do not undo previous progress

• This means that the steps are $A$-conjugate, instead of merely being orthogonal

• The dimension of the Krylow-subspace increases with each step, so we require at most $n$ steps

This approach gives us the exact solution after at most $n$ steps, if we have perfect arithmetic. While we cannot assume this in practice, much fewer steps can give a good approximation.

1.5. The Conjugate Gradient Algorithm (Pseudo-code)¶

Now let’s look at the algorithm we want to implement. The goal is to find $x$ that minimizes $\frac{1}{2}x^T A x - x^T b$.

Given the system $Ax=b$, an initial guess $x_0$, a tolerance tol, and a maximum number of iterations max_iter:

1. Initialize

   • $r_0 = b - Ax_0$

   • $p_0 = r_0$

   • $k = 0$

2. Iterate while $k < \text{maxiter}$ and $\|r_k\| > \text{tol}$:

   • $\alpha_k = \frac{r_k^T r_k}{p_k^T A p_k}$

   • $x_{k+1} = x_k + \alpha_k p_k$

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

   • $\beta_k = \frac{r_{k+1}^T r_{k+1}}{r_k^T r_k}$

   • $p_{k+1} = r_{k+1} + \beta_k p_k$

   • $k = k + 1$

3. Return $x_k$

Key variables:

• x: The solution vector, which we are improving in each step.

• r: The residual vector ($r = b - Ax$), which we can use to compute the current error.

• p: The search direction vector. Note that we always move along the residual, here we choose special “$A$-orthogonal” directions.

• alpha: The step size.

• beta: The improvement factor for the next search direction.

2. Exercise: CG Implementation¶

2.1. Implementing the Algorithm¶

Now it’s your turn! Fill in the missing lines in the conjugate_gradient function below. Follow the pseudo-code from the previous section.

Hint: The term $r_k^T r_k$ appears multiple times. It’s efficient to compute it once per iteration and reuse the value.

2.2. Test your Implementation¶

Let’s run your implementation, check that it works, and see how it performs compared to the direct solver.

Don’t worry if the code does not seem to converge, it is quite possible that you will have errors around e-2, even if you implemented everything correctly.

What is happening?

• Is the result correct? (Compare error_cg to error_naive from earlier).

• How does the computation time compare?

• How many iterations did it take to converge? Did it converge at all?

Let’s plot the convergence history to see how the residual norm decreases over time.

You will probably notice some things with your implementation:

1. The direct solver is usually way faster:

• For one, CG’s strength lies in very large and/or sparse matrices.

• Even more importantly, the CG method often converges in much fewer than the $n$ steps in practice

• As long as our implementation does not converge at most of the time, it is nowhere near as efficient as it could be!

2. For any non-trivial example ($n > 10$) and reasonable error tolerances, the algorithm may not converge: You can see in the plot, that the method has a hard time getting below an error tolerance of 1.e-2.

3. The algorithm is not guaranteed to converge after $n$ steps (you can check this by setting max-iter to a value higher than n): Again, this is due to the lack of perfectly accurate arithmetic.

As you can see, the method is a lot worse than we would have hoped. In the second part of this exercise, we will take a look at preconditioners to speed up convergence and improve these results.

Bonus exercise¶

If you are finished before the break, have a think about how we could verify our assumption, before implementing preconditioning. Are there some easy ways to modify the construction of $A$ (i.e. our generate_spd_matrix method) to make sure it is well-conditioned?

Answer: By increasing the factor in M += np.eye(n) to e.g. M += 100 * np.eye(n), we can make the matrix diagonally dominant. This gives us a well-conditioned matrix, that our simple algorithm should be able to handle well. However, this is only helpful for demonstration purposes, if we really want to solve an ill-conditioned system, we need to go about it in a different way.

3. Preconditioning¶

3.1. The Condition Number¶

The number of iterations CG takes to converge is closely related to the condition number of the matrix $A$, denoted $\kappa(A)$.

$\kappa(A) = \|A\| \|A^{-1}\| = \frac{\lambda_{\max}}{\lambda_{\min}}$ (for symmetric positive-definite matrices)

A high condition number means the matrix is “ill-conditioned”. This has two major negative effects:

1. Error Amplification: Small errors in the input b can lead to large errors in the output x.

2. Slow Convergence: Iterative methods like CG converge very slowly. The number of iterations is roughly proportional to $\sqrt{\kappa(A)}$.

Let’s check the condition number of our matrix A.

That’s a pretty large number! This explains why our simple CG implementation might have taken many iterations to converge. To fix this, we can use preconditioning.

3.2. Preconditioning Theory¶

The idea is to find a preconditioner matrix that is:

1. Also symmetric and positive-definite.

2. A good approximation of $A$.

3. Easy to invert.

Instead of solving $Ax=b$, we solve a modified system. The most common form is left preconditioning, where we solve $P^{-1}Ax = P^{-1}b$.

Note: Even if $A$ and $P$ are symmetric, the new system matrix $M = P^{-1}A$ is generally not symmetric! This means we cannot apply the standard CG algorithm directly to it. The Preconditioned Conjugate Gradient method is a clever variant that is mathematically equivalent to running CG on a transformed, symmetric system, without having to explicitly form $M$.

Our goal is to find a $P$ such that $\kappa(P^{-1}A) \ll \kappa(A)$.

The simplest (but often effective) preconditioner is the Jacobi Preconditioner (or Diagonal Preconditioner). We choose $P$ to be the diagonal of $A$.

$P = \text{diag}(A_{11}, A_{22}, ..., A_{nn})$

Since $A$ is symmetric positive-definite, all its diagonal entries $A_{ii}$ are positive, which means $P$ is also positive-definite. Inverting this $P$ is also trivial, as $P^{-1}$ is just a diagonal matrix with entries $1/A_{ii}$.

3.3. Implementing the Jacobi Preconditioner¶

For this part, we will first modify our matrix to generate a new, ill-conditioned matrix that work well with Jacobi preconditioning.

As for the obvious question: Why bother? The Jacobi preconditioner works well for ill-conditioned, but diagonally dominant matrices, where the main source of issues comes from the diagonal. So this modification is only for demonstrative purposes, but we will get back to the original matrix later.

Let’s now create our preconditioner. We will actually compute its inverse, since that’s what the algorithm needs.

As you can see, the condition number has been significantly reduced! This should lead to much faster convergence.

3.4. The Preconditioned Conjugate Gradients Algorithm¶

The algorithm is slightly modified to include the preconditioner. The key step is to solve the system $Pz_k = r_k$ for $z_k$ (which is easy since we have $P^{-1}$).

1. Initialize

   • $r_0 = b - Ax_0$

   • Solve $Pz_0 = r_0$ (i.e., $z_0 = P^{-1}r_0$)

   • $p_0 = z_0$

   • $k = 0$

2. Iterate:

   • $\alpha_k = \frac{r_k^T z_k}{p_k^T A p_k}$

   • $x_{k+1} = x_k + \alpha_k p_k$

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

   • Check for convergence on $|r_{k+1}|$

   • Solve $Pz_{k+1} = r_{k+1}$

   • $\beta_k = \frac{r_{k+1}^T z_{k+1}}{r_k^T z_k}$

   • $p_{k+1} = z_{k+1} + \beta_k p_k$

Now, implement the PCG algorithm below.

3.5. Testing the PCG Algorithm¶

Let’s run our new PCG implementation and see the difference.

As we can see, the improvement is dramatic and the new algorithm converges after only a handful of steps. Of course, this particular example is designed to show drastic effects, and it will typically not go quite as well in practice, but you get the idea.

Open-ended exercise: Explore the limitations of the methods we implemented¶

• Just like the first part, take some time to try different values or matrices, and get a feel for the limitations of the algorithm

• What happens, if you go back to the first type of matrices that we used, and apply the Jacobi preconditioner?

  • Does the result improve?

  • What may be the reason for its behaviour?

  • What might a suitable preconditioner for this type of matrix look like?

• What happens if you try to apply the CG method to a fully random matrix?

• You could also try to come up with other kinds of matrix and see where the Jacobi preconditioner starts to fail

Outlook¶

As you’ve seen in this exercise, the CG method can be a very elegant and efficient approach to solving many issues, if they have the right properties. But just as there were many cases, where it brought about a significant performance benefit, there are also many situations, where it is either unsuitable, or the preconditioner needs to be carefully chosen.

Always keep in mind: Just because you have a hammer, not every problem is a nail.

If you are interested in further reading, you can move on to more advanced topics, such as:

• GMRES: The Generalized Minimal Residual applies the ideas of CG in a more general setting. It is on of the standard approaches if you are faced with a non-symmetric matrix $A$. It comes with some drawbacks such as increased memory requirements, though.

• Other preconditioners: There are a multitude of preconditioners, and the best choice always depends on the problem you work with. Typical ones include

  • Gauss-Seidel: It uses the lower triangular part of $A$, and is relatively difficult to parallelize

  • (Symmetric) Succesive Over-Relaxation (SORR): This method is based on Gauss-Seidel, but slightly more involved.