## Is it possible to speed-up solving Ax=b for multiple b's by pre-computing the Cholesky factorization of A?

### David DeVries (view profile)

on 21 Sep 2017
Latest activity Edited by Cedric Wannaz

on 22 Sep 2017

### Christine Tobler (view profile)

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?

Cedric Wannaz

on 22 Sep 2017

### Christine Tobler (view profile)

on 22 Sep 2017

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

### Walter Roberson (view profile)

on 22 Sep 2017
And +1 for the up to date information on the new facility!
David DeVries

on 22 Sep 2017