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

A technical blog post by an author writing a paper on the Chebyshev-Jacobi method argues that, contrary to common belief, Gauss-Seidel is not always faster than Jacobi, and presents a modern convex optimization framework showing Jacobi as gradient descent and symmetric Gauss-Seidel as coordinate descent. The post derives the Chebyshev-Jacobi method as the optimal descent method with the best worst-case convergence rate, requiring only the inverse of a diagonal matrix for parallelization.

Glad you asked I am actually in the process of writing a short paper where I re-derive the Chebyshev-Jacobi method from a modern convex optimization point of view. I think adopting this point of view is nice as it enables making link with popular techniques being used elsewhere looking at you the ML community . So here is the very condensed version. I’ll skip all of the proofs for clarity. Suppose you want to solve Pz = q where P is a symmetric positive definite matrix. It turns out that this system of equation is nothing but the optimality condition of the following strictly convex unconstrained quadratic program Why is it the optimality condition ? Because the minimizer is a point z \star where f z \star is minimal and thus its gradient \nabla f z \star \equiv P z \star - q = 0. Hence P z \star = q, i.e. the problem we started with. How to solve this problem? – Most optimization solvers use some sort of descent method incrementally improving an initial guess until some convergence criterion is reached. A large class of these descent methods can be written as Here, z t is our estimate of the solution at iteration t. The term \sum {i=0}^{t-1} \alpha i^{ t } z {i+1} - z i essentially is a momentum term, and the last term \alpha t^{ t } \nabla f z t is just the descent direction given by the gradient. What differentiate most optimizers is the choice of the non-stationary weights \left\{ \alpha i^{ t } \right\} {i=0, \cdots, t}. For instance, standard gradient descent can be recovered if you set \alpha i^{ t } = 0 \forall i = 0, \cdots, t-1 and \alpha t^{ t } = -\eta where \eta is the so-called step size . How does the Jacobi method fit in this framework ? – Suppose you want to solve Ax = b with A a symmetric positive definite matrix. The Jacobi update rule reads It vaguely looks like the generic update rule for z but not quite since D^{-\frac12}A is not symmetric unless D is a constant diagonal matrix d I . However, if you multiply the Jacobi update rule by D^{\frac12} from the left and let z t = D^{\frac12} x t, it can be rewritten as So, up to an invertible diagonal transform, the Jacobi method essentially find the minimizer of the quadratic form using a gradient descent with a constant unit step size. Likewise if I recall correctly , you can show that the symmetric Gauss-Seidel method essentially is a form coordinate-wise descent. So what is the Chebyshev-Jacobi method ? – Basically, to derive the Chebyshev-Jacobi method, you ask the following question : among all descent methods of the form how does I choose the set of weights \left\{ \alpha i^{ t } \right\} such that the method has the best worst-case convergence rate? The proof is not very complicated but a bit tedious so I’ll skip it. Essentially, it has to do with matrix polynomials, and in particular the Chebyshev polynomials, hence the name Chebyshev-Jacobi. At the end of the day, you can show that the descent method having the best worst-case convergence rate involves only z t, the gradient \nabla f z t = Pz t - q and the first difference z t - z {t-1}. For our original problem in terms of x , this leads to the non-stationary update rule with \omega 0 = 2 I think , x 0 is whatever your initial guess is and with L and \mu being the largest and smallest eigenvalues of D^{-\frac12} A D^{-\frac12}. For our test problem, L and \mu can actually be computed analytically. In general, you’d need to have some estimates \hat{L} \geq L 0 and 0 < \hat{\mu} \leq \mu. Obviously, the better your estimates, the closer you are to the optimal method. Except for the extra scalar update for \omega t and buffer to store x {t-1}, you end up with a method that keeps the simplicity of the original Jacobi method, in particular its straightforward parallelization as it requires the inverse of a diagonal matrix a purely local operation while having a significantly better convergence rate. As @rsci https://fortran-lang.discourse.group/u/rsci mentioned, it also has the same worst-case scaling as Conjugate Gradient also, which is nice. So why use Chebyshev-Jacobi instead of CG? – Mainly a matter of personal preferences. In practice though, CG may perform better in terms of iteration count . One key advantage of Chebyshev-Jacobi however is that it requires no inner product while CG does at least 2 or 3 if I recall correctly . For serial computations, it probably doesn’t change things much, but for large-scale parallel ones it does as inner product require a synchronization of all the processes. In contrast, for the Chebyshev-Jacobi method you can perform inner products only once in a while to compute the norm of the first difference x {t+1} - x {t} if you use the classical 2-norm, or even avoid them entirely if using for instance the 1-norm or \infty-norm although it does require some form of all reduce nonetheless . So basically, I think Chebyshev-Jacobi may have better scaling capabilities in terms of parallelization. It is not a silver bullet though, most notably because you really do need good estimates of the extremal eigenvalues which CG does not require . And that’s the only bit of information you use about the spectrum. In contrast, CG is a Krylov method so it essentially adapts on-the-fly to the actual distribution of eigenvalues and tends to have somewhat better convergence properties even though the worst case is the same . The picture also changes quite drastically if you consider preconditioning. I don’t think Chebyshev-Jacobi can outperform preconditioned CG if you use a suitable preconditioner.