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

2 views (last 30 days)
Alec Jacobson
Alec Jacobson on 24 Jun 2020
Commented: Christine Tobler on 30 Jun 2020
If I issue:
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;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:
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)
  • 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

Answers (1)

Christine Tobler
Christine Tobler on 25 Jun 2020
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.
Christine Tobler
Christine Tobler on 30 Jun 2020
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.

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!