Conjugate gradient performance on sparse SPD matrices

What this research found

Conjugate gradient is the standard iterative solver for large sparse symmetric positive-definite linear systems, and its speed depends almost entirely on preconditioning. K-Dense wrote its own conjugate gradient iteration and benchmarked three variants — unpreconditioned, Jacobi, and incomplete Cholesky — on two matrices from the SuiteSparse collection chosen to span very different size and conditioning regimes. Incomplete Cholesky cut iteration counts by roughly 30 to 45 times, and on the badly conditioned matrix it also turned a wrong answer into a correct one: the unpreconditioned solve met its residual tolerance while returning a solution with an error of order one.

• On the ill-conditioned structural-stiffness matrix bcsstk16 (4,884 unknowns, estimated condition number 4.94 x 10^9), incomplete Cholesky converged in 7 iterations against 197 for Jacobi and 313 unpreconditioned — roughly a 45-fold reduction — and was also the fastest in wall-clock terms at 0.054 seconds versus 0.176.
• The unpreconditioned solve on that matrix declared convergence at a relative residual of 9.41 x 10^-9, comfortably under the 10^-8 tolerance, yet its solution differed from the known exact answer by 1.00 in the infinity norm. With a condition number near 5 x 10^9, that residual permits a relative error of order 50, so a small residual is no certificate of a small error.
• Preconditioning restored accuracy as well as speed. Jacobi alone cut the solution error by more than six orders of magnitude to 5.33 x 10^-7, and incomplete Cholesky brought it to 1.94 x 10^-8, an accuracy matching the requested tolerance.
• On the half-million-unknown parabolic finite-element matrix (525,825 unknowns, 3,674,625 nonzeros, condition number about 2.10 x 10^5), incomplete Cholesky needed 66 iterations and 12.9 seconds against 1,321 iterations and 18.8 seconds for Jacobi and 1,951 iterations and 36.8 seconds unpreconditioned.
• Iteration savings do not translate one-for-one into time savings. The 30-fold iteration reduction on the large matrix became only a 2.9-fold solve-time reduction, because each preconditioned step applies triangular solves with factors 6.2 times denser than the matrix, and the factorization itself costs about 11 seconds up front.
• Where conditioning is moderate, residual and error stay coupled: all three variants on the finite-element matrix reached a solution error near 2 x 10^-10 regardless of iteration count, so preconditioning bought speed rather than correctness.

How it was done

Two symmetric positive-definite matrices were drawn from the SuiteSparse Matrix Collection and independently verified rather than trusted from metadata — symmetry checked exactly, positive-definiteness confirmed by dense Cholesky for the small matrix and by shift-invert Lanczos for the large one. Each right-hand side was formed as the matrix times the all-ones vector, so the exact solution is known and the true error is observable at every iteration alongside the residual. A hand-written preconditioned conjugate gradient iteration, recording the relative residual at each step and recomputing it explicitly at termination, was run to a relative residual below 10^-8 with a cap of 20,000 iterations that no solve reached. The incomplete-Cholesky preconditioner was built as an exactly symmetric operator from the incomplete lower factor and its pivots — verified symmetric to machine precision — under a symmetric minimum-degree fill-reducing ordering, and condition numbers were estimated from extreme eigenvalues computed by Lanczos with shift-invert for the smallest.

Data sources

• SuiteSparse Matrix Collection — bcsstk16, Harwell-Boeing group, 4,884 unknowns, 290,378 nonzeros
• SuiteSparse Matrix Collection — parabolic_fem, Wissgott group, 525,825 unknowns, 3,674,625 nonzeros
• Hestenes & Stiefel, Journal of Research of the NBS 49:409 (1952) — the conjugate gradient method
• Meijerink & van der Vorst, Mathematics of Computation 31:148 (1977) and Kershaw, Journal of Computational Physics 26:43 (1978) — incomplete Cholesky preconditioning