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

16 visualizaciones (últimos 30 días)
Jianqiao Huang
Jianqiao Huang el 12 de Ag. de 2019
Editada: LEE HUANG el 16 de Ag. de 2019
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

Respuestas (3)

John D'Errico
John D'Errico el 13 de Ag. de 2019
Editada: John D'Errico el 13 de Ag. de 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 comentarios
LEE HUANG
LEE HUANG el 14 de Ag. de 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 el 16 de Ag. de 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

Iniciar sesión para comentar.


Walter Roberson
Walter Roberson el 13 de Ag. de 2019
Since your A is staying the same, you might be able to take advantage of decomposition(A')
  1 comentario
LEE HUANG
LEE HUANG el 13 de Ag. de 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

Iniciar sesión para comentar.


Christine Tobler
Christine Tobler el 13 de Ag. de 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 comentarios
Christine Tobler
Christine Tobler el 15 de Ag. de 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 el 16 de Ag. de 2019
Editada: LEE HUANG el 16 de Ag. de 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

Iniciar sesión para comentar.

Categorías

Más información sobre Elementary Math en Help Center y File Exchange.

Community Treasure Hunt

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

Start Hunting!

Translated by