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

7 次查看（过去 30 天）
David DeVries
David DeVries 2017-9-21
编辑： Cedric 2017-9-22
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?

请先登录，再回答此问题。

采纳的回答

Christine Tobler
Christine Tobler 2017-9-22
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;
  2 个评论
David DeVries
David DeVries 2017-9-22
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".

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Creating and Concatenating Matrices 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by