GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-04-20 08:39:37 +08:00
Code Issues Packages Projects Releases Wiki Activity

Add true determinant to QR and it's variants

Browse Source
This commit is contained in:
sjusju 2022-07-29 18:24:14 +00:00 committed by Rasmus Munk Larsen
parent b7668c0371
commit ef4654bae7
7 changed files with 167 additions and 20 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -219,6 +219,21 @@ template<typename MatrixType_> 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<typename MatrixType_> 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<typename MatrixType_> 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<typename MatrixType_> 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>
typename MatrixType::Scalar ColPivHouseholderQR<MatrixType>::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<HCoeffsType, Scalar, NumTraits<Scalar>::IsComplex>::run(m_hCoeffs, detQ);
return m_qr.diagonal().prod() * detQ * Scalar(m_det_p);
}
template<typename MatrixType>
typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::absDeterminant() const
{
@ -574,7 +599,7 @@ void ColPivHouseholderQR<MatrixType>::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;
}

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

@ -197,6 +197,21 @@ template <typename MatrixType_> 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 <typename MatrixType_> 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 <typename MatrixType_> 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 <typename MatrixType_> class CompleteOrthogonalDecomposition
RowVectorType m_temp;
};
template <typename MatrixType>
typename MatrixType::Scalar
CompleteOrthogonalDecomposition<MatrixType>::determinant() const {
return m_cpqr.determinant();
}
template <typename MatrixType>
typename MatrixType::RealScalar
CompleteOrthogonalDecomposition<MatrixType>::absDeterminant() const {

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

@ -210,6 +210,21 @@ template<typename MatrixType_> 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<typename MatrixType_> 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<typename MatrixType_> 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<typename MatrixType_> 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>
typename MatrixType::Scalar FullPivHouseholderQR<MatrixType>::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<HCoeffsType, Scalar, NumTraits<Scalar>::IsComplex>::run(m_hCoeffs, detQ);
return m_qr.diagonal().prod() * detQ * Scalar(m_det_p);
}
template<typename MatrixType>
typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType>::absDeterminant() const
{
@ -531,7 +556,7 @@ void FullPivHouseholderQR<MatrixType>::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;
}

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

@ -184,6 +184,21 @@ template<typename MatrixType_> 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<typename MatrixType_> 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<typename MatrixType_> 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<typename MatrixType_> class HouseholderQR
bool m_isInitialized;
};
namespace internal {
/** \internal */
template<typename HCoeffs, typename Scalar, bool IsComplex>
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<typename HCoeffs, typename Scalar>
struct householder_determinant<HCoeffs, Scalar, false>
{
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>
typename MatrixType::Scalar HouseholderQR<MatrixType>::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<HCoeffsType, Scalar, NumTraits<Scalar>::IsComplex>::run(m_hCoeffs, detQ);
return m_qr.diagonal().prod() * detQ;
}
template<typename MatrixType>
typename MatrixType::RealScalar HouseholderQR<MatrixType>::absDeterminant() const
{

11
test/qr.cpp
View File

@ -75,15 +75,17 @@ template<typename MatrixType> void qr_invertible()
// now construct a matrix with prescribed determinant
m1.setZero();
for(int i = 0; i < size; i++) m1(i,i) = internal::random<Scalar>();
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<RealScalar>(abs(absdet),abs(qr.absDeterminant()))) );
RealScalar tol = numext::maxi(RealScalar(pow(0.5,size)), numext::maxi<RealScalar>(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<typename MatrixType> void qr_verify_assert()
@ -96,6 +98,7 @@ template<typename MatrixType> 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())
}

8
test/qr_colpivoting.cpp
View File

@ -273,10 +273,12 @@ template<typename MatrixType> void qr_invertible()
// now construct a matrix with prescribed determinant
m1.setZero();
for(int i = 0; i < size; i++) m1(i,i) = internal::random<Scalar>();
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<typename MatrixType> 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<typename MatrixType> 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())
}

7
test/qr_fullpivoting.cpp
View File

@ -98,10 +98,12 @@ template<typename MatrixType> void qr_invertible()
// now construct a matrix with prescribed determinant
m1.setZero();
for(int i = 0; i < size; i++) m1(i,i) = internal::random<Scalar>();
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<typename MatrixType> 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())
}