40
$\begingroup$

I was comparing a few of my codes to "stock" MATLAB codes. I am surprised at the results.

I ran a sample code (Sparse Matrix)

n = 5000;
a = diag(rand(n,1));
b = rand(n,1);
disp('For a\b');
tic;a\b;toc;
disp('For LU');
tic;LULU;toc;
disp('For Conj Grad');
tic;conjgrad(a,b,1e-8);toc;
disp('Inv(A)*B');
tic;inv(a)*b;toc;

Results :

    For a\b
    Elapsed time is 0.052838 seconds.

    For LU
    Elapsed time is 7.441331 seconds.

    For Conj Grad
    Elapsed time is 3.819182 seconds.

    Inv(A)*B
    Elapsed time is 38.511110 seconds.

For Dense Matrix:

n = 2000;
a = rand(n,n);
b = rand(n,1);
disp('For a\b');
tic;a\b;toc;
disp('For LU');
tic;LULU;toc;
disp('For Conj Grad');
tic;conjgrad(a,b,1e-8);toc;
disp('For INV(A)*B');
tic;inv(a)*b;toc;

Results:

For a\b
Elapsed time is 0.575926 seconds.

For LU
Elapsed time is 0.654287 seconds.

For Conj Grad
Elapsed time is 9.875896 seconds.

Inv(A)*B
Elapsed time is 1.648074 seconds.

How the heck is a\b so awesome?

Improve this question
$\endgroup$
8
  • 1
    $\begingroup$ The built-in backslash of MATLAB, in other words the direct solver for a linear equations system, uses Multifrontal method for sparse matrix, that is why A\B is so awesome. $\endgroup$
    – Shuhao Cao
    Jan 25, 2012 at 17:31
  • 1
    $\begingroup$ It uses Tim Davis's code available at cise.ufl.edu/research/sparse. Also the awesomeness goes away when you have a non-trivial problem $\endgroup$
    – stali
    Jan 25, 2012 at 18:45
  • 1
    $\begingroup$ That link deals with Sparse matrices. But, repeatedly, I have noticed that the fastest mode of solving dense matrices as well is a\b. Any directions for that? $\endgroup$
    – Inquest
    Jan 25, 2012 at 19:25
  • 1
    $\begingroup$ What is "LULU"? Why do you think it is a good implementation of an LU factorization and subsequent direct solve? $\endgroup$
    – Jed Brown
    Jan 26, 2012 at 0:26
  • 3
    $\begingroup$ @Nunoxic What implementation? Did you write it yourself? High-performance dense linear algebra, while usually well-understood algorithmically, is not easy to implement efficiently on modern hardware. The best BLAS/Lapack implementations should get close to peak for a matrix of that size. Also, from your comments, I'm getting the impression that you think LU and Gaussian Elimination are different algorithms. $\endgroup$
    – Jed Brown
    Jan 26, 2012 at 13:16

4 Answers 4

Reset to default
41
$\begingroup$

In Matlab, the ‘\’ command invokes an algorithm which depends upon the structure of the matrix A and includes checks (small overhead) on properties of A.

  1. If A is sparse and banded, employ a banded solver.
  2. If A is an upper or lower triangular matrix, employ a backward substitution algorithm.
  3. If A is symmetric and has real positive diagonal elements, attempt a Cholesky factorization. If A is sparse, employ reordering first to minimize fill-in.
  4. If none of criteria above is fulfilled, do a general triangular factorization using Gaussian elimination with partial pivoting.
  5. If A is sparse, then employ the UMFPACK library.
  6. If A is not square, employ algorithms based on QR factorization for undetermined systems.

To reduce overhead it is possible to use the linsolve command in Matlab and select a suitable solver among these options yourself.

Improve this answer
$\endgroup$
3
  • $\begingroup$ Assuming I am dealing with a 10000x10000 unstructured dense matrix with all elements nonzero (High level of density), what would be my best bet? I want to isolate that 1 algorithm which works for dense matrices. Is it LU, QR or Gaussian Elimination? $\endgroup$
    – Inquest
    Jan 25, 2012 at 20:14
  • 1
    $\begingroup$ Sounds like a Step 4 where Gaussian Elimination is invoked which corresponds to the most general case where no structure of A can be exploited to boost performance. So, basically this is a LU factorization and subsequent one forward followed by a backward substitution step. $\endgroup$
    – Allan P. Engsig-Karup
    Jan 25, 2012 at 20:18
  • $\begingroup$ Thanks! I think that gives me a direction to think. Currently, Gaussian Elimination is the best we have for solving such unstructured problems, is that correct? $\endgroup$
    – Inquest
    Jan 25, 2012 at 20:25
39
$\begingroup$

If you want to see what a\b does for your particular matrix you can set spparms('spumoni',1) and figure exactly what algorithm you were impressed by. For example:

spparms('spumoni',1);
A = delsq(numgrid('B',256));
b = rand(size(A,2),1);
mldivide(A,b);  % another way to write A\b

will output

sp\: bandwidth = 254+1+254.
sp\: is A diagonal? no.
sp\: is band density (0.01) > bandden (0.50) to try banded solver? no.
sp\: is A triangular? no.
sp\: is A morally triangular? no.
sp\: is A a candidate for Cholesky (symmetric, real positive diagonal)? yes.
sp\: is CHOLMOD's symbolic Cholesky factorization (with automatic reordering) successful? yes.
sp\: is CHOLMOD's numeric Cholesky factorization successful? yes.
sp\: is CHOLMOD's triangular solve successful? yes.

