C# DenseMatrix.Inverse方法代码示例

        public void Optimize(Dictionary<double, double> values)
        {
            var n = _f.Functions.Count();
              var xs = values.Select(v => v.Key).ToList();
              var ys = values.Select(v => v.Value).ToList();
              var fs = new List<List<double>>(n);

              for (var i = 0; i < n; i++)
              {
            fs[i] = _f.Functions[i].Evaluate(xs);
              }

              var matrix = new DenseMatrix(n, n);
              var vector = new DenseVector(n);
              for (var i = 0; i < n; i++)
              {
            for (var j = 0; j < n; j++)
            {
              matrix[i, j] = fs[i].ScalarProduct(fs[j]);
            }
            vector[i] = ys.ScalarProduct(fs[i]);
              }

              var matrixInverse = matrix.Inverse();

              var result = matrixInverse * vector;

              for (var i = 0; i < n; i++)
              {
            _f.LinearParameters[i].Value = result[i];
              }
        }

 public FedKF(KF[] filters, DenseMatrix dcm, DenseMatrix covariances)
 {
     _filters = filters;
     _filteredSignals = new DenseMatrix(_filters.Length, 1);
     _dcm = dcm;
     _dcmt = dcm.Transpose();
     _cInv = covariances.Inverse();
 }

        public void Optimize(Dictionary<double, double> values)
        {
            var ys = new List<double>();
              var fs = new List<double>();
              var gs = new List<double>();
              var hs = new List<double>();

              foreach (var value in values)
              {
            ys.Add(value.Value);
            fs.Add(_lppl.F(value.Key));
            gs.Add(_lppl.G(value.Key));
            hs.Add(_lppl.H(value.Key));
              }

              var N = values.Count;
              var sfs = fs.Sum();
              var sgs = gs.Sum();
              var shs = hs.Sum();
              var sfsfs = fs.ScalarProduct(fs);
              var sgsgs = gs.ScalarProduct(gs);
              var shshs = hs.ScalarProduct(hs);
              var sfsgs = fs.ScalarProduct(gs);
              var sfshs = fs.ScalarProduct(hs);
              var sgshs = gs.ScalarProduct(hs);

              var sys = ys.Sum();
              var sysfs = ys.ScalarProduct(fs);
              var sysgs = ys.ScalarProduct(gs);
              var syshs = ys.ScalarProduct(hs);

              var matrixArray = new[]
                          {
                            N, sfs, sgs, shs,
                            sfs, sfsfs, sfsgs, sfshs,
                            sgs, sfsgs, sgsgs, sgshs,
                            shs, sfshs, sgshs, shshs
                          };

              var matrix = new DenseMatrix(4, 4, matrixArray);

              var matrixInverse = matrix.Inverse();

              var vector = new DenseVector(new[] { sys, sysfs, sysgs, syshs });

              var result = matrixInverse * vector;
              _lppl.A = result[0];
              _lppl.B = result[1];
              _lppl.C1 = result[2];
              _lppl.C2 = result[3];
        }

        public static Matrix<double> WorldToImagePoints(this Matrix<double> worldPoints, DenseMatrix cameraCalibration, DenseVector posePosition, DenseMatrix poseRotation)
        {
            return
            cameraCalibration.Multiply(
                        poseRotation.Inverse()
                            .Multiply(worldPoints)
                            .Translate(-posePosition)
                            .ProjectCamCenteredWorldToHomogImagePoints());

            /*This version is consistent with my DPixelToWorld but different from openCV ProjectPoints
             * cameraCalibration.Multiply(
                poseRotation.Inverse()
                .Multiply(worldPoints.Translate(-posePosition))
                .ProjectCamCenteredWorldToHomogImagePoints());*/
        }

        public void Optimize(Dictionary<double, double> values)
        {
            var ys = values.Select(v => v.Value).ToList();
              var fs = _f.Functions.Select(f => values.Select(v => NaNToZero(f(v.Key))).ToList()).ToList();

              var matrix = new DenseMatrix(_f.LinearCount, _f.LinearCount);
              var vector = new DenseVector(_f.LinearCount);
              for (var i = 0; i < _f.LinearCount; i++)
              {
            for (var j = 0; j < _f.LinearCount; j++)
            {
              matrix[i, j] = fs[i].ScalarProduct(fs[j]);
            }
            vector[i] = ys.ScalarProduct(fs[i]);
              }

              var matrixInverse = matrix.Inverse();
              var result = matrixInverse * vector;

              for (var i = 0; i < _f.LinearCount; i++)
              {
            _f.SetA(i, result[i]);
              }
        }

        /// <summary>Optimizes the specified data</summary>
        /// <param name="data">data</param>
        /// <param name="W">W</param>
        /// <param name="H">H</param>
        protected virtual void Optimize(IBooleanMatrix data, Matrix<float> W, Matrix<float> H)
        {
            var HH          = new Matrix<double>(num_factors, num_factors);
            var HC_minus_IH = new Matrix<double>(num_factors, num_factors);
            var HCp         = new double[num_factors];

            var m = new DenseMatrix(num_factors, num_factors);

            // source code comments are in terms of computing the user factors
            // works the same with users and items exchanged

            // (1) create HH in O(f^2|Items|)
            // HH is symmetric
            for (int f_1 = 0; f_1 < num_factors; f_1++)
                for (int f_2 = 0; f_2 < num_factors; f_2++)
                {
                    double d = 0;
                    for (int i = 0; i < H.dim1; i++)
                        d += H[i, f_1] * H[i, f_2];
                    HH[f_1, f_2] = d;
                }
            // (2) optimize all U
            // HC_minus_IH is symmetric
            for (int u = 0; u < W.dim1; u++)
            {
                var row = data.GetEntriesByRow(u);
                // create HC_minus_IH in O(f^2|S_u|)
                for (int f_1 = 0; f_1 < num_factors; f_1++)
                    for (int f_2 = 0; f_2 < num_factors; f_2++)
                    {
                        double d = 0;
                        foreach (int i in row)
                            //d += H[i, f_1] * H[i, f_2] * (c_pos - 1);
                            d += H[i, f_1] * H[i, f_2] * c_pos;
                        HC_minus_IH[f_1, f_2] = d;
                    }
                // create HCp in O(f|S_u|)
                for (int f = 0; f < num_factors; f++)
                {
                    double d = 0;
                    foreach (int i in row)
                        //d += H[i, f] * c_pos;
                        d += H[i, f] * (1 + c_pos);
                    HCp[f] = d;
                }
                // create m = HH + HC_minus_IH + reg*I
                // m is symmetric
                // the inverse m_inv is symmetric
                for (int f_1 = 0; f_1 < num_factors; f_1++)
                    for (int f_2 = 0; f_2 < num_factors; f_2++)
                    {
                        double d = HH[f_1, f_2] + HC_minus_IH[f_1, f_2];
                        if (f_1 == f_2)
                            d += regularization;
                        m[f_1, f_2] = d;
                    }
                var m_inv = m.Inverse();
                // write back optimal W
                for (int f = 0; f < num_factors; f++)
                {
                    double d = 0;
                    for (int f_2 = 0; f_2 < num_factors; f_2++)
                        d += m_inv[f, f_2] * HCp[f_2];
                    W[u, f] = (float) d;
                }
            }
        }

        private void optimize(DenseMatrix coefficients, DenseVector objFunValues, bool artifical)
        {
            //for calculations on the optimal solution row
            int cCounter,
                width = coefficients.ColumnCount;
            DenseVector cBVect = new DenseVector(basics.Count);

            //Sets up the b matrix
            DenseMatrix b = new DenseMatrix(basics.Count, 1);

            //basics will have values greater than coefficients.ColumnCount - 1 if there are still artificial variables
            //or if Nathan is bad and didn't get rid of them correctly
            foreach (int index in basics)
            {
                b = (DenseMatrix)b.Append(DenseVector.OfVector(coefficients.Column(index)).ToColumnMatrix());
            }
            // removes the first column
            b = (DenseMatrix)b.SubMatrix(0, b.RowCount, 1, b.ColumnCount - 1);

            double[] cPrimes = new double[width];
            double[] rhsOverPPrime;
            DenseMatrix[] pPrimes = new DenseMatrix[width];
            DenseMatrix bInverse;

            int newEntering, exitingRow;

            bool optimal = false;

            if(artifical)
            {
                rhsOverPPrime = new double[numConstraints + 1];
            }
            else
            {
                rhsOverPPrime = new double[numConstraints];
            }

            while (!optimal)
            {
                //calculates the inverse of b for this iteration
                bInverse = (DenseMatrix)b.Inverse();

                //updates the C vector with the most recent basic variables
                cCounter = 0;
                foreach (int index in basics)
                {
                    cBVect[cCounter++] = objFunValues.At(index);
                }

                //calculates the pPrimes and cPrimes
                for (int i = 0; i < coefficients.ColumnCount; i++)
                {
                    if (!basics.Contains(i))
                    {
                        pPrimes[i] = (DenseMatrix)bInverse.Multiply((DenseMatrix)coefficients.Column(i).ToColumnMatrix());

                        //c' = objFunVals - cB * P'n
                        //At(0) to turn it into a double
                        cPrimes[i] = objFunValues.At(i) - (pPrimes[i].LeftMultiply(cBVect)).At(0);
                    }
                    else
                    {
                        pPrimes[i] = null;
                    }
                }

                //RHS'
                xPrime = (DenseMatrix)bInverse.Multiply((DenseMatrix)rhsValues.ToColumnMatrix());

                //Starts newEntering as the first nonbasic
                newEntering = -1;
                int iter = 0;
                while(newEntering == -1)
                {
                    if(!basics.Contains(iter))
                    {
                        newEntering = iter;
                    }

                    iter++;
                }

                //new entering becomes the small cPrime that corresponds to a non-basic value
                for (int i = 0; i < cPrimes.Length; i++)
                {
                    if (cPrimes[i] < cPrimes[newEntering] && !basics.Contains(i))
                    {
                        newEntering = i;
                    }
                }

                //if the smallest cPrime is >= 0, ie they are all positive
                if (cPrimes[newEntering] >= 0)
                {
                    optimal = true;
                }
                else
                {
                    //fix me to deal with if all these values are negative
                    exitingRow = 0;
//.........这里部分代码省略.........

        private void CalcBaseAandC()
        {
            BaseA = new DenseMatrix(CurrBaseJ.Count);
            BaseC = new DenseVector(CurrBaseJ.Count);

            int i = 0;
            foreach (var j in CurrBaseJ)
            {
                BaseA.SetColumn(i, A.Column(j));
                BaseC[i] = c[j];
                i++;
            }

            InvBaseA = (DenseMatrix)BaseA.Inverse();

            CurrNonBaseJ = new List<int>(Enumerable.Range(0, A.ColumnCount).Except(CurrBaseJ));
        }

        public void placeCompensatorPoles(double value)
        {
            // Check for dimensions
            if (A.ColumnCount != A.RowCount || A.ColumnCount <= 0) {
                // TODO: Throw exception
            }

            int n = A.ColumnCount;

            // Calculate controllability matrix
            Matrix<double> controllabilityMatrix = new DenseMatrix(n, n);
            for (int i = 0; i < n; i++) {
                Vector<double> vec = B.Column(0);

                for (int j = 0; j < i; j++) {
                    vec = A * vec;
                }

                controllabilityMatrix.SetColumn(i, vec);
            }

            // Unity vector
            Matrix<double> unityVector = new DenseMatrix(1,n);
            for (int i = 0; i < n; i++) {
                // Set 1 at last index
                unityVector.At(0, i,
                    (i == n-1 ? 1 : 0)
                    );
            }

            // Coefficients matrix
            Matrix<double> preparedMatrix = A.Clone();
            for (int i = 0; i < n; i++) {
                // Substract value from diagonal
                preparedMatrix.At(i, i,
                    preparedMatrix.At(i, i) - value
                    );
            }
            Matrix<double> coefficientsMatrix = preparedMatrix.Clone();
            for (int i = 0; i < n-1; i++) {
                // Multiply n-1 times
                coefficientsMatrix = preparedMatrix * coefficientsMatrix;
            }

            // Calculate new K using Ackermann's formula
            K = unityVector * controllabilityMatrix.Inverse() * coefficientsMatrix;
        }

        /// <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();
        }

        private double[] linearApproxym(double[] x, double[] y)
        {
            int n = x.Length;
            double sumOfX = 0;
            double sumOfXSquare = 0;
            double sumOfY = 0;
            double sumOfXY = 0;

            for (int i = 0; i < n; i++)
            {
                sumOfX += x[i];
                sumOfXSquare += x[i] * x[i];
                sumOfY += y[i];
                sumOfXY += x[i] * y[i];
            }

            var matrixX = new DenseMatrix(new double[,] { { n, sumOfX }, { sumOfX, sumOfXSquare } });
            var matrixY = new DenseMatrix(2, 1, new double[] { sumOfY, sumOfXY });

            var matrixP = matrixX.Inverse().Multiply(matrixY);

            double[] p = new double[2];

            p[0] = matrixP[1, 0];
            p[1] = matrixP[0, 0];

            return p;
        }

        private void CalculateCoefficients(int points)
        {
            var c = new double[points][,];

            // For ever possible center given the number of points, compute ever possible coefficient for all possible orders.
            for (int center = 0; center < points; center++)
            {
                // Deltas matrix for center located at 'center'.
                var A = new DenseMatrix(points);
                var l = points - center - 1;
                for (int row = points - 1; row >= 0; row--)
                {
                    A[row, 0] = 1.0;
                    for (int col = 1; col < points; col++)
                    {
                        A[row, col] = A[row, col - 1] * l / col;
                    }
                    l -= 1;
                }

                c[center] = A.Inverse().ToArray();

                // "Polish" results by rounding.
                var fac = SpecialFunctions.Factorial(points);
                for (int j = 0; j < points; j++)
                    for (int k = 0; k < points; k++)
#if PORTABLE
                        c[center][j, k] = (Math.Round(c[center][j, k] * fac)) / fac;
#else
                        c[center][j, k] = (Math.Round(c[center][j, k] * fac, MidpointRounding.AwayFromZero)) / fac;
#endif
            }

            _coefficients = c;
        }

        //Given an observation, mu and sigma. What is N(O | mu, sigma)?
        public double queryGuassian(DenseVector observation, DenseVector mean, DenseMatrix covar)
        {
            double scaletemp, scale, exponent, prob;
            DenseMatrix v1, v2; //Temp matrices, for multiplying

            scaletemp = (Math.Pow(2 * Math.PI, mean.Count / 2));
            scale = (Math.Sqrt(covar.Determinant()));
            scale *= scaletemp;
            scale = 1 / scale;

            v1 = (DenseMatrix)(observation - mean).ToRowMatrix();
            v2 = (DenseMatrix)(observation - mean).ToColumnMatrix();
            v2 = (DenseMatrix)(covar.Inverse()) * v2;

            exponent = (-0.5) *  ((v1 * v2).ToArray()[0,0]);
            prob = scale * Math.Pow(Math.E, exponent);

            return prob;
        }

		private void Optimize(int u, IBooleanMatrix data, Matrix<float> W, Matrix<float> H, Matrix<double> HH)
		{
			var row = data.GetEntriesByRow(u);
			// HC_minus_IH is symmetric
			// create HC_minus_IH in O(f^2|S_u|)
			var HC_minus_IH = new Matrix<double>(num_factors, num_factors);
			for (int f_1 = 0; f_1 < num_factors; f_1++)
				for (int f_2 = f_1; f_2 < num_factors; f_2++)
				{
					double d = 0;
					foreach (int i in row)
						d += H[i, f_1] * H[i, f_2];
					HC_minus_IH[f_1, f_2] = d * alpha;
					HC_minus_IH[f_2, f_1] = d * alpha;
				}
			// create HCp in O(f|S_u|)
			var HCp = new double[num_factors];
			for (int f = 0; f < num_factors; f++)
			{
				double d = 0;
				foreach (int i in row)
					d += H[i, f];
				HCp[f] = d * (1 + alpha);
			}
			// create m = HH + HC_minus_IH + reg*I
			// m is symmetric
			// the inverse m_inv is symmetric
			var m = new DenseMatrix(num_factors, num_factors);
			for (int f_1 = 0; f_1 < num_factors; f_1++)
				for (int f_2 = f_1; f_2 < num_factors; f_2++)
				{
					double d = HH[f_1, f_2] + HC_minus_IH[f_1, f_2];
					if (f_1 == f_2)
						d += regularization;
					m[f_1, f_2] = d;
					m[f_2, f_1] = d;
				}
			var m_inv = m.Inverse();
			// write back optimal W
			for (int f = 0; f < num_factors; f++)
			{
				double d = 0;
				for (int f_2 = 0; f_2 < num_factors; f_2++)
					d += m_inv[f, f_2] * HCp[f_2];
				W[u, f] = (float) d;
			}
		}

 public Vector<double> Solve(DenseMatrix systemMatrix, Vector<double> freeTerms)
 {
     return systemMatrix.Inverse().Multiply(freeTerms);
 }