Take the 2-minute tour ×
Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. It's 100% free, no registration required.

I am using MATLAB to solve a problem that involves solving $A\mathbf{x}=\mathbf{b}$ at every timestep ($\mathbf{b}$ changes with time). Right now I am accomplishing this by MATLAB's mldivide:

x = A\b

I have the flexibility to make as many precomputations as needed, so I am wondering if there is a faster and/or more accurate method than mldivide. What is typically done here? Thanks all!

share|improve this question
    
Do you have specific knowledge about the structure of $A$? For instance, is it symmetric? Positive definite? Tridiagonal? Orthogonal? –  Dominique Feb 16 '13 at 3:44
    
The matrix $A$ is a dense square matrix. –  Doubt Feb 16 '13 at 22:23
2  
If you have no other knowledge on $A$, the $LU$ factorization as described in the answer below is your best bet. –  Dominique Feb 16 '13 at 22:26
add comment

5 Answers

up vote 8 down vote accepted

The most obvious thing you can do is to precompute

[L,U] = lu(A) ~ O(n^3)

Then you just compute

x = U \ (L \ b) ~ O(2 n^2)

This would reduce the cost enormously and make it faster. Accuracy would be the same.

share|improve this answer
1  
Note, from the documentation, L isn't necessarily lower triangular. This answer would likely be faster than a direct solve, however I would be careful to make sure that the L \ b command is smart enough to know to solve L in the correct order (it likely is, but it doesn't say for certain in the documentation). –  Godric Seer Feb 16 '13 at 2:49
    
Yeah you're right, L is the product of a lower-triangular matrix and a permutation matrix. But I'll be damned if it doesn't recognize that all it has to do is backward substitution with L\b. Because I have seen this exact line being used in high-performance code by those I consider experts. –  Milind R Feb 16 '13 at 3:14
6  
mldivide recognizes permuted triangular matrices and does the right $O(n^{2})$ thing in solving such a system. However, in my experiments, this seems to make the solution process slower by a factor of 10 or so for matrix of size 2000 by 2000 to 10000 by 10000. Thus you'd be better off keeping track of the permutation explicitly by using [L,U,P]=lu(P). –  Brian Borchers Feb 16 '13 at 3:26
1  
Also, if your matrix is sparse, you should take advantage of the sparsity in solving the system. The simplest way to do this is to make sure that $A$ is stored in sparse format by using A=sparse(A) before computing the LU factorization. You can also try to permute the rows of A so as to reduce fill-in during the LU factorization. –  Brian Borchers Feb 16 '13 at 3:28
3  
@BrianBorcher As far as I know, the best way to keep track of permutation is [L,U,p] = lu(A,'vector'); x = U\(L\b(p)); See example 3 in the lu docs. –  Stefano M Feb 18 '13 at 14:05
show 1 more comment

We did some extensive computer labs in our scientific computing courses about this topic. For the "small" calculations we did there, Matlab's backslash operator was always faster than anything else, even after we had optimized our code as much as possible and re-orded all matrices beforehand (for instance with Reverse Cuthill McKee ordering for sparse matrices).

You can check out one of our lab instructions. The answer to your question is covered (shortly) on page 4.

A good book on the topic is written for instance by Cheney.

share|improve this answer
add comment

Suppose $A$ is a $n\times n$ dense matrix and you have to solve $Ax_i = b_i$, $i=1\dots m$. If $m$ is big enough then there is nothing wrong in

V = inv(A);
...
x = V*b;

Flops are $O(n^3)$ for inv(A) and $O(n^2)$ for V*b, therefore in order to determine the break-even value for $m$ some experimentation is needed...

>> n = 5000;
>> A = randn(n,n);
>> x = randn(n,1);
>> b = A*x;
>> rcond(A)
ans =
   1.3837e-06
