C# Matrix.Svd方法代码示例

        public Molecule(double[,] inputAtoms, bool getDirections)
        {
            numAtoms = inputAtoms.GetLength(0);
            atoms = DenseMatrix.OfArray(inputAtoms);
            centroid = atoms.ColumnSums() / (double)numAtoms;

            if (getDirections)
            {
                for (int i = 0; i < 3; i++)
                {
                    atoms.SetColumn(i, atoms.Column(i) - centroid[i]);
                }
                direction = atoms.Svd(true).VT.Row(0);
                if ((atoms.Row(atoms.RowCount - 1) - atoms.Row(0)).DotProduct(direction) < 0)
                {
                    direction *= -1;
                }
            }
        }

        public static Vector<float> Solve(Matrix<float> A, Vector<float> b)
        {
            MathNet.Numerics.LinearAlgebra.Factorization.Svd<float> svd = A.Svd();

            double dmax = svd.S[0];
            Vector<float> y = new MathNet.Numerics.LinearAlgebra.Single.DenseVector(A.ColumnCount);
            Vector<float> bp = svd.U.Transpose() * b;

            // TODO: use the fact that S is sorted ( while S[i]/dmax < e then assume its greater )
            int minSize = Math.Min(A.ColumnCount, A.RowCount);
            for(int i = 0; i < minSize; ++i)
            {
                if(svd.S[i] / dmax < 1e-16)
                {
                    y[i] = 0.0f;
                }
                else
                {
                    y[i] = bp[i] / svd.S[i];
                }
            }

            return svd.VT.Transpose() * y;
        }

        public void CreateEpiGeometry_Normalized()
        {
            var ChL = new DenseVector(4);
            _C_L.CopySubVectorTo(ChL, 0, 0, 3);
            ChL[3] = 1.0;
            var ChR = new DenseVector(4);
            _C_R.CopySubVectorTo(ChR, 0, 0, 3);
            ChR[3] = 1.0;

            // Find e_R = P_R*C_L, e_L = P_L*C_R
            _epi_R = _CMN_R * (_Nr * ChL);
            _epi_L = _CMN_L * (_Nr * ChR);
            _epi_R.DivideThis(_epi_R[2]);
            _epi_L.DivideThis(_epi_L[2]);

            _epiX_L = new DenseMatrix(3, 3);
            _epiX_L[0, 0] = 0.0;
            _epiX_L[1, 0] = _epi_L[2];
            _epiX_L[2, 0] = -_epi_L[1];
            _epiX_L[0, 1] = -_epi_L[2];
            _epiX_L[1, 1] = 0.0;
            _epiX_L[2, 1] = _epi_L[0];
            _epiX_L[0, 2] = _epi_L[1];
            _epiX_L[1, 2] = -_epi_L[0];
            _epiX_L[2, 2] = 0.0;

            _epiX_R = new DenseMatrix(3, 3);
            _epiX_R[0, 0] = 0.0;
            _epiX_R[1, 0] = _epi_R[2];
            _epiX_R[2, 0] = -_epi_R[1];
            _epiX_R[0, 1] = -_epi_R[2];
            _epiX_R[1, 1] = 0.0;
            _epiX_R[2, 1] = _epi_R[0];
            _epiX_R[0, 2] = _epi_R[1];
            _epiX_R[1, 2] = -_epi_R[0];
            _epiX_R[2, 2] = 0.0;

            // F = [er]x * Pr * pseudoinv(Pl)
            _F = _epiX_R * (_CMN_R * _CMN_L.PseudoInverse());
            int rank = _F.Rank();
            if(rank == 3)
            {
                // Need to ensure rank 2, so set smallest singular value to 0
                var svd = _F.Svd();
                var E = svd.W;
                E[2, 2] = 0;
                var oldF = _F;
                _F = svd.U * E * svd.VT;
                var diff = _F - oldF; // Difference should be very small if all is correct
            }

            // Scale F, so that F33 = 1
            _F = _F.Divide(_F[2, 2]);
        }

        private void ComputeEssentialFundamental()
        {
            if(!IsCamLeftCalibrated || !IsCamRightCalibrated)
                return;

            // Find e_R = P_R*C_L, e_L = P_L*C_R
            _epiPoleRight = _camRight * new DenseVector(new double[] { _translationLeft.At(0), _translationLeft.At(1), _translationLeft.At(2), 1.0 });
            _epiPoleLeft = _camLeft * new DenseVector(new double[] { _translationRight.At(0), _translationRight.At(1), _translationRight.At(2), 1.0 });

            _epiLeftInInfinity = false;
            _epiRightInInfinity = false;

            // Check if any epipole is in infinity:
            if(Math.Abs(_epiPoleLeft.At(2)) < 1e-9)
            {
                _epiLeftInInfinity = true;
                // Normalize epipole not by dividing by w, but so |e| = 1
                _epiPoleLeft = _epiPoleLeft.Normalize(2);
            }
            if(Math.Abs(_epiPoleRight.At(2)) < 1e-9)
            {
                _epiRightInInfinity = true;
                // Normalize epipole not by dividing by w, but so |e| = 1
                _epiPoleRight = _epiPoleRight.Normalize(2);
            }

            //             |  0 -e3  e2|
            // Find [e]x = | e3   0 -e1|
            //             |-e2  e1   0|

            _epiCrossRight = new DenseMatrix(3, 3);
            _epiCrossRight[0, 0] = 0.0;
            _epiCrossRight[1, 0] = _epiPoleRight[2];
            _epiCrossRight[2, 0] = -_epiPoleRight[1];
            _epiCrossRight[0, 1] = -_epiPoleRight[2];
            _epiCrossRight[1, 1] = 0.0;
            _epiCrossRight[2, 1] = _epiPoleRight[0];
            _epiCrossRight[0, 2] = _epiPoleRight[1];
            _epiCrossRight[1, 2] = -_epiPoleRight[0];
            _epiCrossRight[2, 2] = 0.0;

            _epiCrossLeft = new DenseMatrix(3, 3);
            _epiCrossLeft[0, 0] = 0.0;
            _epiCrossLeft[1, 0] = _epiPoleLeft[2];
            _epiCrossLeft[2, 0] = -_epiPoleLeft[1];
            _epiCrossLeft[0, 1] = -_epiPoleLeft[2];
            _epiCrossLeft[1, 1] = 0.0;
            _epiCrossLeft[2, 1] = _epiPoleLeft[0];
            _epiCrossLeft[0, 2] = _epiPoleLeft[1];
            _epiCrossLeft[1, 2] = -_epiPoleLeft[0];
            _epiCrossLeft[2, 2] = 0.0;

            // F = [er]x * Pr * pseudoinv(Pl)
            _fundamental = _epiCrossRight * _camRight * (_camLeft.PseudoInverse());
            int rank = _fundamental.Rank();
            if(rank == 3)
            {
                // Need to ensure rank 2, so set smallest singular value to 0
                var svd = _fundamental.Svd();
                var E = svd.W;
                E[2, 2] = 0;
                var oldF = _fundamental;
                _fundamental = svd.U * E * svd.VT;
                var diff = _fundamental - oldF; // Difference should be very small if all is correct

                // NOPE
                // Get new SVD -> last singular value should be 0
                // Vectors corresponding to this value are right/left epipoles (change them to ensure e'^TFe = 0)
               // var evd = _fundamental.Evd();
               // svd = _fundamental.Svd();
            }

            // Scale F, so that F33 = 1
            _fundamental = _fundamental.Divide(_fundamental[2, 2]);

            // E = Kr^T F Kl
            _essential = _calibrationRight.Transpose() * _fundamental * _calibrationLeft;
        }

