Matrix decomposition, a technique that breaks down a square numeric matrix into two different square matrices, is the basis for efficiently solving a system of equations, which in turn is the basis for inverting a matrix. And inverting a matrix is a part of many important algorithms. This article presents and explains C# code that performs matrix decomposition, matrix inversion, a system of equations solution and related operations.

Admittedly, matrix decomposition isn’t a flashy topic, but a collection of matrix methods can be an important addition to your personal code library. The methods are explained so you can modify the source code to meet your own needs. Additionally, some of the techniques used in the matrix methods can be reused in other coding scenarios.

The best way for you to get a feel for the kind of information presented in this article is to take a look at the screenshot in **Figure 1**. The demo program begins by creating a 4 x 4 square matrix and displaying its values. Next, the matrix is decomposed into what’s called an LUP matrix. L stands for lower and U stands for upper. The P part (P stands for permutation) is an array with values {3,1,2,0} and indicates that rows 0 and 3 were exchanged during the decomposition process. The decomposition also generated a toggle value of -1, indicating that an odd number of row exchanges occurred. The demo program displays the decomposition in two ways: first as a combined LU matrix and then as separate L and U matrices. Next, the program computes and displays the inverse of the original matrix, using the LUP matrix behind the scenes. The demo program computes the determinant of the original matrix, again using the decomposition. It then uses the inverse of the matrix to solve a system of linear equations and concludes by combining the L and U matrices back into the original matrix.

**Figure 1 Matrix Decomposition Demo**

But why go to all the trouble of creating a custom matrix decomposition method and a library of related methods? Although there are many standalone matrix tools available, they can sometimes be difficult to integrate into an application or system. And in spite of the fundamental importance of matrix decomposition, there are few free, non-copyrighted .NET code implementations available; those that do exist are often not explained in enough detail for you to modify the source code to suit your coding scenarios.

This article assumes you have intermediate C# programming skills and at least a basic understanding of matrix operations and terminology. All of the key C# code is presented in this article. The code is also available from the MSDN code download site at msdn.microsoft.com/magazine/msdnmag1212.

## Matrix Definition

There are several ways to implement a matrix in C#. The traditional approach, and the one used in this article, is to use an array of arrays, sometimes called a jagged array. For example, this code defines a matrix with three rows and two columns:

```
double[][] m = new double[3][];
m[0] = new double[2];
m[1] = new double[2];
m[2] = new double[2];
m[2][1] = 5.0; // set row 2, col 1
```

Unlike most programming languages, C# has a built-in multidimensional array type, which provides an alternative approach. For example:

```
double[,] m = new double[3,2];
m[2,1] = 5.0;
```

A third approach to implementing matrices in C# is to use a single array combined with array index manipulation, like this:

```
int rows = 3;
int cols = 2;
double[] m = new double[rows * cols];
int i = 2;
int j = 1;
m[i * cols + j] = 5.0;
```

Regardless of the storage scheme used, matrices can be implemented using either an OOP or a static-method approach. For example, an OOP approach could resemble:

```
public class MyMatrix
{
int m; // number rows
int n; // number columns
data[][]; // the values
...
}
```

There’s no single best choice for matrix design; the best design depends on the particular coding scenario you’re operating in and on your personal coding preferences. This article uses a static-method approach because it’s the easiest to understand and refactor.

When using an array-of-arrays design for matrices, because each row must be allocated separately, it’s often convenient to define a helper method to perform memory allocation. For example:

```
static double[][] MatrixCreate(int rows, int cols)
{
// creates a matrix initialized to all 0.0s
// do error checking here?
double[][] result = new double[rows][];
for (int i = 0; i < rows; ++i)
result[i] = new double[cols]; // auto init to 0.0
return result;
}
```

The method can be called like so:

```
double[][] m = MatrixCreate(3,2);
m[2][1] = 5.0;
```

This method demonstrates one advantage of creating your own library of matrix methods: If you want to improve performance, you can omit error-checking the input parameters—at the expense of increasing the risk of the calling code causing an exception. (To keep this article short, most error checking has been removed.) Another advantage is that you can customize your library to optimize for your exact scenario. The main disadvantage of creating your own library is that it can take longer than using an existing library.

