Cholesky/LU/SparseCholesky: Add absDeterminant, logAbsDeterminant and signDeterminant libeigen/eigen!2991 Closes #757 Co-authored-by: Rasmus Munk Larsen <rmlarsen@gmail.com>
diff --git a/Eigen/src/Cholesky/BunchKaufman.h b/Eigen/src/Cholesky/BunchKaufman.h index 2f9e37f..46a270a 100644 --- a/Eigen/src/Cholesky/BunchKaufman.h +++ b/Eigen/src/Cholesky/BunchKaufman.h
@@ -262,6 +262,61 @@ MatrixType reconstructedMatrix() const; + /** \returns the determinant of the matrix of which *this is the Bunch-Kaufman decomposition. + * + * It has only linear complexity (that is, O(n) where n is the dimension of the square matrix) + * as the decomposition has already been computed. + * + * \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() and signDeterminant() instead. + * Also, do not rely on the determinant being exactly zero for testing + * singularity or rank-deficiency. + * + * \sa absDeterminant(), logAbsDeterminant(), signDeterminant(), MatrixBase::determinant() + */ + Scalar determinant() const; + + /** \returns the absolute value of the determinant of the matrix of which *this is the Bunch-Kaufman + * decomposition. + * + * It has only linear complexity (that is, O(n) where n is the dimension of the square matrix) + * as the decomposition has already been computed. + * + * \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 determinant(), logAbsDeterminant(), signDeterminant(), MatrixBase::determinant() + */ + RealScalar absDeterminant() const; + + /** \returns the natural log of the absolute value of the determinant of the matrix of which *this is the + * Bunch-Kaufman decomposition. + * + * It has only linear complexity (that is, O(n) where n is the dimension of the square matrix) + * as the decomposition has already been computed. + * + * \note This method is useful to work around the risk of overflow/underflow that's inherent + * to determinant computation. The 2x2 blocks of D are scaled by their off-diagonal entry before their + * determinant is formed, so a block whose own determinant is out of range still contributes a + * finite term. + * + * \sa determinant(), absDeterminant(), signDeterminant(), MatrixBase::determinant() + */ + RealScalar logAbsDeterminant() const; + + /** \returns the sign of the determinant of the matrix of which *this is the Bunch-Kaufman decomposition, + * that is, \c 1, \c -1, or \c 0 if the matrix is singular. + * + * The sign is read off the inertia rather than computed from a product, so it is exact and cannot + * overflow: \f$ \mathrm{sign}(\det A) = (-1)^{n_-} \f$, with \f$ n_- \f$ the number of negative + * eigenvalues, which congruence leaves invariant (Sylvester's law of inertia). + * + * \sa determinant(), absDeterminant(), logAbsDeterminant(), MatrixBase::determinant() + */ + Scalar signDeterminant() const; + /** \returns the adjoint of \c *this, that is, a const reference to the decomposition itself as the underlying matrix * is self-adjoint. * @@ -306,6 +361,16 @@ /** \internal Compute the inertia (counts of positive / negative / zero eigenvalues) from D. */ void computeInertia(); + /** \internal \returns \f$ \det(D_k)/|d_{21}|^2 \f$ for a 2x2 block of D, which is real and shares the + * sign of \f$ \det(D_k) \f$ since \f$ |d_{21}| > 0 \f$ there. */ + static RealScalar scaledBlockDeterminant(const RealScalar& d11, const RealScalar& d22, const RealScalar& d21); + + /** \internal \returns \f$ \det(D) \f$, the product over the 1x1 and 2x2 diagonal blocks of D. D is + * Hermitian, so a block determinant is real: \f$ d_{11} \f$ for a 1x1 block and + * \f$ d_{11} d_{22} - |d_{21}|^2 \f$ for a 2x2 one. Accumulated as a mantissa and a power of two, + * so that only the result has to be representable. */ + RealScalar determinantD() const; + MatrixType m_matrix; RealScalar m_l1_norm; TranspositionType m_transpositions; @@ -802,6 +867,18 @@ } } +// det(D_k) itself over- or underflows on an extreme-scaled 2x2 block, so it is only ever formed scaled by +// |d21|^2. Divide by |d21| twice rather than multiply by its reciprocal, which overflows once |d21| is +// subnormal. The pivot criterion bounds |q| by a^2 < 1 with a = (1+sqrt(17))/8, so q >= 1 or NaN is an +// artifact of d22/d21 overflowing -- the criterion bounds |d22| only against its own row -- and +// det(D_k) = -|d21|^2 to within that same bound there. +template <typename MatrixType, int UpLo_> +typename BunchKaufman<MatrixType, UpLo_>::RealScalar BunchKaufman<MatrixType, UpLo_>::scaledBlockDeterminant( + const RealScalar& d11, const RealScalar& d22, const RealScalar& d21) { + const RealScalar q = (d11 / d21) * (d22 / d21); + return q < RealScalar(1) ? q - RealScalar(1) : RealScalar(-1); +} + template <typename MatrixType, int UpLo_> void BunchKaufman<MatrixType, UpLo_>::computeInertia() { const Index n = m_matrix.rows(); @@ -811,11 +888,8 @@ if (k + 1 < n && !numext::is_exactly_zero(m_subdiag.coeff(k))) { const RealScalar d11 = numext::real(m_matrix.coeff(k, k)); const RealScalar d22 = numext::real(m_matrix.coeff(k + 1, k + 1)); - const Scalar d21 = m_subdiag.coeff(k); - // Scaled determinant denom = det/|d21|^2 (|d21|^2 > 0 for a 2x2 block), so sign(denom) == sign(det); - // avoids forming det = d11*d22 - |d21|^2, which over/underflows on extreme-scaled 2x2 blocks. - const Scalar id = Scalar(1) / d21; - const RealScalar denom = numext::real((d22 * id) * (d11 * numext::conj(id))) - RealScalar(1); + const RealScalar d21 = numext::abs(m_subdiag.coeff(k)); + const RealScalar denom = scaledBlockDeterminant(d11, d22, d21); if (denom < RealScalar(0)) { // Indefinite 2x2 block: one positive and one negative eigenvalue. ++m_n_pos; @@ -917,6 +991,87 @@ return true; } +template <typename MatrixType, int UpLo_> +typename BunchKaufman<MatrixType, UpLo_>::RealScalar BunchKaufman<MatrixType, UpLo_>::determinantD() const { + const Index n = m_matrix.rows(); + // det(D) is accumulated as mantissa * 2^exponent, the mantissa renormalized to [1/2, 1) after each + // block. Both a single 2x2 block determinant and the running product can leave the representable range + // while det(D) itself stays in it, and an overflowed block meeting an underflowed one gives NaN. + RealScalar mantissa(1); + Index exponent = 0; + Index k = 0; + while (k < n) { + RealScalar blockM, blockE; + if (k + 1 < n && !numext::is_exactly_zero(m_subdiag.coeff(k))) { + const RealScalar d11 = numext::real(m_matrix.coeff(k, k)); + const RealScalar d22 = numext::real(m_matrix.coeff(k + 1, k + 1)); + const RealScalar d21 = numext::abs(m_subdiag.coeff(k)); + RealScalar e21; + const RealScalar m21 = internal::pfrexp<RealScalar>(d21, e21); + blockM = m21 * m21 * scaledBlockDeterminant(d11, d22, d21); + blockE = RealScalar(2) * e21; + k += 2; + } else { + blockM = internal::pfrexp<RealScalar>(numext::real(m_matrix.coeff(k, k)), blockE); + k += 1; + } + // |blockM| < 2 and, unless the block is singular, above 1/16, so the product below stays normal and + // its rounding is the only error this step adds. + RealScalar renorm; + mantissa = internal::pfrexp<RealScalar>(mantissa * blockM, renorm); + exponent += Index(blockE) + Index(renorm); + } + // ldexp() saturates to zero or infinity but takes an int exponent; past this magnitude it has already + // saturated, so clamping first cannot change the result. + const Index limit = Index(NumTraits<RealScalar>::max_exponent()) - Index(NumTraits<RealScalar>::min_exponent()) + + Index(NumTraits<RealScalar>::digits()); + return numext::ldexp(mantissa, int(numext::mini(numext::maxi(exponent, -limit), limit))); +} + +// A = P^T L D L^* P with L unit triangular, so det(A) = det(D) = prod over the 1x1 and 2x2 blocks of D. + +template <typename MatrixType, int UpLo_> +typename BunchKaufman<MatrixType, UpLo_>::Scalar BunchKaufman<MatrixType, UpLo_>::determinant() const { + eigen_assert(m_isInitialized && "BunchKaufman is not initialized."); + return Scalar(determinantD()); +} + +template <typename MatrixType, int UpLo_> +typename BunchKaufman<MatrixType, UpLo_>::RealScalar BunchKaufman<MatrixType, UpLo_>::absDeterminant() const { + eigen_assert(m_isInitialized && "BunchKaufman is not initialized."); + return numext::abs(determinantD()); +} + +template <typename MatrixType, int UpLo_> +typename BunchKaufman<MatrixType, UpLo_>::RealScalar BunchKaufman<MatrixType, UpLo_>::logAbsDeterminant() const { + eigen_assert(m_isInitialized && "BunchKaufman is not initialized."); + const Index n = m_matrix.rows(); + RealScalar result(0); + Index k = 0; + while (k < n) { + if (k + 1 < n && !numext::is_exactly_zero(m_subdiag.coeff(k))) { + // log|det(D_k)| = 2 log|d21| + log|det(D_k)/|d21|^2|. + const RealScalar d11 = numext::real(m_matrix.coeff(k, k)); + const RealScalar d22 = numext::real(m_matrix.coeff(k + 1, k + 1)); + const RealScalar d21 = numext::abs(m_subdiag.coeff(k)); + const RealScalar scaled = scaledBlockDeterminant(d11, d22, d21); + result += RealScalar(2) * numext::log(d21) + numext::log(numext::abs(scaled)); + k += 2; + } else { + result += numext::log(numext::abs(numext::real(m_matrix.coeff(k, k)))); + k += 1; + } + } + return result; +} + +template <typename MatrixType, int UpLo_> +typename BunchKaufman<MatrixType, UpLo_>::Scalar BunchKaufman<MatrixType, UpLo_>::signDeterminant() const { + eigen_assert(m_isInitialized && "BunchKaufman is not initialized."); + if (m_n_zero > 0) return Scalar(0); + return Scalar((m_n_neg % 2 == 0) ? 1 : -1); +} + /** \returns the matrix represented by the decomposition, i.e., the product \f$ P^T L D L^* P \f$. * This function is provided for debug purposes. */ template <typename MatrixType, int UpLo_>
diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h index aca8eb3..5ecbda8 100644 --- a/Eigen/src/Cholesky/LDLT.h +++ b/Eigen/src/Cholesky/LDLT.h
@@ -227,6 +227,65 @@ template <typename Derived> LDLT& rankUpdate(const MatrixBase<Derived>& w, const RealScalar& alpha = 1); + /** \returns the determinant of the matrix of which *this is the Cholesky decomposition. + * + * It has only linear complexity (that is, O(n) where n is the dimension of the square matrix) + * as the Cholesky decomposition has already been computed. + * + * \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() and signDeterminant() instead. + * Also, do not rely on the determinant being exactly zero for testing + * singularity or rank-deficiency. + * + * \pre info() returns \c Success. A failed factorization does not represent the input matrix. + * + * \sa absDeterminant(), logAbsDeterminant(), signDeterminant(), MatrixBase::determinant() + */ + Scalar determinant() const; + + /** \returns the absolute value of the determinant of the matrix of which *this is the Cholesky decomposition. + * + * It has only linear complexity (that is, O(n) where n is the dimension of the square matrix) + * as the Cholesky decomposition has already been computed. + * + * \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. + * + * \pre info() returns \c Success. A failed factorization does not represent the input matrix. + * + * \sa determinant(), logAbsDeterminant(), signDeterminant(), MatrixBase::determinant() + */ + RealScalar absDeterminant() const; + + /** \returns the natural log of the absolute value of the determinant of the matrix of which *this is the Cholesky + * decomposition. + * + * It has only linear complexity (that is, O(n) where n is the dimension of the square matrix) + * as the Cholesky decomposition has already been computed. + * + * \note This method is useful to work around the risk of overflow/underflow that's inherent + * to determinant computation. + * + * \pre info() returns \c Success. A failed factorization does not represent the input matrix. + * + * \sa determinant(), absDeterminant(), signDeterminant(), MatrixBase::determinant() + */ + RealScalar logAbsDeterminant() const; + + /** \returns the sign of the determinant of the matrix of which *this is the Cholesky decomposition, + * that is, \c 1, \c -1, or \c 0 if the matrix is singular. + * + * It has only linear complexity (that is, O(n) where n is the dimension of the square matrix) + * as the Cholesky decomposition has already been computed. + * + * \pre info() returns \c Success. A failed factorization does not represent the input matrix. + * + * \sa determinant(), absDeterminant(), logAbsDeterminant(), MatrixBase::determinant() + */ + Scalar signDeterminant() const; + /** \returns the internal LDLT decomposition matrix * * TODO: document the storage layout. @@ -559,6 +618,36 @@ return *this; } +// A = P^T L D L^* P with L unit lower triangular and D real diagonal, so det(A) = prod(D_ii). + +template <typename MatrixType_, int UpLo_> +typename LDLT<MatrixType_, UpLo_>::Scalar LDLT<MatrixType_, UpLo_>::determinant() const { + eigen_assert(m_isInitialized && "LDLT is not initialized."); + eigen_assert(m_info == Success && "LDLT failed because of a zero pivot."); + return Scalar(vectorD().real().prod()); +} + +template <typename MatrixType_, int UpLo_> +typename LDLT<MatrixType_, UpLo_>::RealScalar LDLT<MatrixType_, UpLo_>::absDeterminant() const { + eigen_assert(m_isInitialized && "LDLT is not initialized."); + eigen_assert(m_info == Success && "LDLT failed because of a zero pivot."); + return numext::abs(vectorD().real().prod()); +} + +template <typename MatrixType_, int UpLo_> +typename LDLT<MatrixType_, UpLo_>::RealScalar LDLT<MatrixType_, UpLo_>::logAbsDeterminant() const { + eigen_assert(m_isInitialized && "LDLT is not initialized."); + eigen_assert(m_info == Success && "LDLT failed because of a zero pivot."); + return vectorD().real().cwiseAbs().array().log().sum(); +} + +template <typename MatrixType_, int UpLo_> +typename LDLT<MatrixType_, UpLo_>::Scalar LDLT<MatrixType_, UpLo_>::signDeterminant() const { + eigen_assert(m_isInitialized && "LDLT is not initialized."); + eigen_assert(m_info == Success && "LDLT failed because of a zero pivot."); + return Scalar(vectorD().real().array().sign().prod()); +} + #ifndef EIGEN_PARSED_BY_DOXYGEN template <typename MatrixType_, int UpLo_> template <typename RhsType, typename DstType>
diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h index 1a1bb48..5f3d186 100644 --- a/Eigen/src/Cholesky/LLT.h +++ b/Eigen/src/Cholesky/LLT.h
@@ -194,6 +194,64 @@ template <typename VectorType> LLT& rankUpdate(const VectorType& vec, const RealScalar& sigma = 1); + /** \returns the determinant of the matrix of which *this is the Cholesky decomposition. + * + * It has only linear complexity (that is, O(n) where n is the dimension of the square matrix) + * as the Cholesky decomposition has already been computed. + * + * \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. + * + * \pre info() returns \c Success. A failed factorization does not represent the input matrix. + * + * \sa absDeterminant(), logAbsDeterminant(), signDeterminant(), MatrixBase::determinant() + */ + Scalar determinant() const; + + /** \returns the absolute value of the determinant of the matrix of which *this is the Cholesky decomposition. + * + * It has only linear complexity (that is, O(n) where n is the dimension of the square matrix) + * as the Cholesky decomposition has already been computed. + * + * \note The decomposed matrix is positive definite, so this is the determinant itself. + * + * \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. + * + * \pre info() returns \c Success. A failed factorization does not represent the input matrix. + * + * \sa determinant(), logAbsDeterminant(), signDeterminant(), MatrixBase::determinant() + */ + RealScalar absDeterminant() const; + + /** \returns the natural log of the absolute value of the determinant of the matrix of which *this is the Cholesky + * decomposition. + * + * It has only linear complexity (that is, O(n) where n is the dimension of the square matrix) + * as the Cholesky decomposition has already been computed. + * + * \note This method is useful to work around the risk of overflow/underflow that's inherent + * to determinant computation. + * + * \pre info() returns \c Success. A failed factorization does not represent the input matrix. + * + * \sa determinant(), absDeterminant(), signDeterminant(), MatrixBase::determinant() + */ + RealScalar logAbsDeterminant() const; + + /** \returns the sign of the determinant of the matrix of which *this is the Cholesky decomposition, + * which is \c 1 since that matrix is positive definite. + * + * This method is provided for compatibility with the other decompositions, thus enabling generic code. + * + * \pre info() returns \c Success. A failed factorization does not represent the input matrix. + * + * \sa determinant(), absDeterminant(), logAbsDeterminant(), MatrixBase::determinant() + */ + Scalar signDeterminant() const; + #ifndef EIGEN_PARSED_BY_DOXYGEN template <typename RhsType, typename DstType> void _solve_impl(const RhsType& rhs, DstType& dst) const; @@ -435,6 +493,34 @@ return *this; } +// A = L L^*, with L real and positive on the diagonal, so det(A) = prod(L_ii)^2 > 0. + +template <typename MatrixType_, int UpLo_> +typename LLT<MatrixType_, UpLo_>::Scalar LLT<MatrixType_, UpLo_>::determinant() const { + return Scalar(absDeterminant()); +} + +template <typename MatrixType_, int UpLo_> +typename LLT<MatrixType_, UpLo_>::RealScalar LLT<MatrixType_, UpLo_>::absDeterminant() const { + eigen_assert(m_isInitialized && "LLT is not initialized."); + eigen_assert(m_info == Success && "LLT failed because matrix appears to be negative"); + return numext::abs2(m_matrix.diagonal().real().prod()); +} + +template <typename MatrixType_, int UpLo_> +typename LLT<MatrixType_, UpLo_>::RealScalar LLT<MatrixType_, UpLo_>::logAbsDeterminant() const { + eigen_assert(m_isInitialized && "LLT is not initialized."); + eigen_assert(m_info == Success && "LLT failed because matrix appears to be negative"); + return RealScalar(2) * m_matrix.diagonal().real().array().log().sum(); +} + +template <typename MatrixType_, int UpLo_> +typename LLT<MatrixType_, UpLo_>::Scalar LLT<MatrixType_, UpLo_>::signDeterminant() const { + eigen_assert(m_isInitialized && "LLT is not initialized."); + eigen_assert(m_info == Success && "LLT failed because matrix appears to be negative"); + return Scalar(1); +} + #ifndef EIGEN_PARSED_BY_DOXYGEN template <typename MatrixType_, int UpLo_> template <typename RhsType, typename DstType>
diff --git a/Eigen/src/LU/FullPivLU.h b/Eigen/src/LU/FullPivLU.h index 13dca4d..f2d406a 100644 --- a/Eigen/src/LU/FullPivLU.h +++ b/Eigen/src/LU/FullPivLU.h
@@ -269,11 +269,66 @@ * * \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() and signDeterminant() instead. + * Also, do not rely on the determinant being exactly zero for testing + * singularity or rank-deficiency. * - * \sa MatrixBase::determinant() + * \sa absDeterminant(), logAbsDeterminant(), signDeterminant(), MatrixBase::determinant() */ typename internal::traits<MatrixType>::Scalar determinant() const; + /** \returns the absolute value of the determinant of the matrix of which + * *this is the LU decomposition. It has only linear complexity + * (that is, O(n) where n is the dimension of the square matrix) + * as the LU 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. + * + * \note Returns exactly zero when rank() finds the decomposition rank-deficient, as the + * rank-revealing QR decompositions do. determinant() is not gated that way: it returns the + * product of the pivots whatever the rank. + * + * \sa determinant(), logAbsDeterminant(), signDeterminant(), MatrixBase::determinant() + */ + RealScalar absDeterminant() const; + + /** \returns the natural log of the absolute value of the determinant of the matrix of which + * *this is the LU decomposition. It has only linear complexity + * (that is, O(n) where n is the dimension of the square matrix) + * as the LU 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. + * + * \note Returns \c -infinity when rank() finds the decomposition rank-deficient. + * + * \sa determinant(), absDeterminant(), signDeterminant(), MatrixBase::determinant() + */ + RealScalar logAbsDeterminant() const; + + /** \returns the sign of the determinant of the matrix of which + * *this is the LU decomposition. It has only linear complexity + * (that is, O(n) where n is the dimension of the square matrix) + * as the LU 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. + * + * \note Returns zero when rank() finds the decomposition rank-deficient, matching the documented + * sign of a singular matrix. + * + * \sa determinant(), absDeterminant(), logAbsDeterminant(), MatrixBase::determinant() + */ + Scalar signDeterminant() const; + /** \returns the absolute value of the i-th pivot coefficient (for RankRevealingBase). */ RealScalar pivotCoeff(Index i) const { using std::abs; @@ -453,6 +508,30 @@ return Scalar(m_det_pq) * Scalar(m_lu.diagonal().prod()); } +template <typename MatrixType, typename PermutationIndex> +typename FullPivLU<MatrixType, PermutationIndex>::RealScalar FullPivLU<MatrixType, PermutationIndex>::absDeterminant() + const { + eigen_assert(m_isInitialized && "LU is not initialized."); + eigen_assert(m_lu.rows() == m_lu.cols() && "You can't take the determinant of a non-square matrix!"); + return isInjective() ? numext::abs(m_lu.diagonal().prod()) : RealScalar(0); +} + +template <typename MatrixType, typename PermutationIndex> +typename FullPivLU<MatrixType, PermutationIndex>::RealScalar +FullPivLU<MatrixType, PermutationIndex>::logAbsDeterminant() const { + eigen_assert(m_isInitialized && "LU is not initialized."); + eigen_assert(m_lu.rows() == m_lu.cols() && "You can't take the determinant of a non-square matrix!"); + return isInjective() ? m_lu.diagonal().cwiseAbs().array().log().sum() : -NumTraits<RealScalar>::infinity(); +} + +template <typename MatrixType, typename PermutationIndex> +typename FullPivLU<MatrixType, PermutationIndex>::Scalar FullPivLU<MatrixType, PermutationIndex>::signDeterminant() + const { + eigen_assert(m_isInitialized && "LU is not initialized."); + eigen_assert(m_lu.rows() == m_lu.cols() && "You can't take the determinant of a non-square matrix!"); + return isInjective() ? Scalar(m_det_pq) * m_lu.diagonal().array().sign().prod() : Scalar(0); +} + /** \returns the matrix represented by the decomposition, * i.e., it returns the product: \f$ P^{-1} L U Q^{-1} \f$. * This function is provided for debug purposes. */
diff --git a/Eigen/src/LU/PartialPivLU.h b/Eigen/src/LU/PartialPivLU.h index 851bb49..ff92663 100644 --- a/Eigen/src/LU/PartialPivLU.h +++ b/Eigen/src/LU/PartialPivLU.h
@@ -204,11 +204,59 @@ * * \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() and signDeterminant() instead. + * Also, do not rely on the determinant being exactly zero for testing + * singularity or rank-deficiency. * - * \sa MatrixBase::determinant() + * \sa absDeterminant(), logAbsDeterminant(), signDeterminant(), MatrixBase::determinant() */ Scalar determinant() const; + /** \returns the absolute value of the determinant of the matrix of which + * *this is the LU decomposition. It has only linear complexity + * (that is, O(n) where n is the dimension of the square matrix) + * as the LU 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. + * Also, do not rely on the determinant being exactly zero for testing + * singularity or rank-deficiency. + * + * \sa determinant(), logAbsDeterminant(), signDeterminant(), MatrixBase::determinant() + */ + RealScalar absDeterminant() const; + + /** \returns the natural log of the absolute value of the determinant of the matrix of which + * *this is the LU decomposition. It has only linear complexity + * (that is, O(n) where n is the dimension of the square matrix) + * as the LU 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(), signDeterminant(), MatrixBase::determinant() + */ + RealScalar logAbsDeterminant() const; + + /** \returns the sign of the determinant of the matrix of which + * *this is the LU decomposition. It has only linear complexity + * (that is, O(n) where n is the dimension of the square matrix) + * as the LU 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() + */ + Scalar signDeterminant() const; + MatrixType reconstructedMatrix() const; constexpr Index rows() const noexcept { return m_lu.rows(); } @@ -540,6 +588,27 @@ return Scalar(m_det_p) * m_lu.diagonal().prod(); } +template <typename MatrixType, typename PermutationIndex> +typename PartialPivLU<MatrixType, PermutationIndex>::RealScalar +PartialPivLU<MatrixType, PermutationIndex>::absDeterminant() const { + eigen_assert(m_isInitialized && "PartialPivLU is not initialized."); + return numext::abs(m_lu.diagonal().prod()); +} + +template <typename MatrixType, typename PermutationIndex> +typename PartialPivLU<MatrixType, PermutationIndex>::RealScalar +PartialPivLU<MatrixType, PermutationIndex>::logAbsDeterminant() const { + eigen_assert(m_isInitialized && "PartialPivLU is not initialized."); + return m_lu.diagonal().cwiseAbs().array().log().sum(); +} + +template <typename MatrixType, typename PermutationIndex> +typename PartialPivLU<MatrixType, PermutationIndex>::Scalar +PartialPivLU<MatrixType, PermutationIndex>::signDeterminant() const { + eigen_assert(m_isInitialized && "PartialPivLU is not initialized."); + return Scalar(m_det_p) * m_lu.diagonal().array().sign().prod(); +} + /** \returns the matrix represented by the decomposition, * i.e., it returns the product: P^{-1} L U. * This function is provided for debug purpose. */
diff --git a/Eigen/src/SparseCholesky/SimplicialCholesky.h b/Eigen/src/SparseCholesky/SimplicialCholesky.h index 6f7db68..fbe209c 100644 --- a/Eigen/src/SparseCholesky/SimplicialCholesky.h +++ b/Eigen/src/SparseCholesky/SimplicialCholesky.h
@@ -264,6 +264,8 @@ }; mutable ComputationInfo m_info; + // Set once factorize() has run, success or not: factorize_preordered() breaks out on a bad pivot, leaving + // the tails of m_diag and of m_matrix's diagonal unwritten. Readers of those also need m_info == Success. bool m_factorizationIsOk; bool m_analysisIsOk; @@ -455,9 +457,40 @@ /** \returns the determinant of the underlying matrix from the current factorization */ Scalar determinant() const { + eigen_assert(Base::m_factorizationIsOk && Base::m_info == Success && + "Simplicial LLT is not factorized, or its factorization failed"); Scalar detL = Base::m_matrix.diagonal().prod(); return numext::abs2(detL); } + + /** \returns the absolute value of the determinant of the underlying matrix from the current factorization */ + RealScalar absDeterminant() const { + eigen_assert(Base::m_factorizationIsOk && Base::m_info == Success && + "Simplicial LLT is not factorized, or its factorization failed"); + return numext::abs2(Base::m_matrix.diagonal().prod()); + } + + /** \returns the natural log of the absolute value of the determinant of the underlying matrix from the current + * factorization. + * + * Unlike determinant(), this stays finite for the large factorizations where a determinant overflows or underflows. + */ + RealScalar logAbsDeterminant() const { + eigen_assert(Base::m_factorizationIsOk && Base::m_info == Success && + "Simplicial LLT is not factorized, or its factorization failed"); + return RealScalar(2) * Base::m_matrix.diagonal().cwiseAbs().array().log().sum(); + } + + /** \returns the sign of the determinant of the underlying matrix, which is \c 1 since that matrix is positive + * definite. + * + * This method is provided for compatibility with the other decompositions, thus enabling generic code. + */ + Scalar signDeterminant() const { + eigen_assert(Base::m_factorizationIsOk && Base::m_info == Success && + "Simplicial LLT is not factorized, or its factorization failed"); + return Scalar(1); + } }; /** \ingroup SparseCholesky_Module @@ -543,7 +576,36 @@ void factorize(const MatrixType& a) { Base::template factorize<true, false>(a); } /** \returns the determinant of the underlying matrix from the current factorization */ - Scalar determinant() const { return Base::m_diag.prod(); } + Scalar determinant() const { + eigen_assert(Base::m_factorizationIsOk && Base::m_info == Success && + "Simplicial LDLT is not factorized, or its factorization failed"); + return Base::m_diag.prod(); + } + + /** \returns the absolute value of the determinant of the underlying matrix from the current factorization */ + RealScalar absDeterminant() const { + eigen_assert(Base::m_factorizationIsOk && Base::m_info == Success && + "Simplicial LDLT is not factorized, or its factorization failed"); + return numext::abs(Base::m_diag.real().prod()); + } + + /** \returns the natural log of the absolute value of the determinant of the underlying matrix from the current + * factorization. + * + * Unlike determinant(), this stays finite for the large factorizations where a determinant overflows or underflows. + */ + RealScalar logAbsDeterminant() const { + eigen_assert(Base::m_factorizationIsOk && Base::m_info == Success && + "Simplicial LDLT is not factorized, or its factorization failed"); + return Base::m_diag.real().cwiseAbs().array().log().sum(); + } + + /** \returns the sign of the determinant of the underlying matrix from the current factorization */ + Scalar signDeterminant() const { + eigen_assert(Base::m_factorizationIsOk && Base::m_info == Success && + "Simplicial LDLT is not factorized, or its factorization failed"); + return Scalar(Base::m_diag.real().array().sign().prod()); + } }; /** \ingroup SparseCholesky_Module @@ -626,9 +688,37 @@ /** \returns the determinant of the underlying matrix from the current factorization */ Scalar determinant() const { + eigen_assert(Base::m_factorizationIsOk && Base::m_info == Success && + "Simplicial LLT is not factorized, or its factorization failed"); Scalar detL = Base::m_matrix.diagonal().prod(); return detL * detL; } + + /** \returns the absolute value of the determinant of the underlying matrix from the current factorization */ + RealScalar absDeterminant() const { + eigen_assert(Base::m_factorizationIsOk && Base::m_info == Success && + "Simplicial LLT is not factorized, or its factorization failed"); + return numext::abs2(Base::m_matrix.diagonal().prod()); + } + + /** \returns the natural log of the absolute value of the determinant of the underlying matrix from the current + * factorization. + * + * Unlike determinant(), this stays finite for the large factorizations where a determinant overflows or underflows. + */ + RealScalar logAbsDeterminant() const { + eigen_assert(Base::m_factorizationIsOk && Base::m_info == Success && + "Simplicial LLT is not factorized, or its factorization failed"); + return RealScalar(2) * Base::m_matrix.diagonal().cwiseAbs().array().log().sum(); + } + + /** \returns the sign of the determinant of the underlying matrix from the current factorization */ + Scalar signDeterminant() const { + eigen_assert(Base::m_factorizationIsOk && Base::m_info == Success && + "Simplicial LLT is not factorized, or its factorization failed"); + Scalar signL = Base::m_matrix.diagonal().array().sign().prod(); + return signL * signL; + } }; /** \ingroup SparseCholesky_Module @@ -715,7 +805,36 @@ void factorize(const MatrixType& a) { Base::template factorize<true, true>(a); } /** \returns the determinant of the underlying matrix from the current factorization */ - Scalar determinant() const { return Base::m_diag.prod(); } + Scalar determinant() const { + eigen_assert(Base::m_factorizationIsOk && Base::m_info == Success && + "Simplicial LDLT is not factorized, or its factorization failed"); + return Base::m_diag.prod(); + } + + /** \returns the absolute value of the determinant of the underlying matrix from the current factorization */ + RealScalar absDeterminant() const { + eigen_assert(Base::m_factorizationIsOk && Base::m_info == Success && + "Simplicial LDLT is not factorized, or its factorization failed"); + return numext::abs(Base::m_diag.prod()); + } + + /** \returns the natural log of the absolute value of the determinant of the underlying matrix from the current + * factorization. + * + * Unlike determinant(), this stays finite for the large factorizations where a determinant overflows or underflows. + */ + RealScalar logAbsDeterminant() const { + eigen_assert(Base::m_factorizationIsOk && Base::m_info == Success && + "Simplicial LDLT is not factorized, or its factorization failed"); + return Base::m_diag.cwiseAbs().array().log().sum(); + } + + /** \returns the sign of the determinant of the underlying matrix from the current factorization */ + Scalar signDeterminant() const { + eigen_assert(Base::m_factorizationIsOk && Base::m_info == Success && + "Simplicial LDLT is not factorized, or its factorization failed"); + return Base::m_diag.array().sign().prod(); + } }; /** \deprecated use SimplicialLDLT or class SimplicialLLT @@ -846,6 +965,8 @@ } Scalar determinant() const { + eigen_assert(Base::m_factorizationIsOk && Base::m_info == Success && + "Simplicial Cholesky is not factorized, or its factorization failed"); if (m_LDLT) { return Base::m_diag.prod(); } else {
diff --git a/test/bunchkaufman.cpp b/test/bunchkaufman.cpp index 9d4d4a1..8bb679f 100644 --- a/test/bunchkaufman.cpp +++ b/test/bunchkaufman.cpp
@@ -17,6 +17,7 @@ #include <Eigen/Cholesky> #include <Eigen/QR> #include <Eigen/Eigenvalues> +#include "fp_control.h" #include "solverbase.h" template <typename MatrixType, int UpLo> @@ -272,6 +273,172 @@ } } +// A = Q D Q^*, with Q unitary and D real, is Hermitian with det(A) = prod(D_ii). Mixed signs make A +// indefinite, so the factorization mixes 1x1 and 2x2 blocks of D; drawing the |D_ii| from an annulus keeps +// A well conditioned, hence the inertia -- and with it signDeterminant() -- unambiguous. +template <typename MatrixType> +void bunchkaufman_determinant(Index size) { + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits<Scalar>::Real RealScalar; + typedef Matrix<RealScalar, Dynamic, 1> RealVectorType; + + MatrixType q = MatrixType::Random(size, size).householderQr().householderQ(); + RealVectorType d(size); + for (Index i = 0; i < size; ++i) { + d(i) = internal::random<RealScalar>(RealScalar(1.2), RealScalar(2.8)); + if (internal::random<bool>()) d(i) = -d(i); + } + const MatrixType a = q * d.template cast<Scalar>().asDiagonal() * q.adjoint(); + + const RealScalar det = d.prod(); + const RealScalar logabsdet = d.array().abs().log().sum(); + + BunchKaufman<MatrixType, Lower> bklo(a); + VERIFY(bklo.info() == Success); + check_determinant(bklo, Scalar(det), logabsdet); + // Unlike the other decompositions, this sign is read off the inertia rather than accumulated from a + // product of signs, so it is exact. + VERIFY_IS_EQUAL(bklo.signDeterminant(), Scalar(numext::sign(det))); + + BunchKaufman<MatrixType, Upper> bkup(a); + VERIFY(bkup.info() == Success); + check_determinant(bkup, Scalar(det), logabsdet); + // Unlike the other decompositions, this sign is read off the inertia rather than accumulated from a + // product of signs, so it is exact. + VERIFY_IS_EQUAL(bkup.signDeterminant(), Scalar(numext::sign(det))); +} + +// The determinant of an empty matrix is the empty product, 1. +template <typename MatrixType> +void bunchkaufman_determinant_empty() { + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits<Scalar>::Real RealScalar; + + BunchKaufman<MatrixType> bk{MatrixType(0, 0)}; + VERIFY_IS_EQUAL(bk.determinant(), Scalar(1)); + VERIFY_IS_EQUAL(bk.absDeterminant(), RealScalar(1)); + VERIFY_IS_EQUAL(bk.logAbsDeterminant(), RealScalar(0)); + VERIFY_IS_EQUAL(bk.signDeterminant(), Scalar(1)); +} + +// det(D) can be representable while the block determinants it is built from are not. The blocks below are +// [[s/2, s], [s, s/2]] with det = -3s^2/4 at s = 2^600 and s = 2^-600, one overflowing and one +// underflowing, whose product is exactly 9/16. Multiplying the blocks directly gives inf * 0 = NaN. +void bunchkaufman_determinant_mixed_scale() { + MatrixXd a = MatrixXd::Zero(4, 4); + for (int b = 0; b < 2; ++b) { + const double s = numext::ldexp(1.0, b == 0 ? 600 : -600); + a(2 * b, 2 * b) = a(2 * b + 1, 2 * b + 1) = 0.5 * s; + a(2 * b + 1, 2 * b) = a(2 * b, 2 * b + 1) = s; + } + + // (-3/4)^2 (2^600 2^-600)^2 = 9/16, exactly, both factors being powers of two. + const double det = 0.5625; + + BunchKaufman<MatrixXd, Lower> bklo(a); + VERIFY(bklo.info() == Success); + check_determinant(bklo, det, numext::log(det)); + + BunchKaufman<MatrixXd, Upper> bkup(a); + VERIFY(bkup.info() == Success); + check_determinant(bkup, det, numext::log(det)); +} + +// Bunch-Kaufman selects a 2x2 block only where |d11 d22| <= alpha^2 |d21|^2, alpha = (1+sqrt(17))/8 < 1, so +// det D_k = d11 d22 - |d21|^2 < 0 always. Subnormal |d21| underflows that determinant to zero, leaving the +// log and the sign as the only accessors that can report it; both need det D_k / |d21|^2 = O(1), which +// 1/|d21| overflows on, and the resulting infinity flips the block's inertia. +template <typename MatrixType> +void bunchkaufman_determinant_subnormal_block() { + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits<Scalar>::Real RealScalar; + + if (!subnormalDivisionIsExact<RealScalar>()) { + const char* reason = ScopedFlushToZero::hardwareFlushesSubnormalInputs() ? "the hardware flushes subnormal inputs" + : "the compiler relaxed the division"; + std::cout << "SKIP: bunchkaufman_determinant_subnormal_block needs an environment that divides by subnormals " + "per IEEE 754 (" + << reason << ")." << std::endl; + return; + } + + // Entries are exact integer multiples of the smallest subnormal u, so det A = (k11 k22 - k21^2) u^2 with + // the bracket an exact int, and log|det A| = log|k11 k22 - k21^2| + 2 log u -- a reference that shares no + // expression with the code under test. The off-diagonal stays real to keep the entries exact; the complex + // instantiation still reaches numext::abs() on a subnormal. + const RealScalar u = (std::numeric_limits<RealScalar>::denorm_min)(); + const int k11[] = {0, 1, 2, 100}; + const int k22[] = {0, 1, 2, 100}; + const int k21[] = {1, 4, 8, 300}; + + for (int c = 0; c < 4; ++c) { + const int idet = k11[c] * k22[c] - k21[c] * k21[c]; + VERIFY(idet < 0); + const RealScalar logabsdet = numext::log(RealScalar(-idet)) + RealScalar(2) * numext::log(u); + + MatrixType a(2, 2); + a << Scalar(RealScalar(k11[c]) * u), Scalar(RealScalar(k21[c]) * u), Scalar(RealScalar(k21[c]) * u), + Scalar(RealScalar(k22[c]) * u); + // Flushing subnormal results to zero empties the matrix even where the division probe passed. + if (numext::is_exactly_zero(numext::abs(a.coeff(1, 0)))) return; + + BunchKaufman<MatrixType, Lower> bk(a); + // c == 0 has d11 = d22 = 0, where unblocked()'s own scaled determinant is 0*inf and it reports + // NumericalIssue; the accessors below are exact either way. + if (k11[c] != 0) VERIFY(bk.info() == Success); + + VERIFY_IS_EQUAL(bk.signDeterminant(), Scalar(-1)); + VERIFY_IS_APPROX(bk.logAbsDeterminant(), logabsdet); + // det A is of order u^2, so zero is all these two can report. + VERIFY_IS_EQUAL(bk.absDeterminant(), RealScalar(0)); + VERIFY_IS_EQUAL(bk.determinant(), Scalar(0)); + VERIFY(!bk.isPositive()); + VERIFY(!bk.isNegative()); + } +} + +// The criterion bounds |d11| against alpha|d21| but bounds |d22| only against the largest entry of its own +// row, which can dwarf |d21|. So d22/d21 can overflow, and the scaled determinant comes out +inf, or 0*inf +// = NaN where d11 is zero; either one counts the block as definite with the sign of its trace, and the +// trace is positive here. The criterion puts the true value below 1, which is what rejects both. +void bunchkaufman_inertia_wide_2x2_block() { + // The 2x2 block is (d11, d21, d22) = (a00, 1e-10, 1e299), determinant a00*1e299 - 1e-20 < 0 in both + // rows below, so it contributes one eigenvalue of each sign. info() is NumericalIssue: unblocked() + // forms the same product through 1/d21 and NaNs the trailing update, so the inertia is the only part + // of the factorization that is meaningful for these. + for (double a00 : {0.0, 1e-320}) { + MatrixXd a = MatrixXd::Zero(4, 4); + a(0, 0) = a00; + a(1, 0) = a(0, 1) = 1e-10; + a(1, 1) = 1e299; + a(3, 1) = a(1, 3) = 2e299; + + BunchKaufman<MatrixXd> bk(a); + VERIFY(!bk.isPositive()); + VERIFY(!bk.isNegative()); + } +} + +// The scaling above divides by numext::abs(d21), so it relies on that magnitude not underflowing: for +// |3u + 4iu| Eigen's hypot scales by the larger component instead of summing squares, giving 5u exactly +// where u^2 would be zero. Not template code, because building the off-diagonal needs a complex literal. +void bunchkaufman_determinant_subnormal_block_complex() { + typedef std::complex<double> Scalar; + + if (!subnormalDivisionIsExact<double>()) return; + + const double u = (std::numeric_limits<double>::denorm_min)(); + MatrixXcd a(2, 2); + a << Scalar(u, 0), Scalar(3 * u, 4 * u), Scalar(3 * u, -4 * u), Scalar(u, 0); + BunchKaufman<MatrixXcd, Lower> bk(a); + + // det = u^2 - |3u + 4iu|^2 = -24 u^2. + VERIFY_IS_EQUAL(bk.signDeterminant(), Scalar(-1)); + VERIFY_IS_APPROX(bk.logAbsDeterminant(), numext::log(24.0) + 2.0 * numext::log(u)); + VERIFY(!bk.isPositive()); + VERIFY(!bk.isNegative()); +} + template <typename MatrixType> void bunchkaufman_verify_assert() { MatrixType tmp; @@ -286,6 +453,10 @@ VERIFY_RAISES_ASSERT(bk.matrixLDLT()) VERIFY_RAISES_ASSERT(bk.reconstructedMatrix()) VERIFY_RAISES_ASSERT(bk.solve(tmp)) + VERIFY_RAISES_ASSERT(bk.determinant()) + VERIFY_RAISES_ASSERT(bk.absDeterminant()) + VERIFY_RAISES_ASSERT(bk.logAbsDeterminant()) + VERIFY_RAISES_ASSERT(bk.signDeterminant()) } // Build a random Hermitian (real symmetric) indefinite matrix of the same type/size as `m`. @@ -397,6 +568,12 @@ CALL_SUBTEST_6(bunchkaufman(MatrixXcd(s, s))); TEST_SET_BUT_UNUSED_VARIABLE(s); + // Bounded so that the determinant itself, not just its logarithm, stays in range. + s = internal::random<int>(1, 30); + CALL_SUBTEST_5(bunchkaufman_determinant<MatrixXd>(s)); + CALL_SUBTEST_6(bunchkaufman_determinant<MatrixXcd>(s)); + TEST_SET_BUT_UNUSED_VARIABLE(s); + s = internal::random<int>(2, EIGEN_TEST_MAX_SIZE); CALL_SUBTEST_5(bunchkaufman_inertia_and_conditioning<MatrixXd>(s)); s = internal::random<int>(2, EIGEN_TEST_MAX_SIZE / 2); @@ -417,6 +594,17 @@ // Empty-matrix edge case. CALL_SUBTEST_5(bunchkaufman(MatrixXd(0, 0))); + CALL_SUBTEST_5(bunchkaufman_determinant_empty<MatrixXd>()); + + // Subnormal 2x2 block: the determinant underflows, its log and sign do not. + CALL_SUBTEST_5(bunchkaufman_determinant_subnormal_block<MatrixXd>()); + CALL_SUBTEST_6(bunchkaufman_determinant_subnormal_block<MatrixXcd>()); + CALL_SUBTEST_6(bunchkaufman_determinant_subnormal_block_complex()); + CALL_SUBTEST_8(bunchkaufman_determinant_subnormal_block<MatrixXf>()); + CALL_SUBTEST_5(bunchkaufman_inertia_wide_2x2_block()); + + // Mixed-scale 2x2 blocks: the block determinants leave the representable range, their product does not. + CALL_SUBTEST_5(bunchkaufman_determinant_mixed_scale()); // Problem-size constructors. CALL_SUBTEST_8(BunchKaufman<MatrixXf>(10));
diff --git a/test/cholesky.cpp b/test/cholesky.cpp index 205a305..08fbb6c 100644 --- a/test/cholesky.cpp +++ b/test/cholesky.cpp
@@ -604,6 +604,128 @@ } } +// A = Q D Q^*, with Q unitary and D real, is Hermitian with det(A) = prod(D_ii). +// Drawing the |D_ii| from an annulus keeps A well conditioned and log|det(A)| away from zero. +template <typename MatrixType> +void cholesky_determinant(Index size) { + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits<Scalar>::Real RealScalar; + typedef Matrix<RealScalar, Dynamic, 1> RealVectorType; + + MatrixType q = MatrixType::Random(size, size).householderQr().householderQ(); + RealVectorType d(size); + for (Index i = 0; i < size; ++i) d(i) = internal::random<RealScalar>(RealScalar(1.2), RealScalar(2.8)); + + const RealScalar absdet = d.prod(); + const RealScalar logabsdet = d.array().log().sum(); + const MatrixType spd = q * d.template cast<Scalar>().asDiagonal() * q.adjoint(); + + LLT<MatrixType, Lower> lltlo(spd); + VERIFY(lltlo.info() == Success); + check_determinant(lltlo, Scalar(absdet), logabsdet); + + LLT<MatrixType, Upper> lltup(spd); + VERIFY(lltup.info() == Success); + check_determinant(lltup, Scalar(absdet), logabsdet); + + // Negating D leaves |det(A)| alone and makes A negative definite, so det(A) picks up a factor (-1)^n. + // LDLT's D is purely diagonal and so does not cover the indefinite case; BunchKaufman does, and + // test/bunchkaufman.cpp checks the same identities there. + d = -d; + const RealScalar det = d.prod(); + const MatrixType negdef = q * d.template cast<Scalar>().asDiagonal() * q.adjoint(); + + LDLT<MatrixType, Lower> ldltlo(negdef); + VERIFY(ldltlo.info() == Success); + check_determinant(ldltlo, Scalar(det), logabsdet); + + LDLT<MatrixType, Upper> ldltup(negdef); + VERIFY(ldltup.info() == Success); + check_determinant(ldltup, Scalar(det), logabsdet); +} + +// The determinant of an empty matrix is the empty product, 1. +template <typename MatrixType> +void cholesky_determinant_empty() { + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits<Scalar>::Real RealScalar; + + MatrixType empty(0, 0); + + LLT<MatrixType> llt(empty); + VERIFY_IS_EQUAL(llt.determinant(), Scalar(1)); + VERIFY_IS_EQUAL(llt.absDeterminant(), RealScalar(1)); + VERIFY_IS_EQUAL(llt.logAbsDeterminant(), RealScalar(0)); + VERIFY_IS_EQUAL(llt.signDeterminant(), Scalar(1)); + + LDLT<MatrixType> ldlt(empty); + VERIFY_IS_EQUAL(ldlt.determinant(), Scalar(1)); + VERIFY_IS_EQUAL(ldlt.absDeterminant(), RealScalar(1)); + VERIFY_IS_EQUAL(ldlt.logAbsDeterminant(), RealScalar(0)); + VERIFY_IS_EQUAL(ldlt.signDeterminant(), Scalar(1)); +} + +// logAbsDeterminant() exists to survive the range where the determinant itself does not: with n = 200 and +// a diagonal of 10^4, det = 10^800 overflows every supported float type while log|det| = 800 log 10 does not. +template <typename MatrixType> +void cholesky_determinant_overflow() { + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits<Scalar>::Real RealScalar; + + const Index size = 200; + for (bool overflow : {true, false}) { + const RealScalar scale = overflow ? RealScalar(1e4) : RealScalar(1e-4); + const MatrixType a = MatrixType::Identity(size, size) * Scalar(scale); + const RealScalar logabsdet = RealScalar(size) * numext::log(scale); + + LLT<MatrixType> llt(a); + VERIFY(llt.info() == Success); + VERIFY(determinant_out_of_range(llt.absDeterminant(), overflow)); + VERIFY_IS_APPROX(llt.logAbsDeterminant(), logabsdet); + VERIFY_IS_EQUAL(llt.signDeterminant(), Scalar(1)); + + LDLT<MatrixType> ldlt(a); + VERIFY(ldlt.info() == Success); + VERIFY(determinant_out_of_range(ldlt.absDeterminant(), overflow)); + VERIFY_IS_APPROX(ldlt.logAbsDeterminant(), logabsdet); + VERIFY_IS_EQUAL(ldlt.signDeterminant(), Scalar(1)); + } +} + +// A failed factorization does not represent the input, so the four accessors assert rather than answer from +// it: for [[0,1],[1,0]] LDLT reports NumericalIssue with D = 0, where det = -1. m_isInitialized alone does +// not catch that, since compute() sets it either way. +template <typename MatrixType> +void cholesky_determinant_failed_factorization(Index size) { + eigen_assert(size >= 2); + + // A zero diagonal with non-zero off-diagonal entries makes the first pivot invalid while the matrix is + // not; both factorizations give up on it. + const MatrixType indefinite = MatrixType::Ones(size, size) - MatrixType::Identity(size, size); + + LDLT<MatrixType, Lower> ldltlo(indefinite); + VERIFY(ldltlo.info() == NumericalIssue); + VERIFY(!ldltlo.reconstructedMatrix().isApprox(indefinite)); + VERIFY_RAISES_ASSERT(ldltlo.determinant()) + VERIFY_RAISES_ASSERT(ldltlo.absDeterminant()) + VERIFY_RAISES_ASSERT(ldltlo.logAbsDeterminant()) + VERIFY_RAISES_ASSERT(ldltlo.signDeterminant()) + + LDLT<MatrixType, Upper> ldltup(indefinite); + VERIFY(ldltup.info() == NumericalIssue); + VERIFY_RAISES_ASSERT(ldltup.determinant()) + VERIFY_RAISES_ASSERT(ldltup.absDeterminant()) + VERIFY_RAISES_ASSERT(ldltup.logAbsDeterminant()) + VERIFY_RAISES_ASSERT(ldltup.signDeterminant()) + + LLT<MatrixType, Lower> lltlo(indefinite); + VERIFY(lltlo.info() == NumericalIssue); + VERIFY_RAISES_ASSERT(lltlo.determinant()) + VERIFY_RAISES_ASSERT(lltlo.absDeterminant()) + VERIFY_RAISES_ASSERT(lltlo.logAbsDeterminant()) + VERIFY_RAISES_ASSERT(lltlo.signDeterminant()) +} + template <typename MatrixType> void cholesky_verify_assert() { MatrixType tmp; @@ -615,6 +737,10 @@ VERIFY_RAISES_ASSERT(llt.transpose().solve(tmp)) VERIFY_RAISES_ASSERT(llt.adjoint().solve(tmp)) VERIFY_RAISES_ASSERT(llt.solveInPlace(tmp)) + VERIFY_RAISES_ASSERT(llt.determinant()) + VERIFY_RAISES_ASSERT(llt.absDeterminant()) + VERIFY_RAISES_ASSERT(llt.logAbsDeterminant()) + VERIFY_RAISES_ASSERT(llt.signDeterminant()) LDLT<MatrixType> ldlt; VERIFY_RAISES_ASSERT(ldlt.matrixL()) @@ -626,6 +752,10 @@ VERIFY_RAISES_ASSERT(ldlt.transpose().solve(tmp)) VERIFY_RAISES_ASSERT(ldlt.adjoint().solve(tmp)) VERIFY_RAISES_ASSERT(ldlt.solveInPlace(tmp)) + VERIFY_RAISES_ASSERT(ldlt.determinant()) + VERIFY_RAISES_ASSERT(ldlt.absDeterminant()) + VERIFY_RAISES_ASSERT(ldlt.logAbsDeterminant()) + VERIFY_RAISES_ASSERT(ldlt.signDeterminant()) } // Test Cholesky decomposition at blocking and vectorization boundaries. @@ -731,9 +861,21 @@ CALL_SUBTEST_2(cholesky_ldlt_rankupdate_zero_components(MatrixXd(s, s))); CALL_SUBTEST_6(cholesky_ldlt_rankupdate_zero_components(MatrixXcd(s, s))); TEST_SET_BUT_UNUSED_VARIABLE(s); + + // Bounded so that the determinant itself, not just its logarithm, stays in range. + s = internal::random<int>(1, 30); + CALL_SUBTEST_2(cholesky_determinant<MatrixXd>(s)); + CALL_SUBTEST_8(cholesky_determinant<MatrixXf>(s)); + CALL_SUBTEST_6(cholesky_determinant<MatrixXcd>(s)); + TEST_SET_BUT_UNUSED_VARIABLE(s); } // empty matrix, regression test for Bug 785: CALL_SUBTEST_2(cholesky(MatrixXd(0, 0))); + CALL_SUBTEST_2(cholesky_determinant_empty<MatrixXd>()); + CALL_SUBTEST_8(cholesky_determinant_overflow<MatrixXf>()); + CALL_SUBTEST_3(cholesky_determinant_failed_factorization<Matrix2d>(2)); + CALL_SUBTEST_2(cholesky_determinant_failed_factorization<MatrixXd>(internal::random<int>(2, 20))); + CALL_SUBTEST_6(cholesky_determinant_failed_factorization<MatrixXcd>(internal::random<int>(2, 20))); // This does not work yet: // CALL_SUBTEST_2( cholesky(Matrix<double,0,0>()) );
diff --git a/test/lu.cpp b/test/lu.cpp index 17675d1..b5e17f0 100644 --- a/test/lu.cpp +++ b/test/lu.cpp
@@ -10,6 +10,7 @@ #include "main.h" #include <Eigen/LU> +#include <Eigen/QR> #include "solverbase.h" using namespace std; @@ -167,16 +168,20 @@ // Regression test: FullPivLU took the maximum of an empty column-sum vector when computing its l1 norm, // so it could not be constructed at all from a matrix with zero columns. PartialPivLU already guarded the -// same reduction. +// same reduction. The determinant of an empty matrix is the empty product, 1. template <typename MatrixType> void lu_empty() { typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits<Scalar>::Real RealScalar; const Index n = 5; FullPivLU<MatrixType> lu{MatrixType(0, 0)}; VERIFY_IS_EQUAL(lu.rank(), Index(0)); VERIFY(lu.isInvertible()); VERIFY_IS_EQUAL(lu.determinant(), Scalar(1)); + VERIFY_IS_EQUAL(lu.absDeterminant(), RealScalar(1)); + VERIFY_IS_EQUAL(lu.logAbsDeterminant(), RealScalar(0)); + VERIFY_IS_EQUAL(lu.signDeterminant(), Scalar(1)); lu.compute(MatrixType(0, n)); VERIFY_IS_EQUAL(lu.rank(), Index(0)); @@ -190,6 +195,98 @@ PartialPivLU<MatrixType> plu{MatrixType(0, 0)}; VERIFY_IS_EQUAL(plu.determinant(), Scalar(1)); + VERIFY_IS_EQUAL(plu.absDeterminant(), RealScalar(1)); + VERIFY_IS_EQUAL(plu.logAbsDeterminant(), RealScalar(0)); + VERIFY_IS_EQUAL(plu.signDeterminant(), Scalar(1)); +} + +// A = Q D Q^*, with Q unitary, has det(A) = det(D) because det(Q) det(Q^*) = |det(Q)|^2 = 1. +template <typename MatrixType> +void lu_determinant(Index size) { + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits<Scalar>::Real RealScalar; + + MatrixType d = MatrixType::Zero(size, size); + setRandomWellConditionedDiagonal(d); + const MatrixType q = MatrixType::Random(size, size).householderQr().householderQ(); + const MatrixType a = q * d * q.adjoint(); + + const Scalar det = d.diagonal().prod(); + const RealScalar logabsdet = d.diagonal().cwiseAbs().array().log().sum(); + + check_determinant(FullPivLU<MatrixType>(a), det, logabsdet); + check_determinant(PartialPivLU<MatrixType>(a), det, logabsdet); +} + +// logAbsDeterminant() exists to survive the range where the determinant itself does not: with n = 200 and +// a diagonal of 10^4, det = 10^800 overflows every supported float type while log|det| = 800 log 10 does not. +template <typename MatrixType> +void lu_determinant_overflow() { + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits<Scalar>::Real RealScalar; + + const Index size = 200; + for (bool overflow : {true, false}) { + const RealScalar scale = overflow ? RealScalar(1e4) : RealScalar(1e-4); + const MatrixType a = MatrixType::Identity(size, size) * Scalar(scale); + const RealScalar logabsdet = RealScalar(size) * numext::log(scale); + + PartialPivLU<MatrixType> plu(a); + VERIFY(determinant_out_of_range(plu.absDeterminant(), overflow)); + VERIFY_IS_APPROX(plu.logAbsDeterminant(), logabsdet); + VERIFY_IS_EQUAL(plu.signDeterminant(), Scalar(1)); + + FullPivLU<MatrixType> lu(a); + VERIFY(determinant_out_of_range(lu.absDeterminant(), overflow)); + VERIFY_IS_APPROX(lu.logAbsDeterminant(), logabsdet); + VERIFY_IS_EQUAL(lu.signDeterminant(), Scalar(1)); + } +} + +// A rank-deficient decomposition has no determinant worth reporting, and FullPivLU is the one LU that +// knows it. The three accessors that can express that gate on rank, so they agree with the rank-revealing +// QR decompositions; determinant() keeps its documented behaviour of returning the pivot product. +template <typename MatrixType> +void lu_determinant_rank_deficient(Index size) { + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits<Scalar>::Real RealScalar; + + const Index rank = internal::random<Index>(1, size - 1); + MatrixType a(size, size); + createRandomPIMatrixOfRank(rank, size, size, a); + + // The generated singular values are 0 or 1, so any threshold well inside that gap recovers the rank. + const RealScalar threshold(0.01); + + FullPivLU<MatrixType> lu; + lu.setThreshold(threshold); + lu.compute(a); + VERIFY_IS_EQUAL(lu.rank(), rank); + VERIFY(!lu.isInvertible()); + VERIFY_IS_EQUAL(lu.absDeterminant(), RealScalar(0)); + VERIFY_IS_EQUAL(lu.logAbsDeterminant(), -NumTraits<RealScalar>::infinity()); + VERIFY_IS_EQUAL(lu.signDeterminant(), Scalar(0)); + + ColPivHouseholderQR<MatrixType> qr; + qr.setThreshold(threshold); + qr.compute(a); + VERIFY_IS_EQUAL(qr.rank(), rank); + VERIFY_IS_EQUAL(qr.absDeterminant(), lu.absDeterminant()); + VERIFY_IS_EQUAL(qr.logAbsDeterminant(), lu.logAbsDeterminant()); + VERIFY_IS_EQUAL(qr.signDeterminant(), lu.signDeterminant()); + + // determinant() is deliberately not gated. Pin that on a matrix whose smallest pivot is below the + // threshold but nonzero, where the two answers are visibly different rather than both roundoff. + MatrixType b = MatrixType::Identity(size, size); + b(size - 1, size - 1) = Scalar(RealScalar(1e-30)); + FullPivLU<MatrixType> blu; + blu.setThreshold(RealScalar(1e-3)); + blu.compute(b); + VERIFY_IS_EQUAL(blu.rank(), size - 1); + VERIFY_IS_EQUAL(blu.absDeterminant(), RealScalar(0)); + VERIFY_IS_EQUAL(blu.logAbsDeterminant(), -NumTraits<RealScalar>::infinity()); + VERIFY_IS_EQUAL(blu.signDeterminant(), Scalar(0)); + VERIFY_IS_APPROX(numext::abs(blu.determinant()), RealScalar(1e-30)); } template <typename MatrixType> @@ -206,6 +303,9 @@ VERIFY_RAISES_ASSERT(lu.transpose().solve(tmp)) VERIFY_RAISES_ASSERT(lu.adjoint().solve(tmp)) VERIFY_RAISES_ASSERT(lu.determinant()) + VERIFY_RAISES_ASSERT(lu.absDeterminant()) + VERIFY_RAISES_ASSERT(lu.logAbsDeterminant()) + VERIFY_RAISES_ASSERT(lu.signDeterminant()) VERIFY_RAISES_ASSERT(lu.rank()) VERIFY_RAISES_ASSERT(lu.dimensionOfKernel()) VERIFY_RAISES_ASSERT(lu.isInjective()) @@ -220,6 +320,9 @@ VERIFY_RAISES_ASSERT(plu.transpose().solve(tmp)) VERIFY_RAISES_ASSERT(plu.adjoint().solve(tmp)) VERIFY_RAISES_ASSERT(plu.determinant()) + VERIFY_RAISES_ASSERT(plu.absDeterminant()) + VERIFY_RAISES_ASSERT(plu.logAbsDeterminant()) + VERIFY_RAISES_ASSERT(plu.signDeterminant()) VERIFY_RAISES_ASSERT(plu.inverse()) } @@ -324,22 +427,29 @@ CALL_SUBTEST_3(lu_non_invertible<MatrixXf>()); CALL_SUBTEST_3(lu_invertible<MatrixXf>()); CALL_SUBTEST_3(lu_verify_assert<MatrixXf>()); + CALL_SUBTEST_3(lu_determinant<MatrixXf>(internal::random<int>(1, 30))); CALL_SUBTEST_4(lu_non_invertible<MatrixXd>()); CALL_SUBTEST_4(lu_invertible<MatrixXd>()); CALL_SUBTEST_4(lu_partial_piv<MatrixXd>(internal::random<int>(1, EIGEN_TEST_MAX_SIZE))); CALL_SUBTEST_4(lu_verify_assert<MatrixXd>()); CALL_SUBTEST_4(lu_empty<MatrixXd>()); + CALL_SUBTEST_4(lu_determinant<MatrixXd>(internal::random<int>(1, 30))); + CALL_SUBTEST_4(lu_determinant_rank_deficient<MatrixXd>(internal::random<int>(2, 30))); + CALL_SUBTEST_3(lu_determinant_overflow<MatrixXf>()); CALL_SUBTEST_5(lu_non_invertible<MatrixXcf>()); CALL_SUBTEST_5(lu_invertible<MatrixXcf>()); CALL_SUBTEST_5(lu_verify_assert<MatrixXcf>()); + CALL_SUBTEST_5(lu_determinant<MatrixXcf>(internal::random<int>(1, 30))); CALL_SUBTEST_6(lu_non_invertible<MatrixXcd>()); CALL_SUBTEST_6(lu_invertible<MatrixXcd>()); CALL_SUBTEST_6(lu_partial_piv<MatrixXcd>(internal::random<int>(1, EIGEN_TEST_MAX_SIZE))); CALL_SUBTEST_6(lu_verify_assert<MatrixXcd>()); CALL_SUBTEST_6(lu_empty<MatrixXcd>()); + CALL_SUBTEST_6(lu_determinant<MatrixXcd>(internal::random<int>(1, 30))); + CALL_SUBTEST_6(lu_determinant_rank_deficient<MatrixXcd>(internal::random<int>(2, 30))); CALL_SUBTEST_7((lu_non_invertible<Matrix<float, Dynamic, 16> >()));
diff --git a/test/simplicial_cholesky.cpp b/test/simplicial_cholesky.cpp index a98631f..0abd429 100644 --- a/test/simplicial_cholesky.cpp +++ b/test/simplicial_cholesky.cpp
@@ -10,6 +10,61 @@ #include "sparse_solver.h" +template <typename Solver> +void verify_determinant_asserts(const Solver& solver) { + VERIFY_RAISES_ASSERT(solver.determinant()); + VERIFY_RAISES_ASSERT(solver.absDeterminant()); + VERIFY_RAISES_ASSERT(solver.logAbsDeterminant()); + VERIFY_RAISES_ASSERT(solver.signDeterminant()); +} + +// factorize_preordered() breaks out of its column loop on a bad pivot but still sets m_factorizationIsOk, +// leaving the tail of D -- and of L's diagonal -- unwritten. The determinant accessors read exactly those, +// so they check info() too, as the dense LLT ones do. +template <typename SparseMatrixType> +void test_simplicial_cholesky_determinant_asserts() { + typedef typename SparseMatrixType::Scalar Scalar; + const Index n = 6; + + // A zero pivot stops all four classes; SimplicialLLT alone also rejects a negative one. + SparseMatrixType zero_pivot(n, n), negative_pivot(n, n); + for (Index i = 0; i < n; ++i) { + zero_pivot.insert(i, i) = (i == 0) ? Scalar(0) : Scalar(1); + negative_pivot.insert(i, i) = (i == n - 1) ? Scalar(-1) : Scalar(1); + } + zero_pivot.makeCompressed(); + negative_pivot.makeCompressed(); + + SimplicialLLT<SparseMatrixType> llt(zero_pivot); + VERIFY(llt.info() == NumericalIssue); + verify_determinant_asserts(llt); + + SimplicialLLT<SparseMatrixType> indefinite_llt(negative_pivot); + VERIFY(indefinite_llt.info() == NumericalIssue); + verify_determinant_asserts(indefinite_llt); + + SimplicialLDLT<SparseMatrixType> ldlt(zero_pivot); + VERIFY(ldlt.info() == NumericalIssue); + verify_determinant_asserts(ldlt); + + SimplicialNonHermitianLLT<SparseMatrixType> nhllt(zero_pivot); + VERIFY(nhllt.info() == NumericalIssue); + verify_determinant_asserts(nhllt); + + SimplicialNonHermitianLDLT<SparseMatrixType> nhldlt(zero_pivot); + VERIFY(nhldlt.info() == NumericalIssue); + verify_determinant_asserts(nhldlt); + + // The deprecated SimplicialCholesky exposes determinant() alone, and reads the same tail. + SimplicialCholesky<SparseMatrixType> deprecated(zero_pivot); + VERIFY(deprecated.info() == NumericalIssue); + VERIFY_RAISES_ASSERT(deprecated.determinant()); + + // info() alone would not reject an unfactorized decomposition: it is constructed Success. + SimplicialLDLT<SparseMatrixType> unfactorized; + verify_determinant_asserts(unfactorized); +} + template <typename T, typename I_, int flag> void test_simplicial_cholesky_T() { typedef SparseMatrix<T, flag, I_> SparseMatrixType; @@ -50,6 +105,17 @@ check_sparse_nonhermitian_determinant(nhldlt_colmajor_lower_amd); check_sparse_nonhermitian_determinant(nhldlt_colmajor_upper_amd); + check_sparse_spd_log_abs_determinant(llt_colmajor_lower_amd); + check_sparse_spd_log_abs_determinant(llt_colmajor_upper_amd); + check_sparse_spd_log_abs_determinant(ldlt_colmajor_lower_amd); + check_sparse_spd_log_abs_determinant(ldlt_colmajor_upper_amd); + check_sparse_nonhermitian_log_abs_determinant(nhllt_colmajor_lower_amd); + check_sparse_nonhermitian_log_abs_determinant(nhllt_colmajor_upper_amd); + check_sparse_nonhermitian_log_abs_determinant(nhldlt_colmajor_lower_amd); + check_sparse_nonhermitian_log_abs_determinant(nhldlt_colmajor_upper_amd); + + test_simplicial_cholesky_determinant_asserts<SparseMatrixType>(); + check_sparse_spd_solving(ldlt_colmajor_lower_nat, (std::min)(300, EIGEN_TEST_MAX_SIZE), 1000); check_sparse_spd_solving(ldlt_colmajor_upper_nat, (std::min)(300, EIGEN_TEST_MAX_SIZE), 1000); check_sparse_nonhermitian_solving(nhldlt_colmajor_lower_nat, (std::min)(300, EIGEN_TEST_MAX_SIZE), 1000);
diff --git a/test/solverbase.h b/test/solverbase.h index 4429311..dddc3fc 100644 --- a/test/solverbase.h +++ b/test/solverbase.h
@@ -39,4 +39,25 @@ VERIFY_IS_APPROX(matrix * m2, matrix * solver_solution); } +// Checks the four determinant accessors of a decomposition against a reference determinant \a det and a +// reference \a logabsdet = log|det| formed independently of it. Callers must keep \a det itself in range. +template <typename SolverType, typename Scalar> +void check_determinant(const SolverType& solver, const Scalar& det, const typename NumTraits<Scalar>::Real& logabsdet) { + typedef typename NumTraits<Scalar>::Real RealScalar; + VERIFY_IS_APPROX(solver.determinant(), det); + VERIFY_IS_APPROX(solver.absDeterminant(), numext::abs(det)); + // log|det| passes through zero, where a relative comparison says nothing; bound the error absolutely. + VERIFY_IS_MUCH_SMALLER_THAN(solver.logAbsDeterminant() - logabsdet, RealScalar(1)); + VERIFY_IS_APPROX(solver.signDeterminant(), numext::sign(det)); +} + +// True when |det| has left the representable range in the direction it was expected to: exactly infinity +// where the determinant overflowed, exactly zero where it underflowed. logAbsDeterminant() is still +// meaningful there; determinant() and absDeterminant() are not. A NaN is neither, and testing !isfinite() +// would accept one -- with it any inf - inf or 0 * inf artifact. +template <typename RealScalar> +bool determinant_out_of_range(const RealScalar& absdet, bool overflow) { + return overflow ? numext::equal_strict(absdet, NumTraits<RealScalar>::infinity()) : numext::is_exactly_zero(absdet); +} + #endif // TEST_SOLVERBASE_H
diff --git a/test/sparse_solver.h b/test/sparse_solver.h index 212be7b..7290266 100644 --- a/test/sparse_solver.h +++ b/test/sparse_solver.h
@@ -343,6 +343,28 @@ } template <typename Solver, typename DenseMat> +void check_sparse_log_abs_determinant(Solver& solver, const typename Solver::MatrixType& A, const DenseMat& dA) { + typedef typename Solver::MatrixType Mat; + typedef typename Mat::Scalar Scalar; + typedef typename NumTraits<Scalar>::Real RealScalar; + + solver.compute(A); + if (solver.info() != Success) { + std::cerr << "WARNING | sparse solver testing: factorization failed (check_sparse_log_abs_determinant)\n"; + return; + } + + Scalar refDet = dA.determinant(); + RealScalar refAbsDet = numext::abs(refDet); + VERIFY_IS_APPROX(refAbsDet, solver.absDeterminant()); + VERIFY_IS_APPROX(numext::sign(refDet), solver.signDeterminant()); + // log|det| crosses zero, so the meaningful bound here is absolute once |log|det|| drops below one. + RealScalar refLogAbsDet = numext::log(refAbsDet); + VERIFY_IS_MUCH_SMALLER_THAN(solver.logAbsDeterminant() - refLogAbsDet, + numext::maxi(RealScalar(1), numext::abs(refLogAbsDet))); +} + +template <typename Solver, typename DenseMat> int generate_sparse_spd_problem(Solver&, typename Solver::MatrixType& A, typename Solver::MatrixType& halfA, DenseMat& dA, int maxSize = 300) { typedef typename Solver::MatrixType Mat; @@ -492,6 +514,23 @@ } } +template <typename Solver> +void check_sparse_spd_log_abs_determinant(Solver& solver) { + typedef typename Solver::MatrixType Mat; + typedef typename Mat::Scalar Scalar; + typedef Matrix<Scalar, Dynamic, Dynamic> DenseMatrix; + + // generate the problem + Mat A, halfA; + DenseMatrix dA; + generate_sparse_spd_problem(solver, A, halfA, dA, 30); + + for (int i = 0; i < g_repeat; i++) { + check_sparse_log_abs_determinant(solver, A, dA); + check_sparse_log_abs_determinant(solver, halfA, dA); + } +} + template <typename Solver, typename DenseMat> int generate_sparse_nonhermitian_problem(Solver&, typename Solver::MatrixType& A, typename Solver::MatrixType& halfA, DenseMat& dA, int maxSize = 300) { @@ -583,6 +622,23 @@ } template <typename Solver> +void check_sparse_nonhermitian_log_abs_determinant(Solver& solver) { + typedef typename Solver::MatrixType Mat; + typedef typename Mat::Scalar Scalar; + typedef Matrix<Scalar, Dynamic, Dynamic> DenseMatrix; + + // generate the problem + Mat A, halfA; + DenseMatrix dA; + generate_sparse_nonhermitian_problem(solver, A, halfA, dA, 30); + + for (int i = 0; i < g_repeat; i++) { + check_sparse_log_abs_determinant(solver, A, dA); + check_sparse_log_abs_determinant(solver, halfA, dA); + } +} + +template <typename Solver> void check_sparse_zero_matrix(Solver& solver) { typedef typename Solver::MatrixType Mat;