LU Matrix Decomposition in C++
To download:
1 | bzr branch http://bazaar.enseed.com/lib/Generic/ look under Geometry/LUDecomposition.h |
A LU decomposition decomposes a matrix M into a Lower and Upper matrix pair such that L×U=M. This representation is useful in many ways – for example, it’s easy to solve (L×U)x = b. (for more info, see LU Decomposition on Wikipedia or LU Decomposition on MathWorld.
This particular implementation uses the Matrix class and the same philosophy. Since our matrix class has the size specified as a template parameter, the LU decomposition knows the size of the matrix at compile time. Like most LU decomposition algorithms, it will use a recursive method (For a NxN matrix, first solve for a (N-1)xN matrix…) except that this is done here using recursive templates. The compiler can thus generate the entire unrolled decomposition.
Decomposing a Matrix
Decomposing a matrix in its LU pair is done like this:
1 2 3 4 | Matrix4d m; Matrix4d l, u; LUDecomposition::decompose(m, &l, &u); // l is now a lower matrix, u is an upper matrix, l*u = m |
Solving Mx=b
You once you have a LU pair, you can solve Mx = b = (L×U)x = b:
1 2 3 | Matrix xb; // set this to b initially, you'll get x after LUDecomposition::solve(l, u, &xb); // xb now holds x |
Matrix Inverse
In the previous example, x and b were vectors (or 1×4 matrices), but they could be nx4 matrices. For example, say we have MX =B where X = M-1 and B = I, solving for X would be solving for the inverse of our matrix M:
1 2 3 | Matrix4d inverse = Matrix4d::identityMatrix(); // set this to b initially, you'll get x after LUDecomposition::solve(l, u, &xb); // xb now holds x, which is the inverse of m (or l*u) |
Being More Robust
Sometimes, matrices won’t be well formed and the LU decomposition will generate large numerical errors. One way to prevent this is to solve for LU = PM where P is a pivot matrix (it changes the orders of the rows in M. When the ultimate goal is to use the decomposition to solve a system, we can live with this solution and obtain better precision:
1 2 3 4 5 | Matrix4d m; Matrix4d l, u; Point4i p; // pivot LUDecomposition::decompose(m, &l, &u, &p); // we now have l*u = p*m. |
Now, instead of solving Mx = b, we’ll have to solve for PMx = Pb = LUx. In other words, we have to solve for Pb instead of just b. Say we want to inverse our matrix like in the previous example:
1 2 3 4 5 6 7 8 9 | Matrix4d pivotMatrix = Matrix4d::zeroMatrix(); for (int i = 0; i < 4; ++i) pivotMatrix[i][p[i]] = 1; Matrix4d xb = Matrix4d::itentityMatrix(); // this is b xb = pivotMatirx * xb; // this is Pb LUDecomposition::solve(l, u, &xb); // xb is now x, which is the inverse of our matrix |