Computational complexity of direct algorithms for solving linear systems of equations is O(n^{3}). The only way to reduce enormous effect of cubic growth is to take advantage of various matrix properties and use specialized solvers^{1}.

Toolbox follows this strategy by relying on poly-solver which automatically detects matrix properties and applies the best matching algorithm for the particular case. In this post we outline our solver architecture for full/dense matrices.

**Poly-solver for dense input matrices** *(click on flowchart to see it in higher resolution)*:

Toolbox analyses the structure of input matrix, computes its bandwidth, checks symmetry/Hermitianity, determines permutation to convert matrix to trivial cases (diagonal or triangular) if possible. Then solver selects the most appropriate algorithm depending on matrix properties. Special attention is paid to n-diagonal matrices (banded), frequently encountered in solution of PDE.

Solver relies on following decompositions (from slowest to fastest): rank-revealing QR, LU with partial pivoting, LDL^{T} and Cholesky^{2}. Each decomposition is implemented in several variants, tailored for particular matrix type (tridiagonal^{3}, n-diagonal, full) and enabled with parallelism – thus performance scales with number of CPU cores in all cases.

The need for such diversity of algorithms can be seen from timings comparison^{4}:

Solving `2000x2000`

matrix by LU under `25`

seconds is notable fact on its own. However timing can be reduced further to `17`

seconds (LDL^{T}) if matrix is symmetric, or even to `13`

seconds if matrix is positive definite (Cholesky).

In case of n-diagonal matrices, specialized solver delivers significantly better timings as well:

Performance gain is around `20`

times for `2000x2000`

matrices. But most importantly, speed-up increases with matrix size (since complexity growth of banded LU decomposition is sub-cubical).

Similarly to quadruple (34-digits), specialized solvers works well in higher precision:

Still we see substantial speed-up (`>4`

times) but it is lower since higher precision needs more memory and algorithms slowly become memory bound in addition to large constant in front of `n`

.^{3}

^{1}Iterative solvers use another strategy by relying only on matrix-vector operations of O(n

^{2}).

^{2}The LDL^{T} decomposition is considered not stable in some cases, but we rely on Bunch-Kaufman diagonal pivoting which is as stable as LU with pivoting.

^{3}Computation of reciprocal condition number might take up to `50%`

of total solver time, thus in trivial cases (including tridiagonal) it is not computed.

^{4}We used newest version of toolbox (`4.2.2.11893`

) and `MATLAB R2016a`

on `Core i7 990x`

to built the plots above.

{ 2 comments… read them below or add one }

For linear system solve for real symmetric positive definite matrix which is not-banded, flowchart shows Cholesky solver if Cholesky succeeds, and LDL’ solver if Cholesky fails. Is there a way to find out whether Cholesky succeeded without making a separate call [R,p] = chol(A) ?

Consider a case in which A is extremely ill-conditioned. Presumably, for a given A, there will be some minimum precision at which Cholesky succeeds.

@”Is there a way to find out whether Cholesky succeeded without making a separate call [R,p] = chol(A) ?”

Currently MLDIVIDE does not provide this information.

To avoid extra call for “chol” and to get precise control you can use something like:

This mimics the flowchart branch for symmetric matrices, but it might be slower than MLDIVIDE itself (since we do not store factors explicitly and there is no data exchange with MATLAB).