C# Matrix.Inverse方法代码示例

 public void TestInverse()
 {
     double[] data = new double[] { 1, 3, 2, 5};
     Matrix m = new Matrix(data,2,2);
     Assert.AreEqual(m[0, 0], 1);
     Matrix inv = m.Inverse();
     Assert.AreEqual(inv[0, 0], -5);
 }

        public void ExceptionZeroDeterminant()
        {
            var matrix = new Matrix(3, 3);

            // [ 4,  4,  4 ]
            // [ 4,  4,  4 ]
            // [ 4,  4,  4 ]
            // Determinant = 0

            matrix[0, 0] = 4;
            matrix[0, 1] = 4;
            matrix[0, 2] = 4;

            matrix[1, 0] = 4;
            matrix[1, 1] = 4;
            matrix[1, 2] = 4;

            matrix.Inverse();
        }

 public double[] AddCovariance(double[,] matrix, double[] vector)
 {
     var eigen = new Eigen(matrix);
     var p = new Matrix();
     for (int i = 0; i < 3; i++)
     {
         for (int j = 0; j < 3; j++)
         {
             p[i, j] = eigen.Eigenvectors[i, j].Real;
         }
     }
     var y = Matrix.Multiply(Matrix.Multiply(p.Inverse(), matrix), p);
     var diagonalizing = new Matrix();
     for (int i = 0; i < 3; i++)
     {
         diagonalizing[i, i] = Math.Sqrt(y[i, i]);
     }
     return Matrix.Multiply(Matrix.Multiply(p, diagonalizing), vector);
 }

        public void InverseInterface()
        {
            // [  1, -1,  3 ] -1
            // [  2,  1,  2 ]
            // [ -2, -2,  1 ]
            //
            // =
            //
            // [  1,        -1,   -1 ]
            // [  -(6/5),  7/5,  4/5 ]
            // [  -(2/5),  4/5,  3/5 ]

            IMathematicalMatrix matrix = new Matrix(3, 3);
            matrix[0, 0] = 1;
            matrix[0, 1] = -1;
            matrix[0, 2] = 3;

            matrix[1, 0] = 2;
            matrix[1, 1] = 1;
            matrix[1, 2] = 2;

            matrix[2, 0] = -2;
            matrix[2, 1] = -2;
            matrix[2, 2] = 1;

            var a = matrix.Inverse();

            Assert.AreEqual(a[0, 0], 1, 0.000000001);
            Assert.AreEqual(a[0, 1], -1, 0.000000001);
            Assert.AreEqual(a[0, 2], -1, 0.000000001);

            Assert.AreEqual(a[1, 0], (6F / 5F) * -1F, 0.00001);
            Assert.AreEqual(a[1, 1], 7F / 5F, 0.00001);
            Assert.AreEqual(a[1, 2], 4F / 5F, 0.00001);

            Assert.AreEqual(a[2, 0], (2F / 5F) * -1F, 0.00001);
            Assert.AreEqual(a[2, 1], 4F / 5F, 0.00001);
            Assert.AreEqual(a[2, 2], 3F / 5F, 0.00001);
        }

        /// <summary>
        /// Samples a matrix normal distributed random variable.
        /// </summary>
        /// <param name="rnd">The random number generator to use.</param>
        /// <param name="m">The mean of the matrix normal.</param>
        /// <param name="v">The covariance matrix for the rows.</param>
        /// <param name="k">The covariance matrix for the columns.</param>
        /// <exception cref="ArgumentOutOfRangeException">If the dimensions of the mean and two covariance matrices don't match.</exception>
        /// <returns>a sequence of samples from the distribution.</returns>
        public static Matrix<double> Sample(Random rnd, Matrix<double> m, Matrix<double> v, Matrix<double> k)
        {
            if (Control.CheckDistributionParameters && !IsValidParameterSet(m, v, k))
            {
                throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
            }

            var n = m.RowCount;
            var p = m.ColumnCount;

            // Compute the Kronecker product of V and K, this is the covariance matrix for the stacked matrix.
            var vki = v.KroneckerProduct(k.Inverse());

            // Sample a vector valued random variable with VKi as the covariance.
            var vector = SampleVectorNormal(rnd, new DenseVector(n * p, 0.0), vki);

            // Unstack the vector v and add the mean.
            var r = m.Clone();
            for (var i = 0; i < n; i++)
            {
                for (var j = 0; j < p; j++)
                {
                    r[i, j] += vector[(j * n) + i];
                }
            }

            return r;
        }

        /// <summary>
        /// Samples an inverse Wishart distributed random variable by sampling
        /// a Wishart random variable and inverting the matrix.
        /// </summary>
        /// <param name="rnd">The random number generator to use.</param>
        /// <param name="nu">The degrees of freedom.</param>
        /// <param name="s">The scale matrix.</param>
        /// <returns>a sample from the distribution.</returns>
        public static Matrix<double> Sample(Random rnd, double nu, Matrix<double> s)
        {
            if (Control.CheckDistributionParameters && !IsValidParameterSet(nu, s))
            {
                throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
            }

            var r = Wishart.Sample(rnd, nu, s.Inverse());
            return r.Inverse();
        }

        public void Convertor()
        {
            var eigen1 = new Eigen(matrix1);

            var p1 = new Matrix();
            for (int i = 0; i < 3; i++)
            {
                for (int j = 0; j < 3; j++)
                {
                    p1[i, j] = eigen1.Eigenvectors[i, j].Real;
                }
            }
            var y1 = Matrix.Multiply(Matrix.Multiply(p1.Inverse(), matrix1), p1);
            var diagonalizing = new Matrix();
            for (int i = 0; i < 3; i++)
            {
                diagonalizing[i, i] = 1 / Math.Sqrt(y1[i, i]);
            }
            var y2 = Matrix.Multiply(p1.Transpose(),Matrix.Multiply(matrix2,p1));
            var z2 = Matrix.Multiply(Matrix.Multiply(diagonalizing.Transpose(), y2),diagonalizing);

            var eigen2 = new Eigen(z2);
            var eigen_vectors2 = new Matrix();
            for (int i = 0; i < 3; i++)
            {
                for (int j = 0; j < 3; j++)
                {
                    eigen_vectors2[i, j] = eigen2.Eigenvectors[i, j].Real;
                }
            }
            var v2 = Matrix.Multiply(Matrix.Multiply(eigen_vectors2.Transpose(),z2), eigen_vectors2);
            var new_means1 = Matrix.Multiply(eigen_vectors2,Matrix.Multiply(diagonalizing, Matrix.Multiply(p1.Transpose(), means1)));
            var new_means2 = Matrix.Multiply(eigen_vectors2,Matrix.Multiply(Matrix.Multiply(diagonalizing, p1.Transpose()), means2));
            PrintMatrix(y1, "y1");
            PrintMatrix(y2, "y2");
            PrintMatrix(z2, "z2");
            PrintMatrix(v2, "v2");

            eigen2 = new Eigen(matrix2);
            var p2 = new Matrix();
            for (int i = 0; i < 3; i++)
            {
                for (int j = 0; j < 3; j++)
                {
                    p2[i, j] = eigen2.Eigenvectors[i, j].Real;
                }
            }
            var y = Matrix.Multiply(Matrix.Multiply(p2.Inverse(), matrix2), p2);
            var diagonalizing2 = new Matrix();
            for (int i = 0; i < 3; i++)
            {
                diagonalizing2[i, i] = Math.Sqrt(y[i, i]);
            }
            var diagonalizing3 = new Matrix();
            for (int i = 0; i < 3; i++)
            {
                diagonalizing3[i, i] = Math.Sqrt(y1[i, i]);
            }
            for (int i = 0; i < 200; i++)
            {
                var penBlue = new Pen(Color.Blue);
                var penRed = new Pen(Color.Red);
                var vector1 = GenerateGaussianVector(means1);
                vector1 = Matrix.Multiply(p1, Matrix.Multiply(diagonalizing3, vector1));
                var vector2 = GenerateGaussianVector(means2);
                vector2 = Matrix.Multiply(p2, Matrix.Multiply(diagonalizing2, vector2));
                // x1-x2 before
                var g = pn_12before.CreateGraphics();
                var point = new Point((int)(vector1[0] * 6 + 150), (int)(vector1[1] * 6 + 130));
                var s = new System.Drawing.Size(2,2);
                var circle = new Rectangle(point,s);
                g.DrawRectangle(penBlue, circle);

                point = new Point((int)(vector2[0] * 6 + 150), (int)(vector2[1] * 6 + 130));
                circle = new Rectangle(point, s);
                g.DrawRectangle(penRed, circle);
                // x1-x3 before
                g = pn_13before.CreateGraphics();
                point = new Point((int)(vector1[0] * 6 + 150), (int)(vector1[2] * 6 + 130));
                circle = new Rectangle(point, s);
                g.DrawRectangle(penBlue, circle);

                point = new Point((int)(vector2[0] * 6 + 150), (int)(vector2[2] * 6 + 130));
                circle = new Rectangle(point, s);
                g.DrawRectangle(penRed, circle);
                // after
                vector1 = GenerateGaussianVector(new_means1);
                vector2 = GenerateGaussianVector(new_means2);
                var eigen = new Eigen(v2);
                var m = new Matrix();
                for (int k = 0; k < 3; k++)
                {
                    for (int j = 0; j < 3; j++)
                    {
                        m[k, j] = eigen.Eigenvectors[k, j].Real;
                    }
                }
                var result1 = Matrix.Multiply(Matrix.Multiply(m.Inverse(), v2), m);
                var diagonal = new Matrix();
                for (int k = 0; k < 3; k++)
//.........这里部分代码省略.........

        /// <summary>
        /// Compute finite difference coefficients according to the method provided here:
        /// 
        /// http://en.wikipedia.org/wiki/Finite_difference_coefficients
        ///
        /// </summary>
        /// <returns>An array of the coefficients for FD.</returns>
        public double[] CreateCoefficients()
        {
            var result = new double[_pointCount];

            var delts = new Matrix(_pointCount, _pointCount);
            double[][] t = delts.Data;

            for (int j = 0; j < _pointCount; j++)
            {
                double delt = (j - _center);
                double x = 1.0;

                for (int k = 0; k < _pointCount; k++)
                {
                    t[j][k] = x/EncogMath.Factorial(k);
                    x *= delt;
                }
            }

            Matrix invMatrix = delts.Inverse();
            double f = EncogMath.Factorial(_pointCount);

            for (int k = 0; k < _pointCount; k++)
            {
                result[k] = (Math.Round(invMatrix.Data[1][k]*f))/f;
            }

            return result;
        }

        public void InvertPivotTest()
        {
            var matrix = new Matrix<double>(new double[,] { { 5, 2, 1 }, { 1, 4, 3 }, { 1, 10, 2 } });

            var res = matrix.Inverse();
            double delta = 0.0001;
            Assert.AreEqual(0.2157, res[1, 1], delta);
            Assert.AreEqual(-0.0588, res[1, 2], delta);
            Assert.AreEqual(-0.0196, res[1, 3], delta);
            Assert.AreEqual(-0.0098, res[2, 1], delta);
            Assert.AreEqual(-0.0882, res[2, 2], delta);
            Assert.AreEqual(0.1373, res[2, 3], delta);
            Assert.AreEqual(-0.0588, res[3, 1], delta);
            Assert.AreEqual(0.4706, res[3, 2], delta);
            Assert.AreEqual(-0.1765, res[3, 3], delta);
            MatrixHelpers.NotNaNOrInfinity(res);
        }

        public void InverseTest()
        {
            var target = new Matrix<double>(new double[,] { { 57, -76, -32 }, { 27, -72, -74 }, { -93, -2, -4 } });
            var expected = new Matrix<double>(new double[,] { { -0.0004556401744, 0.0007810974419, -0.01080518128 }, { -0.02274946300, 0.01042765085, -0.01091583675 }, { 0.02196836555, -0.02337434095, 0.006678383128 } });
            var actual = target.Inverse();

            double delta = 0.00000000001;

            MatrixHelpers.Compare(expected, actual, delta);
        }

        //DLT算法
        public void Adjustment()
        {
            int count=data.oCount;
            Matrix B=new Matrix(count * 2, 11);
            Matrix l=new Matrix(count * 2, 1);

            //求L系数初值
            for (int j = 0; j < data.Count;j++)
            {
                for (int i = 0; i < data.oCount;i++)
                {
                    double x = data.oIPoints[j][i].X,
                        y = data.oIPoints[j][i].Y,
                        X = data.oPoints[i].X,
                        Y = data.oPoints[i].Y,
                        Z = data.oPoints[i].Z;
                    B[2 * i,0] = X;
                    B[2 * i,1] = Y;
                    B[2 * i,2] = Z;
                    B[2 * i,3] = 1;
                    B[2 * i,8] = -x*X;
                    B[2 * i,9] = -x*Y;
                    B[2 * i,10] = -x*Z;
                    B[2 * i + 1,4] = X;
                    B[2 * i + 1,5] = Y;
                    B[2 * i + 1,6] = Z;
                    B[2 * i + 1,7] = 1;
                    B[2 * i + 1,8] = -y*X;
                    B[2 * i + 1,9] = -y*Y;
                    B[2 * i + 1,10] = -y*Z;
                    l[2 * i,0] = x;
                    l[2 * i + 1,0] =  y;
                }
                int xCount = 12;//未知数个数
                Matrix L0 = ((B.T()*B).Inverse()*(B.T()*l));
                Matrix L=new Matrix(xCount, 1);
                for (int i = 0; i < xCount-1; i++)
                    L[i,0] = L0[i,0];
                Matrix M=new Matrix(count * 2, xCount);
                Matrix W=new Matrix(count * 2, 1);
                double f = 9999;

                //迭代求解L系数
                while (abs(f - inE[j].fx) >= 0.01)
                {
                    f = inE[j].fx;
                    double x0 = (L[0,0] * L[8,0] + L[1,0] * L[9,0] + L[2,0] * L[10,0]) /
                        (pow2(L[8,0]) + pow2(L[9,0]) + pow2(L[10,0])),
                        y0 = (L[4,0] * L[8,0] + L[5,0] * L[9,0] + L[6,0] * L[10,0]) /
                        (pow2(L[8,0]) + pow2(L[9,0]) + pow2(L[10,0]));
                    for (int i = 0; i < data.oCount; i++)
                    {
                        double x = data.oIPoints[j][i].X,
                               y = data.oIPoints[j][i].Y,
                               X = data.oPoints[i].X,
                               Y = data.oPoints[i].Y,
                               Z = data.oPoints[i].Z;
                        double A = X*L[8,0] +Y*L[9,0] +Z*L[10,0] + 1;
                        double r_2 = (x - x0)*(x - x0) + (y - y0)*(y - y0);
                        M[2 * i,0] = X / A;
                        M[2 * i,1] = Y / A;
                        M[2 * i,2] = Z / A;
                        M[2 * i,3] = 1 / A;
                        M[2 * i,8] = -X*x / A;
                        M[2 * i,9] = -Y*x / A;
                        M[2 * i,10] = -Z*x / A;
                        M[2 * i,11] = -(x - x0)*r_2;
                        M[2 * i + 1,4] = X / A;
                        M[2 * i + 1,5] = Y / A;
                        M[2 * i + 1,6] = Z / A;
                        M[2 * i + 1,7] = 1 / A;
                        M[2 * i + 1,8] = -X*y / A;
                        M[2 * i + 1,9] = -Y*y / A;
                        M[2 * i + 1,10] = -Z*y / A;
                        M[2 * i + 1,11] = -(y - y0)*r_2;
                        W[2 * i,0] = x / A;
                        W[2 * i + 1,0] = y / A;
                    }
                    Matrix nL = (M.T()*M).Inverse()*M.T()*W;
                    double dbeta, ds, fx, fy,Xs,Ys,Zs,a3,b3,c3,b1,b2;
                    double gama3 = 1 / sqrt(pow2(nL[8,0]) + pow2(nL[9,0]) + pow2(nL[10,0]));
                    x0 = (nL[0,0] * nL[8,0] + nL[1,0] * nL[9,0] + nL[2,0] * nL[10,0]) /
                        (pow2(nL[8,0]) + pow2(nL[9,0]) + pow2(nL[10,0]));
                    y0 = (nL[4,0] * nL[8,0] + nL[5,0] * nL[9,0] + nL[6,0] * nL[10,0]) /
                        (pow2(nL[8,0]) + pow2(nL[9,0]) + pow2(nL[10,0]));
                    double At = gama3*gama3*(pow2(nL[0,0]) + pow2(nL[1,0]) + pow2(nL[2,0])) - x0*x0,
                        Bt = gama3*gama3*(pow2(nL[4,0]) + pow2(nL[5,0]) + pow2(nL[6,0])) - y0*y0,
                        Ct = gama3*gama3*(nL[0,0] * nL[4,0] + nL[1,0] * nL[5,0] + nL[2,0] * nL[6,0]) - x0*y0;
                    if (Ct >= 0)
                    {
                        dbeta = -asin(sqrt(Ct*Ct / At / Bt));
                    }
                    else
                    {
                        dbeta = asin(sqrt(Ct*Ct / At / Bt));
                    }
                    ds = sqrt(At / Bt) - 1;
                    fx = sqrt((At*Bt - Ct*Ct) / Bt);
                    fy = sqrt((At*Bt - Ct*Ct) / At);
//.........这里部分代码省略.........

        public void Simple2()
        {
            // [  19.6, 24.974999999999998,  36.999999999999993 ] -1
            // [  24.974999999999998,  52.125000000000014,  67.5 ]
            // [  36.999999999999993, 67.5,  100.0 ]
            //
            // =
            //
            // [  1,        -1,   -1 ]
            // [  -(6/5),  7/5,  4/5 ]
            // [  -(2/5),  4/5,  3/5 ]

            var matrix = new Matrix(3, 3);
            matrix[0, 0] = 19.6;
            matrix[0, 1] = 24.974999999999998;
            matrix[0, 2] = 36.999999999999993;

            matrix[1, 0] = 24.974999999999998;
            matrix[1, 1] = 52.125000000000014;
            matrix[1, 2] = 67.5;

            matrix[2, 0] = 36.999999999999993;
            matrix[2, 1] = 67.5;
            matrix[2, 2] = 100.0;

            var inv = matrix.Inverse();

            Assert.AreEqual(inv[0, 0], 0.169204737732656, 0.000000001);
            Assert.AreEqual(inv[0, 1], -1.44028932924345E-16, 0.000000001);
            Assert.AreEqual(inv[0, 2], -0.0626057529610827, 0.000000001);

            Assert.AreEqual(inv[1, 0], -7.20806576118106E-17, 0.00001);
            Assert.AreEqual(inv[1, 1], 0.152380952380952, 0.00001);
            Assert.AreEqual(inv[1, 2], -0.102857142857143, 0.00001);

            Assert.AreEqual(inv[2, 0], -0.0626057529610828, 0.00001);
            Assert.AreEqual(inv[2, 1], -0.102857142857142, 0.00001);
            Assert.AreEqual(inv[2, 2], 0.102592700024172, 0.00001);
        }

        public void Convertor()
        {
            #region Calculate v2 new means
            var eigen1 = new Eigen(matrix1);
            var p1 = new Matrix();
            for (int i = 0; i < 3; i++)
            {
                for (int j = 0; j < 3; j++)
                {
                    p1[i, j] = eigen1.Eigenvectors[i, j].Real;
                }
            }
            var y1 = Matrix.Multiply(Matrix.Multiply(p1.Inverse(), matrix1), p1);
            var diagonalizing = new Matrix();
            for (int i = 0; i < 3; i++)
            {
                diagonalizing[i, i] = 1 / Math.Sqrt(y1[i, i]);
            }
            var y2 = Matrix.Multiply(p1.Transpose(), Matrix.Multiply(matrix2, p1));
            var z2 = Matrix.Multiply(Matrix.Multiply(diagonalizing.Transpose(), y2), diagonalizing);
            var eigen2 = new Eigen(z2);
            var eigen_vectors2 = new Matrix();
            for (int i = 0; i < 3; i++)
            {
                for (int j = 0; j < 3; j++)
                {
                    eigen_vectors2[i, j] = eigen2.Eigenvectors[i, j].Real;
                }
            }
            var v2 = Matrix.Multiply(Matrix.Multiply(eigen_vectors2.Transpose(), z2), eigen_vectors2);
            new_means1 = Matrix.Multiply(eigen_vectors2, Matrix.Multiply(diagonalizing, Matrix.Multiply(p1.Transpose(), means1)));
            new_means2 = Matrix.Multiply(eigen_vectors2, Matrix.Multiply(Matrix.Multiply(diagonalizing, p1.Transpose()), means2));
            #endregion

            #region generate 200 points
            for (int i = 0; i < 200; i++)
            {
                var orig_vector1 = GenerateGaussianVector(means1);
                var orig_vector2 = GenerateGaussianVector(means2);
                var vector1 = AddCovariance(matrix1, orig_vector1);
                var vector2 = AddCovariance(matrix2, orig_vector2);
                //for (int j = 0; i < 3;j++ )
                //{
                //    orig_vectors1[i, j] = orig_vector1[j];
                //    orig_vectors2[i, j] = orig_vector2[j];
                //}
                try
                {
                    Draw(pn_12before, Color.Blue, vector1[0], vector1[1]);
                    Draw(pn_13before, Color.Blue, vector1[0], vector1[2]);
                    Draw(pn_12before, Color.Red, vector2[0], vector2[1]);
                    Draw(pn_13before, Color.Red, vector2[0], vector2[2]);

                    // after
                    //vector1 = ModifyVector(vector1, means1);
                    //vector2 = ModifyVector(vector2, means2);
                    vector1 = GenerateGaussianVector(new_means1);
                    vector2 = GenerateGaussianVector(new_means2);

                    vector2 = AddCovariance(z2, vector2);

                    Draw(pn_12after, Color.Blue, vector1[0], vector1[1]);
                    Draw(pn_13after, Color.Blue, vector1[0], vector1[2]);
                    Draw(pn_12after, Color.Red, vector2[0], vector2[1]);
                    Draw(pn_13after, Color.Red, vector2[0], vector2[2]);
                }
                catch(Exception)
                { }
            }
            #endregion

            #region Calculate discreminant function
            var E = new double[,]{{1,0,0},{0,1,0},{0,0,1}};
            var identity = new Matrix(E);

            DrawDiscreminant(pn_12before, matrix1, matrix2, means1, means2, 0, 1);
            DrawDiscreminant(pn_13before, matrix1, matrix2, means1, means2, 0, 2);
            DrawDiscreminant(pn_12after, identity, z2, new_means1, new_means2, 0, 1);
            DrawDiscreminant(pn_13after, identity, z2, new_means1, new_means2, 0, 2);
            #endregion

            #region Test
            var correct = 0;
            var wrong = 0;
            // Before
            for(int i=0;i<200;i++)
            {
                var vector1 = GenerateGaussianVector(means1);
                var vector2 = GenerateGaussianVector(means2);
                vector1 = AddCovariance(matrix1, vector1);
                vector2 = AddCovariance(matrix2, vector2);

                if (1 == Test(matrix1, matrix2, means1, means2, vector1))
                {
                    correct++;
                }
                else
                {
                    wrong++;
                }
//.........这里部分代码省略.........

 public int Test(Matrix matrix1, Matrix matrix2, double[] means1, double[] means2, double[] vector)
 {
     var m1 = new Matrix(matrix1);
     var m2 = new Matrix(matrix2);
     var matrix_means1 = new Vector(means1);
     var matrix_means2 = new Vector(means2);
     var matrix_vector = new Vector(vector);
     // Calculate A B C
     var A = m2.Inverse() - m1.Inverse();
     var B = 2 * (Matrix.Multiply(matrix_means1.ToMatrix().Transpose(), m1.Inverse()) - Matrix.Multiply(matrix_means2.ToMatrix().Transpose(), m2.Inverse()));
     var C = (Matrix.Multiply(Matrix.Multiply(matrix_means2.ToMatrix().Transpose(), m2.Inverse()), matrix_means2.ToMatrix())).Determinant() - (Matrix.Multiply(Matrix.Multiply(matrix_means1.ToMatrix().Transpose(), m1.Inverse()), matrix_means1.ToMatrix()).Determinant() - 2 * (Math.Log10(m1.Determinant() / m2.Determinant())));
     var result = Matrix.Multiply(Matrix.Multiply(matrix_vector.ToMatrix().Transpose(), A), matrix_vector.ToMatrix()).Determinant() + B[0, 0] * matrix_vector[0] + B[0, 2] * matrix_vector[2] + B[0, 2] * matrix_vector[2] + C;
     if (result > 0)
     {
         return 1;
     }
     return 2;
 }

        public void DrawDiscreminant(Panel p, Matrix matrix1, Matrix matrix2, double[] means1, double[] means2, int x, int y)
        {
            var m1 = new Matrix(matrix1);
            var m2 = new Matrix(matrix2);
            var matrix_means1 = new Vector(means1);
            var matrix_means2 = new Vector(means2);
            // Calculate A B C
            var A = m2.Inverse() - m1.Inverse();
            var B = 2 * (Matrix.Multiply(matrix_means1.ToMatrix().Transpose(), m1.Inverse()) - Matrix.Multiply(matrix_means2.ToMatrix().Transpose(), m2.Inverse()));
            var C = (Matrix.Multiply(Matrix.Multiply(matrix_means2.ToMatrix().Transpose(), m2.Inverse()), matrix_means2.ToMatrix())).Determinant() - (Matrix.Multiply(Matrix.Multiply(matrix_means1.ToMatrix().Transpose(), m1.Inverse()), matrix_means1.ToMatrix()).Determinant() - 2 * (Math.Log10(m1.Determinant() / m2.Determinant())));
            var a11 = A[x, x];
            var a12 = A[x, y];
            var a21 = A[y, x];
            var a22 = A[y, y];
            var b1 = B[0, x];
            var b2 = B[0, y];
            for (double i = -20; i < 20; i = i + 0.1)
            {
                var quad = a22;
                var line = a12 * i + a21 * i + b2;
                var cons = C + b1 * i + a11 * i * i;

                if ((line * line - 4 * quad * cons) >= 0)
                {
                    try
                    {
                        if (quad != 0)
                        {
                            var j = (-line + Math.Sqrt(line * line - 4 * quad * cons)) / (2 * quad);
                            Draw(p, Color.Black, i, j);

                            j = (-line - Math.Sqrt(line * line - 4 * quad * cons)) / (2 * quad);
                            Draw(p, Color.Black, i, j);
                        }
                        else if (line != 0)
                        {
                            var j = -cons / line;
                            Draw(p, Color.Black, i, j);
                        }
                    }
                    catch(Exception)
                    { }
                }
            }
        }