Inverse of Lower Triangular Matrix

[Jes1x1]

Dear folks!
Please help me this problem.
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.
Please help me
Thanks very much
Sincerely yours

 

[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.

 

[Jes1x1]

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.

Thanks for your answer. You are very nice and enthusiastic.
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.
Please help me. Which algorithm can I use, which book can I read?
Thanks so much
Sincerely yours

 

[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.

 

[Jes1x1]

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.

Thanks for your help.

Could you tell me more detail about your ideas.
My matrix is very special.
Thanks very much
Sincerely yours

 

[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.

 

[Jes1x1]

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.

Thanks for your help.

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

Anyway, I am grateful you.

Thank you.

Sincerely yours

 

[chiro]

Thanks for your help.

Could you tell me more detail about your ideas.
My matrix is very special.
Thanks very much
Sincerely yours

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 n³. 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.

 

[Jes1x1]

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

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 n³. 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.

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.

 

[robertsj]

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

 

[Jes1x1]

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

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 n³. 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.

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]];

]

Please help me

Thanks you very much

Sincerely yours

 

[Jes1x1]

Dear folks!

Please help me in this problem. It is very important to me.

Thanks so much

Sincerely yours

 

[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.

 

[Jes1x1]

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 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.

Thank for your help.

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.

If you have any ideas, please help me

Thanks so much

Sincerely yours

 

[chiro]

Thank for your help.

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.

If you have any ideas, please help me

Thanks so much

Sincerely yours

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).

 

[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.

 

[chiro]

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.

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?

 

[Jes1x1]

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?

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]];

]

Please help me

Thanks you very much

Sincerely yours

 

[Jes1x1]

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?

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.

Please help me

Thanks very much

Sincerely yours

 

[Jes1x1]

Dear folks!

Please help me. I need solve this problem.

Thank you very much

Sincerely yours!

 

[chiro]

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.

Please help me

Thanks very much

Sincerely yours

I never said there was going to be a zero determinant. I just said that the determinant can be done in a one line calculation (multiply all the diagonal entries).

I think what you are doing in your code is you are trying to implement a general solution, which is not what I was trying to advocate.

What I was instead advocating is that you get some pen and paper and you expand the minor determinants out for any that are non-zero and then you hard-code these calculations into your routine.

What this will do is that you will set your output matrix to all zero's and then you will modify only the entries that have a non-zero minor determinant (also called a cofactor). If your n is fixed, you can hardcode your routine to initialize your result to zero and only modify the entries that need to be modified thus removing any un-necessary calculations and creating the fastest possible routine.

I'll give you an example: let's say we have a matrix that is 2x2 of the form [a 0; b c] where a <> 0 and c <> 0. Then your routine would look something like this (pseudo-code)

res = zeromatrix(2);
invdet = 1/(a*c);
res(0,0) = invdet*c;
res(1,1) = invdet*a;
res(1,0) = -invdet*b;
return res;

This is the most optimized routine you can have for a lower-triangular matrix for normal non-vector based computation in a normal procedural description. The invdet is used because if you calculate (1/det) then this costs extra computation-time.

Now all I did was use Cramers-Rule and get rid of any un-necessary computations for a lower triangular 2x2 matrix. If you have a fixed n, then this is what you will do. If you have a variable n, then you use the same idea but reduce it down to as low as possible.

Does this make sense?

 

[Jes1x1]

I never said there was going to be a zero determinant. I just said that the determinant can be done in a one line calculation (multiply all the diagonal entries).

I think what you are doing in your code is you are trying to implement a general solution, which is not what I was trying to advocate.

What I was instead advocating is that you get some pen and paper and you expand the minor determinants out for any that are non-zero and then you hard-code these calculations into your routine.

What this will do is that you will set your output matrix to all zero's and then you will modify only the entries that have a non-zero minor determinant (also called a cofactor). If your n is fixed, you can hardcode your routine to initialize your result to zero and only modify the entries that need to be modified thus removing any un-necessary calculations and creating the fastest possible routine.

I'll give you an example: let's say we have a matrix that is 2x2 of the form [a 0; b c] where a <> 0 and c <> 0. Then your routine would look something like this (pseudo-code)

res = zeromatrix(2);
invdet = 1/(a*c);
res(0,0) = invdet*c;
res(1,1) = invdet*a;
res(1,0) = -invdet*b;
return res;

This is the most optimized routine you can have for a lower-triangular matrix for normal non-vector based computation in a normal procedural description. The invdet is used because if you calculate (1/det) then this costs extra computation-time.

Now all I did was use Cramers-Rule and get rid of any un-necessary computations for a lower triangular 2x2 matrix. If you have a fixed n, then this is what you will do. If you have a variable n, then you use the same idea but reduce it down to as low as possible.

Does this make sense?

Thank you, Chiro!

You are very enthusiastic and very nice. I own you a lot.

I will try your hints. Hopefully I can gain the results. Then I will ask you late.

Sincerely yours