# Inverse of Lower Triangular Matrix

1. May 4, 2012

### Jes1x1

Dear folks!
I want to calculate inverse of lower triangular matrix in a fastest way.
I use Mathematica and function LinearSolve. It did well.
But because the matrix is lower triangular matrix. It is very special. So, I want to make a function called NewLinearSolve which operate like LinearSolve but faster because the NewLinearSolve only focuses on the characteristics of Lower Triangular Matrix, that means it is only focus on calculating elements in the lower part of matrix.
I have tried Forward/Back Substitution algorithm but it is very slow comparing to the function LinearSolve.
I don't know why the function LinearSolve is faster than Forward/Back Substitution.
Thanks very much
Sincerely yours

2. May 4, 2012

### robertsj

Under the hood, Mathematica probably already knows the matrix is lower diagonal and uses some specialized routine to perform the operation. My guess is it uses something like a highly tuned BLAS/LAPACK implementation.

As for implementing your own solver, in almost all cases, it will be slower than the built in functions. This has to do with performance details of linear algebra implementations (caches, etc.) as well as overhead introduced by programming in the higher level language (in this case Mathematica).

Basic take away: rarely expect to beat built in functions except for very, very special cases.

3. May 4, 2012

### Jes1x1

In fact, my matrix quite special. It is a Lower Triangular Matrix which has its first 2 columns is different. Others elements in the remain columns (columns 3 to n) have the same elements with the elements in second columns.
[a 0 0 0 0]
[b c 0 0 0]
[d e c 0 0]
[f g e c 0]
[h k g e c]
My algorithm is tried to calculate only elements in the first 2 columns. Then I copy all elements in the second column to others remain columns.
I did it successfully. But the problem is my algorithm is still slower than LinearSolve.
I have tested and checked my algorithm. It is succeed. But run time is too high compare to LinearSolve.
I have another algorithm is I use LinearSolve only for the first 2 columns. Then I copy elements of the second columns to the remain columns. It is also succeed. And, fortunately, its run time is faster than LinearSolve when n is large (n is capacity of matrix).
But I do not want to use LinearSolve. I want to make a algorithm which named NewLinearSolve to solve Lower Triangular Matrix directly.
Thanks so much
Sincerely yours

4. May 4, 2012

### chiro

Hey Jes1x1 and welcome to the forums.

For this kind of problem, you will use the techniques employed in typical LU decomposition problems.

If you apply Cramers method for matrix inversion, you will notice that because we have a triangular matrix we can calculate the determinant very easily (multiply all diagonal elements and that is your determinant). You can also calculate the cofactors very quickly and optimize the routine for a fixed nxn matrix so that you can calculate very very quickly.

So if you read the cofactor method, that will give you the hints required to do what you need to do and when you convert that to code, you can test it to make sure you get the right output.

5. May 4, 2012

### Jes1x1

My matrix is very special.
Thanks very much
Sincerely yours

6. May 4, 2012

### robertsj

Jes1x1: I'll repeat what I noted above. Even though you could write down the exact solution by manually doing back substitution, my guess is that it would not be faster due to overhead. From what I see of your matrix, it doesn't look like you can make any significant decrease in floating point operations by fortuitous cancelation, so your operations will still be proportional to n^2. (I could be wrong, though.)

I checked out the Mathematica site, and they do in fact use a BLAS/LAPACK implementation. That suggests you will never get the performance they achieve because your implementation (as specialized as it might be) will never outperform the highly tuned *compiled* library.

7. May 4, 2012

### Jes1x1

I know it is difficult. But I must do it.

Anyway, I am grateful you.

Thank you.

Sincerely yours

8. May 4, 2012

### chiro

The idea is simple and it's based on this method:

http://en.wikipedia.org/wiki/Cramer's_rule.

The thing though is that you will get a lot of optimization because you can a) calculate the determinant with one formula (multiply all diagonal entries and return the result) and b) calculate the cofactors very quickly.

Now you will only have to calculate a limited number of cofactors and for many of these the formulas can be reduced to simple equations that can be executed quite fast.

What you have to do is take your lower triangular matrix of nxn, use formula for determinant (multiply all diagonal elements and return result) and then for a specific n for your nxn matrix, on paper figure out the formulas for the cofactors that are non-zero and plug these into your code that calculates the matrix. Then use the fact that your inverse is 1/Det(A) x CofactorMatrix(A) (as given by Cramers rule) in the most optimized way and return your matrix.

The number of operations should be very small and the determinant can be calculated in one formula with n multiplications (1 calculation), while the cofactor calculatiion should be something like (n+1)n/2 calculations for the non-zero cofactor terms. This is a lot less than n3. If it is for a fixed n, then you can create really fast code.

I gaurantee that if you write this code in a compiled library, especially for large n with good code then it should run very quick.

9. May 4, 2012

### Jes1x1

Thank you very much. I really thank you.

I will try coding. Hopefully, I will get the results. I only care about the Run time compare with LinearSolve.

You are very nice and enthusiastic. I own you. Thank you

I hope I can solve this problem.

Sincerely yours.

10. May 4, 2012

### robertsj

Just to clarify, backsubstitution costs n^2 flops. Full elimination costs O(n^3), but the underlying library should recognize the triangular system.

11. May 4, 2012

### Jes1x1

Dear Chiro.

Please help me. I have programmed but it is not faster, even if it is slower than my algorithm.
My code is:

IdenMatrix=IdentityMatrix[n];
t=Det[MyLowerTriangularMatrix];

