You don't need matrix inversion.

Kind of, sometimes anyway.

Math

Tricks

C++

Speed up inversion of definite matrices using the code in (6).

I suppose that, if you are reading this, I don't need to explain why you might want to speed up matrix inversion (or avoid it entirely) wherever possible. There is a trick you can use to speed up matrix inversion by something between

20\%-70\%

, depending on your use case.

I don't claim to have discovered this — in fact a specialized version of this alogrithm is commonly used to calculate the squared Mahalanobis distance. Another version, specialized for the inversion of positive definite Hermitian matrices, is probably used here. I just haven't seen this generalization anywhere.

Limitations

This trick only works for definite matrices and you need to know the definiteness. However, I think this is commonly the case (think covariances etc.).

Generalization

Let's assume you need to invert a matrix

M > 0

(for

M < 0

, just use

-M

) in order to calculate the result of some multiplication with another matrix:

\begin{split} A &= M^{-1}B \\ M &\in \text{GL}(n, \mathbb{R}), M > 0, B\in\mathbb{R}^{n\times k} \end{split}

(1)

This is a generalization of three common problems:

Inverting a Matrix: $B=I \implies A=M^{-1}$
Calculating a matrix-vector product of the inverse: $B = v \in \mathbb{R}^k \implies A = M^{-1}v$
Calculating a matrix-matrix of the inverse product as in (1)

You could implement a solution to (1) like this:

A = M.inverse() * B;

(2)

But should you?

Equations aren't instructions

A pitfall when implementing equations in code is that there is a difference between '=' and '

=

'. The former is an assignment; the latter tells us that its left side is equal to its right side. This allows us to transform (1) into

MA = B

(3)

which we can now solve for

A

. Intuition may tell you that this implies going back to (1), which leads to the code in (2) — not helpful.

Now, the trick: Let your computer solve for

A

. First, perform a Cholesky decomposition on

M

\begin{split} M &= LL^\top \implies A = (LL^\top)^{-1}B = L^{-\top}L^{-1}B \end{split}

(4)

So I have made the problem worse

In addition to the decomposition, there are now more matrix multiplications and the number of inverses hasn't changed. Let's avoid any matrix inversion by applying the same logic from above to solve:

\begin{split} C = L^{-1}B \iff LC = B \\ A = L^{-\top}C \iff L^\top A = C \end{split}

(5)

Isn't this the exact problem that we attempted to solve in the beginning… but twice?

Actually, no

The matrix

L

has a triangular shape. This means (5) can be trivially solved using forward substitution. An implementation might look like this:

LLT = M.llt();
C = LLT.matrixL().solve(B);
A = LLT.matrixL().transpose().solve(C);

(6)

Is this faster?

Before going into the benchmarks, keep in mind that classic matrix inversion and the technique involving Cholesky decomposition both operate in

O(n^3)

I ran some benchmarks on these problems:

Computing the inverse $M^{-1}$ using $B=I$
Computing $M^{-1}v$ where $v$ is a vector
Computing $M^{-1}B_k$ where $B \in \mathbb{R}^{n\times k}$

Benchmarks were implemented in C++ using Eigen and can be found here. I chose random positive definite matrices

M

, random values for

B

and measured the computation time relative to the respective implementation using the .inverse() method.

Here are the relative computation times:

	$M^{-1}$	$M^{-1}B_{k=n}$	$M^{-1}B_{k=4}$	$M^{-1}v$
$n=8$	$70\%$	$54\%$	$50\%$	$45\%$
$n=20$	$77\%$	$55\%$	$38\%$	$30\%$

As you can see from the benchmarks, the greatest benefit can be gained when calculating a matrix-matrix product, especially as the number of columns

k

B

decrease.