[prev in list] [next in list] [prev in thread] [next in thread]
List: kde-commits
Subject: kdesupport/eigen2
From: Benoît Jacob <jacob.benoit.1 () gmail ! com>
Date: 2009-04-01 0:21:17
Message-ID: 1238545277.305746.12683.nullmailer () svn ! kde ! org
[Download RAW message or body]
SVN commit 947641 by bjacob:
More Cholesky fixes.
* Cholesky decs are NOT rank revealing so remove all the rank/isPositiveDefinite etc \
stuff.
* fix bug in LLT: s/return/continue/
* introduce machine_epsilon constants, they are actually needed for Higman's formula \
determining the cutoff in Cholesky. Btw fix the page reference to his book (chat \
with Keir).
* solve methods always return true, since this isn't a rank revealing dec. \
Actually... they already did always return true!! Now it's explicit.
* updated dox and unit-test
M +10 -30 Eigen/src/Cholesky/LDLT.h
M +17 -19 Eigen/src/Cholesky/LLT.h
M +8 -0 Eigen/src/Core/MathFunctions.h
M +0 -50 test/cholesky.cpp
--- trunk/kdesupport/eigen2/Eigen/src/Cholesky/LDLT.h #947640:947641
@@ -43,6 +43,9 @@
* zeros in the bottom right rank(A) - n submatrix. Avoiding the square root
* on D also stabilizes the computation.
*
+ * Remember that Cholesky decompositions are not rank-revealing. Also, do not use \
a Cholesky decomposition to determine + * whether a system of equations has a \
solution. + *
* \sa MatrixBase::ldlt(), class LLT
*/
/* THIS PART OF THE DOX IS CURRENTLY DISABLED BECAUSE INACCURATE BECAUSE OF BUG IN \
THE DECOMPOSITION CODE @@ -88,25 +91,6 @@
/** \returns true if the matrix is negative (semidefinite) */
inline bool isNegative(void) const { return m_sign == -1; }
- /** \returns true if the matrix is invertible */
- inline bool isInvertible(void) const { return m_rank == m_matrix.rows(); }
-
- /** \returns true if the matrix is positive definite */
- inline bool isPositiveDefinite(void) const { return isPositive() && \
isInvertible(); }
-
- /** \returns true if the matrix is negative definite */
- inline bool isNegativeDefinite(void) const { return isNegative() && \
isInvertible(); }
-
- /** \returns the rank of the matrix of which *this is the LDLT decomposition.
- *
- * \note This is computed at the time of the construction of the LDLT \
decomposition. This
- * method does not perform any further computation.
- */
- inline int rank() const
- {
- return m_rank;
- }
-
template<typename RhsDerived, typename ResDerived>
bool solve(const MatrixBase<RhsDerived> &b, MatrixBase<ResDerived> *result) \
const;
@@ -125,7 +109,7 @@
MatrixType m_matrix;
IntColVectorType m_p;
IntColVectorType m_transpositions;
- int m_rank, m_sign;
+ int m_sign;
};
/** Compute / recompute the LDLT decomposition A = L D L^* = U^* D U of \a matrix
@@ -135,7 +119,6 @@
{
ei_assert(a.rows()==a.cols());
const int size = a.rows();
- m_rank = size;
m_matrix = a;
@@ -168,8 +151,8 @@
// to the largest overall, the algorithm bails. This cutoff is suggested
// in "Analysis of the Cholesky Decomposition of a Semi-definite Matrix" by
// Nicholas J. Higham. Also see "Accuracy and Stability of Numerical
- // Algorithms" page 208, also by Higham.
- cutoff = ei_abs(precision<RealScalar>() * size * biggest_in_corner);
+ // Algorithms" page 217, also by Higham.
+ cutoff = ei_abs(machine_epsilon<Scalar>() * size * biggest_in_corner);
m_sign = ei_real(m_matrix.diagonal().coeff(index_of_biggest_in_corner)) > 0 ? \
1 : -1; }
@@ -178,7 +161,6 @@
if(biggest_in_corner < cutoff)
{
for(int i = j; i < size; i++) m_transpositions.coeffRef(i) = i;
- m_rank = j;
break;
}
@@ -200,11 +182,9 @@
m_matrix.coeffRef(j,j) = Djj;
// Finish early if the matrix is not full rank.
- if(ei_abs(Djj) < cutoff) // i made experiments, this is better than \
isMuchSmallerThan(biggest_in_corner), and of course
- // much better than plain sign comparison as used to be \
done before. + if(ei_abs(Djj) < cutoff)
{
for(int i = j; i < size; i++) m_transpositions.coeffRef(i) = i;
- m_rank = j;
break;
}
@@ -230,7 +210,7 @@
/** Computes the solution x of \f$ A x = b \f$ using the current decomposition of A.
* The result is stored in \a result
*
- * \returns true in case of success, false otherwise.
+ * \returns true always! If you need to check for existence of solutions, use \
another decomposition like LU, QR, or SVD.
*
* In other words, it computes \f$ b = A^{-1} b \f$ with
* \f$ P^T{L^{*}}^{-1} D^{-1} L^{-1} P b \f$ from right to left.
@@ -252,6 +232,8 @@
*
* \param bAndX represents both the right-hand side matrix b and result x.
*
+ * \returns true always! If you need to check for existence of solutions, use \
another decomposition like LU, QR, or SVD. + *
* This version avoids a copy when the right hand side matrix b is not
* needed anymore.
*
@@ -264,8 +246,6 @@
const int size = m_matrix.rows();
ei_assert(size == bAndX.rows());
- if (m_rank != size) return false;
-
// z = P b
for(int i = 0; i < size; ++i) \
bAndX.row(m_transpositions.coeff(i)).swap(bAndX.row(i));
--- trunk/kdesupport/eigen2/Eigen/src/Cholesky/LLT.h #947640:947641
@@ -41,6 +41,10 @@
* and even faster. Nevertheless, this standard Cholesky decomposition remains \
useful in many other
* situations like generalised eigen problems with hermitian matrices.
*
+ * Remember that Cholesky decompositions are not rank-revealing. This LLT \
decomposition is only stable on positive definite matrices, + * use LDLT instead for \
the semidefinite case. Also, do not use a Cholesky decomposition to determine whether \
a system of equations + * has a solution.
+ *
* \sa MatrixBase::llt(), class LDLT
*/
/* HEY THIS DOX IS DISABLED BECAUSE THERE's A BUG EITHER HERE OR IN LDLT ABOUT THAT \
(OR BOTH) @@ -70,12 +74,6 @@
/** \returns the lower triangular matrix L */
inline Part<MatrixType, LowerTriangular> matrixL(void) const { return m_matrix; \
}
- /** \returns true if the matrix is positive definite */
- inline bool isPositiveDefinite(void) const { return m_isPositiveDefinite; }
-
- /** \returns true if the matrix is invertible, which in this context is \
equivalent to positive definite */
- inline bool isInvertible(void) const { return m_isPositiveDefinite; }
-
template<typename RhsDerived, typename ResDerived>
bool solve(const MatrixBase<RhsDerived> &b, MatrixBase<ResDerived> *result) \
const;
@@ -90,7 +88,6 @@
* The strict upper part is not used and even not initialized.
*/
MatrixType m_matrix;
- bool m_isPositiveDefinite;
};
/** Computes / recomputes the Cholesky decomposition A = LL^* = U^*U of \a matrix
@@ -101,23 +98,24 @@
assert(a.rows()==a.cols());
const int size = a.rows();
m_matrix.resize(size, size);
- const RealScalar reference = size * a.diagonal().cwise().abs().maxCoeff();
+ // The biggest overall is the point of reference to which further diagonals
+ // are compared; if any diagonal is negligible compared
+ // to the largest overall, the algorithm bails. This cutoff is suggested
+ // in "Analysis of the Cholesky Decomposition of a Semi-definite Matrix" by
+ // Nicholas J. Higham. Also see "Accuracy and Stability of Numerical
+ // Algorithms" page 217, also by Higham.
+ const RealScalar cutoff = machine_epsilon<Scalar>() * size * \
a.diagonal().cwise().abs().maxCoeff(); RealScalar x;
x = ei_real(a.coeff(0,0));
- m_isPositiveDefinite = !ei_isMuchSmallerThan(x, reference) && \
ei_isMuchSmallerThan(ei_imag(a.coeff(0,0)), reference); m_matrix.coeffRef(0,0) = \
ei_sqrt(x); if(size==1)
return;
m_matrix.col(0).end(size-1) = a.row(0).end(size-1).adjoint() / \
ei_real(m_matrix.coeff(0,0)); for (int j = 1; j < size; ++j)
{
- Scalar tmp = ei_real(a.coeff(j,j)) - m_matrix.row(j).start(j).squaredNorm();
- x = ei_real(tmp);
- if (ei_isMuchSmallerThan(x, reference) || (!ei_isMuchSmallerThan(ei_imag(tmp), \
reference)))
- {
- m_isPositiveDefinite = false;
- return;
- }
+ x = ei_real(a.coeff(j,j)) - m_matrix.row(j).start(j).squaredNorm();
+ if (ei_abs(x) < cutoff) continue;
+
m_matrix.coeffRef(j,j) = x = ei_sqrt(x);
int endSize = size-j-1;
@@ -137,7 +135,7 @@
/** Computes the solution x of \f$ A x = b \f$ using the current decomposition of A.
* The result is stored in \a result
*
- * \returns true in case of success, false otherwise.
+ * \returns true always! If you need to check for existence of solutions, use \
another decomposition like LU, QR, or SVD.
*
* In other words, it computes \f$ b = A^{-1} b \f$ with
* \f$ {L^{*}}^{-1} L^{-1} b \f$ from right to left.
@@ -160,6 +158,8 @@
*
* \param bAndX represents both the right-hand side matrix b and result x.
*
+ * \returns true always! If you need to check for existence of solutions, use \
another decomposition like LU, QR, or SVD. + *
* This version avoids a copy when the right hand side matrix b is not
* needed anymore.
*
@@ -171,8 +171,6 @@
{
const int size = m_matrix.rows();
ei_assert(size==bAndX.rows());
- if (!m_isPositiveDefinite)
- return false;
matrixL().solveTriangularInPlace(bAndX);
m_matrix.adjoint().template part<UpperTriangular>().solveTriangularInPlace(bAndX);
return true;
--- trunk/kdesupport/eigen2/Eigen/src/Core/MathFunctions.h #947640:947641
@@ -26,6 +26,7 @@
#define EIGEN_MATHFUNCTIONS_H
template<typename T> inline typename NumTraits<T>::Real precision();
+template<typename T> inline typename NumTraits<T>::Real machine_epsilon();
template<typename T> inline T ei_random(T a, T b);
template<typename T> inline T ei_random();
template<typename T> inline T ei_random_amplitude()
@@ -50,6 +51,7 @@
**************/
template<> inline int precision<int>() { return 0; }
+template<> inline int machine_epsilon<int>() { return 0; }
inline int ei_real(int x) { return x; }
inline int ei_imag(int) { return 0; }
inline int ei_conj(int x) { return x; }
@@ -102,6 +104,7 @@
**************/
template<> inline float precision<float>() { return 1e-5f; }
+template<> inline float machine_epsilon<float>() { return 1.192e-07f; }
inline float ei_real(float x) { return x; }
inline float ei_imag(float) { return 0.f; }
inline float ei_conj(float x) { return x; }
@@ -147,6 +150,8 @@
**************/
template<> inline double precision<double>() { return 1e-11; }
+template<> inline double machine_epsilon<double>() { return 2.220e-16; }
+
inline double ei_real(double x) { return x; }
inline double ei_imag(double) { return 0.; }
inline double ei_conj(double x) { return x; }
@@ -192,6 +197,7 @@
*********************/
template<> inline float precision<std::complex<float> >() { return \
precision<float>(); } +template<> inline float machine_epsilon<std::complex<float> \
>() { return machine_epsilon<float>(); } inline float ei_real(const \
> std::complex<float>& x) { return std::real(x); }
inline float ei_imag(const std::complex<float>& x) { return std::imag(x); }
inline std::complex<float> ei_conj(const std::complex<float>& x) { return \
std::conj(x); } @@ -225,6 +231,7 @@
**********************/
template<> inline double precision<std::complex<double> >() { return \
precision<double>(); } +template<> inline double machine_epsilon<std::complex<double> \
>() { return machine_epsilon<double>(); } inline double ei_real(const \
> std::complex<double>& x) { return std::real(x); }
inline double ei_imag(const std::complex<double>& x) { return std::imag(x); }
inline std::complex<double> ei_conj(const std::complex<double>& x) { return \
std::conj(x); } @@ -259,6 +266,7 @@
******************/
template<> inline long double precision<long double>() { return precision<double>(); \
} +template<> inline long double machine_epsilon<long double>() { return 1.084e-19l; \
} inline long double ei_real(long double x) { return x; }
inline long double ei_imag(long double) { return 0.; }
inline long double ei_conj(long double x) { return x; }
--- trunk/kdesupport/eigen2/test/cholesky.cpp #947640:947641
@@ -86,7 +86,6 @@
{
LLT<SquareMatrixType> chol(symm);
- VERIFY(chol.isPositiveDefinite());
VERIFY_IS_APPROX(symm, chol.matrixL() * chol.matrixL().adjoint());
chol.solve(vecB, &vecX);
VERIFY_IS_APPROX(symm * vecX, vecB);
@@ -103,18 +102,6 @@
{
LDLT<SquareMatrixType> ldlt(symm);
- VERIFY(ldlt.isInvertible());
- if(sign == 1)
- {
- VERIFY(ldlt.isPositive());
- VERIFY(ldlt.isPositiveDefinite());
- }
- if(sign == -1)
- {
- VERIFY(ldlt.isNegative());
- VERIFY(ldlt.isNegativeDefinite());
- }
-
// TODO(keir): This doesn't make sense now that LDLT pivots.
//VERIFY_IS_APPROX(symm, ldlt.matrixL() * ldlt.vectorD().asDiagonal() * \
ldlt.matrixL().adjoint()); ldlt.solve(vecB, &vecX);
@@ -123,15 +110,6 @@
VERIFY_IS_APPROX(symm * matX, matB);
}
- // test isPositiveDefinite on non definite matrix
- if (rows>4)
- {
- SquareMatrixType symm = a0.block(0,0,rows,cols-4) * \
a0.block(0,0,rows,cols-4).adjoint();
- LLT<SquareMatrixType> chol(symm);
- VERIFY(!chol.isPositiveDefinite());
- LDLT<SquareMatrixType> cholnosqrt(symm);
- VERIFY(!cholnosqrt.isPositiveDefinite());
- }
}
template<typename Derived>
@@ -156,30 +134,7 @@
}
}
-template<typename MatrixType> void ldlt_rank()
-{
- // NOTE there seems to be a problem with too small sizes -- could easily lie in \
the doSomeRankPreservingOperations function
- int rows = ei_random<int>(50,200);
- int rank = ei_random<int>(40, rows-1);
-
- // generate a random positive matrix a of given rank
- MatrixType m = MatrixType::Random(rows,rows);
- QR<MatrixType> qr(m);
- typedef typename MatrixType::Scalar Scalar;
- typedef typename NumTraits<Scalar>::Real RealScalar;
- typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> DiagVectorType;
- DiagVectorType d(rows);
- d.setZero();
- for(int i = 0; i < rank; i++) d(i)=RealScalar(1);
- MatrixType a = qr.matrixQ() * d.asDiagonal() * qr.matrixQ().adjoint();
-
- LDLT<MatrixType> ldlt(a);
-
- VERIFY( ei_abs(ldlt.rank() - rank) <= rank / 20 ); // yes, LDLT::rank is a bit \
inaccurate...
-}
-
-
void test_cholesky()
{
for(int i = 0; i < g_repeat; i++) {
@@ -191,9 +146,4 @@
CALL_SUBTEST( cholesky(MatrixXd(17,17)) );
CALL_SUBTEST( cholesky(MatrixXf(200,200)) );
}
- for(int i = 0; i < g_repeat/3; i++) {
- CALL_SUBTEST( ldlt_rank<MatrixXd>() );
- CALL_SUBTEST( ldlt_rank<MatrixXf>() );
- CALL_SUBTEST( ldlt_rank<MatrixXcd>() );
- }
}
[prev in list] [next in list] [prev in thread] [next in thread]
Configure |
About |
News |
Add a list |
Sponsored by KoreLogic