Is it possible to recover *the* decomposition that's ultimately employed by mldivide (backslash)?

If I issue:

spparms('spumoni',1);
A\b;

Then I can see the decision made by mldivide on how to solve against A. In my case it chooses MA57:

sp\: bandwidth = 26633+1+26633.
sp\: is A diagonal? no.
sp\: is band density (0.00299541) > bandden (0.5) 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 or negative diagonal)? no.
sp\: But A is real symmetric; try MA57.
    

I would like to cache the decomposition (to conduct more solves against the same matrix later). If I call:

    D = decomposition(A)

I see

D = 

  decomposition with properties:

    MatrixSize: [27045 27045]
          Type: 'ldl'

Sifting through the documentation I can match up ldl to MA57, so I thought this meant I'm using the same decomposition as above.

However a few things hint otherwise:

• \ is much faster

tic;A\b;toc
tic;D = decomposition(A);toc
tic;D = decomposition(A,'ldl');toc

on my MacPro with 1.5TB memory and 28 cores produces:

Elapsed time is 1.226560 seconds.
Elapsed time is 2.343294 seconds.
Elapsed time is 2.440572 seconds.

• solving with the decomposition gives me a warning that \ doesn't:

D\b;

alerts me:

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND =  8.766761e-23. 
> In decomposition/warnIfIllConditioned (line 714)
  In  \  (line 404) 
  

Questions:

• Do decomposition and \ make the same decisions on the choice of decomposition to use?
• If not, is there way to get at whatever decomposition \ is actually using? (in my case this "faster" MA57 decomposition)

I'm imagining/wish for something like:

  [x,D] = A\b;

where D would contain the reusable decomposition created and used during \.

Here's A and b data

Answer 1

[Christine Tobler]

The same decomposition is used, the difference is in the estimation of the condition number. MA57 provides an estimate of the condition number which is dependent on the right-hand side vector used. Since this wasn't practical to use with the decomposition object, we call another estimator which unfortunately is more expensive. It's value will also not usually match the estimate from MA57 (which changes significantly based on the right-hand side vector), which would be why you're not getting that warning consistently.

You could turn off computing the condition number using decomposition(A, 'CheckCondition', false), this should show a noteable speedup in constructing the object.

[Alec Jacobson]

This is great!

tic;D = decomposition(A,'ldl','CheckCondition',false);D\b;toc

indeed now matches \ 's timing:

Elapsed time is 1.227668 seconds.

One quick followup, is there a complete set of flags/options that I can pass to decomposition to be sure that I always get the same decomposition that would have been used by / ?

[Christine Tobler]

It's all about compromise here: By using 'CheckCondition', false, every solver will be sped up, at the cost of never getting a warning about a matrix being badly conditioned. the decomposition that's used is the same, just the extra computation of the condition number is more expensive for MA57 than for other solvers.

[Alec Jacobson]

I guess I was hoping to make the experience identical (regardless of tradeoff/warnings). Is it documented or somehow discoverable how to set the flags of decomposition to mimic perfectly the behavior of \ ? Or is this a guess-and-check scenario?

[Christine Tobler]

Unfortunately it wasn't quite possible to make the behavior exactly identical. The problem is that there are several cases in mldivide's behavior where it's not possible to separate the factorization step and the solution step as neatly as decomposition requires. These are mostly things that can't be fixed by choosing different options, for that probably MA57 is the only one.

Some examples:

• If A is double but B is single, A\B will cast A to single and then factorize - but decomposition has already factorized A as double, so instead it casts B to double, solves and casts the result to single. That is A\B is the same as single(A)\B, but decomposition(A)\B is the same as single(A \ double(B)).
• For the sparse QR solver, mldivide uses a method that computes and discards the Q matrix on-the-fly, which is faster. This is not possible for decomposition, where the Q factor must be stored (implicitly in an efficient format). That means that decomposition with sparse QR is slower and does not match results from mldivide exactly.

In terms of performance, using 'CheckCondition', false will speed up all solvers to some extent, although it's most noticeable for sparse LDL.