Features

LUP Matrix Decomposition Using Templates

Mike Rumore

Solving linear equations is not for the amateur, or the faint of heart. Templates, and some good algorithms, can be a real help.

Introduction

Solving a set of linear equations is a very common problem in engineering. Two examples are the sets of equations that arise from analyzing an electrical circuit or from statics. In either case, the set of linear equations has the following form:

A11X1 + A12X2 +…+ A1NXN = Y1
A21X1 + A22X2 +…+ A2NXN = Y2
 .
 .
 .
AN1X1 = AN2X2 +…+ ANNXN = YN

where Aij is the i,j element of an NxN matrix, and Xi and Yi are vectors of N dimensions. In general, you know A and Y, and you have to solve for X.

There are several methods used to solve for the vector X in the previous of set of equations. This article describes an implementation of one of these methods, LUP decomposition. LUP decomposition is the method of choice for matrices that are dense (Aij mostly non-zero) and for N less than a few hundred [1]. For matrices that are sparse (Aij mostly zero) and for N greater than a few hundred, a set of techniques known as iterative methods is used to solve for X [2, 3].

The LUP Decomposition Algorithm

The set of equations listed before can be rewritten more concisely using matrix notation as follows:

A*X = Y,

where A is an NxN matrix, and X and Y are vectors of N dimensions. Again, A and Y are known, and X must be determined.

The LUP decomposition method assumes that the matrix A can be rewritten as the product of two matrices,

A = L*U,

where the matrix L is assumed to be lower triangular (all Lij above the main diagonal are zero) and the matrix U is assumed to be upper triangular (all Uij below the main diagonal are zero). Substituting into the original equation we get:

A = L*U 
A*X = Y 
L*U*X = Y

Now suppose that U*X = Z, where Z is a new vector of N dimensions. Then we can rewrite the original problem A*X = Y as follows:

U*X = Z 
L*Z = Y

What has this gained us? It turns out that the two new sets of equations can be solved very efficiently, since L and U are triangular matrices. First we solve L*Z = Y for Z using forward substitution, then we solve U*X = Z for X, the original problem, using backward substitution [4].

There are two common factorizations for the L and U matrices, Crout's and Doolittle's. Crout's factorization results in a U matrix which has 1s along the main diagonal, and Doolittle's factorization results in an L matrix with ones along the the main diagonal. I have chosen to use Gaussian elimination with pivoting and scaling to generate a Doolittle factorization of the matrix A [4]. Since pivoting is required to maintain stability during Gaussian elimination, it is necessary to track how rows are switched. A permutation vector (the P in the name LUP) is generated by the LUP decomposition function to record the row switches.

Implementation

The LUP code is displayed in Listings 1, 2, 3, 4, 5, and 6. It is based on the algorithms described in reference [4]. All the code was tested using the SUN Solaris 4.2 C++ compiler (without the use of #pragmas). All the code is also available from the CUJ ftp site (see p. 3 for downloading instructions).

Listings 1 and 2 display the headers of the vector and matrix classes used in the LUP code. The LUP code requires the Matrix class to have a constructor of the form Matrix<T>(unsigned int number_of_rows, unsigned int number_of_columns) and an access operator of the form T &operator()(unsigned int row, unsigned int col). Purists may argue that the correct method to access a matrix is with the form, matrix[i][j], but I went with the simpler solution. The vector class must support a constructor of form Vector<T>(unsigned int dimension) and an access operator of the form T &operator[](unsigned int).

Listings 4 and 5 display the header and source file containing the LUP code. Listing 5 lists the following functions.

Listing 5 displays a template function calcEpsilon, which calculates the smallest difference that is representable for a given data type. These values are used to determine if a matrix is singular. If the values are too small, then they can be scaled. On a Sun workstation the epsilon calculated for float and double match the values stored in header float.h, FLT_EPSILON and DBL_EPSILON, respectively.

An Example

Listing 6 displays a sample program which calculates the LUP decomposition for a matrix, then uses the results to calculate the X vector for a given Y vector, and calculates the inverse matrix. The solution vector X and the inverse matrix are printed.

Conclusion

The original version of this code was written in Fortran 77 over ten years ago. Recently I rewrote the code for a numerical methods class that I teach at a local two-year college. Finally, I completely rewrote and updated the algorithm using C++. Since then its main use has been for least-squares fits and optimization.

Acknowledgment

I would like to thank K. B. Williams for his comments and review of the LUP code. His input is gratefully acknowledged.

References

[1] William H. Press, Brian P. Flannery, Saul A. Teukolsky, and William T. Vetterling. Numerical Methods: The Art of Scientific Computing (Cambridge University Press, 1987).

[2] Terrence J. Akai. Applied Numerical Methods for Engineers (John Wiley, 1994).

[3] Lee W. Johnson and R. Dean Riess. Numerical Methods (Addison-Wesley, 1982).

[4] David Kincaid and Ward Cheney. Numerical Methods: Mathematics of Scientific Computing (Brooks/Cole, 1991).

Mike Rumore has a Ph.D. in Nuclear Physics from the University of Colorado, Boulder. He has been working for Lucent Technologies (formerly AT&T) Bell Labs for 11 years as a software engineer. He has been using C++ for over six years for tool generation. His current interests are data structures and numerical algorithms (linear inversion methods applied to radiation spectra).