Further improvement of matrix inversion

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!

[John Chilleri]

Not sure if this is the place to ask for help on such a problem, but here's a hint on number 1.

Ironically Linking to MATLAB Answers

Answer 1

[James Tursa]

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.

[acivil acivil]

Thank you for your answers. I want to clarify that for the second question I mean that I need as an output many different vectors as well. Is it faster if I store the inverse matrix and use it a few times for different vectors b and different output vectors?

[John D'Errico]

If you have many vectors, then just supply b as a matrix, with one b vector in each column. A\b will automatically solve the whole set at once.

Or, if you must solve the problems one at a time, so each b vector is created on the fly, then just factorize A using chol.

Since we know that L as returned from chol has the property that

A = L*L'

then solve each individual problem as

x_b = L'\(L\b);

Even better, use linsolve to do each of those solves, telling linsolve that the matrix is upper or lower triangular. That will be even more efficient.

There is never a reason to use inv, nor do you really want or need to do so.

[Stephen23]

"Is it faster if I store the inverse matrix and use it a few times for different vectors b and different output vectors"

Probably not. And even if it was faster, you also have to consider that calculating the inverse is notoriously numerically inaccurate compared to solving the system using \. Use \ because it gives a much better solution, exactly as the inv documentation recommends (and has been discussed on this forum many times as well): "A better way, from the standpoint of both execution time and numerical accuracy, is to use the matrix backslash operator x = A\b"

[Walter Roberson]

"Is it faster if I store the inverse matrix and use it a few times for different vectors b and different output vectors?"

No, not typically. Possibly if you were doing enough of them. Also, you lose accuracy / stability if you do this.

"I need as an output many different vectors as well"

From an efficiency standpoint, compute all of those different vectors first, and place each of them as a column vector in a b matrix and do one A\b and then extract the columns of the result.

Working iteratively would be for the case where you need one result in order to compute the b vector for the next A\b

[acivil acivil]

I am indeed referring to an iterative procedure! So as I can conclude, I should use cholesky factorization and then multiple times linsolve?

[John D'Errico]

Call linsolve with L and L', using the OPTS argument to indicate if the matrix is lower or upper triangular.

Linsolve is indeed faster than simple use of backslash, since it does not need to determine if the matrix is really lower or upper triangular. Not needing to check will give a little improve speed boost.

    The possible field names in OPTS and their corresponding matrix
    properties are:

    Field Name : Matrix Property
    ------------------------------------------------
    LT         : Lower Triangular
    UT         : Upper Triangular

Or, you can use this nice alternative from Tim Davis: factorize .

[acivil acivil]

Thank you very much for your answers!