Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 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 | |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 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 | { |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 91 | SquareMatrixType symmUp = symm.template triangularView<Upper>(); |
| 92 | SquareMatrixType symmLo = symm.template triangularView<Lower>(); |
| 93 | |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 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())); |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 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)); |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 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; |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 183 | |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 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 | } |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 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 | |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 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 | { |
| 324 | mat << 1, 0, 0, -1; |
| 325 | LDLT<MatrixType> ldlt(mat); |
| 326 | VERIFY(!ldlt.isNegative()); |
| 327 | VERIFY(!ldlt.isPositive()); |
| 328 | } |
| 329 | { |
| 330 | mat << 1, 2, 2, 1; |
| 331 | LDLT<MatrixType> ldlt(mat); |
| 332 | VERIFY(!ldlt.isNegative()); |
| 333 | VERIFY(!ldlt.isPositive()); |
| 334 | } |
| 335 | { |
| 336 | mat << 0, 0, 0, 0; |
| 337 | LDLT<MatrixType> ldlt(mat); |
| 338 | VERIFY(ldlt.isNegative()); |
| 339 | VERIFY(ldlt.isPositive()); |
| 340 | } |
| 341 | { |
| 342 | mat << 0, 0, 0, 1; |
| 343 | LDLT<MatrixType> ldlt(mat); |
| 344 | VERIFY(!ldlt.isNegative()); |
| 345 | VERIFY(ldlt.isPositive()); |
| 346 | } |
| 347 | { |
| 348 | mat << -1, 0, 0, 0; |
| 349 | LDLT<MatrixType> ldlt(mat); |
| 350 | VERIFY(ldlt.isNegative()); |
| 351 | VERIFY(!ldlt.isPositive()); |
| 352 | } |
| 353 | } |
| 354 | |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 355 | template<typename MatrixType> void cholesky_verify_assert() |
| 356 | { |
| 357 | MatrixType tmp; |
| 358 | |
| 359 | LLT<MatrixType> llt; |
| 360 | VERIFY_RAISES_ASSERT(llt.matrixL()) |
| 361 | VERIFY_RAISES_ASSERT(llt.matrixU()) |
| 362 | VERIFY_RAISES_ASSERT(llt.solve(tmp)) |
| 363 | VERIFY_RAISES_ASSERT(llt.solveInPlace(&tmp)) |
| 364 | |
| 365 | LDLT<MatrixType> ldlt; |
| 366 | VERIFY_RAISES_ASSERT(ldlt.matrixL()) |
| 367 | VERIFY_RAISES_ASSERT(ldlt.permutationP()) |
| 368 | VERIFY_RAISES_ASSERT(ldlt.vectorD()) |
| 369 | VERIFY_RAISES_ASSERT(ldlt.isPositive()) |
| 370 | VERIFY_RAISES_ASSERT(ldlt.isNegative()) |
| 371 | VERIFY_RAISES_ASSERT(ldlt.solve(tmp)) |
| 372 | VERIFY_RAISES_ASSERT(ldlt.solveInPlace(&tmp)) |
| 373 | } |
| 374 | |
| 375 | void test_cholesky() |
| 376 | { |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 377 | int s = 0; |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 378 | for(int i = 0; i < g_repeat; i++) { |
| 379 | CALL_SUBTEST_1( cholesky(Matrix<double,1,1>()) ); |
| 380 | CALL_SUBTEST_3( cholesky(Matrix2d()) ); |
| 381 | CALL_SUBTEST_3( cholesky_bug241(Matrix2d()) ); |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 382 | CALL_SUBTEST_3( cholesky_definiteness(Matrix2d()) ); |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 383 | CALL_SUBTEST_4( cholesky(Matrix3f()) ); |
| 384 | CALL_SUBTEST_5( cholesky(Matrix4d()) ); |
| 385 | s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE); |
| 386 | CALL_SUBTEST_2( cholesky(MatrixXd(s,s)) ); |
| 387 | s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE/2); |
| 388 | CALL_SUBTEST_6( cholesky_cplx(MatrixXcd(s,s)) ); |
| 389 | } |
| 390 | |
| 391 | CALL_SUBTEST_4( cholesky_verify_assert<Matrix3f>() ); |
| 392 | CALL_SUBTEST_7( cholesky_verify_assert<Matrix3d>() ); |
| 393 | CALL_SUBTEST_8( cholesky_verify_assert<MatrixXf>() ); |
| 394 | CALL_SUBTEST_2( cholesky_verify_assert<MatrixXd>() ); |
| 395 | |
| 396 | // Test problem size constructors |
| 397 | CALL_SUBTEST_9( LLT<MatrixXf>(10) ); |
| 398 | CALL_SUBTEST_9( LDLT<MatrixXf>(10) ); |
| 399 | |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 400 | TEST_SET_BUT_UNUSED_VARIABLE(s) |
| 401 | TEST_SET_BUT_UNUSED_VARIABLE(nb_temporaries) |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 402 | } |