MATLAB: Further improvement of matrix inversion

efficiencymatrix inversionmex filespeed

Using A\b instead of A^(-1) can be used to speed up a code. 1) How can we further speed up this inversion if we know from the beginning that Cholesky decomposition can apply in our matrix? 2)If we are about to use the inverse matrix to multiply it for more than one b vectors then A\b is still faster or we should save the inverse A^(-1) and then only do the multiplications? Is there a general rule for the number of inversions that one way is faster than the other? 3)If we create a mex file from a function that only contains the A\b command can we speed up the code execution or A\b is fully optimized?

Thank you in advance!

Best Answer

How can we further speed up this inversion if we know from the beginning that Cholesky decomposition can apply in our matrix?

You could possibly speed things up by a small fraction by skipping some of the up-front checking, but I don't know if the effort would be worth it.

If we are about to use the inverse matrix to multiply it for more than one b vectors then A\b is still faster or we should save the inverse A^(-1) and then only do the multiplications? Is there a general rule for the number of inversions that one way is faster than the other?

It would be faster to simply combine all of your b vectors into one matrix and then do A\bcombined.

If we create a mex file from a function that only contains the A\b command can we speed up the code execution or A\b is fully optimized?

For speed & accuracy the mex file would likely simply call the exact same BLAS and LAPACK library routines that MATLAB is already calling. So no speed improvement would be expected.

Related Solutions

MATLAB: Updating inverse of a lower triangular matrix

NO. Computing the inverse of A, if it is lower triangular is a BAD idea.

Solving the problem x = A\b is a forward substitution, so fast as hell. On the order of a matrix*vector multiply in terms of the computational load. I.e., essentially an O(n^2) operation.

Wasting the time to do do an update of an inverse matrix, so that you can then do a matrix multiply is just silly.

A = tril(rand(10000));
b = rand(10000,1);
timeit(@() A*b)
ans =
     0.067386
timeit(@() A\b)
ans =
     0.090263

So the backslash took only slightly more time than a multiply. The update of a matrix inverse would take way more time than that. And for no good reason at all, only to still need to do a matrix*vector multiply.

For example, one can compute the change to a matrix inverse under a rank 1 update. (which is exactly what you want to do.) The classic solution is due to Sherman & Morrison. (I remember Woodbury in the name when I learned about it a zillion years ago, but Sherman-Morrison is a special case.)

https://en.wikipedia.org/wiki/Sherman–Morrison_formula

Now, look carefully at the computations in the formula shown. While you could do them, you need to recognize that they would involve more work than a simple matrix*vector multiply would entail. And since the forward substitution itself is trivial, you cannot possibly gain here.

MATLAB: Matrix inverse using Cholesky decomposition

If you have the Cholesky factor, then just use back substitution (twice), applied to an identity matrix. Actually, it will be a forward substitution, then a back subs.

A = L*L'

where L is a lower triangular matrix. So, the problem is to compute the inverse, Ainv here:

L*L'*Ainv = eye(n,n)

Think of it as first solving the problem

L*u = eye(n,n)

L is lower triangular. The right hand side is known, the identity. Apply forward substitution to each column of the RHS.

Then solve a second problem.

L'*Ainv = u

Again, u is now completely known. L' is a upper triangular matrix. Use back substitution to solve for the square matrix Ainv.

Easy, peasy. If you got as far as the Cholesky factor, then this part will be easy for you.

When you are done, VERIFY YOUYR RESULT! Test that

A*Ainv = eye(n,n)

Don't worry if there is some floating point trash in the off diagonal elements of the product if they are on the order of 1e-16. That is to be expected.

Best Answer

Related Solutions

MATLAB: Updating inverse of a lower triangular matrix

MATLAB: Matrix inverse using Cholesky decomposition

Related Question