so I can see that "\" ended up using "CHOLMOD" in this case.

Improve this answer
$\endgroup$
2
  • 4
    $\begingroup$ +1 for new MATLAB settings I'd never heard of. Well played, sir. $\endgroup$
    – Geoff Oxberry
    Feb 28, 2012 at 1:35
  • 2
    $\begingroup$ Hey thanks! It's in help mldivide. $\endgroup$
    – dranxo
    Feb 28, 2012 at 3:43
17
$\begingroup$

For sparse matrices, Matlab uses UMFPACK for the "\" operation, which, in your example, basically runs through the values of a, inverts them, and multiplies them with the values of b. For this example, though, you should use b./diag(a), which is a lot faster.

For dense systems, the backslash-operator is a bit more complicated. A brief description of what is done when is given here. According to that description, in your example, Matlab would solve a\b using backward substitution. For general square matrices, an LU-decomposition is used.

Improve this answer
$\endgroup$
4
  • $\begingroup$ Regd. Sparsity, the inv of a diag matrix would simply be the reciprocal of the diagonal elements so b./diag(a) would work but a\b works awesomely for general sparse matrices as well. Why isn't linsolve or LULU (My optimized version of LU) not faster that a\b for dense matrices in that case. $\endgroup$
    – Inquest
    Jan 25, 2012 at 18:14
  • $\begingroup$ @Nunoxic Does your LULU do anything to detect diagonality or triangularity of the dense matrix? If not, it will take just as long with every matrix, irregardless of its contents or structure. $\endgroup$
    – Pedro
    Jan 25, 2012 at 18:34
  • $\begingroup$ Somewhat. But, with the linsolve OPT flags, I defined everything there is to define about the structure. Yet, a\b is faster. $\endgroup$
    – Inquest
    Jan 25, 2012 at 19:26
  • 1
    $\begingroup$ @Nunoxic, every time you call a user function, you cause overheads. Matlab does everything for backslash internally, e.g. the decomposition and subsequent application of the right hand side, with very little overhead, and will thus be faster. Also, in your tests, you should use more than just one call to obtain reliable timings, e.g. tic; for k=1:100, a\b; end; toc. $\endgroup$
    – Pedro
    Jan 25, 2012 at 22:02
5
$\begingroup$

As a rule of thumb, if you have a sparse matrix of reasonable complexity (i.e., it doesn't have to be the 5-point stencil but can in fact be a discretization of the Stokes equations for which the number of nonzeros per row is much larger than 5), then a sparse direct solver such as UMFPACK typically beats an iterative Krylov solver if the problem is no larger than around maybe 100,000 unknowns.

In other words, for most sparse matrices resulting from 2d discretizations, a direct solver is the fastest alternative. Only for 3d problems where you quickly get above 100,000 unknowns does it become imperative to use an iterative solver.

Improve this answer
$\endgroup$
5
  • 3
    $\begingroup$ It's not clear to me how this answers the question, but I also take issue with the premise. It's true that direct solvers tend to work well for modest size 2D problems, but direct solvers tend to suffer in 3D well before 100k unknowns (the vertex separators are much larger than in 2D). Furthermore, I claim that in most cases (e.g. elliptic operators), iterative solvers can be made faster than direct solvers even for moderate sized 2D problems, but it may take significant effort to do so (e.g. using the right ingredients to make multigrid perform). $\endgroup$
    – Jed Brown
    Jan 26, 2012 at 0:31
  • 1
    $\begingroup$ If your problems are reasonably small and you don't want to think about designing implicit solvers, or if your problems are very difficult (like high-frequency Maxwell) and you don't want to devote your career to designing good solvers, then I agree that sparse direct solvers are a great choice. $\endgroup$
    – Jed Brown
    Jan 26, 2012 at 0:33
  • $\begingroup$ Actually my problem is to design a good dense solver. I don't have a particular application as such (Hence, no structure). I wanted to tweak a few of my codes and make them competitive with what is currently used. I was just curious about a\b and how it always beats the crap out of my code. $\endgroup$
    – Inquest
    Jan 26, 2012 at 5:59
  • $\begingroup$ @JedBrown: Yes, maybe with a significant amount of effort one can get iterative solvers to beat a direct solver even for small 2d problems. But why do it? For problems with <100k unknowns, sparse direct solvers in 2d are plenty fast enough. What I wanted to say is: for such small problems, don't spend your time tinkering and figuring out the best combination of parameters -- just take the blackbox. $\endgroup$
    – Wolfgang Bangerth
    Feb 10, 2012 at 9:05
  • $\begingroup$ There is already an order of magnitude run-time difference between a sparse Cholesky and (matrix-based) geometric multigrid for "easy" 2D problems with 100k dofs using a 5-point stencil (~0.5 second compared to ~0.05 seconds). If your stencil uses second neighbors (e.g. fourth order discretization, Newton for some choices of nonlinear rheology, stabilization, etc), the size of the vertex separators roughly doubles, so the cost of the direct solve goes up by about 8x (cost is more problem-dependent for MG). With many time steps or optimization/UQ loops, these differences are significant. $\endgroup$
    – Jed Brown
    Feb 10, 2012 at 16:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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