不穩定計算誤差

Question

我需要計算矩陣：（X ^（T） * X）^（ - 1）。 傳說的代碼&評論：

x is double[,] array; 
xT - transposed matrix 
^(-1) - inverted matrix 

每次我產生新的隨機矩陣，它的工作，我發現該程序是非常不穩定的，因爲它不與任何輸入數據的正常工作。我很確定，因爲如果一切正常，我需要最終得到標識矩陣，但有時候我得到一個完全可怕的Ineverted矩陣，所以我沒有得到一個標識矩陣。我很失望，因爲我總是使用相同類型的數據，不會轉換任何東西。編譯器是MVS 2010.希望你能幫助我。

這裏是我的的Program.cs：

static void Main(string[] args) 
    { 

     Matrix x = new Matrix(5, 4); 
     //Matrix temp = new Matrix(x.Row, x.Col); 
     //double[] y = new double[x.Row]; 
     //double[] b = new double[x.Row]; 

     //this data isn't calculated correctly. used for debugging 
     x.MatrixX[0, 0] = 7; x.MatrixX[0, 1] = 6; x.MatrixX[0, 2] = 5; x.MatrixX[0, 3] = 8; 
     x.MatrixX[1, 0] = 7; x.MatrixX[1, 1] = 5; x.MatrixX[1, 2] = 8; x.MatrixX[1, 3] = 5; 
     x.MatrixX[2, 0] = 6; x.MatrixX[2, 1] = 8; x.MatrixX[2, 2] = 6; x.MatrixX[2, 3] = 8; 
     x.MatrixX[3, 0] = 8; x.MatrixX[3, 1] = 5; x.MatrixX[3, 2] = 8; x.MatrixX[3, 3] = 7; 
     x.MatrixX[4, 0] = 8; x.MatrixX[4, 1] = 5; x.MatrixX[4, 2] = 6; x.MatrixX[4, 3] = 7; 
     /* 
     7,00000 6,00000 5,00000 8,00000 
     7,00000 5,00000 8,00000 5,00000 
     6,00000 8,00000 6,00000 8,00000 
     8,00000 5,00000 8,00000 7,00000 
     8,00000 5,00000 6,00000 7,00000 
     */ 


     //random matrix generation 
     /* 
     Random rnd = new Random(); 
     for (int i = 0; i < x.Row; i++) 
      for (int j = 0; j < x.Col; j++) 
       x.MatrixX[i, j] = rnd.Next(5, 10); 
     */ 


     /*i'm going to calculate: ( X^(T) * X )^(-1) 
     * 1. transpose X 
     * 2. multiply X and (1) 
     * 3. invert matrix (2) 
     * +4. i wanna check the results: Multilate of (2) and (3) = Identity_matrix. 
     * */ 
     Matrix.Display(x); 

     //1 
     Matrix xt = Matrix.Transpose(x); 
     Matrix.Display(xt); 

     //2 
     Matrix xxt = Matrix.Multiply(x, xt); 
     Matrix.Display(xxt); 

     //3 
     Matrix xxtinv = Matrix.Invert(Matrix.Multiply(x, xt)); 
     Matrix.Display(xxtinv); 

     //4 
     Console.WriteLine("Invert(xxt) * xxt. IdentityMatrix:"); 
     Matrix IdentityMatrix = Matrix.Multiply(xxtinv, xxt); 
     Matrix.Display(IdentityMatrix); 

     Console.ReadKey(); 
    } 

這裏是Matrix.cs所有功能：

