deteminant N * N矩阵

Question

我想实现代码在这里找到&算法：

deteminant of matrix 这里： How to calculate matrix determinant? n*n or just 5*5

但我还是坚持了下来。

我的第一个问题实际上是什么统治这个算法使用（因为有明显的数学由某人可以计算行列式的几个规则） - 所以我想，如果算法是正确应用上首先要检查。

我的第二个问题是我做错了（我指的是执行）或什么是错的算法本身，因为它看起来像3×3和正常工作的4x4，但5x5的它给人NaN的。用几个在线矩阵行列式计算器检查结果，除了5x5以外，结果都很好。

这是我的代码：

using System; 

public class Matrix 
{ 
    private int row_matrix; //number of rows for matrix 
    private int column_matrix; //number of colums for matrix 
    private double[,] matrix; //holds values of matrix itself 

    //create r*c matrix and fill it with data passed to this constructor 
    public Matrix(double[,] double_array) 
    { 
     matrix = double_array; 
     row_matrix = matrix.GetLength(0); 
     column_matrix = matrix.GetLength(1); 
     Console.WriteLine("Contructor which sets matrix size {0}*{1} and fill it with initial data executed.", row_matrix, column_matrix); 
    } 

    //returns total number of rows 
    public int countRows() 
    { 
     return row_matrix; 
    } 

    //returns total number of columns 
    public int countColumns() 
    { 
     return column_matrix; 
    } 

    //returns value of an element for a given row and column of matrix 
    public double readElement(int row, int column) 
    { 
     return matrix[row, column]; 
    } 

    //sets value of an element for a given row and column of matrix 
    public void setElement(double value, int row, int column) 
    { 
     matrix[row, column] = value; 
    } 

    public double deterMatrix() 
    { 
     double det = 0; 
     double value = 0; 
     int i, j, k; 

     i = row_matrix; 
     j = column_matrix; 
     int n = i; 

     if (i != j) 
     { 
      Console.WriteLine("determinant can be calculated only for sqaure matrix!"); 
      return det; 
     } 

     for (i = 0; i < n - 1; i++) 
     { 
      for (j = i + 1; j < n; j++) 
      { 
       det = (this.readElement(j, i)/this.readElement(i, i)); 

       //Console.WriteLine("readElement(j, i): " + this.readElement(j, i)); 
       //Console.WriteLine("readElement(i, i): " + this.readElement(i, i)); 
       //Console.WriteLine("det is" + det); 
       for (k = i; k < n; k++) 
       { 
        value = this.readElement(j, k) - det * this.readElement(i, k); 

        //Console.WriteLine("Set value is:" + value); 
        this.setElement(value, j, k); 
       } 
      } 
     } 
     det = 1; 
     for (i = 0; i < n; i++) 
      det = det * this.readElement(i, i); 

     return det; 
    } 
} 

internal class Program 
{ 
    private static void Main(string[] args) 
    { 
     Matrix mat03 = new Matrix(new[,] 
     { 
      {1.0, 2.0, -1.0}, 
      {-2.0, -5.0, -1.0}, 
      {1.0, -1.0, -2.0}, 
     }); 

     Matrix mat04 = new Matrix(new[,] 
     { 
      {1.0, 2.0, 1.0, 3.0}, 
      {-2.0, -5.0, -2.0, 1.0}, 
      {1.0, -1.0, -3.0, 2.0}, 
      {4.0, -1.0, -3.0, 1.0}, 
     }); 

     Matrix mat05 = new Matrix(new[,] 
     { 
      {1.0, 2.0, 1.0, 2.0, 3.0}, 
      {2.0, 1.0, 2.0, 2.0, 1.0}, 
      {3.0, 1.0, 3.0, 1.0, 2.0}, 
      {1.0, 2.0, 4.0, 3.0, 2.0}, 
      {2.0, 2.0, 1.0, 2.0, 1.0}, 
     }); 

     double determinant = mat03.deterMatrix(); 
     Console.WriteLine("determinant is: {0}", determinant); 

     determinant = mat04.deterMatrix(); 
     Console.WriteLine("determinant is: {0}", determinant); 

     determinant = mat05.deterMatrix(); 
     Console.WriteLine("determinant is: {0}", determinant); 
    } 
} 

