# 矩陣解構

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

## 總結

**博士。**詹姆斯 · 麥卡弗裡為伏資訊科學 Inc.，他在管理工作在微軟雷德蒙德，華盛頓州，校園的軟體工程師的技術培訓工作。他曾經參與過多項 Microsoft 產品的研發，包括 Internet Explorer 和 MSN Search。麥卡弗裡是.NET 測試自動化食譜 （Apress，2006年） 的作者。他可以在達成 jammc@microsoft.com