Gauss-Seidel is faster than Jacobi. Or is it?

Uncategorized
1

Hej,

It’s been a long time! Hope everyone enjoyed their summer break. It has been a long time also since I had time to write a new blog post so here it comes: Gauss-Seidel is faster than Jacobi. Or is it?

It is a continuation of my previous post (almost a year ago already) on the Jacobi method in Fortran. I wanted to keep things at a high level but still explain why, on modern CPUs, the Jacobi method could actually be faster than Gauss-Seidel despite the maths saying its convergence rate is worse. I guess most of the material will be pretty obvious for our members, but who knows, it may be useful for students :slight_smile:

9 Likes
2

I would find it interesting to compare your Jacobi implementation used as a pre-conditioner for GMRES or CG with ILU with parallelization in mind. Jacobi is still attractive as a pre-conditioner in some applications due to the fact that it can be parallelized when other (faster converging) methods can’t.

3

That is part of the series i have in mind :slight_smile:

Eventually, I want to address red/black Gauss-Seidel, Chebyshev-Jacobi, Conjugate gradient, Krylov methods, and preconditioning. While I’d love to be able to stick to 1 post/month, I ain’t sure I’ll manage. Too many other things to do and too little time.

4

Great work! One thing I was wondering about is the roughly 10× difference in cost per iteration between the two methods.

If vectorization is the main reason for the difference, then for double precision and AVX-512 I would have expected something closer to an 8× improvement in the ideal case. Could the remaining difference perhaps be related to the loop-carried dependency in lexicographic Gauss–Seidel? Besides preventing straightforward vectorization, this dependency might also limit instruction-level parallelism and make the loop more sensitive to instruction latency, whereas Jacobi has many independent updates that the CPU can overlap.

It might be interesting to repeat the test with vectorization explicitly disabled for Jacobi (novec or the corresponding compiler option). Comparing optimized Jacobi with and without vectorization could help separate the SIMD contribution from other effects. If Jacobi remains noticeably faster than Gauss–Seidel in the non-vectorized case, that might give some indication of how important these other effects are.

Also, I think there is a small typo in the Gauss–Seidel formula in the post: in the last term on the RHS the superscript should be (t+1), rather than (t).

5

Typo fixed. Thanks a lot.

Yes, there are many things that could be explored using the Jacobi and Gauss-Seidel kernels, both mathematically and code/compiler-wise. I came to realize that, even though they can be massively outperformed by other techniques and may seem outdated, these two kernels are simple enough that any undergraduate student can understand the gist of it. Yet, despite this mathematical (relative) simplicity, they offer a superb platform to explore good practice in scientific computing (and programming in general I guess) as well as to understand subtle concepts such as vectorization and how a compiler may leverage fundamental properties to make the code run faster.

I do have in mind a set of more computer science/compiler oriented posts in the future, but I ain’t quite there yet.

6

Very interesting post!

Since you use Poisson equation in the post, could I ask if these basic iterative methods (like Jacobi, Gauss-Seidel, SOR) work well on this problem? Or Krylov-based methods are way better than these basic iterative methods for Poisson equation?

7

Very roughly, for the standard 2D Poisson problem on an N \times N grid, Jacobi needs O(N^2) iterations for a fixed error reduction, while Gauss-Seidel has the same scaling but usually takes about half as many iterations. Optimal parameter SOR and Chebyshev iteration are O(N), and CG has the same worst-case scaling.

In practice, CG is often noticeably better than that bound, since its convergence depends on the actual eigenvalue distribution and on which spectral components are present in the error.

It is also worth mentioning multigrid methods, whose iteration count can be essentially independent of the problem size.

1 Like
8

Let’s be specific here: the Poisson equation I consider as test case is the 2D Poisson equation on the unit square (or rectangle) discretized with second-order accurate central finite differences on a uniform mesh.

Do basic iterative methods work well? (Jacobi, Gauss-Seidel, SOR) – They work, for sure. Do they work well ? It kind of depends. I would never use Jacobi or Gauss-Seidel as actual solvers outside of a pedagogical objective. While they are very simple to implement, they are simply too slow. On the blog post, you can see that for a 512² grid, the fastest Jacobi takes 16 seconds to solve the system. Once. This is not practical in any way if you have to solve this equation repeatedly like in an incompressible Navier-Stokes simulation.

SOR, that is a different story. For this particular discretized equation, I think it has one of the best trade-off between code complexity and computational performance. Code-wise, it is basically a minor variation of Gauss-Seidel, so very easy to implement. Performance-wise, it is way more efficient. Even more so for this particular test case since you actually have a closed-form expression for the optimal relaxation weight. There are some caveats though:

  • Just like Gauss-Seidel, SOR suffers from the same loop-carried dependency. It’ll thus be harder for the compiler to optimize the code. But at the same time, the method requires far fewer iteration to converge. So, trade off again.
  • If you run a 3D case (or even a big 2D case) and have to implement the linear solver yourself, SOR can be a pain to parallelize.

If you want to stay in Jacobi-land in terms of solver complexity, but still have pretty good performances and relatively straightforward parallelization, Chebyshev-Jacobi is in my opinion your best alternative. You get convergence rate somewhat similar to SOR but retain the simplicity of Jacobi. And it might actually be more computationally efficient than SOR (even though it requires more iterations) precisely because the compiler can vectorize the code easily. Another important feature is that, except for periodic computation of the norm of the correction vector, no inner product is required and thus no all_reduce in parallel which would force a synchronization of the different processes.

How about Krylov methods ? (Conjugate Gradient, GMRES) – I’ll consider only CG since the test case is a symmetric positive definite operator. You’ll have to be slightly more careful when implementing it, but it ain’t that complicated. And it’ll converge a lot faster than Jacobi or Gauss-Seidel. It also turns out to be very closely related to the Chebyshev-Jacobi method (the two can essentially be derived in the same way).

I’ve just ran a quick test case with LightKrylov, and basically it converges in roughly 1000 iterations. That is 1000 matrix-vector products, compared to the 74 000 / 138 000 needed for Gauss-Seidel/Jacobi. So even if a single iteration may be slightly more costly (in particular because it involves a handful of inner products), time-to-solution is well under a second. With a block-Jacobi preconditionner, the number of iterations drops around 500, so even faster. And in general, having a good preconditionner is key for maximum computational performance but is very case-dependent.