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
\mathrm{minimize} \quad f(z) \equiv \dfrac12 z^\top P z - z^\top q
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
z_{t+1} = z_t + \sum_{i=0}^{t-1} \alpha_i^{(t)} \left( z_{i+1} - z_i \right) + \alpha_t^{(t)} \nabla f(z_t).
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. Gradient descent with optimal step size can be obtained by letting \alpha_t^{(t)} = \mathrm{argmin} \ f \left( x_t + \alpha_t^{(t)} \nabla f (x_t) \right), and Polyak’s heavy-ball is yet another example.
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
x_{t+1} = x_t - D^{-1} \left( Ax_t - b \right).
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
z_{t+1} = z_t - \left( D^{-\frac12} AD^{-\frac12} z_t - D^{-\frac12} b \right).
So, up to an invertible diagonal transform, the Jacobi method essentially finds the minimizer of the quadratic form
f(z) \equiv \dfrac12 z^\top D^{-\frac12} A D^{-\frac12} z - z^\top D^{-\frac12}b
using a gradient descent with a constant unit step size ^\dagger. Likewise (if I recall correctly), you can show that the symmetric Gauss-Seidel method basically is a form coordinate-wise descent.
So what is the Chebyshev-Jacobi method ? – To derive the Chebyshev-Jacobi method, you ask the following question : among all descent methods of the form
z_{t+1} = z_t + \sum_{i=0}^{t-1} \alpha_i^{(t)} \left( z_{i+1} - z_i \right) + \alpha_t^{(t)} \nabla f(z_t).
how do 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. It mostly 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
\left\{
\begin{aligned}
\omega_t & = \left( 1 - \dfrac{1}{4a^2} \omega_{t-1} \right)^{-1}, \\
x_{t+1} & = \omega_t \left( I - \dfrac{2}{L + \mu} D^{-1} A \right) x_t + \dfrac{2}{L+\mu} \omega_t D^{-1}b + \left( \omega_{t} - 1 \right) x_{t-1},
\end{aligned}
\right.
with \omega_0 = 2 (I think), x_0 is whatever your initial guess is and
a = \dfrac{L - \mu}{L + \mu}
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 mentioned, it also has the same worst-case convergence rate 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.
^\dagger : Since it uses a constant step size no matter what the Lipchitz constant of the quadratic form is, I believe this is another way to prove that the Jacobi method converges only for a subset of symmetric positive definite matrices. I may explore that idea some time.
I think I do not have typos, but ain’t entirely sure. I’ll double check with the latest draft whenever I have time.
I compared Jacoby, GS, SOR and preconditioned CG for two small problems. Preconditioners were a tridiagonal solver (which mimics an ILU), Jacobi, and a polynomial preconditioner that I never could get to converge. The two matrices were a periodic tri-diagonal system, and solution of a 2D Laplace equation Finite Difference system. SOR if you can find something close to an optimal acceleration parameter is hard to beat in terms of convergence but Jacobi-CG is usually a little better.