Another convenient method to add to your matrix library is one that can be used to display a matrix as a string. Here’s one possibility:

```
static string MatrixAsString(double[][] matrix)
{
string s = "";
for (int i = 0; i < matrix.Length; ++i)
{
for (int j = 0; j < matrix[i].Length; ++j)
s += matrix[i][j].ToString("F3").PadLeft(8) + " ";
s += Environment.NewLine;
}
return s;
}
```

You may want to parameterize the number of decimals to display, or the column-width padding, or both.

## Matrix Multiplication

The Microsoft .NET Framework 4 and later provide a neat way to significantly improve the performance of a matrix multiplication method. Matrix multiplication is illustrated in **Figure 2**.

**Figure 2 Matrix Multiplication**

Note that the computation of each cell value of the result does not depend on any other cell values in the result, so each computation is independent and they can potentially be performed in parallel on a machine with multiple processors. Here’s a standard approach to matrix multiplication:

```
static double[][] MatrixProduct(double[][] matrixA,
double[][] matrixB)
{
int aRows = matrixA.Length; int aCols = matrixA[0].Length;
int bRows = matrixB.Length; int bCols = matrixB[0].Length;
if (aCols != bRows)
throw new Exception("Non-conformable matrices in MatrixProduct");
double[][] result = MatrixCreate(aRows, bCols);
for (int i = 0; i < aRows; ++i) // each row of A
for (int j = 0; j < bCols; ++j) // each col of B
for (int k = 0; k < aCols; ++k)
result[i][j] += matrixA[i][k] * matrixB[k][j];
return result;
}
```

Because matrix multiplication is an O(n^3) operation, performance can be an issue. For example, if matrix A has size 300 x 200 and matrix B has size 200 x 400, computing the product of A and B requires 300 * 200 * 400 = 24,000,000 multiplications. The Task Parallel Library (TPL) in the System.Threading.Tasks namespace in the .NET Framework 4 and later makes it easy to code a simple parallelized version of matrix multiplication. One possibility is:

```
static double[][] MatrixProduct(double[][] matrixA,
double[][] matrixB)
{
// error check, compute aRows, aCols, bCols
double[][] result = MatrixCreate(aRows, bCols);
Parallel.For(0, aRows, i =>
{
for (int j = 0; j < bCols; ++j)
for (int k = 0; k < aCols; ++k)
result[i][j] += matrixA[i][k] * matrixB[k][j];
}
);
return result;
}
```

This version chops up the computations by rows. Behind the scenes, the TPL generates all the complex synchronization plumbing code to perform the computations on multiple processors.

## Consistency Testing

An interesting aspect of a library of methods that are related to each other is that it’s often possible to test them by checking to see if they produce consistent results. For example, suppose you have a method that creates a random matrix:

```
static double[][] MatrixRandom(int rows, int cols,
double minVal, double maxVal, int seed)
{
// return matrix with values between minVal and maxVal
Random ran = new Random(seed);
double[][] result = MatrixCreate(rows, cols);
for (int i = 0; i < rows; ++i)
for (int j = 0; j < cols; ++j)
result[i][j] = (maxVal - minVal) * ran.NextDouble() + minVal;
return result;
}
```

Additionally, suppose you have a matrix that creates the identity matrix—that is, a square matrix with 1.0s on the main diagonal, 0.0s elsewhere:

```
static double[][] MatrixIdentity(int n)
{
double[][] result = MatrixCreate(n, n);
for (int i = 0; i < n; ++i)
result[i][i] = 1.0;
return result;
}
```

And suppose you have a method that compares two matrices for equality:

```
static bool MatrixAreEqual(double[][] matrixA,
double[][] matrixB, double epsilon)
{
// true if all values in A == corresponding values in B
int aRows = matrixA.Length;
int bCols = matrixB[0].Length;
for (int i = 0; i < aRows; ++i) // each row of A and B
for (int j = 0; j < aCols; ++j) // each col of A and B
if (Math.Abs(matrixA[i][j] - matrixB[i][j]) > epsilon)
return false;
return true;
}
```

Notice the MatrixAreEqual method does not compare cell values for exact equality because the values are type double. Instead, the method checks to see if all the cell values are very close (within epsilon) to each other.

