From ef4654bae77b7a8d77d2b9bc93ef10ffe0863def Mon Sep 17 00:00:00 2001 From: sjusju <1124wwqq@gmail.com> Date: Fri, 29 Jul 2022 18:24:14 +0000 Subject: [PATCH] Add true determinant to QR and it's variants --- Eigen/src/QR/ColPivHouseholderQR.h | 33 +++++++-- .../src/QR/CompleteOrthogonalDecomposition.h | 25 ++++++- Eigen/src/QR/FullPivHouseholderQR.h | 33 +++++++-- Eigen/src/QR/HouseholderQR.h | 70 ++++++++++++++++++- test/qr.cpp | 11 +-- test/qr_colpivoting.cpp | 8 ++- test/qr_fullpivoting.cpp | 7 +- 7 files changed, 167 insertions(+), 20 deletions(-) diff --git a/Eigen/src/QR/ColPivHouseholderQR.h b/Eigen/src/QR/ColPivHouseholderQR.h index fec53213a..c906997fc 100644 --- a/Eigen/src/QR/ColPivHouseholderQR.h +++ b/Eigen/src/QR/ColPivHouseholderQR.h @@ -219,6 +219,21 @@ template class ColPivHouseholderQR return m_colsPermutation; } + /** \returns 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. + * + * \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. + * + * \sa absDeterminant(), logAbsDeterminant(), MatrixBase::determinant() + */ + typename MatrixType::Scalar determinant() const; + /** \returns the absolute value 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) @@ -230,7 +245,7 @@ template class ColPivHouseholderQR * of large enough dimension, there is a risk of overflow/underflow. * One way to work around that is to use logAbsDeterminant() instead. * - * \sa logAbsDeterminant(), MatrixBase::determinant() + * \sa determinant(), logAbsDeterminant(), MatrixBase::determinant() */ typename MatrixType::RealScalar absDeterminant() const; @@ -244,7 +259,7 @@ template class ColPivHouseholderQR * \note This method is useful to work around the risk of overflow/underflow that's inherent * to determinant computation. * - * \sa absDeterminant(), MatrixBase::determinant() + * \sa determinant(), absDeterminant(), MatrixBase::determinant() */ typename MatrixType::RealScalar logAbsDeterminant() const; @@ -442,9 +457,19 @@ template class ColPivHouseholderQR bool m_isInitialized, m_usePrescribedThreshold; RealScalar m_prescribedThreshold, m_maxpivot; Index m_nonzero_pivots; - Index m_det_pq; + Index m_det_p; }; +template +typename MatrixType::Scalar ColPivHouseholderQR::determinant() 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 m_qr.diagonal().prod() * detQ * Scalar(m_det_p); +} + template typename MatrixType::RealScalar ColPivHouseholderQR::absDeterminant() const { @@ -574,7 +599,7 @@ void ColPivHouseholderQR::computeInPlace() for(PermIndexType k = 0; k < size/*m_nonzero_pivots*/; ++k) m_colsPermutation.applyTranspositionOnTheRight(k, PermIndexType(m_colsTranspositions.coeff(k))); - m_det_pq = (number_of_transpositions%2) ? -1 : 1; + m_det_p = (number_of_transpositions%2) ? -1 : 1; m_isInitialized = true; } diff --git a/Eigen/src/QR/CompleteOrthogonalDecomposition.h b/Eigen/src/QR/CompleteOrthogonalDecomposition.h index bfaaa93ef..02583a2fd 100644 --- a/Eigen/src/QR/CompleteOrthogonalDecomposition.h +++ b/Eigen/src/QR/CompleteOrthogonalDecomposition.h @@ -197,6 +197,21 @@ template class CompleteOrthogonalDecomposition return m_cpqr.colsPermutation(); } + /** \returns 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. + * + * \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. + * + * \sa absDeterminant(), logAbsDeterminant(), MatrixBase::determinant() + */ + typename MatrixType::Scalar determinant() const; + /** \returns the absolute value 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) @@ -208,7 +223,7 @@ template class CompleteOrthogonalDecomposition * of large enough dimension, there is a risk of overflow/underflow. * One way to work around that is to use logAbsDeterminant() instead. * - * \sa logAbsDeterminant(), MatrixBase::determinant() + * \sa determinant(), logAbsDeterminant(), MatrixBase::determinant() */ typename MatrixType::RealScalar absDeterminant() const; @@ -223,7 +238,7 @@ template class CompleteOrthogonalDecomposition * \note This method is useful to work around the risk of overflow/underflow * that's inherent to determinant computation. * - * \sa absDeterminant(), MatrixBase::determinant() + * \sa determinant(), absDeterminant(), MatrixBase::determinant() */ typename MatrixType::RealScalar logAbsDeterminant() const; @@ -407,6 +422,12 @@ template class CompleteOrthogonalDecomposition RowVectorType m_temp; }; +template +typename MatrixType::Scalar +CompleteOrthogonalDecomposition::determinant() const { + return m_cpqr.determinant(); +} + template typename MatrixType::RealScalar CompleteOrthogonalDecomposition::absDeterminant() const { diff --git a/Eigen/src/QR/FullPivHouseholderQR.h b/Eigen/src/QR/FullPivHouseholderQR.h index dcb9e1a96..ec7e19b55 100644 --- a/Eigen/src/QR/FullPivHouseholderQR.h +++ b/Eigen/src/QR/FullPivHouseholderQR.h @@ -210,6 +210,21 @@ template class FullPivHouseholderQR return m_rows_transpositions; } + /** \returns 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. + * + * \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. + * + * \sa absDeterminant(), logAbsDeterminant(), MatrixBase::determinant() + */ + typename MatrixType::Scalar determinant() const; + /** \returns the absolute value 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) @@ -221,7 +236,7 @@ template class FullPivHouseholderQR * of large enough dimension, there is a risk of overflow/underflow. * One way to work around that is to use logAbsDeterminant() instead. * - * \sa logAbsDeterminant(), MatrixBase::determinant() + * \sa determinant(), logAbsDeterminant(), MatrixBase::determinant() */ typename MatrixType::RealScalar absDeterminant() const; @@ -235,7 +250,7 @@ template class FullPivHouseholderQR * \note This method is useful to work around the risk of overflow/underflow that's inherent * to determinant computation. * - * \sa absDeterminant(), MatrixBase::determinant() + * \sa determinant(), absDeterminant(), MatrixBase::determinant() */ typename MatrixType::RealScalar logAbsDeterminant() const; @@ -419,9 +434,19 @@ template class FullPivHouseholderQR RealScalar m_prescribedThreshold, m_maxpivot; Index m_nonzero_pivots; RealScalar m_precision; - Index m_det_pq; + Index m_det_p; }; +template +typename MatrixType::Scalar FullPivHouseholderQR::determinant() 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 m_qr.diagonal().prod() * detQ * Scalar(m_det_p); +} + template typename MatrixType::RealScalar FullPivHouseholderQR::absDeterminant() const { @@ -531,7 +556,7 @@ void FullPivHouseholderQR::computeInPlace() for(Index k = 0; k < size; ++k) m_cols_permutation.applyTranspositionOnTheRight(k, m_cols_transpositions.coeff(k)); - m_det_pq = (number_of_transpositions%2) ? -1 : 1; + m_det_p = (number_of_transpositions%2) ? -1 : 1; m_isInitialized = true; } diff --git a/Eigen/src/QR/HouseholderQR.h b/Eigen/src/QR/HouseholderQR.h index 5e67463c7..abfefd1e1 100644 --- a/Eigen/src/QR/HouseholderQR.h +++ b/Eigen/src/QR/HouseholderQR.h @@ -184,6 +184,21 @@ template class HouseholderQR return *this; } + /** \returns 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. + * + * \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. + * + * \sa absDeterminant(), logAbsDeterminant(), MatrixBase::determinant() + */ + typename MatrixType::Scalar determinant() const; + /** \returns the absolute value 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) @@ -195,7 +210,7 @@ template class HouseholderQR * of large enough dimension, there is a risk of overflow/underflow. * One way to work around that is to use logAbsDeterminant() instead. * - * \sa logAbsDeterminant(), MatrixBase::determinant() + * \sa determinant(), logAbsDeterminant(), MatrixBase::determinant() */ typename MatrixType::RealScalar absDeterminant() const; @@ -209,7 +224,7 @@ template class HouseholderQR * \note This method is useful to work around the risk of overflow/underflow that's inherent * to determinant computation. * - * \sa absDeterminant(), MatrixBase::determinant() + * \sa determinant(), absDeterminant(), MatrixBase::determinant() */ typename MatrixType::RealScalar logAbsDeterminant() const; @@ -242,6 +257,57 @@ template class HouseholderQR bool m_isInitialized; }; +namespace internal { + +/** \internal */ +template +struct householder_determinant +{ + static void run(const HCoeffs& hCoeffs, Scalar& out_det) + { + out_det = Scalar(1); + Index size = hCoeffs.rows(); + for (Index i = 0; i < size; i ++) + { + // For each valid reflection Q_n, + // det(Q_n) = - conj(h_n) / h_n + // where h_n is the Householder coefficient. + if (hCoeffs(i) != Scalar(0)) + out_det *= - numext::conj(hCoeffs(i)) / hCoeffs(i); + } + } +}; + +/** \internal */ +template +struct householder_determinant +{ + static void run(const HCoeffs& hCoeffs, Scalar& out_det) + { + bool negated = false; + Index size = hCoeffs.rows(); + for (Index i = 0; i < size; i ++) + { + // Each valid reflection negates the determinant. + if (hCoeffs(i) != Scalar(0)) + negated ^= true; + } + out_det = negated ? Scalar(-1) : Scalar(1); + } +}; + +} // end namespace internal + +template +typename MatrixType::Scalar HouseholderQR::determinant() 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 m_qr.diagonal().prod() * detQ; +} + template typename MatrixType::RealScalar HouseholderQR::absDeterminant() const { diff --git a/test/qr.cpp b/test/qr.cpp index c38e3439b..36f312105 100644 --- a/test/qr.cpp +++ b/test/qr.cpp @@ -75,15 +75,17 @@ template void qr_invertible() // now construct a matrix with prescribed determinant m1.setZero(); for(int i = 0; i < size; i++) m1(i,i) = internal::random(); - RealScalar absdet = abs(m1.diagonal().prod()); + Scalar det = m1.diagonal().prod(); + RealScalar absdet = abs(det); m3 = qr.householderQ(); // get a unitary - m1 = m3 * m1 * m3; + m1 = m3 * m1 * m3.adjoint(); qr.compute(m1); VERIFY_IS_APPROX(log(absdet), qr.logAbsDeterminant()); // 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 - VERIFY_IS_MUCH_SMALLER_THAN( abs(absdet-qr.absDeterminant()), numext::maxi(RealScalar(pow(0.5,size)),numext::maxi(abs(absdet),abs(qr.absDeterminant()))) ); - + RealScalar tol = numext::maxi(RealScalar(pow(0.5,size)), numext::maxi(abs(absdet), abs(qr.absDeterminant()))); + VERIFY_IS_MUCH_SMALLER_THAN(abs(det - qr.determinant()), tol); + VERIFY_IS_MUCH_SMALLER_THAN(abs(absdet - qr.absDeterminant()), tol); } template void qr_verify_assert() @@ -96,6 +98,7 @@ template void qr_verify_assert() VERIFY_RAISES_ASSERT(qr.transpose().solve(tmp)) VERIFY_RAISES_ASSERT(qr.adjoint().solve(tmp)) VERIFY_RAISES_ASSERT(qr.householderQ()) + VERIFY_RAISES_ASSERT(qr.determinant()) VERIFY_RAISES_ASSERT(qr.absDeterminant()) VERIFY_RAISES_ASSERT(qr.logAbsDeterminant()) } diff --git a/test/qr_colpivoting.cpp b/test/qr_colpivoting.cpp index e25167138..4185f514e 100644 --- a/test/qr_colpivoting.cpp +++ b/test/qr_colpivoting.cpp @@ -273,10 +273,12 @@ template void qr_invertible() // now construct a matrix with prescribed determinant m1.setZero(); for(int i = 0; i < size; i++) m1(i,i) = internal::random(); - RealScalar absdet = abs(m1.diagonal().prod()); + Scalar det = m1.diagonal().prod(); + RealScalar absdet = abs(det); m3 = qr.householderQ(); // get a unitary - m1 = m3 * m1 * m3; + m1 = m3 * m1 * m3.adjoint(); qr.compute(m1); + VERIFY_IS_APPROX(det, qr.determinant()); VERIFY_IS_APPROX(absdet, qr.absDeterminant()); VERIFY_IS_APPROX(log(absdet), qr.logAbsDeterminant()); } @@ -296,6 +298,7 @@ template void qr_verify_assert() VERIFY_RAISES_ASSERT(qr.isSurjective()) VERIFY_RAISES_ASSERT(qr.isInvertible()) VERIFY_RAISES_ASSERT(qr.inverse()) + VERIFY_RAISES_ASSERT(qr.determinant()) VERIFY_RAISES_ASSERT(qr.absDeterminant()) VERIFY_RAISES_ASSERT(qr.logAbsDeterminant()) } @@ -315,6 +318,7 @@ template void cod_verify_assert() VERIFY_RAISES_ASSERT(cod.isSurjective()) VERIFY_RAISES_ASSERT(cod.isInvertible()) VERIFY_RAISES_ASSERT(cod.pseudoInverse()) + VERIFY_RAISES_ASSERT(cod.determinant()) VERIFY_RAISES_ASSERT(cod.absDeterminant()) VERIFY_RAISES_ASSERT(cod.logAbsDeterminant()) } diff --git a/test/qr_fullpivoting.cpp b/test/qr_fullpivoting.cpp index f2d8cb33e..cca9a8c1e 100644 --- a/test/qr_fullpivoting.cpp +++ b/test/qr_fullpivoting.cpp @@ -98,10 +98,12 @@ template void qr_invertible() // now construct a matrix with prescribed determinant m1.setZero(); for(int i = 0; i < size; i++) m1(i,i) = internal::random(); - RealScalar absdet = abs(m1.diagonal().prod()); + Scalar det = m1.diagonal().prod(); + RealScalar absdet = abs(det); m3 = qr.matrixQ(); // get a unitary - m1 = m3 * m1 * m3; + m1 = m3 * m1 * m3.adjoint(); qr.compute(m1); + VERIFY_IS_APPROX(det, qr.determinant()); VERIFY_IS_APPROX(absdet, qr.absDeterminant()); VERIFY_IS_APPROX(log(absdet), qr.logAbsDeterminant()); } @@ -121,6 +123,7 @@ template void qr_verify_assert() VERIFY_RAISES_ASSERT(qr.isSurjective()) VERIFY_RAISES_ASSERT(qr.isInvertible()) VERIFY_RAISES_ASSERT(qr.inverse()) + VERIFY_RAISES_ASSERT(qr.determinant()) VERIFY_RAISES_ASSERT(qr.absDeterminant()) VERIFY_RAISES_ASSERT(qr.logAbsDeterminant()) }