//.........这里部分代码省略.........

            if (sampleWeight != null)
            {
                solver = RidgeSolver.DenseCholesky;
            }

            if (solver == RidgeSolver.Lsqr)
            {
                // According to the lsqr documentation, alpha = damp^2.
                double sqrtAlpha = Math.Sqrt(alpha);
                Matrix coefs = new DenseMatrix(y.ColumnCount, x.ColumnCount);
                foreach (var column in y.ColumnEnumerator())
                {
                    Vector<double> c = Lsqr.lsqr(
                        x,
                        column.Item2,
                        damp: sqrtAlpha,
                        atol: tol,
                        btol: tol,
                        iterLim: maxIter).X;

                    coefs.SetRow(column.Item1, c);
                }

                return coefs;
            }

            if (solver == RidgeSolver.DenseCholesky)
            {
                //# normal equations (cholesky) method
                if (nFeatures > nSamples || sampleWeight != null)
                {
                    // kernel ridge
                    // w = X.T * inv(X X^t + alpha*Id) y
                    var k = x.TransposeAndMultiply(x);
                    Vector<double> sw = null;
                    if (sampleWeight != null)
                    {
                        sw = sampleWeight.Sqrt();
                        // We are doing a little danse with the sample weights to
                        // avoid copying the original X, which could be big

                        y = y.MulColumnVector(sw);

                        k = k.PointwiseMultiply(sw.Outer(sw));
                    }

                    k.Add(DenseMatrix.Identity(k.RowCount)*alpha, k);
                    try
                    {
                        var dualCoef = k.Cholesky().Solve(y);
                        if (sampleWeight != null)
                            dualCoef = dualCoef.MulColumnVector(sw);

                        return x.TransposeThisAndMultiply(dualCoef).Transpose();
                    }
                    catch (Exception) //todo:
                    {
                        // use SVD solver if matrix is singular
                        solver = RidgeSolver.Svd;
                    }
                }
                else
                {
                    // ridge
                    // w = inv(X^t X + alpha*Id) * X.T y
                    var a = x.TransposeThisAndMultiply(x);
                    a.Add(DenseMatrix.Identity(a.ColumnCount)*alpha, a);

                    var xy = x.TransposeThisAndMultiply(y);

                    try
                    {
                        return a.Cholesky().Solve(xy).Transpose();
                    }
                    catch (Exception) //todo:
                    {
                        // use SVD solver if matrix is singular
                        solver = RidgeSolver.Svd;
                    }
                }
            }

            if (solver == RidgeSolver.Svd)
            {
                // slower than cholesky but does not break with
                // singular matrices
                var svd = x.Svd(true);
                //U, s, Vt = linalg.svd(X, full_matrices=False)
                int k = Math.Min(x.ColumnCount, x.RowCount);
                var d = svd.S().SubVector(0, k);
                d.MapInplace(v => v > 1e-15 ? v/(v*v + alpha) : 0.0);

                var ud = svd.U().SubMatrix(0, x.RowCount, 0, k).TransposeThisAndMultiply(y).Transpose();
                ud = ud.MulRowVector(d);
                return ud.Multiply(svd.VT().SubMatrix(0, k, 0, x.ColumnCount));
            }

            return null;
        }