Because the product of any square matrix m and the identity matrix of the same dimension is equal to the original matrix m, you can test the matrix product method along the lines of this code:

```
double[][] m = MatrixRandom(4, 4, -9.0, 9.0, 0);
double[][] i = MatrixIdentity(4);
double[][] mi = MatrixProduct(m, i);
if (MatrixAreEqual(m, mi, 0.00000001) == true)
Console.WriteLine("Consistent result");
else
Console.WriteLine("Something is wrong");
Consistency checking lends itself well to random input testing.
```

## Matrix Decomposition

Matrix decomposition takes a square matrix M and computes two new square matrices that when multiplied together give the original matrix M. The idea is similar to ordinary number factoring: the number 6 can be factored into 2 and 3 because 2 * 3 = 6. At first it may seem like there’s little point in decomposing a matrix, but it turns out that matrix decomposition makes the very difficult task of matrix inversion quite a bit simpler. There are many different kinds of matrix decomposition, and each kind can be computed using several different algorithms. The technique presented in this article is called LUP decomposition and uses Doolittle’s method with partial pivoting.

To understand LUP decomposition, it’s helpful to first understand the simpler LU decomposition, which was introduced by the famous mathematician Alan Turing. Suppose you have this 4 x 4 matrix M:

```
9.000 5.000 3.000 4.000
4.000 8.000 2.000 5.000
3.000 5.000 7.000 1.000
2.000 6.000 0.000 8.000
```

One possible LU decomposition of M is L =

```
1.000 0.000 0.000 0.000
0.444 1.000 0.000 0.000
0.333 0.577 1.000 0.000
0.222 0.846 -0.219 1.000
```

And U =

```
9.000 5.000 3.000 4.000
0.000 5.778 0.667 3.222
0.000 0.000 5.615 -2.192
0.000 0.000 0.000 3.904
```

This works because L * U = M. Notice that the lower L matrix has 1.0s on the diagonal and 0.0s in the upper right. In other words, the significant cell values of the lower matrix are in the lower left. Similarly the significant cell values of the upper matrix are on the main diagonal and in the upper right.

Notice as well that there’s no overlap of the locations of the significant cell values in L and U. So, instead of generating two results matrices, L and U, matrix decomposition usually stores both the lower and upper results into a single matrix that holds both L and U to save memory space.

LUP matrix decomposition is a slight but important variation on LU decomposition. LUP decomposition takes a matrix M and produces L and U matrices but also a P array. The product of L and U in LUP is not exactly the original matrix M but instead is a version of M where some of the rows have been rearranged. The way in which the rows have been rearranged is stored into the P array; this information can be used to reconstruct the original matrix M.

A close cousin to the Doolittle decomposition presented in this article is called Crout’s decomposition. The main difference between Doolittle and Crout is that Doolittle places 1.0s on the main diagonal of the L matrix and Crout places 1.0s on the main diagonal of the U matrix.

The reason LUP decomposition is used more often than LU decomposition is subtle. As you’ll see shortly, matrix decomposition is used to compute the inverse of a matrix. However, when matrix decomposition is used as a helper for matrix inversion, it turns out that the inversion will fail if there’s a 0.0 value on the main diagonal of the LU matrix. So in LUP decomposition, when a 0.0 value ends up on the main diagonal, the algorithm exchanges two rows to move the 0.0 value off the diagonal and keeps track of which rows were exchanged in the P array.

**Figure 3** lists a matrix decomposition method.

**Figure 3 A Matrix Decomposition Method**

