diff options
| -rw-r--r-- | Eigen/src/SVD/SVD.h | 33 | ||||
| -rw-r--r-- | test/svd.cpp | 32 |
2 files changed, 59 insertions, 6 deletions
diff --git a/Eigen/src/SVD/SVD.h b/Eigen/src/SVD/SVD.h index 4705855e0..9d269169e 100644 --- a/Eigen/src/SVD/SVD.h +++ b/Eigen/src/SVD/SVD.h @@ -79,6 +79,9 @@ template<typename MatrixType> class SVD void compute(const MatrixType& matrix); SVD& sort(); + void computeUnitaryPositive(MatrixUType *unitary, MatrixType *positive) const; + void computePositiveUnitary(MatrixType *positive, MatrixVType *unitary) const; + protected: /** \internal */ MatrixUType m_matU; @@ -534,6 +537,36 @@ bool SVD<MatrixType>::solve(const MatrixBase<OtherDerived> &b, ResultType* resul return true; } +/** Computes the polar decomposition of the matrix, as a product unitary x positive. + * + * If either pointer is zero, the corresponding computation is skipped. + * + * Only for square matrices. + */ +template<typename MatrixType> +void SVD<MatrixType>::computeUnitaryPositive(typename SVD<MatrixType>::MatrixUType *unitary, + MatrixType *positive) const +{ + ei_assert(m_matU.cols() == m_matV.cols() && "Polar decomposition is only for square matrices"); + if(unitary) *unitary = m_matU * m_matV.adjoint(); + if(positive) *positive = m_matV * m_sigma.asDiagonal() * m_matV.adjoint(); +} + +/** Computes the polar decomposition of the matrix, as a product positive x unitary. + * + * If either pointer is zero, the corresponding computation is skipped. + * + * Only for square matrices. + */ +template<typename MatrixType> +void SVD<MatrixType>::computePositiveUnitary(MatrixType *positive, + typename SVD<MatrixType>::MatrixVType *unitary) const +{ + ei_assert(m_matU.rows() == m_matV.rows() && "Polar decomposition is only for square matrices"); + if(unitary) *unitary = m_matU * m_matV.adjoint(); + if(positive) *positive = m_matU * m_sigma.asDiagonal() * m_matU.adjoint(); +} + /** \svd_module * \returns the SVD decomposition of \c *this */ diff --git a/test/svd.cpp b/test/svd.cpp index 4776ef7e7..688c3f402 100644 --- a/test/svd.cpp +++ b/test/svd.cpp @@ -44,13 +44,15 @@ template<typename MatrixType> void svd(const MatrixType& m) if (ei_is_same_type<RealScalar,float>::ret) largerEps = 1e-3f; - SVD<MatrixType> svd(a); - MatrixType sigma = MatrixType::Zero(rows,cols); - MatrixType matU = MatrixType::Zero(rows,rows); - sigma.block(0,0,cols,cols) = svd.singularValues().asDiagonal(); - matU.block(0,0,rows,cols) = svd.matrixU(); + { + SVD<MatrixType> svd(a); + MatrixType sigma = MatrixType::Zero(rows,cols); + MatrixType matU = MatrixType::Zero(rows,rows); + sigma.block(0,0,cols,cols) = svd.singularValues().asDiagonal(); + matU.block(0,0,rows,cols) = svd.matrixU(); + VERIFY_IS_APPROX(a, matU * sigma * svd.matrixV().transpose()); + } - VERIFY_IS_APPROX(a, matU * sigma * svd.matrixV().transpose()); if (rows==cols) { @@ -63,6 +65,24 @@ template<typename MatrixType> void svd(const MatrixType& m) svd.solve(b, &x); VERIFY_IS_APPROX(a * x,b); } + + + if(rows==cols) + { + SVD<MatrixType> svd(a); + MatrixType unitary, positive; + svd.computeUnitaryPositive(&unitary, &positive); + VERIFY_IS_APPROX(unitary * unitary.adjoint(), MatrixType::Identity(unitary.rows(),unitary.rows())); + VERIFY_IS_APPROX(positive, positive.adjoint()); + for(int i = 0; i < rows; i++) VERIFY(positive.diagonal()[i] >= 0); // cheap necessary (not sufficient) condition for positivity + VERIFY_IS_APPROX(unitary*positive, a); + + svd.computePositiveUnitary(&positive, &unitary); + VERIFY_IS_APPROX(unitary * unitary.adjoint(), MatrixType::Identity(unitary.rows(),unitary.rows())); + VERIFY_IS_APPROX(positive, positive.adjoint()); + for(int i = 0; i < rows; i++) VERIFY(positive.diagonal()[i] >= 0); // cheap necessary (not sufficient) condition for positivity + VERIFY_IS_APPROX(positive*unitary, a); + } } void test_svd() |