项目
07/15/2016

2016 年 7 月

第 31 卷，第 7 期

此文章由机器翻译。

测试运行 - 使用 C# 实现矩阵反转

James 麦卡弗里一种机器学习 (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 矩阵中，使用虚拟的 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 比较 L * U，您将看到 L * U m，但 L 的行是几乎相同 * 您必须已排列（重新排列）。行的排列信息是︰

3  1  2  0

这意味着 m 的行 [0] 行 [3] 的 L * U;m 的行 [1] 位于 L * U; 的行 [1]行 [2] 是 m 的在 L * U; 一行 [2]一行从 [3] 的 m U L * [0] 行。

了解矩阵反转

在正常算术中，数字 z 的逆是在乘以 z 发出 1 的数字。例如，如果 z = 3，逆 z 为 1/3 = 0.33，因为 3 * (1/3) = 1。

矩阵反转扩展了此概念。反函数 nxn （称为"正方形矩阵"，因为的行数等于列数）的矩阵 m 是矩阵 mi 如 m * mi = 我我所在单位矩阵（慢速的对角线，0.0s 在其他位置上）。

并非所有的矩阵具有逆实例。事实证明，没有名为矩阵的行列式的标量值。如果矩阵的行列式为零，则矩阵会计具有逆实例。

请注意，若要全面了解矩阵反转，您必须了解矩阵乘法。更好地示例说明矩阵乘法。看一看中的示例 图 2。结果矩阵单元 [r] [c] 处的值是行的第一个矩阵的 r 中的值和列 c 的第二个矩阵中的值的乘积。

图 2 矩阵乘法

查找矩阵反转时, 您仅使用方形矩阵但矩阵乘法可以应用于具有不同形状的矩阵。在这些情况下矩阵必须调用 conformable。如果矩阵 A 具有形状 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 然后会自动重命名 Program 类对我来说。

在编辑器窗口顶部我删除了所有 using 语句引用不必要的命名空间，而只是对顶级 System 命名空间的一个引用。

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;

在许多情况下，您的源数据将存储在一个文本文件因此您必须编写帮助器方法来从文件加载一个矩阵。该演示使用数组的数组样式矩阵。不同于大多数编程语言，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 方法执行所有工作。帮助器方法 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 参数 perm 是编码方式行具有已排列的整数值的数组。

接下来，演示从组合 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 实际上不需要执行矩阵反转的方法，但用于说明矩阵分解的工作原理。

演示程序最后显示行的排列信息、下限和上限分解矩阵相乘并显示结果︰

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 方法干净便利。如前所述，较低的结果 * 左上原始矩阵中，只是结果的行排列根据 perm 数组中的信息。

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 没有计算矩阵反转的方法。当然，矩阵反转是非常棘手。我测试本文中介绍的随机生成 1 亿个方形矩阵与 2 x 2 到 50 x 50 之间的形状、计算逆方法，并以编程方式验证反转该产品的代码和原始矩阵是适当大小的单位矩阵。

Dr.James McCaffrey 供职于华盛顿地区雷蒙德市沃 Microsoft Research 他参与过多个 Microsoft 产品的工作，包括 Internet Explorer 和 Bing。Scripto可通过 jammc@microsoft.com 与 McCaffrey 取得联系。

衷心感谢以下 Microsoft 技术专家对本文的审阅： Chris Lee 和 Kirk Olynyk

测试运行 - 使用 C# 实现矩阵反转

了解矩阵反转

演示程序

MatrixInverse 方法

矩阵的行列式

总结

其他资源