```
static double[][] MatrixDecompose(double[][] matrix,
out int[] perm, out int toggle)
{
// Doolittle LUP decomposition.
// assumes matrix is square.
int n = matrix.Length; // convenience
double[][] result = MatrixDuplicate(matrix);
perm = new int[n];
for (int i = 0; i < n; ++i) { perm[i] = i; }
toggle = 1;
for (int j = 0; j < n - 1; ++j) // each column
{
double colMax = Math.Abs(result[j][j]); // largest val in col j
int pRow = j;
for (int i = j + 1; i < n; ++i)
{
if (result[i][j] > colMax)
{
colMax = result[i][j];
pRow = i;
}
}
if (pRow != j) // swap rows
{
double[] rowPtr = result[pRow];
result[pRow] = result[j];
result[j] = rowPtr;
int tmp = perm[pRow]; // and swap perm info
perm[pRow] = perm[j];
perm[j] = tmp;
toggle = -toggle; // row-swap toggle
}
if (Math.Abs(result[j][j]) < 1.0E-20)
return null; // consider a throw
for (int i = j + 1; i < n; ++i)
{
result[i][j] /= result[j][j];
for (int k = j + 1; k < n; ++k)
result[i][k] -= result[i][j] * result[j][k];
}
} // main j column loop
return result;
}
```

The method could be called like this:

```
double[][] m = MatrixRandom(4, 4, -9.0, 9.0, 0);
int[] perm;
int toggle;
double luMatrix = MatrixDecompose(m, out perm, out toggle);
```

The MatrixDecompose method accepts as its input a square matrix. The method has three return values. The explicit return is a permuted LU matrix. The method returns two values as out parameters. One is a permutation array that holds information about how the rows were permuted. The second out parameter is a toggle value that’s either +1 or -1 depending on whether the number of row exchanges was even (+1) or odd (-1). The toggle value is not used for matrix inversion but is needed if matrix decomposition is used to compute the determinant of a matrix.

The MatrixDecompose method is fairly tricky, but realistically, there are only a few details you need to understand to modify the code. The version presented here allocates new memory for the LU return matrix using a helper method MatrixDuplicate:

```
static double[][] MatrixDuplicate(double[][] matrix)
{
// assumes matrix is not null.
double[][] result = MatrixCreate(matrix.Length, matrix[0].Length);
for (int i = 0; i < matrix.Length; ++i) // copy the values
for (int j = 0; j < matrix[i].Length; ++j)
result[i][j] = matrix[i][j];
return result;
}
```

An alternative approach is to compute the result into the input matrix in order to save memory. With C# semantics, this would make the matrix parameter a ref parameter because it’s used for both input and output. Using this approach, the method signature would be:

```
static void MatrixDecompose(ref double[][] matrix, out int[] perm,
out int toggle)
```

Or, because the explicit return value has been eliminated, you could use it for the permutation array or the exchanges toggle. For example:

`static int[] MatrixDecompose(ref double[][] matrix, out int toggle)`

You might want to eliminate the toggle parameter to simplify the method signature if you don’t intend to use matrix decomposition to compute a determinant.

Another area of MatrixDecompose you might want to modify is this statement:

```
if (Math.Abs(result[j][j]) < 1.0E-20)
return null;
```

In words, this code means: “If, even after exchanging two rows because there was a 0.0 value on the main diagonal, there’s still a 0.0 there, then return null.” You may want to modify the arbitrary epsilon value for the equality of zero check from 1.0E-20 to some other value based on the precision characteristics of your system. And instead of returning null, you might want to throw an exception; if the method were to be called as part of matrix inversion, the inversion would fail. Finally, if you’re using matrix decomposition for some purpose other than matrix inversion, you may want to eliminate this statement altogether.

## Matrix Inversion

The key to using matrix decomposition to invert a matrix is to write a helper method that solves a system of equations. This key helper method is presented in **Figure 4**.

**Figure 4 The HelperSolve Method**

```
static double[] HelperSolve(double[][] luMatrix,
double[] b)
{
// solve luMatrix * x = b
int n = luMatrix.Length;
double[] x = new double[n];
b.CopyTo(x, 0);
for (int i = 1; i < n; ++i)
{
double sum = x[i];
for (int j = 0; j < i; ++j)
sum -= luMatrix[i][j] * x[j];
x[i] = sum;
}
x[n - 1] /= luMatrix[n - 1][n - 1];
for (int i = n - 2; i >= 0; --i)
{
double sum = x[i];
for (int j = i + 1; j < n; ++j)
sum -= luMatrix[i][j] * x[j];
x[i] = sum / luMatrix[i][i];
}
return x;
}
```

