mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-03-07 18:27:40 +08:00
Code Issues Packages Projects Releases Wiki Activity

Add method signDeterminant() to QR and related decompositions.

Browse Source
This commit is contained in:
Rasmus Munk Larsen 2024-02-20 23:44:28 +00:00
parent db6b9db33b
commit 3859e8d5b2
7 changed files with 138 additions and 7 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

33
Eigen/src/QR/ColPivHouseholderQR.h
View File

@ -238,6 +238,20 @@ class ColPivHouseholderQR : public SolverBase<ColPivHouseholderQR<MatrixType_, P
*/
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.
*
* \sa determinant(), absDeterminant(), logAbsDeterminant(), MatrixBase::determinant()
*/
typename MatrixType::Scalar signDeterminant() const;
/** \returns the rank of the matrix of which *this is the QR decomposition.
*
* \note This method has to determine which pivots should be considered nonzero.
@ -424,26 +438,39 @@ class ColPivHouseholderQR : public SolverBase<ColPivHouseholderQR<MatrixType_, P
template <typename MatrixType, typename PermutationIndex>
typename MatrixType::Scalar ColPivHouseholderQR<MatrixType, PermutationIndex>::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<HCoeffsType, Scalar, NumTraits<Scalar>::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, typename PermutationIndex>
typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType, PermutationIndex>::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, typename PermutationIndex>
typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType, PermutationIndex>::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<RealScalar>::infinity();
}
template <typename MatrixType, typename PermutationIndex>
typename MatrixType::Scalar ColPivHouseholderQR<MatrixType, PermutationIndex>::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<HCoeffsType, Scalar, NumTraits<Scalar>::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

20
Eigen/src/QR/CompleteOrthogonalDecomposition.h
View File

@ -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<MatrixType, Perm
return m_cpqr.logAbsDeterminant();
}
template <typename MatrixType, typename PermutationIndex>
typename MatrixType::Scalar CompleteOrthogonalDecomposition<MatrixType, PermutationIndex>::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.

33
Eigen/src/QR/FullPivHouseholderQR.h
View File

@ -248,6 +248,20 @@ class FullPivHouseholderQR : public SolverBase<FullPivHouseholderQR<MatrixType_,
*/
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.
*
* \sa determinant(), absDeterminant(), logAbsDeterminant(), MatrixBase::determinant()
*/
typename MatrixType::Scalar signDeterminant() const;
/** \returns the rank of the matrix of which *this is the QR decomposition.
*
* \note This method has to determine which pivots should be considered nonzero.
@ -421,26 +435,39 @@ class FullPivHouseholderQR : public SolverBase<FullPivHouseholderQR<MatrixType_,
template <typename MatrixType, typename PermutationIndex>
typename MatrixType::Scalar FullPivHouseholderQR<MatrixType, PermutationIndex>::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<HCoeffsType, Scalar, NumTraits<Scalar>::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, typename PermutationIndex>
typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType, PermutationIndex>::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, typename PermutationIndex>
typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType, PermutationIndex>::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<RealScalar>::infinity();
}
template <typename MatrixType, typename PermutationIndex>
typename MatrixType::Scalar FullPivHouseholderQR<MatrixType, PermutationIndex>::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<HCoeffsType, Scalar, NumTraits<Scalar>::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

33
Eigen/src/QR/HouseholderQR.h
View File

@ -187,6 +187,8 @@ class HouseholderQR : public SolverBase<HouseholderQR<MatrixType_>> {
* \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<HouseholderQR<MatrixType_>> {
* \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<HouseholderQR<MatrixType_>> {
* \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<MatrixType>::logAbsDeterminant() c
return m_qr.diagonal().cwiseAbs().array().log().sum();
}
template <typename MatrixType>
typename MatrixType::Scalar HouseholderQR<MatrixType>::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<HCoeffsType, Scalar, NumTraits<Scalar>::IsComplex>::run(m_hCoeffs, detQ);
return detQ * m_qr.diagonal().array().sign().prod();
}
namespace internal {
/** \internal */

3
test/qr.cpp
View File

@ -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) {

21
test/qr_colpivoting.cpp
View File

@ -21,6 +21,7 @@ void cod() {
Index rank = internal::random<Index>(1, (std::min)(rows, cols) - 1);
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> 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<Index>(2, 20);
matrix.setZero(size, size);
for (int i = 0; i < size; i++) {
matrix(i, i) = internal::random<Scalar>();
}
Scalar det = matrix.diagonal().prod();
RealScalar absdet = abs(det);
CompleteOrthogonalDecomposition<MatrixType> 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 <typename MatrixType, int Cols2>
@ -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 <typename MatrixType>
@ -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 <typename MatrixType>
@ -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) {

2
test/qr_fullpivoting.cpp
View File

@ -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 <typename MatrixType>
@ -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) {