How can I solve a sparse linear system efficiently in MATLAB?

17 views (last 30 days)
Hi,
I need to calculate A*W*A' in each iteration, where W=(B)^(-1), A is a m by n matrix, and W is a n by n matrix. B is updated in each iteration. And n=7500, m=11250.
How can I calculate A*W*A' efficiently?
Now what I do is as follows:
A*B\A.
But it takes more than 20 seconds to compute B\A. I need to run a lot of iterations, so it's time consuming.
Thanks,
Jianqiao Huang

Answers (3)

John D'Errico
John D'Errico on 13 Aug 2019
Edited: John D'Errico on 13 Aug 2019
B changes in each iteration. It is different. We are given no clue as to how it changes, or how much. We don't even know anything about the properties of B. But even if the perturbation is small, the difference in the solution from the previous iteration can still be quite large, especialy so if B is poorly conditioned.
That means the system of equations you need to solve changes in each iteration. Sorry. But you just need to accept that solving large problems takes large time.
For example, you might decide to try using an iterative solver, that would employ the previous solution for the last problem as the start point. But unless the matrix is quite well posed then the solution can be quite different, so you would not gain much from the starting values so provided.
Is B a symmetric, positive definite matrix? If so, then you might gain, by use of a cholesky factorization of B, in the sense that the linsolve tool can be used. Linsolve does not test that the matrix has the indicated property, so the small amount of time saved might save you a little bit.
Does B not change much at all? There are tricks you might decide to try using iterative solvers and a very intelligently chosen preconditioning matrix.
All of the above ideas depend very much upon properties of B that we are not told. They may also require some degree of knowledge of your part to use the necessary tools well. A good choice of preconditioning matrix can make a large difference for an iterative solver, for example.
(Note that a decomposition of A is useless, because B is the matrix that must be factored here.)
  4 Comments
LEE HUANG
LEE HUANG on 14 Aug 2019
Hi Christine Tobler,
You're right.
A is not symmetric positive definite, but B is symmetric positive definite.
Thanks,
Jianqiao Huang
LEE HUANG
LEE HUANG on 16 Aug 2019
Thanks for your comments. I make improvements now.
The computation time of calculating B\A is reduced from 22 seconds to 2.9 seconds.
What I do is as follows.
(1) B=sparse(B);
(2) B=(B+B')/2;
(1) can reduce computation time to 18 secnods, and (2) can reduce computation time to 2.9 seconds.
B is symmetric theretically, but it's not symmetric numerically. Norm(B-B',inf)=1.5*e-5.
Is there any other way to handle numercial issue instead of using B=(B+B')/2?
Thanks,
Jianqiao Huang

Sign in to comment.


Walter Roberson
Walter Roberson on 13 Aug 2019
Since your A is staying the same, you might be able to take advantage of decomposition(A')
  1 Comment
LEE HUANG
LEE HUANG on 13 Aug 2019
Thanks. I don't need to calculate inv(A), so I don't think decomposition(A') works here.
What do you think?
Thanks,
Jianqiao Huang

Sign in to comment.


Christine Tobler
Christine Tobler on 13 Aug 2019
Do you need the whole matrix A*inv(B)*A', or do you just need to apply this matrix to a vector (A*inv(B)*A'*x)? In the second case, you could change the order of operations to be A*( B \ (A'*x) ), which would mean solving a linear system B*y = A'*x with just one vector, instead of solving B*y = A' which has many right-hand side vectors.
I can't recommend much else without additional information: What are you doing with A*inv(B)*A' once you've computed it? Are m and n about the same size, or is one typically larger than the other? Are either A or B sparse / symmetric / any other special structure?
  7 Comments
Christine Tobler
Christine Tobler on 15 Aug 2019
You were saying you need to compute inv(C)*b, where C = (A*W*A').*(A*W*A'), so I rewrote that this way. So since you're not solving a linear system with the matrix C.*C, what operations are you doing with this matrix? In some (rare) cases, it might be possible to avoid computing the matrix explicitly at all.
Also, C = A\B*A' is equivalent to inv(A)*B*A', while your original post mentioned A*B\A, which would be equivalent to A*inv(B)*A.
I'm having a hard time understanding your algorithm, partially because the variable names and exact formulas seem to change in every post. Could you send a complete code (as it pertains to these two matrices)?
LEE HUANG
LEE HUANG on 16 Aug 2019
Edited: LEE HUANG on 16 Aug 2019
Sorry. I did make mistakes in the formula. It's my first time to use MATLAB Central, so I missed some older comments.
Let me clarify two points as follows.
(1) C = (A*B\A').*(A*B\A')
(2) x=C\b=[(A*B\A').*(A*B\A')]\b
Is it possible to avoid calculating C explicitly here?
Moreover, your advice of Cholesky decomposition of B = R*R' is good. In my case, B is symmetrically theretically, but not numercially. It takes very long time to compute B\A, but it's very quickly to compute Cholesky decomposition of B. I use B=(B+B')/2 to make B symmetrical. Is there any other way to handle the numercial issue?
Thanks,
Jianqiao Huang

Sign in to comment.

Categories

Find more on Elementary Math in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!