C# Decompositions.SingularValueDecomposition类代码示例

        public void MahalanobisTest3()
        {
            // Example from Statistical Distance Calculator
            // http://maplepark.com/~drf5n/cgi-bin/dist.cgi

            double[,] cov = 
            {
                { 1.030303, 2.132728, 0.576716 },
                { 2.132728, 4.510515, 1.185771 },
                { 0.576716, 1.185771, 0.398922 }
            };



            double[] x, y;
            double actual, expected;

            var svd = new SingularValueDecomposition(cov, true, true, true);

            var inv = cov.Inverse();
            var pinv = svd.Inverse();
            Assert.IsTrue(inv.IsEqual(pinv, 1e-6));

            x = new double[] { 2, 4, 1 };
            y = new double[] { 0, 0, 0 };

            {
                var bla = cov.Solve(x);
                var blo = svd.Solve(x);
                var ble = inv.Multiply(x);
                var bli = pinv.Multiply(x);

                Assert.IsTrue(bla.IsEqual(blo, 1e-6));
                Assert.IsTrue(bla.IsEqual(ble, 1e-6));
                Assert.IsTrue(bla.IsEqual(bli, 1e-6));
            }

            expected = 2.0773536867741504;
            actual = Distance.Mahalanobis(x, y, inv);
            Assert.AreEqual(expected, actual, 1e-6);

            actual = Distance.Mahalanobis(x, y, svd);
            Assert.AreEqual(expected, actual, 1e-6);


            x = new double[] { 7, 5, 1 };
            y = new double[] { 1, 0.52, -79 };

            expected = 277.8828871106366;
            actual = Distance.Mahalanobis(x, y, inv);
            Assert.AreEqual(expected, actual, 1e-5);
            actual = Distance.Mahalanobis(x, y, svd);
            Assert.AreEqual(expected, actual, 1e-5);
        }

        public static double SquareMahalanobis(this double[] x, double[] y, SingularValueDecomposition covariance)
        {
            double[] d = new double[x.Length];
            for (int i = 0; i < x.Length; i++)
                d[i] = x[i] - y[i];

            double[] z = covariance.Solve(d);

            double r = 0.0;
            for (int i = 0; i < d.Length; i++)
                r += d[i] * z[i];
            return Math.Abs(r);
        }

        /// <summary>
        ///   Constructs a multivariate Gaussian distribution
        ///   with given mean vector and covariance matrix.
        /// </summary>
        public NormalDistribution(double[] mean, double[,] covariance)
            : base(mean.Length)
        {
            int k = mean.Length;

            this.mean = mean;
            this.covariance = covariance;
            variance = covariance.Diagonal();

            double detSqrt = System.Math.Sqrt(System.Math.Abs(covariance.Determinant()));
            constant = 1.0/(System.Math.Pow(2.0*System.Math.PI, k/2.0)*detSqrt);

            chol = new CholeskyDecomposition(covariance, true);

            if (chol.Determinant == 0)
            {
                // The covariance matrix is singular, use pseudo-inverse
                svd = new SingularValueDecomposition(covariance);
            }
        }

        public void InverseTest()
        {
            double[,] value = new double[,]
            { 
                  { 1.0, 1.0 },
                  { 2.0, 2.0 }
            };

            SingularValueDecomposition target = new SingularValueDecomposition(value);

            double[,] expected = new double[,]
            {
               { 0.1, 0.2 },
               { 0.1, 0.2 }
            };

            double[,] actual = target.Solve(Matrix.Identity(2));
            Assert.IsTrue(Matrix.IsEqual(expected, actual, 0.001));

            actual = target.Inverse();
            Assert.IsTrue(Matrix.IsEqual(expected, actual, 0.001));
        }

        public void InverseTest2()
        {
            int n = 5;

            var I = Matrix.Identity(n);

            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j < n; j++)
                {
                    double[,] value = Matrix.Magic(n);

                    var target = new SingularValueDecomposition(value);

                    double[,] solution = target.Solve(I);
                    double[,] inverse = target.Inverse();
                    double[,] reverse = target.Reverse();

                    Assert.IsTrue(Matrix.IsEqual(solution, inverse, 1e-4));
                    Assert.IsTrue(Matrix.IsEqual(value, reverse, 1e-4));
                }
            }
        }

        public void ConstructorTest1()
        {
            double[,] value = 
            {
               {  4, -2 },
               {  3,  1 },
            };

            var svd = new SingularValueDecomposition(value);
            Assert.AreEqual(2, svd.Rank);


            var target = new GramSchmidtOrthogonalization(value);

            var Q = target.OrthogonalFactor;
            var R = target.UpperTriangularFactor;

            double[,] inv = Q.Inverse();
            double[,] transpose = Q.Transpose();
            double[,] m = Matrix.Multiply(Q, transpose);

            Assert.IsTrue(Matrix.IsEqual(m, Matrix.Identity(2), 1e-10));
            Assert.IsTrue(Matrix.IsEqual(inv, transpose, 1e-10));
        }

        public void SingularValueDecompositionConstructorTest7()
        {
            int count = 100;
            double[,] value = new double[count, 3];
            double[] output = new double[count];

            for (int i = 0; i < count; i++)
            {
                double x = i + 1;
                double y = 2 * (i + 1) - 1;
                value[i, 0] = x;
                value[i, 1] = y;
                value[i, 2] = 1;
                output[i] = 4 * x - y + 3;
            }



            SingularValueDecomposition target = new SingularValueDecomposition(value,
                computeLeftSingularVectors: true,
                computeRightSingularVectors: true);

            {
                double[,] expected = value;
                double[,] actual = target.LeftSingularVectors.Multiply(
                    Matrix.Diagonal(target.Diagonal)).Multiply(target.RightSingularVectors.Transpose());

                // Checking the decomposition
                Assert.IsTrue(Matrix.IsEqual(actual, expected, 1e-8));
            }

            {
                double[] solution = target.Solve(output);

                double[] expected= output;
                double[] actual = value.Multiply(solution);

                Assert.IsTrue(Matrix.IsEqual(actual, expected, 1e-8));
            }
        }

        public void SingularValueDecompositionConstructorTest5()
        {
            // Test using SVD assumption auto-correction feature
            // without computing the left singular vectors.

            double[,] value = new double[,]
             { 
                 { 1, 2 },
                 { 3, 4 },
                 { 5, 6 },
                 { 7, 8 }
             }.Transpose(); // value is 2x4, having less rows than columns.

            SingularValueDecomposition target = new SingularValueDecomposition(value, false, true, true);


            // Checking values
            double[,] U =
            {
                { 0.0, 0.0 },
                { 0.0, 0.0 },
            };

            // U should not have been computed
            Assert.IsTrue(Matrix.IsEqual(target.LeftSingularVectors, U));


            double[,] V = // economy svd
            {
                {  0.152483233310201,  0.822647472225661,  },
                {  0.349918371807964,  0.421375287684580,  },
                {  0.547353510305727,  0.0201031031435023, },
                {  0.744788648803490, -0.381169081397574,  },
            };

            // V can be different, but for the economy SVD it is often equal
            Assert.IsTrue(Matrix.IsEqual(target.RightSingularVectors, V, 0.0001));



            double[,] S = 
            {
                { 14.2690954992615, 0.000000000000000 },
                {  0.0000000000000,	0.626828232417543 },
            };

            // The diagonal values should be equal
            Assert.IsTrue(Matrix.IsEqual(target.Diagonal, Matrix.Diagonal(S), 0.001));
        }

        public void SingularValueDecompositionConstructorTest6()
        {
            // Test using SVD assumption auto-correction feature in place

            double[,] value1 =
            {
                { 2.5,  2.4 },
                { 0.5,  0.7 },
                { 2.2,  2.9 },
                { 1.9,  2.2 },
                { 3.1,  3.0 },
                { 2.3,  2.7 },
                { 2.0,  1.6 },
                { 1.0,  1.1 },
                { 1.5,  1.6 },
                { 1.1,  0.9 }
            };

            double[,] value2 = value1.Transpose();

            var target1 = new SingularValueDecomposition(value1, true, true, true, true);
            var target2 = new SingularValueDecomposition(value2, true, true, true, true);

            Assert.IsTrue(target1.LeftSingularVectors.IsEqual(target2.RightSingularVectors));
            Assert.IsTrue(target1.RightSingularVectors.IsEqual(target2.LeftSingularVectors));
            Assert.IsTrue(target1.DiagonalMatrix.IsEqual(target2.DiagonalMatrix));

        }

        public void SingularValueDecompositionConstructorTest2()
        {
            // test for m-x-n matrices where m > n (4 > 2)

            double[,] value = new double[,]
             { 
                 { 1, 2 },
                 { 3, 4 },
                 { 5, 6 },
                 { 7, 8 }
             }; // value is 4x2, thus having more rows than columns


            SingularValueDecomposition target = new SingularValueDecomposition(value, true, true, false);

            double[,] actual = target.LeftSingularVectors.Multiply(
                Matrix.Diagonal(target.Diagonal)).Multiply(target.RightSingularVectors.Transpose());

            // Checking the decomposition
            Assert.IsTrue(Matrix.IsEqual(actual, value, 0.01));

            double[,] U = // economy svd
            {
                {  0.152483233310201,  0.822647472225661,  },
                {  0.349918371807964,  0.421375287684580,  },
                {  0.547353510305727,  0.0201031031435023, },
                {  0.744788648803490, -0.381169081397574,  },
            };

            // U should be equal except for some sign changes
            Assert.IsTrue(Matrix.IsEqual(target.LeftSingularVectors, U, 0.001));



            // Checking values
            double[,] V =
            {
                {  0.641423027995072, -0.767187395072177 },
                {  0.767187395072177,  0.641423027995072 },
            };

            // V should be equal except for some sign changes
            Assert.IsTrue(Matrix.IsEqual(target.RightSingularVectors, V, 0.0001));


            double[,] S = 
            {
                { 14.2690954992615, 0.000000000000000 },
                {  0.0000000000000,	0.626828232417543 },
            };

            // The diagonal values should be equal
            Assert.IsTrue(Matrix.IsEqual(target.Diagonal, Matrix.Diagonal(S), 0.001));
        }

        public void SingularValueDecompositionConstructorTest4()
        {
            // Test using SVD assumption auto-correction feature
            // without computing the right singular vectors.

            double[,] value = new double[,]
             { 
                 { 1, 2 },
                 { 3, 4 },
                 { 5, 6 },
                 { 7, 8 }
             }.Transpose(); // value is 2x4, having less rows than columns.

            SingularValueDecomposition target = new SingularValueDecomposition(value, true, false, true);


            // Checking values
            double[,] U =
            {
                { 0.641423027995072, -0.767187395072177 },
                { 0.767187395072177,  0.641423027995072 },
            };

            // U should be equal despite some sign changes
            Assert.IsTrue(Matrix.IsEqual(target.LeftSingularVectors, U, 0.001));


            // Checking values
            double[,] V =
            {
                {  0.0, 0.0 },
                {  0.0, 0.0 },
                {  0.0, 0.0 },
                {  0.0, 0.0 },
            };

            // V should not have been computed.
            Assert.IsTrue(Matrix.IsEqual(target.RightSingularVectors, V));


            double[,] S = 
            {
                { 14.2690954992615, 0.000000000000000 },
                {  0.0000000000000,	0.626828232417543 },
            };

            // The diagonal values should be equal
            Assert.IsTrue(Matrix.IsEqual(target.Diagonal, Matrix.Diagonal(S), 0.001));
        }

        private double regress(double[][] inputs, double[] outputs, out double[,] designMatrix, bool robust)
        {
            if (inputs.Length != outputs.Length)
                throw new ArgumentException("Number of input and output samples does not match", "outputs");

            int parameters = Inputs;
            int rows = inputs.Length;   // number of instance points
            int cols = Inputs;          // dimension of each point

            for (int i = 0; i < inputs.Length; i++)
            {
                if (inputs[i].Length != parameters)
                {
                    throw new DimensionMismatchException("inputs", String.Format(
                        "Input vectors should have length {0}. The row at index {1} of the" +
                        " inputs matrix has length {2}.", parameters, i, inputs[i].Length));
                }
            }

            ISolverMatrixDecomposition<double> solver;


            // Create the problem's design matrix. If we
            //  have to add an intercept term, add a new
            //  extra column at the end and fill with 1s.

            if (!addIntercept)
            {
                // Just copy values over
                designMatrix = new double[rows, cols];
                for (int i = 0; i < inputs.Length; i++)
                    for (int j = 0; j < inputs[i].Length; j++)
                        designMatrix[i, j] = inputs[i][j];
            }
            else
            {
                // Add an intercept term
                designMatrix = new double[rows, cols + 1];
                for (int i = 0; i < inputs.Length; i++)
                {
                    for (int j = 0; j < inputs[i].Length; j++)
                        designMatrix[i, j] = inputs[i][j];
                    designMatrix[i, cols] = 1;
                }
            }

            // Check if we have an overdetermined or underdetermined
            //  system to select an appropriate matrix solver method.

            if (robust || cols >= rows)
            {
                // We have more variables than equations, an
                // underdetermined system. Solve using a SVD:
                solver = new SingularValueDecomposition(designMatrix,
                    computeLeftSingularVectors: true,
                    computeRightSingularVectors: true,
                    autoTranspose: true);
            }
            else
            {
                // We have more equations than variables, an
                // overdetermined system. Solve using the QR:
                solver = new QrDecomposition(designMatrix);
            }


            // Solve V*C = B to find C (the coefficients)
            coefficients = solver.Solve(outputs);


            // Calculate Sum-Of-Squares error
            double error = 0.0;
            double e;
            for (int i = 0; i < outputs.Length; i++)
            {
                e = outputs[i] - Compute(inputs[i]);
                error += e * e;
            }

            return error;
        }

        public void SingularValueDecompositionConstructorTest3()
        {
            // Test using SVD assumption auto-correction feature.

            // Test for m-x-n matrices where m < n. The available SVD
            // routine was not meant to be used in this case.

            double[,] value = new double[,]
             { 
                 { 1, 2 },
                 { 3, 4 },
                 { 5, 6 },
                 { 7, 8 }
             }.Transpose(); // value is 2x4, having less rows than columns.

            SingularValueDecomposition target = new SingularValueDecomposition(value, true, true, true);

            double[,] actual = target.LeftSingularVectors.Multiply(
                Matrix.Diagonal(target.Diagonal)).Multiply(target.RightSingularVectors.Transpose());

            // Checking the decomposition
            Assert.IsTrue(Matrix.IsEqual(actual, value, 0.01));

            // Checking values
            double[,] U =
            {
                { 0.641423027995072, -0.767187395072177 },
                { 0.767187395072177,  0.641423027995072 },
            };

            // U should be equal despite some sign changes
            Assert.IsTrue(Matrix.IsEqual(target.LeftSingularVectors, U, 0.001));


            double[,] V = // economy svd
            {
                {  0.152483233310201,  0.822647472225661,  },
                {  0.349918371807964,  0.421375287684580,  },
                {  0.547353510305727,  0.0201031031435023, },
                {  0.744788648803490, -0.381169081397574,  },
            };

            // V can be different, but for the economy SVD it is often equal
            Assert.IsTrue(Matrix.IsEqual(target.RightSingularVectors, V, 0.0001));


            double[,] S = 
            {
                { 14.2690954992615, 0.000000000000000 },
                {  0.0000000000000,	0.626828232417543 },
            };

            // The diagonal values should be equal
            Assert.IsTrue(Matrix.IsEqual(target.Diagonal, Matrix.Diagonal(S), 0.001));
        }

        public void InverseTest3x3()
        {
            double[,] value = 
            { 
                { 6.0, 1.0, 2.0 },
                { 0.0, 8.0, 1.0 },  
                { 2.0, 4.0, 5.0 }  
            };

            Assert.IsFalse(value.IsSingular());

            double[,] expected = new SingularValueDecomposition(value).Inverse();

            double[,] actual = Matrix.Inverse(value);

            Assert.IsTrue(Matrix.IsEqual(expected, actual, 1e-6));
        }

        // ------------------------------------------------------------
        /// <summary>
        ///   Computes the whitening transform for the given data, making
        ///   its covariance matrix equals the identity matrix.
        /// </summary>
        /// <param name="value">A matrix where each column represent a
        ///   variable and each row represent a observation.</param>
        /// <param name="transformMatrix">The base matrix used in the
        ///   transformation.</param>
        /// <returns>
        ///   The transformed source data (which now has unit variance).
        /// </returns>
        /// 
        public static double[,] Whitening(double[,] value, out double[,] transformMatrix)
        {
            if (value == null)
                throw new ArgumentNullException("value");

            int cols = value.GetLength(1);

            double[,] cov = Tools.Covariance(value);

            // Diagonalizes the covariance matrix
            SingularValueDecomposition svd = new SingularValueDecomposition(cov,
                true,  // compute left vectors (to become a transformation matrix)
                false, // do not compute right vectors since they aren't necessary
                true,  // transpose if necessary to avoid erroneous assumptions in SVD
                true); // perform operation in-place, reducing memory usage

            // Retrieve the transformation matrix
            transformMatrix = svd.LeftSingularVectors;

            // Perform scaling to have unit variance
            double[] singularValues = svd.Diagonal;
            for (int i = 0; i < cols; i++)
                for (int j = 0; j < singularValues.Length; j++)
                    transformMatrix[i, j] /= Math.Sqrt(singularValues[j]);

            // Return the transformed data
            return value.Multiply(transformMatrix);
        }