ground/gcs/src/libs/eigen/test/cholesky.cpp

   1 // This file is part of Eigen, a lightweight C++ template library
   2 // for linear algebra.
   3 //
   4 // Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
   5 //
   6 // This Source Code Form is subject to the terms of the Mozilla
   7 // Public License v. 2.0. If a copy of the MPL was not distributed
   8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
   9
  10 #ifndef EIGEN_NO_ASSERTION_CHECKING
  11 #define EIGEN_NO_ASSERTION_CHECKING
  12 #endif
  13
  14 static int nb_temporaries;
  15
  16 #define EIGEN_DENSE_STORAGE_CTOR_PLUGIN { if(size!=0) nb_temporaries++; }
  17
  18 #include "main.h"
  19 #include <Eigen/Cholesky>
  20 #include <Eigen/QR>
  21
  22 #define VERIFY_EVALUATION_COUNT(XPR,N) {\
  23     nb_temporaries = 0; \
  24     XPR; \
  25     if(nb_temporaries!=N) std::cerr << "nb_temporaries == " << nb_temporaries << "\n"; \
  26     VERIFY( (#XPR) && nb_temporaries==N ); \
  27   }
  28
  29 template<typename MatrixType,template <typename,int> class CholType> void test_chol_update(const MatrixType& symm)
  30 {
  31   typedef typename MatrixType::Scalar Scalar;
  32   typedef typename MatrixType::RealScalar RealScalar;
  33   typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
  34
  35   MatrixType symmLo = symm.template triangularView<Lower>();
  36   MatrixType symmUp = symm.template triangularView<Upper>();
  37   MatrixType symmCpy = symm;
  38
  39   CholType<MatrixType,Lower> chollo(symmLo);
  40   CholType<MatrixType,Upper> cholup(symmUp);
  41
  42   for (int k=0; k<10; ++k)
  43   {
  44     VectorType vec = VectorType::Random(symm.rows());
  45     RealScalar sigma = internal::random<RealScalar>();
  46     symmCpy += sigma * vec * vec.adjoint();
  47
  48     // we are doing some downdates, so it might be the case that the matrix is not SPD anymore
  49     CholType<MatrixType,Lower> chol(symmCpy);
  50     if(chol.info()!=Success)
  51       break;
  52
  53     chollo.rankUpdate(vec, sigma);
  54     VERIFY_IS_APPROX(symmCpy, chollo.reconstructedMatrix());
  55
  56     cholup.rankUpdate(vec, sigma);
  57     VERIFY_IS_APPROX(symmCpy, cholup.reconstructedMatrix());
  58   }
  59 }
  60
  61 template<typename MatrixType> void cholesky(const MatrixType& m)
  62 {
  63   typedef typename MatrixType::Index Index;
  64   /* this test covers the following files:
  65      LLT.h LDLT.h
  66   */
  67   Index rows = m.rows();
  68   Index cols = m.cols();
  69
  70   typedef typename MatrixType::Scalar Scalar;
  71   typedef typename NumTraits<Scalar>::Real RealScalar;
  72   typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> SquareMatrixType;
  73   typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
  74
  75   MatrixType a0 = MatrixType::Random(rows,cols);
  76   VectorType vecB = VectorType::Random(rows), vecX(rows);
  77   MatrixType matB = MatrixType::Random(rows,cols), matX(rows,cols);
  78   SquareMatrixType symm =  a0 * a0.adjoint();
  79   // let's make sure the matrix is not singular or near singular
  80   for (int k=0; k<3; ++k)
  81   {
  82     MatrixType a1 = MatrixType::Random(rows,cols);
  83     symm += a1 * a1.adjoint();
  84   }
  85
  86   // to test if really Cholesky only uses the upper triangular part, uncomment the following
  87   // FIXME: currently that fails !!
  88   //symm.template part<StrictlyLower>().setZero();
  89
  90   {
  91     SquareMatrixType symmUp = symm.template triangularView<Upper>();
  92     SquareMatrixType symmLo = symm.template triangularView<Lower>();
  93
  94     LLT<SquareMatrixType,Lower> chollo(symmLo);
  95     VERIFY_IS_APPROX(symm, chollo.reconstructedMatrix());
  96     vecX = chollo.solve(vecB);
  97     VERIFY_IS_APPROX(symm * vecX, vecB);
  98     matX = chollo.solve(matB);
  99     VERIFY_IS_APPROX(symm * matX, matB);
 100
 101     // test the upper mode
 102     LLT<SquareMatrixType,Upper> cholup(symmUp);
 103     VERIFY_IS_APPROX(symm, cholup.reconstructedMatrix());
 104     vecX = cholup.solve(vecB);
 105     VERIFY_IS_APPROX(symm * vecX, vecB);
 106     matX = cholup.solve(matB);
 107     VERIFY_IS_APPROX(symm * matX, matB);
 108
 109     MatrixType neg = -symmLo;
 110     chollo.compute(neg);
 111     VERIFY(chollo.info()==NumericalIssue);
 112
 113     VERIFY_IS_APPROX(MatrixType(chollo.matrixL().transpose().conjugate()), MatrixType(chollo.matrixU()));
 114     VERIFY_IS_APPROX(MatrixType(chollo.matrixU().transpose().conjugate()), MatrixType(chollo.matrixL()));
 115     VERIFY_IS_APPROX(MatrixType(cholup.matrixL().transpose().conjugate()), MatrixType(cholup.matrixU()));
 116     VERIFY_IS_APPROX(MatrixType(cholup.matrixU().transpose().conjugate()), MatrixType(cholup.matrixL()));
 117
 118     // test some special use cases of SelfCwiseBinaryOp:
 119     MatrixType m1 = MatrixType::Random(rows,cols), m2(rows,cols);
 120     m2 = m1;
 121     m2 += symmLo.template selfadjointView<Lower>().llt().solve(matB);
 122     VERIFY_IS_APPROX(m2, m1 + symmLo.template selfadjointView<Lower>().llt().solve(matB));
 123     m2 = m1;
 124     m2 -= symmLo.template selfadjointView<Lower>().llt().solve(matB);
 125     VERIFY_IS_APPROX(m2, m1 - symmLo.template selfadjointView<Lower>().llt().solve(matB));
 126     m2 = m1;
 127     m2.noalias() += symmLo.template selfadjointView<Lower>().llt().solve(matB);
 128     VERIFY_IS_APPROX(m2, m1 + symmLo.template selfadjointView<Lower>().llt().solve(matB));
 129     m2 = m1;
 130     m2.noalias() -= symmLo.template selfadjointView<Lower>().llt().solve(matB);
 131     VERIFY_IS_APPROX(m2, m1 - symmLo.template selfadjointView<Lower>().llt().solve(matB));
 132   }
 133
 134   // LDLT
 135   {
 136     int sign = internal::random<int>()%2 ? 1 : -1;
 137
 138     if(sign == -1)
 139     {
 140       symm = -symm; // test a negative matrix
 141     }
 142
 143     SquareMatrixType symmUp = symm.template triangularView<Upper>();
 144     SquareMatrixType symmLo = symm.template triangularView<Lower>();
 145
 146     LDLT<SquareMatrixType,Lower> ldltlo(symmLo);
 147     VERIFY_IS_APPROX(symm, ldltlo.reconstructedMatrix());
 148     vecX = ldltlo.solve(vecB);
 149     VERIFY_IS_APPROX(symm * vecX, vecB);
 150     matX = ldltlo.solve(matB);
 151     VERIFY_IS_APPROX(symm * matX, matB);
 152
 153     LDLT<SquareMatrixType,Upper> ldltup(symmUp);
 154     VERIFY_IS_APPROX(symm, ldltup.reconstructedMatrix());
 155     vecX = ldltup.solve(vecB);
 156     VERIFY_IS_APPROX(symm * vecX, vecB);
 157     matX = ldltup.solve(matB);
 158     VERIFY_IS_APPROX(symm * matX, matB);
 159
 160     VERIFY_IS_APPROX(MatrixType(ldltlo.matrixL().transpose().conjugate()), MatrixType(ldltlo.matrixU()));
 161     VERIFY_IS_APPROX(MatrixType(ldltlo.matrixU().transpose().conjugate()), MatrixType(ldltlo.matrixL()));
 162     VERIFY_IS_APPROX(MatrixType(ldltup.matrixL().transpose().conjugate()), MatrixType(ldltup.matrixU()));
 163     VERIFY_IS_APPROX(MatrixType(ldltup.matrixU().transpose().conjugate()), MatrixType(ldltup.matrixL()));
 164
 165     if(MatrixType::RowsAtCompileTime==Dynamic)
 166     {
 167       // note : each inplace permutation requires a small temporary vector (mask)
 168
 169       // check inplace solve
 170       matX = matB;
 171       VERIFY_EVALUATION_COUNT(matX = ldltlo.solve(matX), 0);
 172       VERIFY_IS_APPROX(matX, ldltlo.solve(matB).eval());
 173
 174
 175       matX = matB;
 176       VERIFY_EVALUATION_COUNT(matX = ldltup.solve(matX), 0);
 177       VERIFY_IS_APPROX(matX, ldltup.solve(matB).eval());
 178     }
 179
 180     // restore
 181     if(sign == -1)
 182       symm = -symm;
 183
 184     // check matrices coming from linear constraints with Lagrange multipliers
 185     if(rows>=3)
 186     {
 187       SquareMatrixType A = symm;
 188       int c = internal::random<int>(0,rows-2);
 189       A.bottomRightCorner(c,c).setZero();
 190       // Make sure a solution exists:
 191       vecX.setRandom();
 192       vecB = A * vecX;
 193       vecX.setZero();
 194       ldltlo.compute(A);
 195       VERIFY_IS_APPROX(A, ldltlo.reconstructedMatrix());
 196       vecX = ldltlo.solve(vecB);
 197       VERIFY_IS_APPROX(A * vecX, vecB);
 198     }
 199
 200     // check non-full rank matrices
 201     if(rows>=3)
 202     {
 203       int r = internal::random<int>(1,rows-1);
 204       Matrix<Scalar,Dynamic,Dynamic> a = Matrix<Scalar,Dynamic,Dynamic>::Random(rows,r);
 205       SquareMatrixType A = a * a.adjoint();
 206       // Make sure a solution exists:
 207       vecX.setRandom();
 208       vecB = A * vecX;
 209       vecX.setZero();
 210       ldltlo.compute(A);
 211       VERIFY_IS_APPROX(A, ldltlo.reconstructedMatrix());
 212       vecX = ldltlo.solve(vecB);
 213       VERIFY_IS_APPROX(A * vecX, vecB);
 214     }
 215
 216     // check matrices with a wide spectrum
 217     if(rows>=3)
 218     {
 219       RealScalar s = (std::min)(16,std::numeric_limits<RealScalar>::max_exponent10/8);
 220       Matrix<Scalar,Dynamic,Dynamic> a = Matrix<Scalar,Dynamic,Dynamic>::Random(rows,rows);
 221       Matrix<RealScalar,Dynamic,1> d =  Matrix<RealScalar,Dynamic,1>::Random(rows);
 222       for(int k=0; k<rows; ++k)
 223         d(k) = d(k)*std::pow(RealScalar(10),internal::random<RealScalar>(-s,s));
 224       SquareMatrixType A = a * d.asDiagonal() * a.adjoint();
 225       // Make sure a solution exists:
 226       vecX.setRandom();
 227       vecB = A * vecX;
 228       vecX.setZero();
 229       ldltlo.compute(A);
 230       VERIFY_IS_APPROX(A, ldltlo.reconstructedMatrix());
 231       vecX = ldltlo.solve(vecB);
 232       VERIFY_IS_APPROX(A * vecX, vecB);
 233     }
 234   }
 235
 236   // update/downdate
 237   CALL_SUBTEST(( test_chol_update<SquareMatrixType,LLT>(symm)  ));
 238   CALL_SUBTEST(( test_chol_update<SquareMatrixType,LDLT>(symm) ));
 239 }
 240
 241 template<typename MatrixType> void cholesky_cplx(const MatrixType& m)
 242 {
 243   // classic test
 244   cholesky(m);
 245
 246   // test mixing real/scalar types
 247
 248   typedef typename MatrixType::Index Index;
 249
 250   Index rows = m.rows();
 251   Index cols = m.cols();
 252
 253   typedef typename MatrixType::Scalar Scalar;
 254   typedef typename NumTraits<Scalar>::Real RealScalar;
 255   typedef Matrix<RealScalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> RealMatrixType;
 256   typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
 257
 258   RealMatrixType a0 = RealMatrixType::Random(rows,cols);
 259   VectorType vecB = VectorType::Random(rows), vecX(rows);
 260   MatrixType matB = MatrixType::Random(rows,cols), matX(rows,cols);
 261   RealMatrixType symm =  a0 * a0.adjoint();
 262   // let's make sure the matrix is not singular or near singular
 263   for (int k=0; k<3; ++k)
 264   {
 265     RealMatrixType a1 = RealMatrixType::Random(rows,cols);
 266     symm += a1 * a1.adjoint();
 267   }
 268
 269   {
 270     RealMatrixType symmLo = symm.template triangularView<Lower>();
 271
 272     LLT<RealMatrixType,Lower> chollo(symmLo);
 273     VERIFY_IS_APPROX(symm, chollo.reconstructedMatrix());
 274     vecX = chollo.solve(vecB);
 275     VERIFY_IS_APPROX(symm * vecX, vecB);
 276 //     matX = chollo.solve(matB);
 277 //     VERIFY_IS_APPROX(symm * matX, matB);
 278   }
 279
 280   // LDLT
 281   {
 282     int sign = internal::random<int>()%2 ? 1 : -1;
 283
 284     if(sign == -1)
 285     {
 286       symm = -symm; // test a negative matrix
 287     }
 288
 289     RealMatrixType symmLo = symm.template triangularView<Lower>();
 290
 291     LDLT<RealMatrixType,Lower> ldltlo(symmLo);
 292     VERIFY_IS_APPROX(symm, ldltlo.reconstructedMatrix());
 293     vecX = ldltlo.solve(vecB);
 294     VERIFY_IS_APPROX(symm * vecX, vecB);
 295 //     matX = ldltlo.solve(matB);
 296 //     VERIFY_IS_APPROX(symm * matX, matB);
 297   }
 298 }
 299
 300 // regression test for bug 241
 301 template<typename MatrixType> void cholesky_bug241(const MatrixType& m)
 302 {
 303   eigen_assert(m.rows() == 2 && m.cols() == 2);
 304
 305   typedef typename MatrixType::Scalar Scalar;
 306   typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
 307
 308   MatrixType matA;
 309   matA << 1, 1, 1, 1;
 310   VectorType vecB;
 311   vecB << 1, 1;
 312   VectorType vecX = matA.ldlt().solve(vecB);
 313   VERIFY_IS_APPROX(matA * vecX, vecB);
 314 }
 315
 316 // LDLT is not guaranteed to work for indefinite matrices, but happens to work fine if matrix is diagonal.
 317 // This test checks that LDLT reports correctly that matrix is indefinite.
 318 // See http://forum.kde.org/viewtopic.php?f=74&t=106942 and bug 736
 319 template<typename MatrixType> void cholesky_definiteness(const MatrixType& m)
 320 {
 321   eigen_assert(m.rows() == 2 && m.cols() == 2);
 322   MatrixType mat;
 323   LDLT<MatrixType> ldlt(2);
 324
 325   {
 326     mat << 1, 0, 0, -1;
 327     ldlt.compute(mat);
 328     VERIFY(!ldlt.isNegative());
 329     VERIFY(!ldlt.isPositive());
 330   }
 331   {
 332     mat << 1, 2, 2, 1;
 333     ldlt.compute(mat);
 334     VERIFY(!ldlt.isNegative());
 335     VERIFY(!ldlt.isPositive());
 336   }
 337   {
 338     mat << 0, 0, 0, 0;
 339     ldlt.compute(mat);
 340     VERIFY(ldlt.isNegative());
 341     VERIFY(ldlt.isPositive());
 342   }
 343   {
 344     mat << 0, 0, 0, 1;
 345     ldlt.compute(mat);
 346     VERIFY(!ldlt.isNegative());
 347     VERIFY(ldlt.isPositive());
 348   }
 349   {
 350     mat << -1, 0, 0, 0;
 351     ldlt.compute(mat);
 352     VERIFY(ldlt.isNegative());
 353     VERIFY(!ldlt.isPositive());
 354   }
 355 }
 356
 357 template<typename MatrixType> void cholesky_verify_assert()
 358 {
 359   MatrixType tmp;
 360
 361   LLT<MatrixType> llt;
 362   VERIFY_RAISES_ASSERT(llt.matrixL())
 363   VERIFY_RAISES_ASSERT(llt.matrixU())
 364   VERIFY_RAISES_ASSERT(llt.solve(tmp))
 365   VERIFY_RAISES_ASSERT(llt.solveInPlace(&tmp))
 366
 367   LDLT<MatrixType> ldlt;
 368   VERIFY_RAISES_ASSERT(ldlt.matrixL())
 369   VERIFY_RAISES_ASSERT(ldlt.permutationP())
 370   VERIFY_RAISES_ASSERT(ldlt.vectorD())
 371   VERIFY_RAISES_ASSERT(ldlt.isPositive())
 372   VERIFY_RAISES_ASSERT(ldlt.isNegative())
 373   VERIFY_RAISES_ASSERT(ldlt.solve(tmp))
 374   VERIFY_RAISES_ASSERT(ldlt.solveInPlace(&tmp))
 375 }
 376
 377 void test_cholesky()
 378 {
 379   int s = 0;
 380   for(int i = 0; i < g_repeat; i++) {
 381     CALL_SUBTEST_1( cholesky(Matrix<double,1,1>()) );
 382     CALL_SUBTEST_3( cholesky(Matrix2d()) );
 383     CALL_SUBTEST_3( cholesky_bug241(Matrix2d()) );
 384     CALL_SUBTEST_3( cholesky_definiteness(Matrix2d()) );
 385     CALL_SUBTEST_4( cholesky(Matrix3f()) );
 386     CALL_SUBTEST_5( cholesky(Matrix4d()) );
 387     s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE);
 388     CALL_SUBTEST_2( cholesky(MatrixXd(s,s)) );
 389     s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE/2);
 390     CALL_SUBTEST_6( cholesky_cplx(MatrixXcd(s,s)) );
 391   }
 392
 393   CALL_SUBTEST_4( cholesky_verify_assert<Matrix3f>() );
 394   CALL_SUBTEST_7( cholesky_verify_assert<Matrix3d>() );
 395   CALL_SUBTEST_8( cholesky_verify_assert<MatrixXf>() );
 396   CALL_SUBTEST_2( cholesky_verify_assert<MatrixXd>() );
 397
 398   // Test problem size constructors
 399   CALL_SUBTEST_9( LLT<MatrixXf>(10) );
 400   CALL_SUBTEST_9( LDLT<MatrixXf>(10) );
 401
 402   TEST_SET_BUT_UNUSED_VARIABLE(s)
 403   TEST_SET_BUT_UNUSED_VARIABLE(nb_temporaries)
 404 }