• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
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 
test_chol_update(const MatrixType & symm)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 
cholesky(const MatrixType & m)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 
cholesky_cplx(const MatrixType & m)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
cholesky_bug241(const MatrixType & m)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
cholesky_definiteness(const MatrixType & m)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 
cholesky_verify_assert()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 
test_cholesky()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 }
405