计算错误不稳定

我需要计算矩阵：(X^(T(*X(^(-1(。代码图例&评论：

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

每次我生成新的随机矩阵来处理它时，我都发现这个程序非常不稳定，因为它不能正确地处理任何输入数据。我确信这一点，因为如果一切都好的话，我最终需要得到恒等矩阵，但有时我会得到一个非常糟糕的内倾矩阵，所以我得不到恒等矩阵。我很失望，因为我总是使用相同类型的数据，不转换任何内容。编译器是MVS 2010。希望你能帮助我。

这是我的程序。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();
    }

这是矩阵.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;
        }
    }