Answer 1

为什么要重新发明轮子？一种众所周知的用于获取确定矩阵反转矩阵的方法是进行LU分解。事实上，MSDN杂志最近在C#这里http://msdn.microsoft.com/en-us/magazine/jj863137.aspx做了一个关于如何做到这一点的帖子。

示例代码是

矩阵行列式

在手的矩阵分解方法，可以很容易地进行编码的方法来计算矩阵的行列式：

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

行列式根据分解矩阵上对角线的乘积求出与一个标志检查。阅读文章了解更多详情。

请注意，它们对矩阵使用锯齿阵列，但是您可以将自己的矩阵存储替换为lum[i][j]到lum[i,j]。

Answer 2

@ ja72 非常感谢您的指导。计算任何n * n行列式的最终解决方案如下所示：

using System; 

internal class MatrixDecompositionProgram 
{ 
    private static void Main(string[] args) 
    { 
     float[,] m = MatrixCreate(4, 4); 
     m[0, 0] = 3.0f; m[0, 1] = 7.0f; m[0, 2] = 2.0f; m[0, 3] = 5.0f; 
     m[1, 0] = 1.0f; m[1, 1] = 8.0f; m[1, 2] = 4.0f; m[1, 3] = 2.0f; 
     m[2, 0] = 2.0f; m[2, 1] = 1.0f; m[2, 2] = 9.0f; m[2, 3] = 3.0f; 
     m[3, 0] = 5.0f; m[3, 1] = 4.0f; m[3, 2] = 7.0f; m[3, 3] = 1.0f; 

     int[] perm; 
     int toggle; 

     float[,] luMatrix = MatrixDecompose(m, out perm, out toggle); 

     float[,] lower = ExtractLower(luMatrix); 
     float[,] upper = ExtractUpper(luMatrix); 

     float det = MatrixDeterminant(m); 

     Console.WriteLine("Determinant of m computed via decomposition = " + det.ToString("F1")); 
    } 

    // -------------------------------------------------------------------------------------------------------------- 
    private static float[,] MatrixCreate(int rows, int cols) 
    { 
     // allocates/creates a matrix initialized to all 0.0. assume rows and cols > 0 
     // do error checking here 
     float[,] result = new float[rows, cols]; 
     return result; 
    } 

