36
$\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?

$\endgroup$
  • 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 '12 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 '12 at 18:45
  • 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 '12 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 '12 at 13:16
  • 1
    $\begingroup$ It calls a Fortran code written using Intel MKL. $\endgroup$ – Inquest Jan 26 '12 at 13:34
37
$\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.

$\endgroup$
  • $\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 '12 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 '12 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 '12 at 20:25
37
$\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.

$\endgroup$
  • 3
    $\begingroup$ +1 for new MATLAB settings I'd never heard of. Well played, sir. $\endgroup$ – Geoff Oxberry Feb 28 '12 at 1:35
  • 2
    $\begingroup$ Hey thanks! It's in help mldivide. $\endgroup$ – dranxo Feb 28 '12 at 3:43
16
$\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.

$\endgroup$
  • $\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 '12 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 '12 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 '12 at 19:26
  • $\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 '12 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.

$\endgroup$
  • 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 '12 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 '12 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 '12 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 '12 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 '12 at 16:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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