## Unscrambling a Lower-Upper Matrix Decomposition

In my last blog entry I gave a C# implementation of a lower-upper matrix decomposition. For example:

double[][] m = (initialize a 4×4 matrix and assign some values)
int[] indx;
double d = 0;
double[][] lum = MatrixDecomposition(m, indx, out d);

But the result of the method is a single matrix which stores both the lower and upper decomposition of a matrix in row-permuted order, where the 1.0s on the diagonal of the lower part are implied, and where the information necessary to put the rows into the correct order is contained in array indx, and where the out d parameter can be used for other purposes. Here’s how to unscramble the rows:

double[][] lower = ExtractLower(lum);
Console.WriteLine(“Just the lower part:”);
Console.WriteLine(MatrixAsString(lower));

double[][] upper = ExtractUpper(lum);
Console.WriteLine(“Just the upper part:”);
Console.WriteLine(MatrixAsString(upper));

double[][] prod = MatrixProduct(lower, upper);
Console.WriteLine(“lower * upper:”);
Console.WriteLine(MatrixAsString(prod));

double[][] unprod = UnScramble(prod, indx);
Console.WriteLine(“unscrambled product:”);
Console.WriteLine(MatrixAsString(unprod));

See the screenshot below. Notice that the result is the same as the original matrix as it should be. The idea is to extract the lower and upper parts of the result of MatrixDecomposition (adding 1.0s on the diagomnal of the lower matrix), then multiply, then unscramble using the information in indx.

Here’s the code:

static double[][] ExtractLower(double[][] matrix)
{
int rows = matrix.Length;
int cols = matrix[0].Length;
double[][] result = MatrixCreate(rows, cols);
for (int i = 0; i < rows; ++i)
{
for (int j = 0; j < cols; ++j)
{
if (i == j)
result[i][j] = 1.0;
else if (i > j)
result[i][j] = matrix[i][j];
}
}
return result;
}

static double[][] ExtractUpper(double[][] matrix)
{
int rows = matrix.Length; int cols = matrix[0].Length;
double[][] result = MatrixCreate(rows, cols);
for (int i = 0; i < rows; ++i)
{
for (int j = 0; j < cols; ++j)
{
if (i <= j)
result[i][j] = matrix[i][j];
}
}
return result;
}

static double[][] UnScramble(double[][] luMatrix, int[] indx)
{
double[][] result = MatrixDuplicate(luMatrix);
for (int k = indx.Length-1; k >= 0; –k) // reverse order!
{
int row1 = k;
int row2 = indx[k];
int cols = result[0].Length;
for (int j = 0; j < cols; ++j)
{
double tmp = result[row1][j];
result[row1][j] = result[row2][j];
result[row2][j] = tmp;
}
}
return result;
}

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) // could use k < bRows
result[i][j] += matrixA[i][k] * matrixB[k][j];

return result;
} // MatrixProduct