
Advertisement
Sign in to follow this
Followers
0
LU Decomposition and computing the Determinant of an NxN Matrix[SOLVED]
By
CodeCriminal
, in Math and Physics
This topic is 2926 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.
If you intended to correct an error in the post then please contact us.
Recommended Posts
Buckeye 10747
I'll give it a try, pretty much parroting the video, and see if it helps.
Given a 3x3 (square) matrix:
Don't know if that helps.
Even if you don't understand the derivation, you can follow it like a recipe.
L_{21} = M_{21}/M_{11} .. and so forth.
[Edited by  Buckeye on February 19, 2010 12:19:43 AM]
Given a 3x3 (square) matrix:
25 5 1
64 8 1
144 12 1
An LU decomposition will result in 2 new matrices, L and U, such that:
M (original matrix) = L * U.
......
To form the upper triangle, you want U_{21} (upper triangle 2nd row, 1st column), U_{31} and U_{32} to be zeros. That is, some matrix that looks like:
[ 25 5 1 ]
[ 0 ? ? ]
[ 0 0 ? ]
using "?" for elements we have to determine.
.....
The lower triangle will have a diagonal of 1's, the upperright 3 elements 0's, and some as yet unknown elements:
[ 1 0 0 ]
[ ? 1 0 ]
[ ? ? 1 ]
.....
Determine the lower triangle value L_{21} (the value in the second row, first column),and L_{31}, and convert the upper triangle values U_{12} and U_{31} to 0.
.....
Determine L_{21}:
.....
Multiply the first row by 64 ( M_{21} ) and divide by 25 ( M_{11} ). Then subtract the first row from the second.
Remember the value 64/25 = 2.56. This will be L_{21}.
.....
First row x64 /25 becomes [ 64 12.8 2.56 ]
.....
The subtraction  Second row minus first row (modified as above)
[ 64 8 1 ] < 2nd row
[ 64 12.8 2.56 ] < 1st row modified
=[ 0 4.8 1.56 ] < this becomes the upper triangle second row
.....
Determine L_{31}
.....
Multiply the first row by 144 and divide by 25. Then subtract the first row from the third row. Remember the value 144/25 = 5.76. This will be L_{31}.
.....
First row x144 /25 becomes [ 144 28.8 5.76 ]
.....
The subtraction  third row minus first row
[ 144 12 1 ]
[ 144 28.8 5.76 ]
=[ 0 16.8 4.76 ] < this temporarily becomes the upper triangle third row
.....
Using the subtraction values for rows 2 and 3, the matrix becomes:
25 5 1
0 4.8 1.56
0 16.8 4.76
.....
However, U_{32}, 16.8, needs to be converted to 0 to complete the upper triangle.
.....
To do that, use the same algorithm on the second row that was used with the first row, and subtract the second row from the third. As follows
.....
Multiply the second row by 16.8 and divide by 4.8, and subtract the second row from the third.
Second row x16.8 /4.8 becomes [ 0 16.8 5.46 ]
Remember the value 16.8/4.8 = 3.5. This will be L_{32}.
.....
The subtraction  third row minus second row
[ 0 16.8 4.76 ]
[ 0 16.8 5.46 ]
=[ 0 0 0.7 ] < this becomes the upper triangle third row
.....
Now you have the upper triangle, the U in LU.
[ 25 5 1 ]
[ 0 4.8 1.56]
[ 0 0 0.7 ]
.....
Form the Lower Triangle, the L in LU.
.....
The lower triangle has 0's in the upper right and 1's along the diagonal.
The other 3 elements are as calculated above.
[ 1 0 0 ]
[ 2.56 1 0 ]
[ 5.76 3.5 1 ]
.....
So, M = L*U, or
.....
[ 25 5 1 ] [ 1 0 0 ] [ 25 5 1 ]
[ 64 8 1 ] = [ 2.56 1 0 ] x [ 0 4.8 1.56 ]
[ 144 12 1 ] [ 5.76 3.5 1 ] [ 0 0 0.7 ]
Don't know if that helps.
Even if you don't understand the derivation, you can follow it like a recipe.
L_{21} = M_{21}/M_{11} .. and so forth.
[Edited by  Buckeye on February 19, 2010 12:19:43 AM]
0
CodeCriminal 290
Thanks, but after your suggestion to watch the videos on LU decomposition i went back for another try and finally constructed a working algorithm (with a little help... im very happy right now :) ).
Here it is (ignore the 4x4 thing, this code will be placed in my template matrix class which allows for nxm matrices):
Would you mind letting me know what you think? are there any possible optimizations etc?..
I tested it using different sized matrices and multiplying L U to see if i would get the origional matrix back. Im assuming thats how you know when its working if you do get it back.
PS: rating you up
[Edited by  CodeCriminal on February 19, 2010 1:03:03 AM]
Here it is (ignore the 4x4 thing, this code will be placed in my template matrix class which allows for nxm matrices):
Matrix4x4::value_t Matrix4x4::Determinant( ) const
{
if (_rows != _cols)
throw ; // non square matrix
const size_t n = _rows; // or _cols since this is a square matrix
Matrix4x4 U(*this);
Matrix4x4 L(Identity);
// compute lower and upper triangular matrices
for (size_t k = 0; k < n1; ++k)
{
for (size_t i = k+1; i < n; ++i)
{
L(i, k) = (U(i, k) / U(k, k));
for (size_t j = 0; j < n; ++j)
U(i, j) = U(k, j) * L(i, k);
}
}
// since L is a matrix with only 1's along its diagonal,
// det(*this) = 1 * det(U) = det(U)
value_t det = U(0, 0);
for (size_t i = 1; i < n; ++i)
det *= U(i, i);
return det;
}
Would you mind letting me know what you think? are there any possible optimizations etc?..
I tested it using different sized matrices and multiplying L U to see if i would get the origional matrix back. Im assuming thats how you know when its working if you do get it back.
PS: rating you up
[Edited by  CodeCriminal on February 19, 2010 1:03:03 AM]
0
Buckeye 10747
Wow. You went from not understanding the video to coding an algorithm. Congrats.
By inspection, it appears to be correct and efficient.
Maybe I'm just being picky but that's really only part of it. In addition to [original matrix] = [L] * , L and U must meet the definition requirements for their elements  zeros above and below the diagonal (respectively) and L has 1's along diagonal. That, after all, is what really makes it a useful algorithm.
Just being picky [headshake] because, and you obviously know this, there are probably an infinite number of matrixes [C], such that [original matrix]=[C].
By inspection, it appears to be correct and efficient.
Quote:
( A = L*U ) Im assuming thats how you know when its working if you do get it back.
Maybe I'm just being picky but that's really only part of it. In addition to [original matrix] = [L] * , L and U must meet the definition requirements for their elements  zeros above and below the diagonal (respectively) and L has 1's along diagonal. That, after all, is what really makes it a useful algorithm.
Just being picky [headshake] because, and you obviously know this, there are probably an infinite number of matrixes [C], such that [original matrix]=[C].
0
Emergent 982
Quote:
Original post by Buckeye
Maybe I'm just being picky but that's really only part of it. In addition to [original matrix] = [L] * , L and U must meet the definition requirements for their elements  zeros above and below the diagonal (respectively) and L has 1's along diagonal. That, after all, is what really makes it a useful algorithm.
Picky, yes, but also right...
I'll throw in that most LU solvers don't actually compute an L and U of that type such that A = L U; rather they also swap rows of the A matrix as they go, so that actually P A = L U, where P is a permutation matrix. If CodeCriminal would like his LU solver to be robust, he should probably look into this. It doesn't greatly complicate the computation of the determinant, though, since det(L)=1, and det(P)=(+/)1 (depending on whether an odd or even number of row swaps was performed), so now det(A)=(+/)det(U)...
0
Sign in to follow this
Followers
0

Advertisement
Again, all i really need to know is how to seperate matrix A into L and U matrices. For some reason I just cannot seem to grasp the teachings of the video tutorial.
Share this post
Link to post
Share on other sites