public class Matrix 
    { 
     private double[,] matrix; 
     private int row; 
     private int col; 

     #region constructors 
     public Matrix(int Row, int Col) 
     { 
      this.row = Row; 
      this.col = Col; 
      matrix = new double[Row, Col]; 
     } 

     public Matrix() 
     { 
      Random rnd = new Random(); 
      Row = rnd.Next(3, 7); 
      Col = rnd.Next(3, 7); 
      matrix = new double[Row, Col]; 
      for (int i = 0; i < Row; i++) 
       for (int j = 0; j < Col; j++) 
        matrix[i, j] = rnd.Next(5, 10); 
     } 

     public Matrix(Matrix a) 
     { 
      this.Col = a.Col; 
      this.Row = a.Row; 
      this.matrix = a.matrix; 
     } 
     #endregion 

     #region properties 
     public int Col 
     { 
      get { return col; } 
      set { col = value; } 
     } 
     public int Row 
     { 
      get { return row; } 
      set { row = value; } 
     } 
     public double[,] MatrixX 
     { 
      get { return matrix; } 
      set { matrix = value; } 
     } 
     #endregion 

     static public Matrix Transpose(Matrix array) 
     { 
      Matrix temp = new Matrix(array.Col, array.Row); 
      for (int i = 0; i < array.Row; i++) 
       for (int j = 0; j < array.Col; j++) 
        temp.matrix[j, i] = array.matrix[i, j]; 
      return temp; 
     } 

     static public void Display(Matrix array) 
     { 
      for (int i = 0; i < array.Row; i++) 
      { 
       for (int j = 0; j < array.Col; j++) 
        Console.Write("{0,5:f2}\t", array.matrix[i, j]); 
       Console.WriteLine(); 
      } 
      Console.WriteLine(); 
     } 

     static public Matrix Multiply(Matrix a, Matrix b) 
     { 
      if (a.Col != b.Row) throw new Exception("multiplication is impossible: a.Col != b.Row"); 

      Matrix r = new Matrix(a.Row, b.Col); 
      for (int i = 0; i < a.Row; i++) 
      { 
       for (int j = 0; j < b.Col; j++) 
       { 
        double sum = 0; 
        for (int k = 0; k < b.Row; k++) 
         sum += a.matrix[i, k] * b.matrix[k, j]; 
        r.matrix[i, j] = sum; 
       } 
      } 
      return r; 
     } 

     static public Matrix Invert(Matrix a) 
     { 
      Matrix E = new Matrix(a.Row, a.Col); 
      double temp = 0; 
      int n = a.Row; 


      for (int i = 0; i < n; i++) 
       for (int j = 0; j < n; j++) 
       { 
        E.matrix[i, j] = 0.0; 

        if (i == j) 
         E.matrix[i, j] = 1.0; 
       } 

      for (int k = 0; k < n; k++) 
      { 
       temp = a.matrix[k, k]; 

       for (int j = 0; j < n; j++) 
       { 
        a.matrix[k, j] /= temp; 
        E.matrix[k, j] /= temp; 
       } 

       for (int i = k + 1; i < n; i++) 
       { 
        temp = a.matrix[i, k]; 

        for (int j = 0; j < n; j++) 
        { 
         a.matrix[i, j] -= a.matrix[k, j] * temp; 
         E.matrix[i, j] -= E.matrix[k, j] * temp; 
        } 
       } 
      } 

      for (int k = n - 1; k > 0; k--) 
      { 
       for (int i = k - 1; i >= 0; i--) 
       { 
        temp = a.matrix[i, k]; 

        for (int j = 0; j < n; j++) 
        { 
         a.matrix[i, j] -= a.matrix[k, j] * temp; 
         E.matrix[i, j] -= E.matrix[k, j] * temp; 
        } 
       } 
      } 


      for (int i = 0; i < n; i++) 
       for (int j = 0; j < n; j++) 
       { 
        a.matrix[i, j] = E.matrix[i, j]; 
       } 

      return a; 
     } 
    } 

Answer 1

在你的榜樣，X *轉的決定（ x）是零。結果，沒有反轉，這可能是你爲什麼會得到奇怪的結果。

我也注意到你的Inverse函數修改傳遞給它的矩陣。這可能應該修改以避免這種情況。