The HelperSolve method finds an array x that when multiplied by an LU matrix gives array b. The method is quite clever, and you can only fully understand it by tracing through a few examples. There are two loops. The first loop uses forward substitution on the lower part of the LU matrix. The second loop uses backward substitution on the upper part of the LU matrix. Some different matrix decomposition implementations call their analogous method something like luBackSub.

Although the code is short, it’s tricky, but there shouldn’t be any scenarios where you’d need to modify the code. Notice that HelperSolve accepts an LU matrix from MatrixDecompose but doesn’t use the perm out-parameter. This implies HelperSolve is in fact a helper method and needs additional wrapper code to solve a system of equations. If you refactor the code in this article to an OOP design, you might want to make HelperSolve a private method.

With the HelperSolve method in place, the matrix inversion method can be implemented, as shown in **Figure 5**.

**Figure 5 The MatrixInverse Method**

```
static double[][] MatrixInverse(double[][] matrix)
{
int n = matrix.Length;
double[][] result = MatrixDuplicate(matrix);
int[] perm;
int toggle;
double[][] lum = MatrixDecompose(matrix, out perm, out toggle);
if (lum == null)
throw new Exception("Unable to compute inverse");
double[] b = new double[n];
for (int i = 0; i < n; ++i)
{
for (int j = 0; j < n; ++j)
{
if (i == perm[j])
b[j] = 1.0;
else
b[j] = 0.0;
}
double[] x = HelperSolve(lum, b);
for (int j = 0; j < n; ++j)
result[j][i] = x[j];
}
return result;
}
```

Again, the code is tricky. The inversion algorithm is based on the idea that the product of a matrix M and its inverse is the identity matrix. Method MatrixInverse essentially solves a system of equations Ax = b where A is an LU matrix decomposition and the b constants are either 1.0 or 0.0 and correspond to the identity matrix. Notice that MatrixInverse uses the perm array from the call to MatrixDecompose.

Calling the MatrixInverse method could look like this:

```
double[][] m = MatrixRandom(4, 4, -9.0, 9.0, 0);
double[][] inverse = MatrixInverse(m);
Console.WriteLine(MatrixAsString(inverse));
```

To summarize, an important matrix operation is matrix inversion, which is quite difficult. One approach is to decompose the original matrix, write a helper solve method that performs forward and backward substitution, and then use the decomposition together with the LU permutation array and the helper solve method to find the inverse. This approach may seem complicated, but it’s usually more efficient and easier than computing a matrix inverse directly.

## Matrix Determinant

With a matrix decomposition method in hand, it’s easy to code a method to compute the determinant of a matrix:

```
static double MatrixDeterminant(double[][] matrix)
{
int[] perm;
int toggle;
double[][] lum = MatrixDecompose(matrix, out perm, out toggle);
if (lum == null)
throw new Exception("Unable to compute MatrixDeterminant");
double result = toggle;
for (int i = 0; i < lum.Length; ++i)
result *= lum[i][i];
return result;
}
```

As it turns out, somewhat surprisingly, the determinant of a matrix is just plus or minus (depending on the toggle value) the product of the values on the main diagonal of the matrix LUP decomposition. Notice that there’s an implicit type conversion of the value of toggle from int to double. In addition to adding error checking to MatrixDeterminant, you might want to add a short-circuit return in situations where the input matrix has size 1 (then return the single value) or size 2 x 2 (then return a*d - b*c). Calling the determinant method could look like this:

```
double[][] m = MatrixRandom(4, 4, -9.0, 9.0, 0);
double det = MatrixDeterminant(m);
Console.WriteLine("Determinant = " + det.ToString("F2"));
```

## Solving Systems of Equations

The HelperSolve method can be easily adapted to solve a system of linear equations:

```
static double[] SystemSolve(double[][] A, double[] b)
{
// Solve Ax = b
int n = A.Length;
int[] perm;
int toggle;
double[][] luMatrix = MatrixDecompose(A,
out perm, out toggle);
if (luMatrix == null)
return null; // or throw
double[] bp = new double[b.Length];
for (int i = 0; i < n; ++i)
bp[i] = b[perm[i]];
double[] x = HelperSolve(luMatrix, bp);
return x;
}
```

Here’s the code that produced the screenshot in **Figure 1** to solve the following system:

