source: pacpussensors/trunk/Vislab/lib3dv/eigen/test/cholesky.cpp@ 140

Last change on this file since 140 was 136, checked in by ldecherf, 8 years ago

Doc
File size: 13.2 KB
Line 
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
14static 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
29template<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
61template<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
241template<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
301template<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
319template<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
357template<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
377void 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}
Note: See TracBrowser for help on using the repository browser.