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 (at) 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 #define TEST_ENABLE_TEMPORARY_TRACKING 15 16 #include "main.h" 17 #include <Eigen/Cholesky> 18 #include <Eigen/QR> 19 20 template<typename MatrixType, int UpLo> 21 typename MatrixType::RealScalar matrix_l1_norm(const MatrixType& m) { 22 MatrixType symm = m.template selfadjointView<UpLo>(); 23 return symm.cwiseAbs().colwise().sum().maxCoeff(); 24 } 25 26 template<typename MatrixType,template <typename,int> class CholType> void test_chol_update(const MatrixType& symm) 27 { 28 typedef typename MatrixType::Scalar Scalar; 29 typedef typename MatrixType::RealScalar RealScalar; 30 typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType; 31 32 MatrixType symmLo = symm.template triangularView<Lower>(); 33 MatrixType symmUp = symm.template triangularView<Upper>(); 34 MatrixType symmCpy = symm; 35 36 CholType<MatrixType,Lower> chollo(symmLo); 37 CholType<MatrixType,Upper> cholup(symmUp); 38 39 for (int k=0; k<10; ++k) 40 { 41 VectorType vec = VectorType::Random(symm.rows()); 42 RealScalar sigma = internal::random<RealScalar>(); 43 symmCpy += sigma * vec * vec.adjoint(); 44 45 // we are doing some downdates, so it might be the case that the matrix is not SPD anymore 46 CholType<MatrixType,Lower> chol(symmCpy); 47 if(chol.info()!=Success) 48 break; 49 50 chollo.rankUpdate(vec, sigma); 51 VERIFY_IS_APPROX(symmCpy, chollo.reconstructedMatrix()); 52 53 cholup.rankUpdate(vec, sigma); 54 VERIFY_IS_APPROX(symmCpy, cholup.reconstructedMatrix()); 55 } 56 } 57 58 template<typename MatrixType> void cholesky(const MatrixType& m) 59 { 60 typedef typename MatrixType::Index Index; 61 /* this test covers the following files: 62 LLT.h LDLT.h 63 */ 64 Index rows = m.rows(); 65 Index cols = m.cols(); 66 67 typedef typename MatrixType::Scalar Scalar; 68 typedef typename NumTraits<Scalar>::Real RealScalar; 69 typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> SquareMatrixType; 70 typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType; 71 72 MatrixType a0 = MatrixType::Random(rows,cols); 73 VectorType vecB = VectorType::Random(rows), vecX(rows); 74 MatrixType matB = MatrixType::Random(rows,cols), matX(rows,cols); 75 SquareMatrixType symm = a0 * a0.adjoint(); 76 // let's make sure the matrix is not singular or near singular 77 for (int k=0; k<3; ++k) 78 { 79 MatrixType a1 = MatrixType::Random(rows,cols); 80 symm += a1 * a1.adjoint(); 81 } 82 83 { 84 SquareMatrixType symmUp = symm.template triangularView<Upper>(); 85 SquareMatrixType symmLo = symm.template triangularView<Lower>(); 86 87 LLT<SquareMatrixType,Lower> chollo(symmLo); 88 VERIFY_IS_APPROX(symm, chollo.reconstructedMatrix()); 89 vecX = chollo.solve(vecB); 90 VERIFY_IS_APPROX(symm * vecX, vecB); 91 matX = chollo.solve(matB); 92 VERIFY_IS_APPROX(symm * matX, matB); 93 94 const MatrixType symmLo_inverse = chollo.solve(MatrixType::Identity(rows,cols)); 95 RealScalar rcond = (RealScalar(1) / matrix_l1_norm<MatrixType, Lower>(symmLo)) / 96 matrix_l1_norm<MatrixType, Lower>(symmLo_inverse); 97 RealScalar rcond_est = chollo.rcond(); 98 // Verify that the estimated condition number is within a factor of 10 of the 99 // truth. 100 VERIFY(rcond_est > rcond / 10 && rcond_est < rcond * 10); 101 102 // test the upper mode 103 LLT<SquareMatrixType,Upper> cholup(symmUp); 104 VERIFY_IS_APPROX(symm, cholup.reconstructedMatrix()); 105 vecX = cholup.solve(vecB); 106 VERIFY_IS_APPROX(symm * vecX, vecB); 107 matX = cholup.solve(matB); 108 VERIFY_IS_APPROX(symm * matX, matB); 109 110 // Verify that the estimated condition number is within a factor of 10 of the 111 // truth. 112 const MatrixType symmUp_inverse = cholup.solve(MatrixType::Identity(rows,cols)); 113 rcond = (RealScalar(1) / matrix_l1_norm<MatrixType, Upper>(symmUp)) / 114 matrix_l1_norm<MatrixType, Upper>(symmUp_inverse); 115 rcond_est = cholup.rcond(); 116 VERIFY(rcond_est > rcond / 10 && rcond_est < rcond * 10); 117 118 119 MatrixType neg = -symmLo; 120 chollo.compute(neg); 121 VERIFY(chollo.info()==NumericalIssue); 122 123 VERIFY_IS_APPROX(MatrixType(chollo.matrixL().transpose().conjugate()), MatrixType(chollo.matrixU())); 124 VERIFY_IS_APPROX(MatrixType(chollo.matrixU().transpose().conjugate()), MatrixType(chollo.matrixL())); 125 VERIFY_IS_APPROX(MatrixType(cholup.matrixL().transpose().conjugate()), MatrixType(cholup.matrixU())); 126 VERIFY_IS_APPROX(MatrixType(cholup.matrixU().transpose().conjugate()), MatrixType(cholup.matrixL())); 127 128 // test some special use cases of SelfCwiseBinaryOp: 129 MatrixType m1 = MatrixType::Random(rows,cols), m2(rows,cols); 130 m2 = m1; 131 m2 += symmLo.template selfadjointView<Lower>().llt().solve(matB); 132 VERIFY_IS_APPROX(m2, m1 + symmLo.template selfadjointView<Lower>().llt().solve(matB)); 133 m2 = m1; 134 m2 -= symmLo.template selfadjointView<Lower>().llt().solve(matB); 135 VERIFY_IS_APPROX(m2, m1 - symmLo.template selfadjointView<Lower>().llt().solve(matB)); 136 m2 = m1; 137 m2.noalias() += symmLo.template selfadjointView<Lower>().llt().solve(matB); 138 VERIFY_IS_APPROX(m2, m1 + symmLo.template selfadjointView<Lower>().llt().solve(matB)); 139 m2 = m1; 140 m2.noalias() -= symmLo.template selfadjointView<Lower>().llt().solve(matB); 141 VERIFY_IS_APPROX(m2, m1 - symmLo.template selfadjointView<Lower>().llt().solve(matB)); 142 } 143 144 // LDLT 145 { 146 int sign = internal::random<int>()%2 ? 1 : -1; 147 148 if(sign == -1) 149 { 150 symm = -symm; // test a negative matrix 151 } 152 153 SquareMatrixType symmUp = symm.template triangularView<Upper>(); 154 SquareMatrixType symmLo = symm.template triangularView<Lower>(); 155 156 LDLT<SquareMatrixType,Lower> ldltlo(symmLo); 157 VERIFY(ldltlo.info()==Success); 158 VERIFY_IS_APPROX(symm, ldltlo.reconstructedMatrix()); 159 vecX = ldltlo.solve(vecB); 160 VERIFY_IS_APPROX(symm * vecX, vecB); 161 matX = ldltlo.solve(matB); 162 VERIFY_IS_APPROX(symm * matX, matB); 163 164 const MatrixType symmLo_inverse = ldltlo.solve(MatrixType::Identity(rows,cols)); 165 RealScalar rcond = (RealScalar(1) / matrix_l1_norm<MatrixType, Lower>(symmLo)) / 166 matrix_l1_norm<MatrixType, Lower>(symmLo_inverse); 167 RealScalar rcond_est = ldltlo.rcond(); 168 // Verify that the estimated condition number is within a factor of 10 of the 169 // truth. 170 VERIFY(rcond_est > rcond / 10 && rcond_est < rcond * 10); 171 172 173 LDLT<SquareMatrixType,Upper> ldltup(symmUp); 174 VERIFY(ldltup.info()==Success); 175 VERIFY_IS_APPROX(symm, ldltup.reconstructedMatrix()); 176 vecX = ldltup.solve(vecB); 177 VERIFY_IS_APPROX(symm * vecX, vecB); 178 matX = ldltup.solve(matB); 179 VERIFY_IS_APPROX(symm * matX, matB); 180 181 // Verify that the estimated condition number is within a factor of 10 of the 182 // truth. 183 const MatrixType symmUp_inverse = ldltup.solve(MatrixType::Identity(rows,cols)); 184 rcond = (RealScalar(1) / matrix_l1_norm<MatrixType, Upper>(symmUp)) / 185 matrix_l1_norm<MatrixType, Upper>(symmUp_inverse); 186 rcond_est = ldltup.rcond(); 187 VERIFY(rcond_est > rcond / 10 && rcond_est < rcond * 10); 188 189 VERIFY_IS_APPROX(MatrixType(ldltlo.matrixL().transpose().conjugate()), MatrixType(ldltlo.matrixU())); 190 VERIFY_IS_APPROX(MatrixType(ldltlo.matrixU().transpose().conjugate()), MatrixType(ldltlo.matrixL())); 191 VERIFY_IS_APPROX(MatrixType(ldltup.matrixL().transpose().conjugate()), MatrixType(ldltup.matrixU())); 192 VERIFY_IS_APPROX(MatrixType(ldltup.matrixU().transpose().conjugate()), MatrixType(ldltup.matrixL())); 193 194 if(MatrixType::RowsAtCompileTime==Dynamic) 195 { 196 // note : each inplace permutation requires a small temporary vector (mask) 197 198 // check inplace solve 199 matX = matB; 200 VERIFY_EVALUATION_COUNT(matX = ldltlo.solve(matX), 0); 201 VERIFY_IS_APPROX(matX, ldltlo.solve(matB).eval()); 202 203 204 matX = matB; 205 VERIFY_EVALUATION_COUNT(matX = ldltup.solve(matX), 0); 206 VERIFY_IS_APPROX(matX, ldltup.solve(matB).eval()); 207 } 208 209 // restore 210 if(sign == -1) 211 symm = -symm; 212 213 // check matrices coming from linear constraints with Lagrange multipliers 214 if(rows>=3) 215 { 216 SquareMatrixType A = symm; 217 Index c = internal::random<Index>(0,rows-2); 218 A.bottomRightCorner(c,c).setZero(); 219 // Make sure a solution exists: 220 vecX.setRandom(); 221 vecB = A * vecX; 222 vecX.setZero(); 223 ldltlo.compute(A); 224 VERIFY_IS_APPROX(A, ldltlo.reconstructedMatrix()); 225 vecX = ldltlo.solve(vecB); 226 VERIFY_IS_APPROX(A * vecX, vecB); 227 } 228 229 // check non-full rank matrices 230 if(rows>=3) 231 { 232 Index r = internal::random<Index>(1,rows-1); 233 Matrix<Scalar,Dynamic,Dynamic> a = Matrix<Scalar,Dynamic,Dynamic>::Random(rows,r); 234 SquareMatrixType A = a * a.adjoint(); 235 // Make sure a solution exists: 236 vecX.setRandom(); 237 vecB = A * vecX; 238 vecX.setZero(); 239 ldltlo.compute(A); 240 VERIFY_IS_APPROX(A, ldltlo.reconstructedMatrix()); 241 vecX = ldltlo.solve(vecB); 242 VERIFY_IS_APPROX(A * vecX, vecB); 243 } 244 245 // check matrices with a wide spectrum 246 if(rows>=3) 247 { 248 using std::pow; 249 using std::sqrt; 250 RealScalar s = (std::min)(16,std::numeric_limits<RealScalar>::max_exponent10/8); 251 Matrix<Scalar,Dynamic,Dynamic> a = Matrix<Scalar,Dynamic,Dynamic>::Random(rows,rows); 252 Matrix<RealScalar,Dynamic,1> d = Matrix<RealScalar,Dynamic,1>::Random(rows); 253 for(Index k=0; k<rows; ++k) 254 d(k) = d(k)*pow(RealScalar(10),internal::random<RealScalar>(-s,s)); 255 SquareMatrixType A = a * d.asDiagonal() * a.adjoint(); 256 // Make sure a solution exists: 257 vecX.setRandom(); 258 vecB = A * vecX; 259 vecX.setZero(); 260 ldltlo.compute(A); 261 VERIFY_IS_APPROX(A, ldltlo.reconstructedMatrix()); 262 vecX = ldltlo.solve(vecB); 263 264 if(ldltlo.vectorD().real().cwiseAbs().minCoeff()>RealScalar(0)) 265 { 266 VERIFY_IS_APPROX(A * vecX,vecB); 267 } 268 else 269 { 270 RealScalar large_tol = sqrt(test_precision<RealScalar>()); 271 VERIFY((A * vecX).isApprox(vecB, large_tol)); 272 273 ++g_test_level; 274 VERIFY_IS_APPROX(A * vecX,vecB); 275 --g_test_level; 276 } 277 } 278 } 279 280 // update/downdate 281 CALL_SUBTEST(( test_chol_update<SquareMatrixType,LLT>(symm) )); 282 CALL_SUBTEST(( test_chol_update<SquareMatrixType,LDLT>(symm) )); 283 } 284 285 template<typename MatrixType> void cholesky_cplx(const MatrixType& m) 286 { 287 // classic test 288 cholesky(m); 289 290 // test mixing real/scalar types 291 292 typedef typename MatrixType::Index Index; 293 294 Index rows = m.rows(); 295 Index cols = m.cols(); 296 297 typedef typename MatrixType::Scalar Scalar; 298 typedef typename NumTraits<Scalar>::Real RealScalar; 299 typedef Matrix<RealScalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> RealMatrixType; 300 typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType; 301 302 RealMatrixType a0 = RealMatrixType::Random(rows,cols); 303 VectorType vecB = VectorType::Random(rows), vecX(rows); 304 MatrixType matB = MatrixType::Random(rows,cols), matX(rows,cols); 305 RealMatrixType symm = a0 * a0.adjoint(); 306 // let's make sure the matrix is not singular or near singular 307 for (int k=0; k<3; ++k) 308 { 309 RealMatrixType a1 = RealMatrixType::Random(rows,cols); 310 symm += a1 * a1.adjoint(); 311 } 312 313 { 314 RealMatrixType symmLo = symm.template triangularView<Lower>(); 315 316 LLT<RealMatrixType,Lower> chollo(symmLo); 317 VERIFY_IS_APPROX(symm, chollo.reconstructedMatrix()); 318 vecX = chollo.solve(vecB); 319 VERIFY_IS_APPROX(symm * vecX, vecB); 320 // matX = chollo.solve(matB); 321 // VERIFY_IS_APPROX(symm * matX, matB); 322 } 323 324 // LDLT 325 { 326 int sign = internal::random<int>()%2 ? 1 : -1; 327 328 if(sign == -1) 329 { 330 symm = -symm; // test a negative matrix 331 } 332 333 RealMatrixType symmLo = symm.template triangularView<Lower>(); 334 335 LDLT<RealMatrixType,Lower> ldltlo(symmLo); 336 VERIFY(ldltlo.info()==Success); 337 VERIFY_IS_APPROX(symm, ldltlo.reconstructedMatrix()); 338 vecX = ldltlo.solve(vecB); 339 VERIFY_IS_APPROX(symm * vecX, vecB); 340 // matX = ldltlo.solve(matB); 341 // VERIFY_IS_APPROX(symm * matX, matB); 342 } 343 } 344 345 // regression test for bug 241 346 template<typename MatrixType> void cholesky_bug241(const MatrixType& m) 347 { 348 eigen_assert(m.rows() == 2 && m.cols() == 2); 349 350 typedef typename MatrixType::Scalar Scalar; 351 typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType; 352 353 MatrixType matA; 354 matA << 1, 1, 1, 1; 355 VectorType vecB; 356 vecB << 1, 1; 357 VectorType vecX = matA.ldlt().solve(vecB); 358 VERIFY_IS_APPROX(matA * vecX, vecB); 359 } 360 361 // LDLT is not guaranteed to work for indefinite matrices, but happens to work fine if matrix is diagonal. 362 // This test checks that LDLT reports correctly that matrix is indefinite. 363 // See http://forum.kde.org/viewtopic.php?f=74&t=106942 and bug 736 364 template<typename MatrixType> void cholesky_definiteness(const MatrixType& m) 365 { 366 eigen_assert(m.rows() == 2 && m.cols() == 2); 367 MatrixType mat; 368 LDLT<MatrixType> ldlt(2); 369 370 { 371 mat << 1, 0, 0, -1; 372 ldlt.compute(mat); 373 VERIFY(ldlt.info()==Success); 374 VERIFY(!ldlt.isNegative()); 375 VERIFY(!ldlt.isPositive()); 376 } 377 { 378 mat << 1, 2, 2, 1; 379 ldlt.compute(mat); 380 VERIFY(ldlt.info()==Success); 381 VERIFY(!ldlt.isNegative()); 382 VERIFY(!ldlt.isPositive()); 383 } 384 { 385 mat << 0, 0, 0, 0; 386 ldlt.compute(mat); 387 VERIFY(ldlt.info()==Success); 388 VERIFY(ldlt.isNegative()); 389 VERIFY(ldlt.isPositive()); 390 } 391 { 392 mat << 0, 0, 0, 1; 393 ldlt.compute(mat); 394 VERIFY(ldlt.info()==Success); 395 VERIFY(!ldlt.isNegative()); 396 VERIFY(ldlt.isPositive()); 397 } 398 { 399 mat << -1, 0, 0, 0; 400 ldlt.compute(mat); 401 VERIFY(ldlt.info()==Success); 402 VERIFY(ldlt.isNegative()); 403 VERIFY(!ldlt.isPositive()); 404 } 405 } 406 407 template<typename> 408 void cholesky_faillure_cases() 409 { 410 MatrixXd mat; 411 LDLT<MatrixXd> ldlt; 412 413 { 414 mat.resize(2,2); 415 mat << 0, 1, 1, 0; 416 ldlt.compute(mat); 417 VERIFY_IS_NOT_APPROX(mat,ldlt.reconstructedMatrix()); 418 VERIFY(ldlt.info()==NumericalIssue); 419 } 420 #if (!EIGEN_ARCH_i386) || defined(EIGEN_VECTORIZE_SSE2) 421 { 422 mat.resize(3,3); 423 mat << -1, -3, 3, 424 -3, -8.9999999999999999999, 1, 425 3, 1, 0; 426 ldlt.compute(mat); 427 VERIFY(ldlt.info()==NumericalIssue); 428 VERIFY_IS_NOT_APPROX(mat,ldlt.reconstructedMatrix()); 429 } 430 #endif 431 { 432 mat.resize(3,3); 433 mat << 1, 2, 3, 434 2, 4, 1, 435 3, 1, 0; 436 ldlt.compute(mat); 437 VERIFY(ldlt.info()==NumericalIssue); 438 VERIFY_IS_NOT_APPROX(mat,ldlt.reconstructedMatrix()); 439 } 440 441 { 442 mat.resize(8,8); 443 mat << 0.1, 0, -0.1, 0, 0, 0, 1, 0, 444 0, 4.24667, 0, 2.00333, 0, 0, 0, 0, 445 -0.1, 0, 0.2, 0, -0.1, 0, 0, 0, 446 0, 2.00333, 0, 8.49333, 0, 2.00333, 0, 0, 447 0, 0, -0.1, 0, 0.1, 0, 0, 1, 448 0, 0, 0, 2.00333, 0, 4.24667, 0, 0, 449 1, 0, 0, 0, 0, 0, 0, 0, 450 0, 0, 0, 0, 1, 0, 0, 0; 451 ldlt.compute(mat); 452 VERIFY(ldlt.info()==NumericalIssue); 453 VERIFY_IS_NOT_APPROX(mat,ldlt.reconstructedMatrix()); 454 } 455 } 456 457 template<typename MatrixType> void cholesky_verify_assert() 458 { 459 MatrixType tmp; 460 461 LLT<MatrixType> llt; 462 VERIFY_RAISES_ASSERT(llt.matrixL()) 463 VERIFY_RAISES_ASSERT(llt.matrixU()) 464 VERIFY_RAISES_ASSERT(llt.solve(tmp)) 465 VERIFY_RAISES_ASSERT(llt.solveInPlace(&tmp)) 466 467 LDLT<MatrixType> ldlt; 468 VERIFY_RAISES_ASSERT(ldlt.matrixL()) 469 VERIFY_RAISES_ASSERT(ldlt.permutationP()) 470 VERIFY_RAISES_ASSERT(ldlt.vectorD()) 471 VERIFY_RAISES_ASSERT(ldlt.isPositive()) 472 VERIFY_RAISES_ASSERT(ldlt.isNegative()) 473 VERIFY_RAISES_ASSERT(ldlt.solve(tmp)) 474 VERIFY_RAISES_ASSERT(ldlt.solveInPlace(&tmp)) 475 } 476 477 void test_cholesky() 478 { 479 int s = 0; 480 for(int i = 0; i < g_repeat; i++) { 481 CALL_SUBTEST_1( cholesky(Matrix<double,1,1>()) ); 482 CALL_SUBTEST_3( cholesky(Matrix2d()) ); 483 CALL_SUBTEST_3( cholesky_bug241(Matrix2d()) ); 484 CALL_SUBTEST_3( cholesky_definiteness(Matrix2d()) ); 485 CALL_SUBTEST_4( cholesky(Matrix3f()) ); 486 CALL_SUBTEST_5( cholesky(Matrix4d()) ); 487 488 s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE); 489 CALL_SUBTEST_2( cholesky(MatrixXd(s,s)) ); 490 TEST_SET_BUT_UNUSED_VARIABLE(s) 491 492 s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE/2); 493 CALL_SUBTEST_6( cholesky_cplx(MatrixXcd(s,s)) ); 494 TEST_SET_BUT_UNUSED_VARIABLE(s) 495 } 496 497 CALL_SUBTEST_4( cholesky_verify_assert<Matrix3f>() ); 498 CALL_SUBTEST_7( cholesky_verify_assert<Matrix3d>() ); 499 CALL_SUBTEST_8( cholesky_verify_assert<MatrixXf>() ); 500 CALL_SUBTEST_2( cholesky_verify_assert<MatrixXd>() ); 501 502 // Test problem size constructors 503 CALL_SUBTEST_9( LLT<MatrixXf>(10) ); 504 CALL_SUBTEST_9( LDLT<MatrixXf>(10) ); 505 506 CALL_SUBTEST_2( cholesky_faillure_cases<void>() ); 507 508 TEST_SET_BUT_UNUSED_VARIABLE(nb_temporaries) 509 } 510
