C# DenseMatrix.QR方法代码示例

        /// <summary>
        /// Linear regression simulation for homework Q5-Q6 of the 
        ///  2nd week of the CS1156x "Learning From Data" at eDX
        /// </summary>
        static void RunQ5Q6Simulation()
        {
            const int EXPERIMENT_COUNT = 1000, N = 100;
              Random rnd = new Random();

              double avgEin = 0, avgEout = 0;
              for (int i = 1; i <= EXPERIMENT_COUNT; i++)
              {
            //pick a random line y = a1 * x + b1
            double x1 = rnd.NextDouble(), y1 = rnd.NextDouble(), x2 = rnd.NextDouble(), y2 = rnd.NextDouble();
            double a = (y1 - y2) / (x1 - x2), b = y1 - a * x1;
            Func<double, double, int> f = (x, y) => a * x + b >= y ? 1 : -1;

            //generate training set of N random points
            var X = new DenseMatrix(N, 3);
            var Y = new DenseVector(N);
            for (int j = 0; j < N; j++)
            {
              X[j, 0] = 1;
              X[j, 1] = rnd.NextDouble() * 2 - 1;
              X[j, 2] = rnd.NextDouble() * 2 - 1;

              Y[j] = f(X[j, 1], X[j, 2]);
            }

            var W = X.QR().Solve(DenseMatrix.Identity(X.RowCount)).Multiply(Y);

            Func<double, double, int> h = (x, y) => W[0] + W[1] * x + W[2] * y >= 0 ? 1 : -1;

            //find Ein
            int count = 0;
            for (int j = 0; j < N; j++) if (h(X[j, 1], X[j, 2]) != Y[j]) count++;
            avgEin += (count + 0.0) / N;

            //find p: f != g
            const int P_SAMPLE_COUNT = 1000;
            count = 0;
            for (int j = 1; j <= P_SAMPLE_COUNT; j++)
            {
              double xx = rnd.NextDouble() * 2 - 1;
              double yy = rnd. NextDouble() * 2 - 1;
              if (f(xx, yy) != h(xx, yy)) count++;
            }

            avgEout += (count + 0.0) / P_SAMPLE_COUNT;
              }

              Console.Out.WriteLine("HW2 Q5:");
              Console.Out.WriteLine("\tEin = {0}", avgEin / EXPERIMENT_COUNT);
              Console.Out.WriteLine("HW2 Q6:");
              Console.Out.WriteLine("\tEout = {0}", avgEout / EXPERIMENT_COUNT);
        }

        /// <summary>
        /// Run example
        /// </summary>
        public void Run()
        {
            // Format matrix output to console
            var formatProvider = (CultureInfo)CultureInfo.InvariantCulture.Clone();
            formatProvider.TextInfo.ListSeparator = " ";

            // Solve next system of linear equations (Ax=b):
            // 5*x + 2*y - 4*z = -7
            // 3*x - 7*y + 6*z = 38
            // 4*x + 1*y + 5*z = 43

            // Create matrix "A" with coefficients
            var matrixA = new DenseMatrix(new[,] { { 5.00, 2.00, -4.00 }, { 3.00, -7.00, 6.00 }, { 4.00, 1.00, 5.00 } });
            Console.WriteLine(@"Matrix 'A' with coefficients");
            Console.WriteLine(matrixA.ToString("#0.00\t", formatProvider));
            Console.WriteLine();

            // Create vector "b" with the constant terms.
            var vectorB = new DenseVector(new[] { -7.0, 38.0, 43.0 });
            Console.WriteLine(@"Vector 'b' with the constant terms");
            Console.WriteLine(vectorB.ToString("#0.00\t", formatProvider));
            Console.WriteLine();

            // 1. Solve linear equations using LU decomposition
            var resultX = matrixA.LU().Solve(vectorB);
            Console.WriteLine(@"1. Solution using LU decomposition");
            Console.WriteLine(resultX.ToString("#0.00\t", formatProvider));
            Console.WriteLine();

            // 2. Solve linear equations using QR decomposition
            resultX = matrixA.QR().Solve(vectorB);
            Console.WriteLine(@"2. Solution using QR decomposition");
            Console.WriteLine(resultX.ToString("#0.00\t", formatProvider));
            Console.WriteLine();

            // 3. Solve linear equations using SVD decomposition
            matrixA.Svd(true).Solve(vectorB, resultX);
            Console.WriteLine(@"3. Solution using SVD decomposition");
            Console.WriteLine(resultX.ToString("#0.00\t", formatProvider));
            Console.WriteLine();

            // 4. Solve linear equations using Gram-Shmidt decomposition
            matrixA.GramSchmidt().Solve(vectorB, resultX);
            Console.WriteLine(@"4. Solution using Gram-Shmidt decomposition");
            Console.WriteLine(resultX.ToString("#0.00\t", formatProvider));
            Console.WriteLine();

            // 5. Verify result. Multiply coefficient matrix "A" by result vector "x"
            var reconstructVecorB = matrixA * resultX;
            Console.WriteLine(@"5. Multiply coefficient matrix 'A' by result vector 'x'");
            Console.WriteLine(reconstructVecorB.ToString("#0.00\t", formatProvider));
            Console.WriteLine();

            // To use Cholesky or Eigenvalue decomposition coefficient matrix must be
            // symmetric (for Evd and Cholesky) and positive definite (for Cholesky)
            // Multipy matrix "A" by its transpose - the result will be symmetric and positive definite matrix
            var newMatrixA = matrixA.TransposeAndMultiply(matrixA);
            Console.WriteLine(@"Symmetric positive definite matrix");
            Console.WriteLine(newMatrixA.ToString("#0.00\t", formatProvider));
            Console.WriteLine();

            // 6. Solve linear equations using Cholesky decomposition
            newMatrixA.Cholesky().Solve(vectorB, resultX);
            Console.WriteLine(@"6. Solution using Cholesky decomposition");
            Console.WriteLine(resultX.ToString("#0.00\t", formatProvider));
            Console.WriteLine();

            // 7. Solve linear equations using eigen value decomposition
            newMatrixA.Evd().Solve(vectorB, resultX);
            Console.WriteLine(@"7. Solution using eigen value decomposition");
            Console.WriteLine(resultX.ToString("#0.00\t", formatProvider));
            Console.WriteLine();

            // 8. Verify result. Multiply new coefficient matrix "A" by result vector "x"
            reconstructVecorB = newMatrixA * resultX;
            Console.WriteLine(@"8. Multiply new coefficient matrix 'A' by result vector 'x'");
            Console.WriteLine(reconstructVecorB.ToString("#0.00\t", formatProvider));
            Console.WriteLine();
        }

        /// <summary>
        /// Linear regressionsimulation with non-separable target function
        /// for homework Q8 of the 2nd week of the CS1156x "Learning From Data" at eDX
        /// </summary>
        static void RunQ8Simulation()
        {
            const int EXPERIMENT_COUNT = 1000, N = 100;
              Random rnd = new Random();

              double avgEin = 0;
              for (int i = 1; i <= EXPERIMENT_COUNT; i++)
              {
            Func<double, double, int> f = (x1, x2) => x1 * x1 + x2 * x2 - 0.6 >= 0 ? 1 : -1;

            //generate training set of N random points
            var X = new DenseMatrix(N, 3);
            var Y = new DenseVector(N);
            for (int j = 0; j < N; j++)
            {
              X[j, 0] = 1;
              X[j, 1] = rnd.NextDouble() * 2 - 1;
              X[j, 2] = rnd.NextDouble() * 2 - 1;

              Y[j] = f(X[j, 1], X[j, 2]);

              //not exactly how it was defined in the problem statement, but shall be good enough
              if (rnd.NextDouble() < 0.1) Y[j] = -Y[j];
            }

            var W = X.QR().Solve(DenseMatrix.Identity(X.RowCount)).Multiply(Y);

            Func<double, double, int> h = (x, y) => W[0] + W[1] * x + W[2] * y >= 0 ? 1 : -1;

            //find Ein
            int count = 0;
            for (int j = 0; j < N; j++) if (h(X[j, 1], X[j, 2]) != Y[j]) count++;
            avgEin += (count + 0.0) / N;
              }

              Console.Out.WriteLine("HW2 Q8:");
              Console.Out.WriteLine("\tEin = {0}", avgEin / EXPERIMENT_COUNT);
        }

        /// <summary>
        /// Non-linear-transformed linear regression simulation for homework Q9, Q10 of the 
        ///  2nd week of the CS1156x "Learning From Data" at eDX
        /// </summary>
        static void RunQ9Q10Simulation()
        {
            const int EXPERIMENT_COUNT = 1000, N = 100;
              Random rnd = new Random();

              double avgEout = 0;
              for (int i = 1; i <= EXPERIMENT_COUNT; i++)
              {
            Func<double, double, int> f = (x1, x2) => x1 * x1 + x2 * x2 - 0.6 >= 0 ? 1 : -1;

            //generate training set of N random points
            var X = new DenseMatrix(N, 3);
            var Y = new DenseVector(N);
            for (int j = 0; j < N; j++)
            {
              X[j, 0] = 1;
              X[j, 1] = rnd.NextDouble() * 2 - 1;
              X[j, 2] = rnd.NextDouble() * 2 - 1;

              Y[j] = f(X[j, 1], X[j, 2]);

              // Just flipping each Y with a 10% chance -
              // not exactly how it was defined in the problem statement, but shall be good enough
              if (rnd.NextDouble() < 0.1) Y[j] = -Y[j];
            }

            var XX = new DenseMatrix(N, 6);
            for (int j = 0; j < N; j++)
            {
              XX[j, 0] = 1;
              XX[j, 1] = X[j, 1];
              XX[j, 2] = X[j, 2];
              XX[j, 3] = X[j, 1] * X[j, 2];
              XX[j, 4] = X[j, 1] * X[j, 1];
              XX[j, 5] = X[j, 2] * X[j, 2];
            }

            var W = XX.QR().Solve(DenseMatrix.Identity(XX.RowCount)).Multiply(Y);

            Func<double, double, int> h = (x, y) => W[0] + W[1] * x + W[2] * y + W[3] * x * y + W[4] * x * x + W[5] * y * y >= 0 ? 1 : -1;

            //find p: f != g
            const int P_SAMPLE_COUNT = 1000;
            int count = 0;
            for (int j = 1; j <= P_SAMPLE_COUNT; j++)
            {
              double xx = rnd.NextDouble() * 2 - 1;
              double yy = rnd.NextDouble() * 2 - 1;
              int ff = f(xx, yy);
              if (rnd.NextDouble() < 0.1) ff = -ff;

              if (ff != h(xx, yy)) count++;
            }

            avgEout += (count + 0.0) / P_SAMPLE_COUNT;
              }

              Console.Out.WriteLine("HW2 Q10:");
              Console.Out.WriteLine("\tEout = {0}", avgEout / EXPERIMENT_COUNT);
        }

        /// <summary>
        /// Run example
        /// </summary>
        /// <seealso cref="http://en.wikipedia.org/wiki/QR_decomposition">QR decomposition</seealso>
        public void Run()
        {
            // Format matrix output to console
            var formatProvider = (CultureInfo)CultureInfo.InvariantCulture.Clone();
            formatProvider.TextInfo.ListSeparator = " ";

            // Create 3 x 2 matrix
            var matrix = new DenseMatrix(new[,] { { 1.0, 2.0 }, { 3.0, 4.0 }, { 5.0, 6.0 } });
            Console.WriteLine(@"Initial 3x2 matrix");
            Console.WriteLine(matrix.ToString("#0.00\t", formatProvider));
            Console.WriteLine();

            // Perform QR decomposition (Householder transformations)
            var qr = matrix.QR();
            Console.WriteLine(@"QR decomposition (Householder transformations)");

            // 1. Orthogonal Q matrix
            Console.WriteLine(@"1. Orthogonal Q matrix");
            Console.WriteLine(qr.Q.ToString("#0.00\t", formatProvider));
            Console.WriteLine();

            // 2. Multiply Q matrix by its transpose gives identity matrix
            Console.WriteLine(@"2. Multiply Q matrix by its transpose gives identity matrix");
            Console.WriteLine(qr.Q.TransposeAndMultiply(qr.Q).ToString("#0.00\t", formatProvider));
            Console.WriteLine();

            // 3. Upper triangular factor R
            Console.WriteLine(@"3. Upper triangular factor R");
            Console.WriteLine(qr.R.ToString("#0.00\t", formatProvider));
            Console.WriteLine();

            // 4. Reconstruct initial matrix: A = Q * R
            var reconstruct = qr.Q * qr.R;
            Console.WriteLine(@"4. Reconstruct initial matrix: A = Q*R");
            Console.WriteLine(reconstruct.ToString("#0.00\t", formatProvider));
            Console.WriteLine();

            // Perform QR decomposition (Gram–Schmidt process)
            var gramSchmidt = matrix.GramSchmidt();
            Console.WriteLine(@"QR decomposition (Gram–Schmidt process)");

            // 5. Orthogonal Q matrix
            Console.WriteLine(@"5. Orthogonal Q matrix");
            Console.WriteLine(gramSchmidt.Q.ToString("#0.00\t", formatProvider));
            Console.WriteLine();

            // 6. Multiply Q matrix by its transpose gives identity matrix
            Console.WriteLine(@"6. Multiply Q matrix by its transpose gives identity matrix");
            Console.WriteLine((gramSchmidt.Q.Transpose() * gramSchmidt.Q).ToString("#0.00\t", formatProvider));
            Console.WriteLine();

            // 7. Upper triangular factor R
            Console.WriteLine(@"7. Upper triangular factor R");
            Console.WriteLine(gramSchmidt.R.ToString("#0.00\t", formatProvider));
            Console.WriteLine();

            // 8. Reconstruct initial matrix: A = Q * R
            reconstruct = gramSchmidt.Q * gramSchmidt.R;
            Console.WriteLine(@"8. Reconstruct initial matrix: A = Q*R");
            Console.WriteLine(reconstruct.ToString("#0.00\t", formatProvider));
            Console.WriteLine();
        }

        /// <summary>
        /// Linear regression/Perceptron simulation for homework Q7 of the 
        ///  2nd week of the CS1156x "Learning From Data" at eDX
        /// </summary>
        static void RunQ7Simulation()
        {
            const int EXPERIMENT_COUNT = 1000, N = 10;
              Random rnd = new Random();

              double avgK = 0;
              for (int i = 1; i <= EXPERIMENT_COUNT; i++)
              {
            //pick a random line y = a1 * x + b1
            double x1 = rnd.NextDouble(), y1 = rnd.NextDouble(), x2 = rnd.NextDouble(), y2 = rnd.NextDouble();
            double a = (y1 - y2) / (x1 - x2), b = y1 - a * x1;
            Func<double, double, int> f = (x, y) => a * x + b >= y ? 1 : -1;

            //generate training set of N random points
            var X = new DenseMatrix(N, 3);
            var Y = new DenseVector(N);
            for (int j = 0; j < N; j++)
            {
              X[j, 0] = 1;
              X[j, 1] = rnd.NextDouble() * 2 - 1;
              X[j, 2] = rnd.NextDouble() * 2 - 1;

              Y[j] = f(X[j, 1], X[j, 2]);
            }

            var W = X.QR().Solve(DenseMatrix.Identity(X.RowCount)).Multiply(Y);

            double w0 = W[0], w1 = W[1], w2 = W[2];

            Func<double, double, int> h = (x, y) => w0 + w1 * x + w2 * y >= 0 ? 1 : -1;

            //run Perceptron
            int k = 1;
            while (Enumerable.Range(0, N).Any(j => f(X[j, 1], X[j, 2]) != h(X[j, 1], X[j, 2])))
            {
              //find all misclasified points
              int[] M = Enumerable.Range(0, N).Where(j => f(X[j, 1], X[j, 2]) != h(X[j, 1], X[j, 2])).ToArray();
              int m = M[rnd.Next(0, M.Length)];

              int sign = f(X[m, 1], X[m, 2]);
              w0 += sign;
              w1 += sign * X[m, 1];
              w2 += sign * X[m, 2];
              k++;
            }

            avgK += k;
              }

              Console.Out.WriteLine("HW2 Q7:");
              Console.Out.WriteLine("\tK = {0}", avgK / EXPERIMENT_COUNT);
        }

        /// <summary>
        /// Run example
        /// </summary>
        /// <seealso cref="http://en.wikipedia.org/wiki/Transpose">Transpose</seealso>
        /// <seealso cref="http://en.wikipedia.org/wiki/Invertible_matrix">Invertible matrix</seealso>
        public void Run()
        {
            // Format matrix output to console
            var formatProvider = (CultureInfo)CultureInfo.InvariantCulture.Clone();
            formatProvider.TextInfo.ListSeparator = " ";

            // Create random square matrix
            var matrix = new DenseMatrix(5);
            var rnd = new Random(1);
            for (var i = 0; i < matrix.RowCount; i++)
            {
                for (var j = 0; j < matrix.ColumnCount; j++)
                {
                    matrix[i, j] = rnd.NextDouble();
                }
            }

            Console.WriteLine(@"Initial matrix");
            Console.WriteLine(matrix.ToString("#0.00\t", formatProvider));
            Console.WriteLine();

            // 1. Get matrix inverse
            var inverse = matrix.Inverse();
            Console.WriteLine(@"1. Matrix inverse");
            Console.WriteLine(inverse.ToString("#0.00\t", formatProvider));
            Console.WriteLine();

            // 2. Matrix multiplied by its inverse gives identity matrix
            var identity = matrix * inverse;
            Console.WriteLine(@"2. Matrix multiplied by its inverse");
            Console.WriteLine(identity.ToString("#0.00\t", formatProvider));
            Console.WriteLine();

            // 3. Get matrix transpose
            var transpose = matrix.Transpose();
            Console.WriteLine(@"3. Matrix transpose");
            Console.WriteLine(transpose.ToString("#0.00\t", formatProvider));
            Console.WriteLine();

            // 4. Get orthogonal  matrix, i.e. do QR decomposition and get matrix Q
            var orthogonal = matrix.QR().Q;
            Console.WriteLine(@"4. Orthogonal  matrix");
            Console.WriteLine(orthogonal.ToString("#0.00\t", formatProvider));
            Console.WriteLine();

            // 5. Transpose and multiply orthogonal matrix by iteslf gives identity matrix
            identity = orthogonal.TransposeAndMultiply(orthogonal);
            Console.WriteLine(@"Transpose and multiply orthogonal matrix by iteslf");
            Console.WriteLine(identity.ToString("#0.00\t", formatProvider));
            Console.WriteLine();
        }

        public static int solve(double[,] A, double[] fitz, CenterArrayNode CenterNode, IRBFPolynomial Poly, bool MathNet)
        {
            var matrixA = new DenseMatrix(A);
               var vectorB = new DenseVector(fitz);
               //Vector<double> resultX = matrixA.LU().Solve(vectorB);
               Vector<double> resultX = matrixA.QR().Solve(vectorB);

               //matrixA.GramSchmidt().Solve(vectorB, resultX);
               List<double> w2 = new List<double>(resultX.ToArray());

               int i = 0;

               w2.ForEach((double weight) =>
               {
                    if (i < CenterNode.Centers.Count)
                         CenterNode[i].w = weight; // set the center's weight
                    else
                         Poly[i - CenterNode.Centers.Count] = weight;//store the polynomial coefficients

                    ++i;
               });

               return 0;
        }

        /// <summary>
        ///     Train.  Single iteration.
        /// </summary>
        public void Iteration()
        {
            int rowCount = _trainingData.Count;
            int inputColCount = _trainingData[0].Input.Length;

            Matrix<double> xMatrix = new DenseMatrix(rowCount, inputColCount + 1);
            Matrix<double> yMatrix = new DenseMatrix(rowCount, 1);

            for (int row = 0; row < _trainingData.Count; row++)
            {
                BasicData dataRow = _trainingData[row];
                int colSize = dataRow.Input.Count();

                xMatrix[row, 0] = 1;
                for (int col = 0; col < colSize; col++)
                {
                    xMatrix[row, col + 1] = dataRow.Input[col];
                }
                yMatrix[row, 0] = dataRow.Ideal[0];
            }

            // Calculate the least squares solution
            QR qr = xMatrix.QR();
            Matrix<double> beta = qr.Solve(yMatrix);

            double sum = 0.0;
            for (int i = 0; i < inputColCount; i++)
                sum += yMatrix[i, 0];
            double mean = sum/inputColCount;

            for (int i = 0; i < inputColCount; i++)
            {
                double dev = yMatrix[i, 0] - mean;
                _sst += dev*dev;
            }

            Matrix<double> residuals = xMatrix.Multiply(beta).Subtract(yMatrix);
            _sse = residuals.L2Norm()*residuals.L2Norm();

            for (int i = 0; i < _algorithm.LongTermMemory.Length; i++)
            {
                _algorithm.LongTermMemory[i] = beta[i, 0];
            }

            // calculate error
            _errorCalculation.Clear();
            foreach (BasicData dataRow in _trainingData)
            {
                double[] output = _algorithm.ComputeRegression(dataRow.Input);
                _errorCalculation.UpdateError(output, dataRow.Ideal, 1.0);
            }
            _error = _errorCalculation.Calculate();
        }

 private double[] Polyfit(double[] x, double[] y, int degree)
 {
     // Vandermonde matrix
     var v = new DenseMatrix(x.Length, degree + 1);
     for (int i = 0; i < v.RowCount; i++)
         for (int j = 0; j <= degree; j++) v[i, j] = Math.Pow(x[i], j);
     var yv = new DenseVector(y).ToColumnMatrix();
     QR<double> qr = v.QR();
     // Math.Net doesn't have an "economy" QR, so:
     // cut R short to square upper triangle, then recompute Q
     var r = qr.R.SubMatrix(0, degree + 1, 0, degree + 1);
     var q = v.Multiply(r.Inverse());
     var p = r.Inverse().Multiply(q.TransposeThisAndMultiply(yv));
     return p.Column(0).ToArray();
 }

            public CurveFit(List<CurvePoint> data)
            {
                var x = new DenseMatrix(data.Count, data.Count);

                for(int i = 0; i < data.Count; i++)
                {
                    for (int j = 0; j < data.Count; j++)
                    {
                        x[i, j] = Pow(data[i].x, data.Count - 1 - j);
                    }
                }

                var y = new DenseVector(data.Select(p => (double)p.y).ToArray());

                var coefficents = x.QR().Solve(y);

                Coefficents = coefficents.Select(c => Round(c)).ToList();
            }