LU Factorization

Under certain conditions the system matrix A of the equation Ax = b can be expressed in the form of a product of a unit lower triangular matrix L with units on the main diagonal and an upper triangular matrix U,and as the result, one has to solve two systems with triangular matrices.

LU Factorization computes matrices L and U such that

	A = L * U

where A is a real, double-precision square matrix, L is a matrix with the same type and dimension as A, and U is an upper triangular matrix with the same type and dimension as A.

LU Factorization

A procedure for decomposing an N X N matrix A into a product of a lower triangular matrix L (has elements only on the diagonal and below) and an upper triangular matrix U (has elements only on the diagonal and above),

L U = A ................... ( I )

Written explicitly for a 4 X 4 matrix A, equation 3.2.4.1 looks like,

................... ( II )

This is same as,

We can use a decomposition such as I to solve the linear set,

...................... ( III )

by first solving for vector y such that,

..................... ( IV )

and then solving,

................. ( V )

Now, equation IV can be solved by forward substitution as follows,

..................... ( VI )

while (3.2.4.5) can be solved by back substitution as follows,

.................... ( VII )

Given a matrix A, we can solve for L and U by first writing the i, j^th component of equation 3.2.4.1 or 3.2.4.2. The component always is a sum beginning with,

The number of terms in the sum depends on whether i or j is the smaller number. There are 3 cases,

.................. ( VIII )

.................. ( IX )

.................. ( X )

Equations ( VIII ) – ( X ) total N² equations for the N² + N unknown ’s and ’s (diagonal represented twice). Since the number of unknowns is greater than the number of equations, N number of unknowns can be specified arbitrarily and used to solve the others. It is always possible to take,

..................... ( XI )

Now, using Crout’s algorithm, all the ’s and ’s can be solved by arranging the equations in a certain order. That order is as follows,

Set (equation XI )
For each do the following 2 procedures:

For each , use equations VIII and IX to solve for

, namely,

....................... ( XII )

For each , use equation X to solve for

, namely,

........................ ( XIII )

By working on a few iterations of the above procedure, it can be seen that ’s and ’s that occur on the right-hand size of equations ( XII ) and ( XIII ) are already determined by the time they are needed. In brief, Crout’s algorithm fills in the combined matrix of ’s and ’s.

............................... ( XIV )

Matrix Determinant using LU Decomposition

Determinant of an LU Decomposed Matrix is the product of its diagonal elements.

det =

Matrix Inverse using LU Decomposition

Using LU Decomposition and back substitution specified above, the inverse of a matrix can be computed column by column in a straightforward manner.

The class os_lufact is a templatized implementation of LU factorization. The following is abrief description of the LU factorization/decomposition technique. The following example code shows how to use various methods supported by the os_lufact class.

Example <ospace/math/examples/lufact.cpp>


  #include <ospace/math.h>


int
main()
  {
  os_math_toolkit math_init;

  double mat_data[ 5 ][ 5 ] =
    {
    { 0, 1, 2, 3, 4 },
    { 5, 6, 7, 8, 9 },
    { 1, 2, 3, 4, 5 },
    { 6, 7, 8, 9, 1 },
    { 2, 3, 4, 5, 6 }
    };

  double vec_data[ 5 ] = { 1, 2, 1, 3, 4 };

  os_num_matrix< double > matrix( 5, 5 );
  os_num_vector< double > vec( 5 );

  for ( int i = 0; i < 5; i++ )
    {
    for ( int j = 0; j < 5; j++ )
      {
      matrix( i, j ) = mat_data[ i ][ j ];
      }
    vec[ i ] = vec_data[ i ];
    }

  os_lufact< double > lufact( matrix );

  cout << "factorized = \n" << lufact.factorized_matrix() << endl;
  cout << "determinant = " << lufact.determinant() << "\n" << endl;
  cout << "inverse = \n" << lufact.inverse() << endl;
  cout << "solve = \n" << lufact.solve( vec ) << endl;


  return 0;
  }