# Gauss-Seidel is faster than Jacobi. Or is it? > Source: > Published: 2026-09-09 14:02:35+00:00 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.