Home > C++ Classes, Programming > LU Matrix Decomposition in C++

LU Matrix Decomposition in C++

August 1st, 2009 Leave a comment Go to comments

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
Categories: C++ Classes, Programming Tags:
  1. sum_pattern
    October 1st, 2011 at 06:47 | #1

    Hello,
    can anybody help me out? I’m a beginner in programming. I’m using MS Visual Studio 2008.
    My purpose (example):
    Solve M * x = b,
    where M is a 7*7 matrix with complex elements (except column 7) and b is a 1*7 matrix with real elements.
    How to use LUDecomposition.h?
    1. copy this header to header files
    2. write a source file (“test.cpp)
    3. run
    Right prosedure? Did I forgot anything? How to include header to test.cpp? How to use complex numbers?
    Any hints for me?
    Sorry for the language, I’m german…

  2. sum_pattern
    October 1st, 2011 at 08:32 | #2

    Never mind…
    I found a solution that works for me on myself.
    Cu

  1. No trackbacks yet.
*

Please leave these two fields as-is:

Protected by Invisible Defender. Showed 403 to 159,832 bad guys.