diff --git a/Eigen/src/QR/ColPivHouseholderQR.h b/Eigen/src/QR/ColPivHouseholderQR.h index f1de6fd0b..bd0849d58 100644 --- a/Eigen/src/QR/ColPivHouseholderQR.h +++ b/Eigen/src/QR/ColPivHouseholderQR.h @@ -238,6 +238,20 @@ class ColPivHouseholderQR : public SolverBase typename MatrixType::Scalar ColPivHouseholderQR::determinant() const { + using Scalar = typename MatrixType::Scalar; eigen_assert(m_isInitialized && "HouseholderQR is not initialized."); eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!"); Scalar detQ; internal::householder_determinant::IsComplex>::run(m_hCoeffs, detQ); - return m_qr.diagonal().prod() * detQ * Scalar(m_det_p); + return isInjective() ? (detQ * Scalar(m_det_p)) * m_qr.diagonal().prod() : Scalar(0); } template typename MatrixType::RealScalar ColPivHouseholderQR::absDeterminant() const { using std::abs; + using RealScalar = typename MatrixType::RealScalar; eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!"); - return abs(m_qr.diagonal().prod()); + return isInjective() ? abs(m_qr.diagonal().prod()) : RealScalar(0); } template typename MatrixType::RealScalar ColPivHouseholderQR::logAbsDeterminant() const { + using RealScalar = typename MatrixType::RealScalar; eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!"); - return m_qr.diagonal().cwiseAbs().array().log().sum(); + return isInjective() ? m_qr.diagonal().cwiseAbs().array().log().sum() : -NumTraits::infinity(); +} + +template +typename MatrixType::Scalar ColPivHouseholderQR::signDeterminant() const { + using Scalar = typename MatrixType::Scalar; + eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); + eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!"); + Scalar detQ; + internal::householder_determinant::IsComplex>::run(m_hCoeffs, detQ); + return isInjective() ? (detQ * Scalar(m_det_p)) * m_qr.diagonal().array().sign().prod() : Scalar(0); } /** Performs the QR factorization of the given matrix \a matrix. The result of diff --git a/Eigen/src/QR/CompleteOrthogonalDecomposition.h b/Eigen/src/QR/CompleteOrthogonalDecomposition.h index 8566e965b..960ccb1e9 100644 --- a/Eigen/src/QR/CompleteOrthogonalDecomposition.h +++ b/Eigen/src/QR/CompleteOrthogonalDecomposition.h @@ -228,6 +228,21 @@ class CompleteOrthogonalDecomposition */ typename MatrixType::RealScalar logAbsDeterminant() const; + /** \returns the sign of the determinant of the + * matrix of which *this is the complete orthogonal decomposition. It has + * only linear complexity (that is, O(n) where n is the dimension of the + * square matrix) as the complete orthogonal decomposition has already been + * computed. + * + * \note This is only for square matrices. + * + * \note This method is useful to work around the risk of overflow/underflow + * that's inherent to determinant computation. + * + * \sa determinant(), absDeterminant(), logAbsDeterminant(), MatrixBase::determinant() + */ + typename MatrixType::Scalar signDeterminant() const; + /** \returns the rank of the matrix of which *this is the complete orthogonal * decomposition. * @@ -424,6 +439,11 @@ typename MatrixType::RealScalar CompleteOrthogonalDecomposition +typename MatrixType::Scalar CompleteOrthogonalDecomposition::signDeterminant() const { + return m_cpqr.signDeterminant(); +} + /** Performs the complete orthogonal decomposition of the given matrix \a * matrix. The result of the factorization is stored into \c *this, and a * reference to \c *this is returned. diff --git a/Eigen/src/QR/FullPivHouseholderQR.h b/Eigen/src/QR/FullPivHouseholderQR.h index d93a5d174..0cadf2317 100644 --- a/Eigen/src/QR/FullPivHouseholderQR.h +++ b/Eigen/src/QR/FullPivHouseholderQR.h @@ -248,6 +248,20 @@ class FullPivHouseholderQR : public SolverBase typename MatrixType::Scalar FullPivHouseholderQR::determinant() const { + using Scalar = typename MatrixType::Scalar; eigen_assert(m_isInitialized && "HouseholderQR is not initialized."); eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!"); Scalar detQ; internal::householder_determinant::IsComplex>::run(m_hCoeffs, detQ); - return m_qr.diagonal().prod() * detQ * Scalar(m_det_p); + return isInjective() ? (detQ * Scalar(m_det_p)) * m_qr.diagonal().prod() : Scalar(0); } template typename MatrixType::RealScalar FullPivHouseholderQR::absDeterminant() const { + using RealScalar = typename MatrixType::RealScalar; using std::abs; eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!"); - return abs(m_qr.diagonal().prod()); + return isInjective() ? abs(m_qr.diagonal().prod()) : RealScalar(0); } template typename MatrixType::RealScalar FullPivHouseholderQR::logAbsDeterminant() const { + using RealScalar = typename MatrixType::RealScalar; eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!"); - return m_qr.diagonal().cwiseAbs().array().log().sum(); + return isInjective() ? m_qr.diagonal().cwiseAbs().array().log().sum() : -NumTraits::infinity(); +} + +template +typename MatrixType::Scalar FullPivHouseholderQR::signDeterminant() const { + using Scalar = typename MatrixType::Scalar; + eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); + eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!"); + Scalar detQ; + internal::householder_determinant::IsComplex>::run(m_hCoeffs, detQ); + return isInjective() ? (detQ * Scalar(m_det_p)) * m_qr.diagonal().array().sign().prod() : Scalar(0); } /** Performs the QR factorization of the given matrix \a matrix. The result of diff --git a/Eigen/src/QR/HouseholderQR.h b/Eigen/src/QR/HouseholderQR.h index 9e736722a..e29737258 100644 --- a/Eigen/src/QR/HouseholderQR.h +++ b/Eigen/src/QR/HouseholderQR.h @@ -187,6 +187,8 @@ class HouseholderQR : public SolverBase> { * \warning a determinant can be very big or small, so for matrices * of large enough dimension, there is a risk of overflow/underflow. * One way to work around that is to use logAbsDeterminant() instead. + * Also, do not rely on the determinant being exactly zero for testing + * singularity or rank-deficiency. * * \sa absDeterminant(), logAbsDeterminant(), MatrixBase::determinant() */ @@ -202,6 +204,8 @@ class HouseholderQR : public SolverBase> { * \warning a determinant can be very big or small, so for matrices * of large enough dimension, there is a risk of overflow/underflow. * One way to work around that is to use logAbsDeterminant() instead. + * Also, do not rely on the determinant being exactly zero for testing + * singularity or rank-deficiency. * * \sa determinant(), logAbsDeterminant(), MatrixBase::determinant() */ @@ -217,10 +221,30 @@ class HouseholderQR : public SolverBase> { * \note This method is useful to work around the risk of overflow/underflow that's inherent * to determinant computation. * + * \warning Do not rely on the determinant being exactly zero for testing + * singularity or rank-deficiency. + * * \sa determinant(), absDeterminant(), MatrixBase::determinant() */ typename MatrixType::RealScalar logAbsDeterminant() const; + /** \returns the sign of the determinant of the matrix of which + * *this is the QR decomposition. It has only linear complexity + * (that is, O(n) where n is the dimension of the square matrix) + * as the QR decomposition has already been computed. + * + * \note This is only for square matrices. + * + * \note This method is useful to work around the risk of overflow/underflow that's inherent + * to determinant computation. + * + * \warning Do not rely on the determinant being exactly zero for testing + * singularity or rank-deficiency. + * + * \sa determinant(), absDeterminant(), MatrixBase::determinant() + */ + typename MatrixType::Scalar signDeterminant() const; + inline Index rows() const { return m_qr.rows(); } inline Index cols() const { return m_qr.cols(); } @@ -306,6 +330,15 @@ typename MatrixType::RealScalar HouseholderQR::logAbsDeterminant() c return m_qr.diagonal().cwiseAbs().array().log().sum(); } +template +typename MatrixType::Scalar HouseholderQR::signDeterminant() const { + eigen_assert(m_isInitialized && "HouseholderQR is not initialized."); + eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!"); + Scalar detQ; + internal::householder_determinant::IsComplex>::run(m_hCoeffs, detQ); + return detQ * m_qr.diagonal().array().sign().prod(); +} + namespace internal { /** \internal */ diff --git a/test/qr.cpp b/test/qr.cpp index de470cae8..f7f6990d9 100644 --- a/test/qr.cpp +++ b/test/qr.cpp @@ -82,6 +82,7 @@ void qr_invertible() { m1 = m3 * m1 * m3.adjoint(); qr.compute(m1); VERIFY_IS_APPROX(log(absdet), qr.logAbsDeterminant()); + VERIFY_IS_APPROX(numext::sign(det), qr.signDeterminant()); // This test is tricky if the determinant becomes too small. // Since we generate random numbers with magnitude range [0,1], the average determinant is 0.5^size RealScalar tol = @@ -102,7 +103,7 @@ void qr_verify_assert() { VERIFY_RAISES_ASSERT(qr.householderQ()) VERIFY_RAISES_ASSERT(qr.determinant()) VERIFY_RAISES_ASSERT(qr.absDeterminant()) - VERIFY_RAISES_ASSERT(qr.logAbsDeterminant()) + VERIFY_RAISES_ASSERT(qr.signDeterminant()) } EIGEN_DECLARE_TEST(qr) { diff --git a/test/qr_colpivoting.cpp b/test/qr_colpivoting.cpp index 4f8711fad..705573c0a 100644 --- a/test/qr_colpivoting.cpp +++ b/test/qr_colpivoting.cpp @@ -21,6 +21,7 @@ void cod() { Index rank = internal::random(1, (std::min)(rows, cols) - 1); typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar RealScalar; typedef Matrix MatrixQType; MatrixType matrix; createRandomPIMatrixOfRank(rank, rows, cols, matrix); @@ -56,6 +57,23 @@ void cod() { MatrixType pinv = cod.pseudoInverse(); VERIFY_IS_APPROX(cod_solution, pinv * rhs); + + // now construct a (square) matrix with prescribed determinant + Index size = internal::random(2, 20); + matrix.setZero(size, size); + for (int i = 0; i < size; i++) { + matrix(i, i) = internal::random(); + } + Scalar det = matrix.diagonal().prod(); + RealScalar absdet = abs(det); + CompleteOrthogonalDecomposition cod2(matrix); + cod2.compute(matrix); + q = cod2.householderQ(); + matrix = q * matrix * q.adjoint(); + VERIFY_IS_APPROX(det, cod2.determinant()); + VERIFY_IS_APPROX(absdet, cod2.absDeterminant()); + VERIFY_IS_APPROX(numext::log(absdet), cod2.logAbsDeterminant()); + VERIFY_IS_APPROX(numext::sign(det), cod2.signDeterminant()); } template @@ -265,6 +283,7 @@ void qr_invertible() { VERIFY_IS_APPROX(det, qr.determinant()); VERIFY_IS_APPROX(absdet, qr.absDeterminant()); VERIFY_IS_APPROX(log(absdet), qr.logAbsDeterminant()); + VERIFY_IS_APPROX(numext::sign(det), qr.signDeterminant()); } template @@ -285,6 +304,7 @@ void qr_verify_assert() { VERIFY_RAISES_ASSERT(qr.determinant()) VERIFY_RAISES_ASSERT(qr.absDeterminant()) VERIFY_RAISES_ASSERT(qr.logAbsDeterminant()) + VERIFY_RAISES_ASSERT(qr.signDeterminant()) } template @@ -305,6 +325,7 @@ void cod_verify_assert() { VERIFY_RAISES_ASSERT(cod.determinant()) VERIFY_RAISES_ASSERT(cod.absDeterminant()) VERIFY_RAISES_ASSERT(cod.logAbsDeterminant()) + VERIFY_RAISES_ASSERT(cod.signDeterminant()) } EIGEN_DECLARE_TEST(qr_colpivoting) { diff --git a/test/qr_fullpivoting.cpp b/test/qr_fullpivoting.cpp index 71f3a51a4..2b6ecc5b8 100644 --- a/test/qr_fullpivoting.cpp +++ b/test/qr_fullpivoting.cpp @@ -105,6 +105,7 @@ void qr_invertible() { VERIFY_IS_APPROX(det, qr.determinant()); VERIFY_IS_APPROX(absdet, qr.absDeterminant()); VERIFY_IS_APPROX(log(absdet), qr.logAbsDeterminant()); + VERIFY_IS_APPROX(numext::sign(det), qr.signDeterminant()); } template @@ -125,6 +126,7 @@ void qr_verify_assert() { VERIFY_RAISES_ASSERT(qr.determinant()) VERIFY_RAISES_ASSERT(qr.absDeterminant()) VERIFY_RAISES_ASSERT(qr.logAbsDeterminant()) + VERIFY_RAISES_ASSERT(qr.signDeterminant()) } EIGEN_DECLARE_TEST(qr_fullpivoting) {