How to solve this Ax=b faster?

Specific Domains Numerics
,
21

Oh, that completely changes the game. Yes, Cholesky factorizations are going to a way faster than any other factorization if the matrix is SPD! I didn’t catch that detail (now that I re-read, there was a quick mention of SPD).

Yes, we should wrap CHOLMOD for this case. Also, for SPD you shouldn’t use GMRES, instead you should use CG (KrylovJL_CG).

2 Likes
22

here I uploaded a JLD2 file and a small script that runs the tests.

following the output of the script

loading variables ... 
A, r, 
size(A) = (162372, 162372)
ishermitian(A) = true
norm(r) = 9027.184755276294

doing u = A\r ...
  3.878 s (52 allocations: 1.39 GiB)
norm(r - A * u) = 5.766109468498555e-12

doing u = lu(A)\r ...
  13.797 s (72 allocations: 2.39 GiB)
norm(r - A * u) = 4.227052608619862e-12

doing LinearSolve.solve(prob, UMFPACKFactorization()) ...
  12.981 s (94 allocations: 2.71 GiB)
norm(r - A * sol.u) = 4.227052608619862e-12

doing LinearSolve.solve(prob, KrylovJL_CG()) ...
  70.360 s (50 allocations: 179.79 MiB)
norm(r - A * sol.u) = 0.00013389304531825536

you can tell KrylovJL_CG() is a conjugate gradient method by looking at how little memory it allocates, I would also have expected a CG method would have been faster on a SPD problem (probably the the matrix is not big enough to outperform a very well optimized direct method)

1 Like
23

@Andrea_Vigliotti, the convergence rate of CG hightly depends on the distribution of the eigenvalues of A.
If you have a good preconditioner, the Krylov method could only requires a few iterations.

2 Likes
24

The \ is probably running cholesky, which is about 3x - 4x faster than lu.

1 Like
25

you are right, u_ch=cholesky(A)\r and u=A\r on this matrix take the same time, do the same number of allocation and find the same solution, norm(u_ch-u) is 0.0