LU Decomposition and computing the Determinant of an NxN Matrix[SOLVED]
Author
Hi all, Ive been trying to find information on LU decomposition after stumbling across this post here unfortunatly im not the greatest mathmatician, and struggle to understand what little information I did find.
I understand how to use the upper and lower trianglular matrices to calculate the determinant of matrix A, my only real problem is obtaining the the L U matrices in the first place.
So, i was wondering if any one could possibly shed some light on the subject or point me in the right direction towards learning how to compute LU decomposition
If there is anything else i should know about computing the determinant of an NxN matrix (perhaps something i should be aware, or something that will make it easier) could you please let me know.
Thanks.
[Edited by - CodeCriminal on February 19, 2010 1:08:57 AM]
I found the following videos informative. They're just Part 1 and Part2 of LU Decomposition. The instructor uses numerical examples, rather than subscipted variables. Don't know whether that'll be good or bad for you. He takes an example 3x3 matrix and forms L and U (each 3x3 of course).
">Part 1
">Part 2
Is that a help?
">Part 1
">Part 2
Is that a help?
Author
Hi, thanks for replying, unfortunatly, I've seen that video a few times before and struggled to understand it. I'll take a look again to see if this time is any different.
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.
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.
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.
L21 = M21/M11 .. and so forth.
[Edited by - Buckeye on February 19, 2010 12:19:43 AM]
Given a 3x3 (square) matrix:
25 5 164 8 1144 12 1An 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 U21 (upper triangle 2nd row, 1st column), U31 and U32 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 upper-right 3 elements 0's, and some as yet unknown elements:[ 1 0 0 ][ ? 1 0 ][ ? ? 1 ].....Determine the lower triangle value L21 (the value in the second row, first column),and L31, and convert the upper triangle values U12 and U31 to 0......Determine L21:.....Multiply the first row by 64 ( M21 ) and divide by 25 ( M11 ). Then subtract the first row from the second.Remember the value 64/25 = 2.56. This will be L21......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 L31.....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 L31......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, U32, -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 x-16.8 /-4.8 becomes [ 0 -16.8 -5.46 ]Remember the value -16.8/-4.8 = 3.5. This will be L32......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.
L21 = M21/M11 .. and so forth.
[Edited by - Buckeye on February 19, 2010 12:19:43 AM]
Author
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 < n-1; ++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]
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].
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)...
This topic is closed to new replies.
Advertisement
Popular Topics
Advertisement
