發行項
02/02/2017

2016 年 7 月

第 31 卷，第 7 期

本文章是由機器翻譯。

測試回合 - 使用 C# 反轉矩陣

James McCaffrey 機器學習 (ML) 的軟體系統中最基本的技術之一是矩陣反轉。不清楚我的理由，Microsoft.NET Framework 似乎沒有矩陣反轉方法（或如果沒有這種方法，很好隱藏）。在本文中我會呈現，並說明會使用稱為 Crout 的 LU 分解演算法矩陣反轉方法的程式碼。

我是第一個承認該矩陣反轉不是很閃亮的主題。但如果您想要建立 ML 系統不需依賴外部程式庫，讓矩陣反轉方法是必要條件，因為矩陣反轉使用數十個重要的 ML 演算法。

請參閱這篇文章向其中的好方法是查看示範程式中 [圖 1。

[圖 1 矩陣反轉示範

示範一開始設定，並且顯示（4 個資料列，4 個資料行） 4 x 4 矩陣 m:

3.0  7.0  2.0  5.0
1.0  8.0  4.0  2.0
2.0  1.0  9.0  3.0
5.0  4.0  7.0  1.0

然後會計算使用程式定義的方法矩陣的反轉，並顯示結果︰

0.097   -0.183   -0.115    0.224
 -0.019    0.146   -0.068    0.010
 -0.087    0.064    0.103   -0.002
  0.204   -0.120    0.123   -0.147

就差不多是這樣。但是如您所見，矩陣反轉非常棘手。接下來，示範驗證計算的反向乘以反向時間原始的矩陣正確︰

1.0  0.0  0.0  0.0
0.0  1.0  0.0  0.0
0.0  0.0  1.0  0.0
0.0  0.0  0.0  1.0

相乘的結果是單位矩陣（對角線上 1.0 的值，在其他地方 0.0 的值）表示反轉結果正確無誤。

在幕後矩陣反轉方法會使用稱為矩陣分解的技術。分解成兩個矩陣，名為 L （下方），而 U （上方），因素矩陣，當一起相乘提供原始的矩陣中，但某些資料列重新排列。這個示範顯示分解︰

5.000    4.000    7.000    1.000
0.200    7.200    2.600    1.800
0.400   -0.083    6.417    2.750
0.600    0.639   -0.602    4.905

分解實際上包含 L 和 U 矩陣。L 矩陣結合 LU 矩陣中，對角線上 dummy 1.0 的值和 0.0 的值，上半部的左下方部分中的值所組成︰

1.000    0.000    0.000    0.000
0.200    1.000    0.000    0.000
0.400   -0.083    1.000    0.000
0.600    0.639   -0.602    1.000

U 矩陣所組成的右上方部分與對角線的結合 LU 矩陣中的值︰

5.000    4.000    7.000    1.000
0.000    7.200    2.600    1.800
0.000    0.000    6.417    2.750
0.000    0.000    0.000    4.905

示範驗證 LU 分解 L，而 U 矩陣相乘，並顯示結果正確︰

5.0  4.0  7.0  1.0
1.0  8.0  4.0  2.0
2.0  1.0  9.0  3.0
3.0  7.0  2.0  5.0

如果比較原始矩陣 m U L *，您會看到 L * U 是幾乎完全相同，m，但 L 的資料列 * 您必須置換（重新整理）。資料列排列資訊為︰

3  1  2  0

這表示資料列 [0] 的 m 是在資料列 [3] 的左 * U;m [1] 列位於資料列 [1] 的左 * U;m [2] 列位於資料列 [2] 的左 * U;m [3] 列位於資料列 [0] 的左 * U。

了解矩陣反轉

在標準的算術、數字 z 反向是乘以 z 提供 1 的數字。例如，如果 z = 3，z 的逆元是 1/3 = 0.33，因為 3 * (1/3) = 1。

矩陣反轉延伸了這個概念。反向 nxn （因為資料列數目等於資料行的數目，稱為「正方形矩陣」）的矩陣 m 是矩陣 mi 滿足下列條件，m * mi = I 我用來為識別矩陣 (斜線，在其他地方 0.0s年 1.0)。

並非所有的矩陣有反轉。其實，沒有純量值稱為矩陣的行列式。如果矩陣的行列式為零，則矩陣 doen't 有反轉。

請注意，若要完全瞭解矩陣反轉，您必須了解矩陣相乘。矩陣相乘最佳範例所說明。看一下範例 [圖 2。結果矩陣的 [r] [c] 儲存格的值是資料列的第一個矩陣的 r 的值和資料行 c 的第二個矩陣中的值的乘積。

[圖 2 矩陣乘法

尋找矩陣的反轉，當您僅使用正方形矩陣，但矩陣相乘可以套用至不同的矩陣。在這些情況下矩陣必須所謂 conformable。如果矩陣的圖形 axn 而且矩陣 B 圖形 nxb，相乘的結果會有圖形 axb。第一個矩陣中的資料行數目必須等於第二個矩陣中的資料列數目。

示範程式會實作與方法 MatrixProduct 和協助程式方法 MatrixCreate，矩陣相乘中所示 [圖 3。示範使用暴力的方法，但矩陣相乘的結果矩陣中的每個資料格計算無關，因為無法執行使用平行 Parallel.For 方法從.NET 工作平行程式庫。

[圖 3 矩陣乘法

static double[][] MatrixCreate(int rows, int cols)
{
  double[][] result = new double[rows][];
  for (int i = 0; i < rows; ++i)
    result[i] = new double[cols];
  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");
  double[][] result = MatrixCreate(aRows, bCols);
  for (int i = 0; i < 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;
}

示範程式

我編寫使用 C# 的範例程式，但是如果您希望，您應該有移植至另一種語言，例如 Visual Basic 或 Python 程式碼不困難。示範程式碼太長而無法完整呈現，而本文所附的下載中取得完整的程式碼。程式碼也會提供在 quaetrix.com/Matrix/code.html （URL 會區分大小寫）。

若要建立示範程式，我啟動 Visual Studio，並建立新 C# 主控台應用程式，名為 MatrixInverse。示範程式有沒有顯著的.NET Framework 相依性，因此任何版本的 Visual Studio 將運作。範本程式碼載入後，請在 [方案總管] 中我以滑鼠右鍵按一下檔案 Program.cs，並重新將它命名為更具描述性的 MatrixInverseProgram.cs 和 Visual Studio 視窗再自動重新命名類別程式我。

在 [編輯器] 視窗頂端我刪除所有 using 陳述式參考不必要的命名空間，但只是一個參考至最上層系統命名空間會保留。

Main 方法一開始會設定要反轉的矩陣︰

Console.WriteLine("Begin matrix inverse demo");
double[][] m = MatrixCreate(4, 4);
m[0][0] = 3; m[0][1] = 7; m[0][2] = 2; m[0][3] = 5;
m[1][0] = 1; m[1][1] = 8; m[1][2] = 4; m[1][3] = 2;
m[2][0] = 2; m[2][1] = 1; m[2][2] = 9; m[2][3] = 3;
m[3][0] = 5; m[3][1] = 4; m[3][2] = 7; m[3][3] = 1;

在許多情況下，來源資料會儲存在文字檔中因此您必須撰寫的 helper 方法，從檔案載入矩陣。範例將使用陣列的陣列樣式矩陣。不同於大部分的程式設計語言，C# 支援，則為 true 的 n 維矩陣型別，但比較偏好使用標準陣列的陣列方法。

接下來，主要顯示矩陣 m，然後計算並顯示反向︰

Console.WriteLine("Original matrix m is ");
Console.WriteLine(MatrixAsString(m));
double[][] inv = MatrixInverse(m);
Console.WriteLine("Inverse matrix inv is ");
Console.WriteLine(MatrixAsString(inv));

所有的工作是由 MatrixInverse 方法執行。Helper 方法 MatrixAsString 傳回的字串表示的矩陣︰

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

以下的小數位數 (3) 和值的寬度 (8) 數目是硬式編碼為求簡化。更常見的方式會將這些值傳遞做為輸入參數。接下來，示範將原始的矩陣和反轉的矩陣相乘才可以確認結果是單位矩陣︰

double[][] prod = MatrixProduct(m, inv);
Console.WriteLine("The product of m * inv is ");
Console.WriteLine(MatrixAsString(prod));

在此版本中，您必須以視覺方式驗證結果是單位矩陣。更複雜的方法是撰寫方法，以接受矩陣，並傳回 true，如果矩陣是單位矩陣，有一些小差異 (1.0 e-5 是典型) 中儲存格的值。

接下來，示範拆解原始矩陣會說明一些工作︰

double[][] lum;
int[] perm;
int toggle = MatrixDecompose(m, out lum, out perm);
Console.WriteLine("The decomposition is");
Console.WriteLine(MatrixAsString(lum));

呼叫的方法簽章 MatrixDecompose 可能會出現您有點不尋常的。明確傳回值為 + 1 或-1，取決於發生的資料列排列的數字 (偶數或奇數，分別)。切換傳回值不會使用所示範的但如果您想要計算行列式的矩陣中，可以讓您知道是否矩陣的反轉存在，我會解釋，則需要。

Out 參數亮度是結合的 LU （上限為較低）分解程序。Out 參數的權限是資料列被置換的方式進行編碼的整數值的陣列。

接下來，示範從結合 LU 矩陣擷取下限和上限矩陣，並加以顯示︰

double[][] lower = ExtractLower(lum);
double[][] upper = ExtractUpper(lum);
Console.WriteLine("The lower part of LUM is");
Console.WriteLine(MatrixAsString(lower));
Console.WriteLine("The upper part of LUM is");
Console.WriteLine(MatrixAsString(upper));

ExtractLower 和 ExtractUpper 的 helper 方法其實不需要執行矩陣反轉，但是用來說明矩陣分解的運作方式。

示範程式結束時，會顯示資料列排列資訊，以相乘的下限和上限分解矩陣，並顯示結果︰

Console.WriteLine("The perm[] array is");
ShowVector(perm);
double[][] lowTimesUp = MatrixProduct(lower, upper);
Console.WriteLine("The product of lower * upper is ");
Console.WriteLine(MatrixAsString(lowTimesUp));
Console.WriteLine("End matrix inverse demo");

程式定義的 helper 方法 ShowVector 是只將全新的 Main 方法的便利性。如前所述，較低的結果 * upper 是原始的矩陣中，不同之處在於結果的資料列會置換根據權限陣列中的資訊。

MatrixInverse 方法

MatrixInverse 方法定義的開頭︰

static double[][] MatrixInverse(double[][] matrix)
{
  int n = matrix.Length;
  double[][] result = MatrixCreate(n, n);
  for (int i = 0; i < n; ++i)
    for (int j = 0; j < n; ++j)
      result[i][j] = matrix[i][j];
...

方法會假設其輸入的參數，事實上，有一個矩陣。這表示您應該先呼叫，順著檢查︰

int d = MatrixDeterminant(m);
if (d == 0.0)
  Console.WriteLine("No inverse");
else
  double[][] mi = MatrixInverse(m);

替代方法是將方法 MatrixInverse，會建立一份輸入矩陣內的這個錯誤檢查程式碼。您也可以執行矩陣反轉位置，節省記憶體，但是會終結原始的矩陣。

接下來，MatrixInverse 打散輸入矩陣的副本︰

double[][] lum; // Combined lower & upper
int[] perm;
int toggle;
toggle = MatrixDecompose(matrix, out lum, out perm);

它可能會奇怪拆解的矩陣，以計算其反向，但是我信任的所有問題，請移，這種方法比直接反轉的矩陣來得容易。

接下來，示範計算反向使用另一個名為協助程式的協助程式方法︰

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 = Helper(lum, b); //
  for (int j = 0; j < n; ++j)
    result[j][i] = x[j];
}

此程式碼是很明顯。幸運的是，您將必須修改這部分的程式碼。

MatrixInverse 方法結束時，會傳回反向︰

...
  return result;
}

它應該很清楚，方法 MatrixInverse 是本質上的包裝函式方法 MatrixDecompose 和協助程式，這是執行大部分的工作項目。這兩種方法的程式碼都在本文所附的下載檔案。

矩陣的行列式

如果矩陣的行列式為零，矩陣就沒有反轉。假設 3x3 矩陣是︰

1.0  4.0  0.0
3.0  2.0  5.0
7.0  8.0  6.0

矩陣的行列式是︰

+1.0 * [(2.0)(6.0) - (5.0)(8.0)]
-4.0 * [(3.0)(6.0) - (5.0)(7.0)]
+0.0 * [(3.0)(8.0) - (2.0)(7.0)]
= +1.0 * (-28.0) -4.0 * (-17.0) = -28.0 + 68.0 = 40.0

每一個正方形矩陣具有行列式。對於與圖形大於 3x3 矩陣計算行列式是非常困難。不過，其實，如果您將矩陣分解您可以使用組合的較低上限結果來計算行列式輕鬆乘以對角線的結果項目。示範程式定義的方法 MatrixDeterminant 為︰

static double MatrixDeterminant(double[][] matrix)
{
  double[][] lum;
  int[] perm;
  int toggle = MatrixDecompose(matrix, out lum,
    out perm);
  double result = toggle;
  for (int i = 0; i < lum.Length; ++i)
    result *= lum[i][i];
  return result;
}

總結

有效率的矩陣反轉的關鍵在於矩陣分解。有數種演算法分解矩陣。示範程式碼會使用稱為 Crout 的演算法的技術。一般替代方式是 Doolittle 的演算法。我使用偏好 Doolittle 的演算法，因為它是比較簡單的非 Crout 的但我現在喜好 Crout 的演算法因為它具有較少失敗的地方。

我永遠不知道為什麼沒有.NET Framework 計算矩陣的反轉的方法。為了安全起見，反轉的矩陣是非常困難。我在測試本文章中的隨機產生圖形 2 x 2 到 50 x 50 100 萬正方形矩陣、計算反向，以程式設計的方式反向產品的呈現的程式碼與原始的矩陣是單位矩陣，適當的大小。

Dr。James McCaffrey 適用於在美國華盛頓州 Redmond 的 Microsoft Research 他曾在包括 Internet Explorer 和 Bing 的數個 Microsoft 產品。Dr.McCaffrey 可以到達 jammc@microsoft.com。

感謝下列 Microsoft 技術專家檢閱這份文件︰ Chris Lee 和 Kirk Olynyk

測試回合 - 使用 C# 反轉矩陣

了解矩陣反轉

示範程式

MatrixInverse 方法

矩陣的行列式

總結

其他資源