Fast calculate the inverse of a lower triangular matrix

Hello, Everyone,

I need to write a code on myself to calculate the inverse of a lower triangular matrix. The matrix is a 3000 x 3000 sparse lower triangular matrix.

I found inv.m took 0.015 seconds to finish the calculation, but my code written using Gaussian elimination takes much more time.

Could anyone tell me how to speed up the time of my code?

Thanks.

Benson

[Steven Lord]

I need to write a code on myself to calculate the inverse of a lower triangular matrix.

Do you really need the inverse of the matrix or do you need to solve a system of equations whose coefficient matrix is that lower triangular matrix? In the latter case, don't invert the matrix. Use the backslash operator \ instead.

Could anyone tell me how to speed up the time of my code?

Without seeing your code, it's highly unlikely anyone could give you any concrete suggestions.

Answer 1

[John D'Errico]

Do not write your own inverse code. Gaussian elimination is great for teaching a student about matrices, but too often they then think they know it all. (This can extend to the doctoral level and beyond.) You don't. And you certainly do not want to write it in MATLAB on a sparse matrix.

Is inv too slow? Well, big problems take time. Your computer is not infinitely fast or infinitely powerful. Not even close.

Most important is to understand why you think you need to compute that inverse! Almost always, you don't. If 0.015 seconds is too slow, then you are doing it a zllion times. So what are you doing with that inverse? Are you multiyplying another matrix with it? WRONG!!!!!!

Use backslash instead. For example...

A = sparse(triu(tril(rand(3000)),-2));

So A is sparse lower triangular. Compare these times...

V = rand(3000,1);
timeit(@() inv(A)*V)

ans = 0.2030

timeit(@() A\V)

ans = 3.5245e-05

DON'T COMPUTE AN INVERSE. You very rarely need to compute an inverse, even if a formula has the inverse written. That inverse is just a nice mathematical way of writing what you need to do, in MATLAB, usually performed best using backslash. If you are post multiplying by the inverse, then forward slash would apply. If you want better help than this, then you need to show exactly what you are doing with that inverse.

Far too often people think just because a formula shows an inverse, that they need to compute that inverse.

[Dyuman Joshi]

I don't believe OP thinks that inv() is slow at all, just the code they wrote for calculating the inverse is very slow.

But yeah, I agree with you (and Steven) that OP does not really need to calculate the inverse (except in the rare case that the task is to specifically find the inverse).

[John D'Errico]

I hear what you say, but I don't come to the same conclusion. It is just a hunch though, and my hunch could be wrong.

If @Benson Gou truly needs to compute an inverse, then inv exists already. There is no need to write inv in that case. This implies they think inv was too slow for their case, so they tried to write it themselves. And that is pretty much always the wrong thing to do.

A big problem is the inverse of a lower tridiagonal matrix is no longer going to be sparse at all. There will be massive fill-in. For example, on a much smaller matrix, but the size is not relevant.

A = sparse(triu(tril(rand(100)),-2));
spy(A)

Ainv = inv(A);
spy(Ainv)

And that will be just bad news, no matter what you are doing.

[Benson Gou]

Thank you all for my question. My problem is that I do need to write code to calculate the inverse of a lower triangular matrix. Because Matlab does not have a specific code to calculate the inverse of lower triangular matrix, I think I can write a specific code to calculate the inverse of a lower triangular matrix. But I found my code is very slow comparing to inv.m. So I want to know why inv.m is too fast and my code is very slow comparing to inv.m.

Could you please tell me how inv.m calculated the inverse? I want to apply the code to a lower triangular matrix.

Thanks a lot for your great help.

Benson

[Dyuman Joshi]

Could you please tell me how inv.m calculated the inverse?

From the Algorithm section of the inv documentation page -

"inv performs an LU decomposition of the input matrix (or an LDL decomposition if the input matrix is Hermitian). It then uses the results to form a linear system whose solution is the matrix inverse inv(X). For sparse inputs, inv(X) creates a sparse identity matrix and uses backslash, X\speye(size(X))."

Note the last sentence.

[John D'Errico]

As I have said, almost always you only think you need to compute the inverse. You have not said why you think that. Just because you have a formula for it, does not mean you need to compute the inverse. I'm sorry, but that is flat out wrong, a complete misconception of many new users of MATLAB, and of numerical computing in general.

As well, you have not said why you cannot just use inv. There is no way in hell you can write code that will run faster than does inv. If you are able to use code like LU, etc., then there is no reason why you cannot use inv anyway, and certainly no reason why you cannot use backslash as I showed.

Next, as I showed, ALWAYS, the inverse of a tridiagonal matrix will no longer be a sparse matrix. And that means your computation of the inverse will be slow, no matter how you write it. HOWEVER, the algorithm used by backslash is NOT slow at all. Again, I showed exactly that.

So you can continue asserting that you know what you need to do, and that you need to do it as you claim. That does not make your blanket assertions valid, made with no cogent arguments for your need,

I'm sorry, but you simply cannot write code by yourself to compute an inverse faster than inv could do it. And since you would apply that code to a sparse matrix, which in turn will be forced to create a full matrix, the result will be slow. The fill-in will be massive. As much as you think you want to do that, you don't, and you can't.

[Steven Lord]

My problem is that I do need to write code to calculate the inverse of a lower triangular matrix.

Why? What problem are you trying to solve that requires inverting the matrix? If you explain more about the problem you're trying to solve in this way, we may be able to suggest an alternate solution that achieves your goal in a more efficient way.

[Benson Gou]

Thank you all for your answers and suggestions.

I actually want to implement inv.m in C++ for a lower triangular sparse matrix. In C++ code, LU decomposition does not take a big CPU time, but the backslash takes a much bigger CPU time. I found an open source C++ code for inverse calculation, I found, for 3000x3000 sparse lower triangluar matrix, LU decomposition takes 0.05 seconds, but to calculate the inverse matrix by solving L*L^-1=I (to solve 3000 linear equations) takes about 4.0 seconds.

However, I calculated the inverse matrix using inv.m on the lower triangular matrix, it took 0.015 seconds in total. This is my quwstion: why inv.m is so fast? In the document of inv.m, inv.m also solves the 3000 linear equations (use backslash operator).

[Matt J]

I found, for 3000x3000 sparse lower triangluar matrix, LU decomposition takes 0.05 seconds,

That seems absurdly slow. For a matrix that is already lower triangular, the LU decomposition should take almost no effort at all.

but to calculate the inverse matrix by solving L*L^-1=I (to solve 3000 linear equations) takes about 4.0 seconds.

Since your 3rd party code is slow at LU decomposition, it's not a big surprise that it is slow at equation solving. Presumably, it is using a very slow equation solving algorithm that doesn't know L is triangular, and perhaps doesn't even know that L is sparse.

[Benson Gou]

Thanks, Matt.

I know it is very slow for the 3rd party code which uses sparse matrix techniques, but it is not the point of my question. My question is why Matlab code (inv.m) is so fast, especially in the linear sparse equation solving?

[Steven Lord]

My question is why Matlab code (inv.m) is so fast, especially in the linear sparse equation solving?

Because we've put in a lot of effort over the decades at making the linear algebra functionality in MATLAB fast and good. [In the project management triangle I'd argue we've chosen "good" and "fast". The amount of design, implementation, and testing work we've done and continue to do pretty much rules out "cheap".]

Since you wrote:

I actually want to implement inv.m in C++ for a lower triangular sparse matrix.

you may be able to generate code using MATLAB Coder for this (the inv documentation page does list "C/C++ Code Generation" as an extended capability for this function and MATLAB Coder does have a documentation page for sparse matrix code generation support).

[Matt J]

My question is why Matlab code (inv.m) is so fast, especially in the linear sparse equation solving?

Again, because it knows and exploits the fact that L is triangular, not just that it is sparse. You haven't confirmed whether your 3rd party code does that.

you may be able to generate code using MATLAB Coder for this (the inv documentation page does list "C/C++ Code Generation" as an extended capability

And for that matter, linsolve is supported by Matlab Coder as well. This would allow you to stipulate directly that the matrix is lower triangular,

opts.LT=true; %notifies linsolve that L will be lower triangular
inverseL=linsolve(L, speye(size(L)),opts);

Answer 2

[Matt J]

Do you really need the inverse of the matrix or do you need to solve a system of equations whose coefficient matrix is that lower triangular matrix? In the latter case, don't invert the matrix. Use the backslash operator \ instead.

Or, perhaps linsolve? This allows you to inform the solver that the matrix is lower triangular.

opts.LT=true;
x=linsolve(A,b,opts);