```
3x0 + 7x1 + 2x2 + 5x3 = 49
x0 + 8x1 + 4x2 + 2x3 = 30
2x0 + x1 + 9x2 + 3x3 = 43
5x0 + 4x1 + 7x2 + x3 = 52
double[][] m = MatrixCreate(4, 4);
m[0][0] = 3.0; m[0][1] = 7.0; m[0][2] = 2.0; m[0][3] = 5.0;
m[1][0] = 1.0; m[1][1] = 8.0; m[1][2] = 4.0; m[1][3] = 2.0;
m[2][0] = 2.0; m[2][1] = 1.0; m[2][2] = 9.0; m[2][3] = 3.0;
m[3][0] = 5.0; m[3][1] = 4.0; m[3][2] = 7.0; m[3][3] = 1.0;
double[] b = new double[] { 49.0, 30.0, 43.0, 52.0 };
double[] x = SystemSolve(m, b);
Console.WriteLine("\nSolution is x = \n" + VectorAsString(x));
```

Notice that SystemSolve rearranges its b input parameter using the perm array from MatrixDecompose before calling HelperSolve.

## Understanding the Permutation Array

The last few lines of output in the screenshot in **Figure 1** indicate that it’s possible to multiply the L and U matrices in such a way as to get the original matrix. Knowing how to do this will not enable you to solve practical matrix problems, but it will help you understand the P part of LUP decomposition. Regenerating an original matrix from its L and U components can also be useful for testing your matrix library methods for consistency.

One way to regenerate an original matrix after LUP decomposition is to multiply L and U, and then permute the rows of the product using the P array:

```
// create matrix m
// call MatrixDecompose, returning lu and perm
// extract the lower part of lu as matrix 'lower'
// extract the upper part of lu as matrix 'upper'
double[][] lu = MatrixProduct(lower, upper);
double[][] orig = UnPermute(lu, perm);
if (MatrixAreEqual(orig, m, 0.000000001) == true)
Console.WriteLine("L and U unpermuted using perm array");
```

The UnPermute method can be coded like this:

```
static double[][] UnPermute(double[][] luProduct, int[] perm)
{
double[][] result = MatrixDuplicate(luProduct);
int[] unperm = new int[perm.Length];
for (int i = 0; i < perm.Length; ++i)
unperm[perm[i]] = i; // create un-perm array
for (int r = 0; r < luProduct.Length; ++r) // each row
result[r] = luProduct[unperm[r]];
return result;
}
```

A second approach for regenerating an original matrix from its LUP decomposition is to convert the perm array to a perm matrix and then multiply the perm matrix and the combined LU matrix:

```
// create matrix m
// call MatrixDecompose, returning lu and perm
// extract the lower part of lu as matrix 'lower'
// extract the upper part of lu as matrix 'upper'
double[][] permMatrix = PermArrayToMatrix(perm);
double[][] orig = MatrixProduct(permMatrix, lu);
if (MatrixAreEqual(orig, m, 0.000000001) == true)
Console.WriteLine("L and U unpermuted using perm matrix");
```

A perm matrix is a square matrix with one 1.0 value in each row and each column. The method that creates a perm matrix from a perm array can be coded like so:

```
static double[][] PermArrayToMatrix(int[] perm)
{
// Doolittle perm array to corresponding matrix
int n = perm.Length;
double[][] result = MatrixCreate(n, n);
for (int i = 0; i < n; ++i)
result[i][perm[i]] = 1.0;
return result;
}
```

## Wrapping Up

There are many algorithms that require solving a system of linear equations, finding the inverse of a matrix or finding the determinant of a matrix. Using matrix decomposition is an effective technique for performing all these operations. The code presented here can be used in situations where you want no external dependencies in your code base or you need the ability to customize the operations in order to improve performance or modify functionality. Take the red pill!

**Dr. James McCaffrey** *works for Volt Information Sciences Inc., where he manages technical training for software engineers working at the Microsoft Redmond, Wash., campus. He has worked on several Microsoft products including Internet Explorer and MSN Search. McCaffrey is the author of .NET Test Automation Recipes (Apress, 2006). He can be reached at jammc@microsoft.com.*

Thanks to the following technical experts for reviewing this article: Paul Koch and Dan Liebling