MATLAB: Is it possible to speed-up solving Ax=b for multiple b’s by pre-computing the Cholesky factorization of A

cholesky factorizationmldividesparse matrices

I'm currently working with large (15000×15000), sparse matrices that are not triangular nor banded, but are real, symmetric, Hermitian, and have all positive values. Following the flow chart for the mldivide function for sparse matrices, it would appear that the Cholesky solver would be used. I solve Ax=b for multiple b's, so I thought that if I pre-computed the Cholesky factorization once, that then solving the equation would be faster, so I decided to do some tests. These were the times I got for each step, and comparison to the usual mldivide time
>> tic; A_factor = chol(A, 'lower'); toc; %use 'lower' for speed boost with sparse matrices
Elapsed time is 51.273020 seconds.
>> tic; x = A_factor\b; toc; % pre-computed solution time
Elapsed time is 0.186439 seconds.
>> tic; x = A\b; toc; % not pre-computed solution time
Elapsed time is 3.275133 seconds.
It is clear that pre-computing the factorization helps (0.18s < 3.28s), but my question is why does the initial pre-computation of the Cholesky factorization take so long (~50s) when the mldivide call without the pre-computation takes only 3s and would be performing the same factorization. Wouldn't we expect the Cholesky factorization to take < 3s? Any ideas?

Best Answer

The main reason the first call is slower is that R = chol(A) for a sparse matrix tends to have a large amount of fill-in (that is R will have many more nonzeros than A). In mldivide, this is avoided by first permuting the elements of A, so that chol(permutedA) has much less fill-in than chol(A) does.
The equivalent of this can be achieved by using
 [R, p, S] = chol(A);
If p is zero (meaning that A is positive definite and CHOL succeeded), the output R'*R is equal to S'*A*S. This can be used to solve the linear system A*x = b as follows
 x = S*(R\(R'\(S'*b)))
By the way, with MATLAB R2017b, the new decomposition object precomputes the factors and solves linear systems more easily, for any matrix type
 dA = decomposition(A);
 x = dA \ b;
Related Question