# 矩阵分解

James McCaffrey

## 矩阵定义

``````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
``````

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

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

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

``````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;
}
``````

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

``````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;
}
``````

## 矩阵乘法

4 及更高版本的 Microsoft.NET 框架提供了一种巧妙的方法，大大改善性能的矩阵乘法运算方法。 矩阵乘法所示图 2

``````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;
}
``````

``````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;
}
``````

## 一致性测试

``````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;
}
``````

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

``````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;
}
``````

``````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.
``````

## 矩阵分解

``````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
``````

``````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
``````

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
``````

``````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;
}
``````

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

MatrixDecompose 方法接受作为其输入一个方阵。 该方法具有三个返回值。 显式的返回是一个排列后的陆矩阵。 该方法作为 out 参数返回两个值。 一个是排列数组，它保存有关 permuted 行是如何的信息。 Out 参数的第二个是 + 1 或-1 根据行交换数量是否甚至 (+ 1) 或 (-1) 为奇数的切换值。 切换值不用于矩阵求逆，但如果用于计算一个矩阵的行列式矩阵分解需要。

MatrixDecompose 方法是相当棘手的事情，但实际上，有只有几个细节，您需要了解修改的代码。 这里提出的版本为陆返回矩阵使用 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;
}
``````

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

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

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

## 矩阵求逆

``````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;
}
``````

HelperSolve 方法查找数组 x 当乘以陆矩阵赋予数组 b。 方法是很聪明的并可以通过跟踪几个例子只有完全理解它。 有两个循环。 第一个循环使用转发替代关于陆矩阵的下半部分。 第二个循环使用落后的替代上陆矩阵的上半部分。 一些不同的矩阵分解实现调用其类似的方法像 luBackSub 一样的东西。

``````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;
}
``````

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

## 矩阵行列式

``````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;
}
``````

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

## 方程组求解

HelperSolve 方法可以轻松地调整，以解决系统的线性方程组：

``````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;
}
``````

``````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));
``````

SystemSolve 重新排列其 b 输入的参数在调用 HelperSolve 之前使用 MatrixDecompose 中的烫发数组的通知。

## 理解排列数组

``````// 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");
``````

UnPermute 方法可以像这样进行编码：

``````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;
}
``````

``````// 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");
``````

``````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;
}
``````

## 总结

**博士。**博士 供职于 Volt Information Sciences, Inc.，在该公司他负责管理对华盛顿州雷蒙德市沃什湾 Microsoft 总部园区的软件工程师进行的技术培训。他参与过多项 Microsoft 产品的研发工作，其中包括 Internet Explorer 和 MSN Search。麦卡弗里是.NET 测试自动化食谱 （Apress，2006年） 的作者。可通过 jammc@microsoft.com 与他联系。