Asked by David DeVries
on 21 Sep 2017

I'm currently working with large (15000x15000), 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?

Answer by Christine Tobler
on 22 Sep 2017

Accepted 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;

Walter Roberson
on 22 Sep 2017

And +1 for the up to date information on the new facility!

David DeVries
on 22 Sep 2017

Thanks Christine! Spot on advice!

For those interested, the permuted Cholesky factorization had a run time of ~4.1s, and then finding solutions was ~0.2s, so a definite speed-up if you're running multiple solutions with more than 1 "b".

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 1 Comment

## Cedric Wannaz (view profile)

## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/357776-is-it-possible-to-speed-up-solving-ax-b-for-multiple-b-s-by-pre-computing-the-cholesky-factorization#comment_486787

Sign in to comment.