For[k=1,k<=n,k++,

A1B=ReplacePart[MyLowerTriangularMatrix,{i_,k} -> IdenMatrix[[i,1]]];

A2B=ReplacePart[MyLowerTriangularMatrix,{i_,k} -> IdenMatrix[[i,2]]];

t1=Det[A1B];

t2=Det[A2B];

InverseMyLowerTriangularMatrix[[k,1]]=t/t1;

InverseMyLowerTriangularMatrix[[k,2]]=t/t2;

];

(*Copy elements from second column to another remains columns*)

For[col=3,col<=n,col++,

InverseMyLowerTriangularMatrix[[col;;n,col]]=InverseMyLowerTriangularMatrix[[2;;n-col+2,2]];

]

Thanks you very much

Sincerely yours

12. May 5, 2012

### Jes1x1

Dear folks!

Thanks so much

Sincerely yours

13. May 5, 2012

### AlephZero

I agree with some other posts that you are probably wasting your time trying to do this.

But you might want to think about first inverting the bottom riight (n-1)x(n-1) matrix, which has the special structure, and then finding the full inverse using that.

The submatrix
[c 0 0 0]
[e c 0 0]
[g e c 0]
[k g e c]
is a special case of a Toeplitz matrix (special because nearly half of it is zero). Toeplitz matrices occur in applications like signal processing where efficient numerical methods are important, so a literature search might find something useful. The references on http://en.wikipedia.org/wiki/Toeplitz_matrix might be a starting point.

You will not succeed by trying to write more efficient code than the professionals who develop programs like mathematica. But you might succeed by finding a more efficent algorithm for this special form of matrix.

14. May 5, 2012

### Jes1x1

I know I can not make a code which is better than professionals who make mathematica's function.

Because my Lower Triangular Matrix is very very special (it has first 2 columns is different.) so I think that I can do something better than LinearSolve.

I will try your hint by Toeplitz matrix.

Thanks so much

Sincerely yours

15. May 5, 2012

### chiro

If you need to work within the mathematica environment, I agree with AlephZero in that you probably won't get a better result than an optimized Mathematica routine.

The only chance I think you would have is if you wrote a routine in C,C++ routine using the Mathematica SDK and then created a compiled DLL which was called by Mathematica. Even then though, Mathematica will have to convert everything to the right format for your DLL and convert back the results after the calculation, so with this in mind, I don't know whether it would be worth it do to this even if you could (but it would be your best option).

16. May 7, 2012

### Jes1x1

My code is

A.X=I (A is Lower Triangular Matrix as I mentioned above, X is Inverse of A, and we must find X, I is Identity Matrix corresponding).

(* Find elements of the first and the second columns*)
X[[1,1]]=1/A[[1,1]];
X[[2,2]]=1/A[[2,2]];
X[[2,1]]=-(A[[2,1]]*X[[1,1]]/A[[2,2]])

For[row=3, row<=m, row++,

X[[row,1]]=-(Sum[A[[row, k]]*X[[k, 1]], {k, 1, row-1}]/A[[2,2]]) (*The first column*)

X[[row,2]]=-(Sum[A[[row, k]]*X[[k, 2]], {k, 2, row-1}]/A[[2,2]]) (*The second column*)

]

(* End Find elements of the first and the second column*)

(* Then copy elements of the second column to other columns*)

For[col=3, col<=m, col++,

X[[col;;m,col]]=X[[2;;m-col+2, 2]];

]

It takes time when I calculate the elements of both the first and the second columns

Please help me improve code to calculate the elements of the first and the second columns. Notice that I want to solve this problem within a large number of m, for example, m=300, 400...

I do not know why my code is slow even when I only calculate elements of the first 2 columns (that means, I only solve to 2*m-1 elements) whereas LinearSolve function calculate all elements of the matrix (that means, LinearSolve solve m^2 elements)

Thank you very much

Sincerely yours.

17. May 7, 2012

### chiro

Did you try using Cramers method to find out what minor determinants are zero and how to simplify the ones that are not zero?

Again your determinant is really easy to calculate which is just the multiple of the diagonal values.

The suggestion was made to you (by me previously) to help you come up with an algorithm that would produce your inverse with minimum operations.

You will have to use some paper and pencils to derive the results for the minor determinants for each cell of the matrix, but you will find that a lot of them will be 0 and thus you don't have to calculate anything for those.

Are you aware of Cramers rule and how to use it to calculate the inverse of a matrix (assuming determinant is non-zero) using pen and paper?

18. May 7, 2012

### Jes1x1

Dear Chiro.

Please help me. I have programmed but it is not faster, even if it is slower than my algorithm.
My code is:

IdenMatrix=IdentityMatrix[n];
t=Det[MyLowerTriangularMatrix];

For[k=1,k<=n,k++,

A1B=ReplacePart[MyLowerTriangularMatrix,{i_,k} -> IdenMatrix[[i,1]]];

A2B=ReplacePart[MyLowerTriangularMatrix,{i_,k} -> IdenMatrix[[i,2]]];

t1=Det[A1B];

t2=Det[A2B];

InverseMyLowerTriangularMatrix[[k,1]]=t/t1;

InverseMyLowerTriangularMatrix[[k,2]]=t/t2;

];

(*Copy elements from second column to another remains columns*)

For[col=3,col<=n,col++,

InverseMyLowerTriangularMatrix[[col;;n,col]]=InverseMyLowerTriangularMatrix[[2;;n-col+2,2]];

]

Thanks you very much

Sincerely yours

19. May 7, 2012

### Jes1x1

Dear Chiro,

Thank for your help. But I am sorry, because I do not think that there is any zero determinant.

The inverse of A is always in this form

[v 0 0 0 0]
[k t 0 0 0]
[s u t 0 0]
[m p u t 0]
[w q p u t]

I focus on calculating the elements in the first 2 columns. Then I will copy the elements of the second column to other columns.

Thanks very much

Sincerely yours

20. May 9, 2012

Dear folks!