C#
矩阵分解
矩阵分解,分解成两个不同的方阵,一个方形数字矩阵的一种技术是有效地解决系统的方程,这反过来逆矩阵的基础的基础。逆矩阵是许多重要算法的一部分。这篇文章介绍并说明执行矩阵分解、 矩阵求逆、 一个方程的求解和相关的业务系统的 C# 代码。
诚然,矩阵分解并不是一个华而不实的话题,但矩阵方法的集合,可以将您的个人代码库的重要补充。方法的解释,以便您可以修改的源代码,以满足您自己的需要。此外,一些在矩阵方法中使用的技术可以在其他编码方案中重用。
您去好好感受一下本文中介绍的信息种类是要看一看的截图中的最佳方法图 1。该演示程序开始通过创建一个 4 x 4 的方阵和显示其值。下一步,该矩阵被分解成了所谓的土地矩阵。L 代表低和 U 站为上限。P 部分 (P 代表置换) {3,120} 的值的数组,指示行 0 和 3 在分解过程中就被交换了。分解也生成切换值为-1,指示发生奇数行交流。该演示程序显示分解两种方式:第一次作为一个组合的 LU 矩阵,然后作为单独的 L 和 U 矩阵。接下来,程序计算并显示原始矩阵,使用在幕后的土地矩阵的逆。该演示程序计算行列式的原始矩阵,请再次使用分解。然后使用矩阵的逆矩阵来解决系统的线性方程组,并通过合并回原始矩阵 L 和 U 矩阵的结论。
图 1 矩阵分解演示
但为什么去创建一个自定义矩阵分解方法和相关方法库的所有的麻烦呢?虽然有许多独立的矩阵工具可用,但他们有时可能很难将集成到应用程序或系统。矩阵分解的根本重要性,尽管有几个免费、 无版权.NET 代码实现可用 ; 那些确实存在往往不是足够为您要修改的源代码,以适合您的编码方案的详细解释。
本文假定您有中间的 C# 编程技巧和矩阵操作和术语的至少一个基本的了解。所有关键的 C# 代码是本文中介绍的。从 MSDN 代码下载站点上的代码是也可用 archve.msdn.microsoft.com/mag201212matrix。
矩阵定义
有几种方法在 C# 中执行一个矩阵。传统的做法和使用在本文中,一是使用一个数组的数组,有时被称为一个交错的数组。例如,此代码定义了一个具有三行和两列的矩阵:
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
与大多数编程语言不同,C# 具有内置的多维数组的类型,其中提供了另一种方法。 例如,
double[,] m = new double[3,2];
m[2,1] = 5.0;
在 C# 中执行矩阵的第三种方法是使用单个数组结合数组索引操作,像这样:
int rows = 3;
int cols = 2;
double[] m = new double[rows * cols];
int i = 2;
int j = 1;
m[i * cols + j] = 5.0;
无论使用的存储方案,可以使用面向对象或静态方法的方法实现矩阵。 例如,OOP 方法可能类似于:
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。
图 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;
}
因为矩阵乘法是 O(n^3) 操作,性能可能是一个问题。 例如,如果矩阵 A 具有大小 300 x 200 和矩阵 B 具有大小 200 x 400,计算产品的 A 和 B 要求 300 * 200 * 400 = 24,000,000 乘法运算。 任务并行库 (TPL) System.Threading.Tasks 命名空间在.NET Framework 4 和更高版本中轻松编写一个简单的矩阵乘法的并行化的版本的代码。 一种可能是:
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;
}
此版本由行印章起计算。 在幕后,TPL 生成所有复杂的同步水暖代码来在多个处理器上执行的计算。
一致性测试
图书馆彼此相关的方法的一个有趣的现象是很多时候可能以测试它们通过检查来看是否他们产生一致的结果。 例如,假设您有一个方法,创建一个随机矩阵:
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;
}
此外,假设您有一个矩阵,创建矩阵的身份 — — 就是一个正方形矩阵与慢速的主对角线,其他地方的 0.0s年上:
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;
}
通知的 MatrixAreEqual 方法并不比单元格的值完全相等的因为都是 double 类型的值。 相反,该方法将查看是否所有单元格的值都非常接近 (epsilon) 内彼此。
因为该产品的任何方阵 m 和同一维度的恒等矩阵等于原始矩阵 m,您可以测试矩阵产品方法沿此代码行:
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.
矩阵分解
矩阵分解需要方阵 M 和计算两个新方阵的相乘时给原始矩阵 M。 这个想法是类似于普通号码保理业务:6 数可计入 2 和 3,因为 2 * 3 = 6。 起初它可能看起来像有小点中分解矩阵,但它原来矩阵分解使得矩阵求逆的非常困难的任务很简单些。 有许多不同种类的矩阵分解,和每一种可以计算使用几种不同的算法。 本文中介绍的技术被称为土地分解,并使用具有部分回转的杜利特尔的方法。
要了解土地分解,最好先了解更简单的 LU 分解,由著名数学家阿兰 · 图灵介绍。 假设您有此 4 × 4 矩阵 m:
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
一个可能的 LU 分解,M 的是 L =
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
这样做是因为 L * U = M。 注意较低的 L 矩阵对角线上的慢速和 0.0s年在右上方。 换句话说,在左下角低矩阵的重要单元格的值。 同样上部矩阵的重要单元格的值是的主对角线上和在右上方。
还注意到有没有重叠的 L 和美国重要的单元格的值的位置 所以,而不是生成两个结果矩阵,L 和 U,矩阵分解通常存储到一个单一的矩阵,持有 L 和 U 以节省内存空间上限和下限结果。
土地矩阵分解是轻微但很重要的变化,对 LU 分解。 土地分解矩阵 M 并生成 L 和 U 的矩阵,但 P 数组。 L 和 U 在图中的产品是不完全的原始矩阵 M,但相反是 M 的版本其中的某些行已经被重新安排。 在这行已经被重新安排的方式存储到 P 的数组 ; 此信息可用于重建原始矩阵 M。
本文中介绍的杜利特尔分解的近亲称为 Crout 的分解。 杜利特尔和 Crout 之间的主要区别是杜利特尔 L 矩阵的主对角线上放置慢速和 Crout U 矩阵的主对角线上放置慢速。
土地分解的原因比 LU 分解是微妙的更频繁的使用。 如您很快就会看到,矩阵分解用于计算矩阵的逆矩阵。 然而,当矩阵分解为佣工用于矩阵求逆,原来反演将失败如果陆矩阵的主对角线上有一个 0.0 的值。 所以在土地分解,当告终的主对角线, 值 0.0 算法会交换移动 0.0 值关闭对角线和跟踪的哪些行被交换 P 数组中的两行。
图 3 列出矩阵分解方法。
图 3 矩阵分解法
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;
}
另一种方法是进入输入矩阵计算出结果,以节省内存。 用 C# 的语义,这将使该矩阵参数的 ref 参数因为它用于输入和输出。 使用此方法,将方法签名:
static void MatrixDecompose(ref double[][] matrix, out int[] perm,
out int toggle)
或者,因为显式的返回值已被消除,就可以使用它为排列数组或交易所切换。 例如,
static int[] MatrixDecompose(ref double[][] matrix, out int toggle)
您可能想要消除要简化方法签名,如果您不打算使用矩阵分解来计算行列式的切换参数。
您可能想要修改的 MatrixDecompose 的另一个领域是此语句:
if (Math.Abs(result[j][j]) < 1.0E-20)
return null;
话说,此代码的含义:"如果,即使后交换两行,因为有一个 0.0 值的主对角线上,仍然是有 0.0,然后返回空值。"您可能想要修改的任意 epsilon 值为 1.0 e 从零复选的平等-20 到一些其他的值基于您的系统的精度特征。 而不是返回 null,您可能想要引发异常 ; 和 如果作为矩阵求逆的一部分被调用的方法,反演将会失败。 最后,如果您使用矩阵分解的一些矩阵求逆以外的目的,您可能想要完全消除这一声明。
矩阵求逆
使用矩阵分解反转矩阵的关键是要编写一个 helper 方法,以解决系统的方程。 此关键的帮助器方法提出在图 4。
图 4 HelperSolve 方法
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 一样的东西。
虽然代码很短,但它是棘手的但就不会有任何情况下,您将需要修改的代码。 通知该帮助器解决接受从 MatrixDecompose LU 矩阵,但不使用该烫发的 out 参数。 这意味着 HelperSolve 其实是要解决系统方程的帮助器方法和需要额外包装代码。 如果您重构到面向对象设计的这篇文章中的代码,您可能想要 HelperSolve 的私有方法。
与地方中的 HelperSolve 方法,矩阵反演方法可以实现,如中所示图 5。
图 5 MatrixInverse 方法
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;
}
再次,代码是棘手的。 反演算法基于一个矩阵 M 和它的逆产品是恒等矩阵的想法。 方法 MatrixInverse 基本上解决了方程组 Ax = 的 b A 是陆矩阵分解和 b 常量是 1.0 或 0.0 和对应单位矩阵。 请注意 MatrixInverse 使用的烫发数组从对 MatrixDecompose 的调用。
调用 MatrixInverse 方法可能看起来像这样:
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;
}
事实证明,有点令人惊讶,矩阵的行列式是只加上或减去 (根据切换值) 的土地分解矩阵的主对角线上的值的乘积。 请注意没有隐式类型转换价值的切换从 int 将翻一番。 除了添加错误检查至 MatrixDeterminant,您可能要添加短路的输入的矩阵有大小 (然后返回单个值) 1 或 2 x 2 的大小的情况下返回 (然后返回 * d-b * c)。 调用的行列式的方法,象这样:
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;
}
这里是制作的截图中的代码图 1 来解决以下系统:
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 中的烫发数组的通知。
理解排列数组
输出中的屏幕快照中的最后几行图 1 表示有可能以一种能获取原始矩阵中的 L 和 U 的矩阵相乘。 了解如何执行此操作,不将使您能够解决实际矩阵的问题,但它将帮助您了解土地分解的 P 部分。 再生原始矩阵从其 L 和 U 的组件也可以用于测试您的矩阵库方法的一致性。
土地分解后重新生成原始矩阵的一种方法是将乘以 L 和 U,然后 permute 产品使用 P 数组中的行:
// 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");
烫发矩阵是一个正方形矩阵与一个 1.0 中的每一行和每一列的值。 从烫发数组创建一个烫发矩阵的方法可以进行编码就像这样:
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 与他联系。
衷心感谢以下技术专家对本文的审阅:保罗 · 科赫和 Dan Liebling