Matrix inverse

Algorithm

ALGLIB contains recursive matrix inversion algorithm. A, given by its triangular factorization (LU-decomposition or Cholesky decomposition) is partitioned into four parts and inverted according to the next formula:

Matrix-matrix products of general form and multiplications by U^-1/L^-1 are calculated with cache-efficient Level 3 ALGLIB BLAS subroutines. Expressions like (LU)^-1 are calculated recursively. Base cases (matrix size is small enough to fit into L1 cache) are handled by another algorithm, which is based on explicit inversion of triangular factors of A. Recursive algorithm has following benefits:

It makes efficient use of CPU cache without explicitly specifying its size. With each iteration matrices become smaller until they fit completely in cache. Such algorithms are called 'cache oblivious'.
Recursive algorithm makes extensive use of ALGLIB BLAS - mostly of the Level 3 BLAS, the most optimized part of ALGLIB package. The matrices being multiplied are usually square (which is more suited for optimization than a low rank updates generated by LAPACK routines).

These features allows to achieve the maximum performance possible with your compiler/language. However, different compilers and different programming languages have different limitations. As for linear algebra algorithms, C++ version of ALGLIB is the most efficient implementation. Implementations in other languages have significantly lower performance.

Subroutines

ALGLIB package contains following matrix inversion subroutines:

rmatrixinverse
spdmatrixinverse
cmatrixinverse
hpdmatrixinverse
rmatrixluinverse
spdmatrixcholeskyinverse
cmatrixluinverse
hpdmatrixcholeskyinverse

First four subroutines are used to invert real (R prefix), symmetric positive definite (SPD prefix), complex (C prefix) or Hermitian positive definite (HPD prefix) matrices. Next four subroutines solve same task, but for matrices given by corresponding triangular factorization (LU decomposition or Cholesky decomposition).

All subroutines check condition number before inversion begins. If matrix is degenerate (one of triangular factors have zero element on its diagonal) or is very ill-conditioned (condition number is so high that there is a risk of overflow during inversion) then subroutine returns -3 error code and zero matrix.