    // -------------------------------------------------------------------------------------------------------------- 
    private static float[,] MatrixDecompose(float[,] matrix, out int[] perm, out int toggle) 
    { 
     // Doolittle LUP decomposition with partial pivoting. 
     // rerturns: result is L (with 1s on diagonal) and U; perm holds row permutations; toggle is +1 or -1 (even or odd) 
     int rows = matrix.GetLength(0); 
     int cols = matrix.GetLength(1); 

     //Check if matrix is square 
     if (rows != cols) 
      throw new Exception("Attempt to MatrixDecompose a non-square mattrix"); 

     float[,] result = MatrixDuplicate(matrix); // make a copy of the input matrix 

     perm = new int[rows]; // set up row permutation result 
     for (int i = 0; i < rows; ++i) { perm[i] = i; } // i are rows counter 

     toggle = 1; // toggle tracks row swaps. +1 -> even, -1 -> odd. used by MatrixDeterminant 

     for (int j = 0; j < rows - 1; ++j) // each column, j is counter for coulmns 
     { 
      float colMax = Math.Abs(result[j, j]); // find largest value in col j 
      int pRow = j; 
      for (int i = j + 1; i < rows; ++i) 
      { 
       if (result[i, j] > colMax) 
       { 
        colMax = result[i, j]; 
        pRow = i; 
       } 
      } 

      if (pRow != j) // if largest value not on pivot, swap rows 
      { 
       float[] rowPtr = new float[result.GetLength(1)]; 

       //in order to preserve value of j new variable k for counter is declared 
       //rowPtr[] is a 1D array that contains all the elements on a single row of the matrix 
       //there has to be a loop over the columns to transfer the values 
       //from the 2D array to the 1D rowPtr array. 
       //----tranfer 2D array to 1D array BEGIN 

       for (int k = 0; k < result.GetLength(1); k++) 
       { 
        rowPtr[k] = result[pRow, k]; 
       } 

       for (int k = 0; k < result.GetLength(1); k++) 
       { 
        result[pRow, k] = result[j, k]; 
       } 

       for (int k = 0; k < result.GetLength(1); k++) 
       { 
        result[j, k] = rowPtr[k]; 
       } 

       //----tranfer 2D array to 1D array END 

       int tmp = perm[pRow]; // and swap perm info 
       perm[pRow] = perm[j]; 
       perm[j] = tmp; 

       toggle = -toggle; // adjust the row-swap toggle 
      } 

      if (Math.Abs(result[j, j]) < 1.0E-20) // if diagonal after swap is zero . . . 
       return null; // consider a throw 

      for (int i = j + 1; i < rows; ++i) 
      { 
       result[i, j] /= result[j, j]; 
       for (int k = j + 1; k < rows; ++k) 
       { 
        result[i, k] -= result[i, j] * result[j, k]; 
       } 
      } 
     } // main j column loop 

     return result; 
    } // MatrixDecompose 

    // -------------------------------------------------------------------------------------------------------------- 
    private static float MatrixDeterminant(float[,] matrix) 
    { 
     int[] perm; 
     int toggle; 
     float[,] lum = MatrixDecompose(matrix, out perm, out toggle); 
     if (lum == null) 
      throw new Exception("Unable to compute MatrixDeterminant"); 
     float result = toggle; 
     for (int i = 0; i < lum.GetLength(0); ++i) 
      result *= lum[i, i]; 

     return result; 
    } 

    // -------------------------------------------------------------------------------------------------------------- 
    private static float[,] MatrixDuplicate(float[,] matrix) 
    { 
     // allocates/creates a duplicate of a matrix. assumes matrix is not null. 
     float[,] result = MatrixCreate(matrix.GetLength(0), matrix.GetLength(1)); 
     for (int i = 0; i < matrix.GetLength(0); ++i) // copy the values 
      for (int j = 0; j < matrix.GetLength(1); ++j) 
       result[i, j] = matrix[i, j]; 
     return result; 
    } 

    // -------------------------------------------------------------------------------------------------------------- 
    private static float[,] ExtractLower(float[,] matrix) 
    { 
     // lower part of a Doolittle decomposition (1.0s on diagonal, 0.0s in upper) 
     int rows = matrix.GetLength(0); int cols = matrix.GetLength(1); 
     float[,] result = MatrixCreate(rows, cols); 
     for (int i = 0; i < rows; ++i) 
     { 
      for (int j = 0; j < cols; ++j) 
      { 
       if (i == j) 
        result[i, j] = 1.0f; 
       else if (i > j) 
        result[i, j] = matrix[i, j]; 
      } 
     } 
     return result; 
    } 

    // -------------------------------------------------------------------------------------------------------------- 
    private static float[,] ExtractUpper(float[,] matrix) 
    { 
     // upper part of a Doolittle decomposition (0.0s in the strictly lower part) 
     int rows = matrix.GetLength(0); int cols = matrix.GetLength(1); 
     float[,] result = MatrixCreate(rows, cols); 
     for (int i = 0; i < rows; ++i) 
     { 
      for (int j = 0; j < cols; ++j) 
      { 
       if (i <= j) 
        result[i, j] = matrix[i, j]; 
      } 
     } 
     return result; 
    } 
}