>> tic, xm = A\b; toc
Elapsed time is 1.907102 seconds.
>> tic, [L,U] = lu(A); toc
Elapsed time is 1.818247 seconds.
>> tic, xl = U\(L\b); toc
Elapsed time is 0.399051 seconds.
>> tic, [L,U,p] = lu(A,'vector'); toc
Elapsed time is 1.581756 seconds.
>> tic, xp = U\(L\b(p)); toc
Elapsed time is 0.060203 seconds.
>> tic, V=inv(A); toc
Elapsed time is 7.614582 seconds.
>> tic, xv = V*b; toc     
Elapsed time is 0.011499 seconds.
>> [norm(xm-x), norm(xp-x), norm(xl-x), norm(xv-x)] ./ norm(x)
ans =
   1.0e-11 *
    0.1912    0.1912    0.1912    0.6183

In this trivial example $A^{-1}$ pre-computation is better than $LU$ forward and backward solution for $m>125$.

Some notes

For stability and error analysis please see the comments to this different answer, especially the one by VictorLiu.

The proposed timings are not "scientific" at all, but are meant to show that the approach proposed in the answer by Milind R, while it makes perfect sense if implemented in C or Fortran by calling relevant LAPACK and BLAS subroutines, may prove not so effective in Matlab, even for $m\ll n$.

Timing were performed with Matlab R2011b on a 12 core computer with a fairly constant UNIX load average of 5; best tic, toc time of three probes.

share|improve this answer
    
Indeed, there is far more parallelism available in a matrix-vector multiply than a triangular solver, so this should be even more apparent if the computations are done in parallel (multicore/GPU/etc...) in any way. –  Aron Ahmadia Feb 18 '13 at 9:24
    
@AronAhmadia I agree: break-even point estimates based only on operation count makes sense only for a serial implementation. –  Stefano M Feb 18 '13 at 9:47
1  
Note that things will be vastly different if the A matrix is sparse- the inverse will typically be quite dense, while the LU factors are typically reasonably sparse, tipping things back in the direction of LU being faster. –  Brian Borchers Feb 18 '13 at 15:34
1  
@Doubt accuracy was already addressed in the comments to this answer. In general there should be only a negligible loss of accuracy, unless you have a pathological $A$ matrix. –  Stefano M Feb 18 '13 at 16:38
1  
@MilindR of course precomputing inv(A) make sense only if you have to repeatedly solve $Ax=b$ for different $b$'s that are not all known at the same time. Usually if you have multiple right-hand sides, just collect them in a matrix $B$ and go with A\B. –  Stefano M Feb 19 '13 at 7:50
show 4 more comments

Take a look at this question, the answers show that mldivide is quite clever, and also gives suggestions as to how to see what Matlab uses to solve A\b. This may give you a hint regarding optimization options.

share|improve this answer
add comment

Using backslash is more or less equivalent to inv(A)*B, if you're coding it freely, the latter might be more intuitive. They're about the same (just different in how the computation is carried out), though you should check Matlab documentation for clarification.

To answer your question, backslash is generally fine, but it does depend on the properties of the mass matrix.

share|improve this answer
1  
Mathematically inv(A)*b is the same as \ however numerically, actually forming the inverse is both less efficient and less accurate. If you are working to learn the linear algebra, this may be acceptable, but I would argue you need a very good reason to form the inverse. –  Godric Seer Feb 16 '13 at 14:49
    
But why would you ever compute inv(A) since that alone is more expensive than A\b? –  Dominique Feb 16 '13 at 18:25
7  
@Godric: There's a recent paper that discusses the "myth" that inv(A)*b is less accurate: on ArXiv. Not saying that there's usually reason to compute the actual inverse, but just sayin'. –  Victor Liu Feb 17 '13 at 4:30
3  
@Dominique: Triangular solves are much less parallelizable than matrix-vector multiplication, and sophisticated preconditioned iterative methods often make use of direct methods on subdomains. It is often useful to explicitly form the inverses of a few modest-sized dense triangular matrices in order to improve the parallelism. –  Jack Poulson Feb 17 '13 at 5:27
    
@VictorLiu: Thank you for the article. I stand corrected on my accuracy statement (atleast for smart implementations of inv(A)). –  Godric Seer Feb 17 '13 at 18:40
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.