1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2008-2014 Gael Guennebaud <gael.guennebaud@inria.fr>
5 // Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
6 //
7 // This Source Code Form is subject to the terms of the Mozilla
8 // Public License v. 2.0. If a copy of the MPL was not distributed
9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
10
11 #ifndef SVD_DEFAULT
12 #error a macro SVD_DEFAULT(MatrixType) must be defined prior to including svd_common.h
13 #endif
14
15 #ifndef SVD_FOR_MIN_NORM
16 #error a macro SVD_FOR_MIN_NORM(MatrixType) must be defined prior to including svd_common.h
17 #endif
18
19 #include "svd_fill.h"
20
21 // Check that the matrix m is properly reconstructed and that the U and V factors are unitary
22 // The SVD must have already been computed.
23 template<typename SvdType, typename MatrixType>
svd_check_full(const MatrixType & m,const SvdType & svd)24 void svd_check_full(const MatrixType& m, const SvdType& svd)
25 {
26 typedef typename MatrixType::Index Index;
27 Index rows = m.rows();
28 Index cols = m.cols();
29
30 enum {
31 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
32 ColsAtCompileTime = MatrixType::ColsAtCompileTime
33 };
34
35 typedef typename MatrixType::Scalar Scalar;
36 typedef typename MatrixType::RealScalar RealScalar;
37 typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime> MatrixUType;
38 typedef Matrix<Scalar, ColsAtCompileTime, ColsAtCompileTime> MatrixVType;
39
40 MatrixType sigma = MatrixType::Zero(rows,cols);
41 sigma.diagonal() = svd.singularValues().template cast<Scalar>();
42 MatrixUType u = svd.matrixU();
43 MatrixVType v = svd.matrixV();
44 RealScalar scaling = m.cwiseAbs().maxCoeff();
45 if(scaling<(std::numeric_limits<RealScalar>::min)())
46 {
47 VERIFY(sigma.cwiseAbs().maxCoeff() <= (std::numeric_limits<RealScalar>::min)());
48 }
49 else
50 {
51 VERIFY_IS_APPROX(m/scaling, u * (sigma/scaling) * v.adjoint());
52 }
53 VERIFY_IS_UNITARY(u);
54 VERIFY_IS_UNITARY(v);
55 }
56
57 // Compare partial SVD defined by computationOptions to a full SVD referenceSvd
58 template<typename SvdType, typename MatrixType>
svd_compare_to_full(const MatrixType & m,unsigned int computationOptions,const SvdType & referenceSvd)59 void svd_compare_to_full(const MatrixType& m,
60 unsigned int computationOptions,
61 const SvdType& referenceSvd)
62 {
63 typedef typename MatrixType::RealScalar RealScalar;
64 Index rows = m.rows();
65 Index cols = m.cols();
66 Index diagSize = (std::min)(rows, cols);
67 RealScalar prec = test_precision<RealScalar>();
68
69 SvdType svd(m, computationOptions);
70
71 VERIFY_IS_APPROX(svd.singularValues(), referenceSvd.singularValues());
72
73 if(computationOptions & (ComputeFullV|ComputeThinV))
74 {
75 VERIFY( (svd.matrixV().adjoint()*svd.matrixV()).isIdentity(prec) );
76 VERIFY_IS_APPROX( svd.matrixV().leftCols(diagSize) * svd.singularValues().asDiagonal() * svd.matrixV().leftCols(diagSize).adjoint(),
77 referenceSvd.matrixV().leftCols(diagSize) * referenceSvd.singularValues().asDiagonal() * referenceSvd.matrixV().leftCols(diagSize).adjoint());
78 }
79
80 if(computationOptions & (ComputeFullU|ComputeThinU))
81 {
82 VERIFY( (svd.matrixU().adjoint()*svd.matrixU()).isIdentity(prec) );
83 VERIFY_IS_APPROX( svd.matrixU().leftCols(diagSize) * svd.singularValues().cwiseAbs2().asDiagonal() * svd.matrixU().leftCols(diagSize).adjoint(),
84 referenceSvd.matrixU().leftCols(diagSize) * referenceSvd.singularValues().cwiseAbs2().asDiagonal() * referenceSvd.matrixU().leftCols(diagSize).adjoint());
85 }
86
87 // The following checks are not critical.
88 // For instance, with Dived&Conquer SVD, if only the factor 'V' is computedt then different matrix-matrix product implementation will be used
89 // and the resulting 'V' factor might be significantly different when the SVD decomposition is not unique, especially with single precision float.
90 ++g_test_level;
91 if(computationOptions & ComputeFullU) VERIFY_IS_APPROX(svd.matrixU(), referenceSvd.matrixU());
92 if(computationOptions & ComputeThinU) VERIFY_IS_APPROX(svd.matrixU(), referenceSvd.matrixU().leftCols(diagSize));
93 if(computationOptions & ComputeFullV) VERIFY_IS_APPROX(svd.matrixV().cwiseAbs(), referenceSvd.matrixV().cwiseAbs());
94 if(computationOptions & ComputeThinV) VERIFY_IS_APPROX(svd.matrixV(), referenceSvd.matrixV().leftCols(diagSize));
95 --g_test_level;
96 }
97
98 //
99 template<typename SvdType, typename MatrixType>
svd_least_square(const MatrixType & m,unsigned int computationOptions)100 void svd_least_square(const MatrixType& m, unsigned int computationOptions)
101 {
102 typedef typename MatrixType::Scalar Scalar;
103 typedef typename MatrixType::RealScalar RealScalar;
104 typedef typename MatrixType::Index Index;
105 Index rows = m.rows();
106 Index cols = m.cols();
107
108 enum {
109 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
110 ColsAtCompileTime = MatrixType::ColsAtCompileTime
111 };
112
113 typedef Matrix<Scalar, RowsAtCompileTime, Dynamic> RhsType;
114 typedef Matrix<Scalar, ColsAtCompileTime, Dynamic> SolutionType;
115
116 RhsType rhs = RhsType::Random(rows, internal::random<Index>(1, cols));
117 SvdType svd(m, computationOptions);
118
119 if(internal::is_same<RealScalar,double>::value) svd.setThreshold(1e-8);
120 else if(internal::is_same<RealScalar,float>::value) svd.setThreshold(2e-4);
121
122 SolutionType x = svd.solve(rhs);
123
124 RealScalar residual = (m*x-rhs).norm();
125 RealScalar rhs_norm = rhs.norm();
126 if(!test_isMuchSmallerThan(residual,rhs.norm()))
127 {
128 // ^^^ If the residual is very small, then we have an exact solution, so we are already good.
129
130 // evaluate normal equation which works also for least-squares solutions
131 if(internal::is_same<RealScalar,double>::value || svd.rank()==m.diagonal().size())
132 {
133 using std::sqrt;
134 // This test is not stable with single precision.
135 // This is probably because squaring m signicantly affects the precision.
136 if(internal::is_same<RealScalar,float>::value) ++g_test_level;
137
138 VERIFY_IS_APPROX(m.adjoint()*(m*x),m.adjoint()*rhs);
139
140 if(internal::is_same<RealScalar,float>::value) --g_test_level;
141 }
142
143 // Check that there is no significantly better solution in the neighborhood of x
144 for(Index k=0;k<x.rows();++k)
145 {
146 using std::abs;
147
148 SolutionType y(x);
149 y.row(k) = (RealScalar(1)+2*NumTraits<RealScalar>::epsilon())*x.row(k);
150 RealScalar residual_y = (m*y-rhs).norm();
151 VERIFY( test_isMuchSmallerThan(abs(residual_y-residual), rhs_norm) || residual < residual_y );
152 if(internal::is_same<RealScalar,float>::value) ++g_test_level;
153 VERIFY( test_isApprox(residual_y,residual) || residual < residual_y );
154 if(internal::is_same<RealScalar,float>::value) --g_test_level;
155
156 y.row(k) = (RealScalar(1)-2*NumTraits<RealScalar>::epsilon())*x.row(k);
157 residual_y = (m*y-rhs).norm();
158 VERIFY( test_isMuchSmallerThan(abs(residual_y-residual), rhs_norm) || residual < residual_y );
159 if(internal::is_same<RealScalar,float>::value) ++g_test_level;
160 VERIFY( test_isApprox(residual_y,residual) || residual < residual_y );
161 if(internal::is_same<RealScalar,float>::value) --g_test_level;
162 }
163 }
164 }
165
166 // check minimal norm solutions, the inoput matrix m is only used to recover problem size
167 template<typename MatrixType>
svd_min_norm(const MatrixType & m,unsigned int computationOptions)168 void svd_min_norm(const MatrixType& m, unsigned int computationOptions)
169 {
170 typedef typename MatrixType::Scalar Scalar;
171 typedef typename MatrixType::Index Index;
172 Index cols = m.cols();
173
174 enum {
175 ColsAtCompileTime = MatrixType::ColsAtCompileTime
176 };
177
178 typedef Matrix<Scalar, ColsAtCompileTime, Dynamic> SolutionType;
179
180 // generate a full-rank m x n problem with m<n
181 enum {
182 RankAtCompileTime2 = ColsAtCompileTime==Dynamic ? Dynamic : (ColsAtCompileTime)/2+1,
183 RowsAtCompileTime3 = ColsAtCompileTime==Dynamic ? Dynamic : ColsAtCompileTime+1
184 };
185 typedef Matrix<Scalar, RankAtCompileTime2, ColsAtCompileTime> MatrixType2;
186 typedef Matrix<Scalar, RankAtCompileTime2, 1> RhsType2;
187 typedef Matrix<Scalar, ColsAtCompileTime, RankAtCompileTime2> MatrixType2T;
188 Index rank = RankAtCompileTime2==Dynamic ? internal::random<Index>(1,cols) : Index(RankAtCompileTime2);
189 MatrixType2 m2(rank,cols);
190 int guard = 0;
191 do {
192 m2.setRandom();
193 } while(SVD_FOR_MIN_NORM(MatrixType2)(m2).setThreshold(test_precision<Scalar>()).rank()!=rank && (++guard)<10);
194 VERIFY(guard<10);
195
196 RhsType2 rhs2 = RhsType2::Random(rank);
197 // use QR to find a reference minimal norm solution
198 HouseholderQR<MatrixType2T> qr(m2.adjoint());
199 Matrix<Scalar,Dynamic,1> tmp = qr.matrixQR().topLeftCorner(rank,rank).template triangularView<Upper>().adjoint().solve(rhs2);
200 tmp.conservativeResize(cols);
201 tmp.tail(cols-rank).setZero();
202 SolutionType x21 = qr.householderQ() * tmp;
203 // now check with SVD
204 SVD_FOR_MIN_NORM(MatrixType2) svd2(m2, computationOptions);
205 SolutionType x22 = svd2.solve(rhs2);
206 VERIFY_IS_APPROX(m2*x21, rhs2);
207 VERIFY_IS_APPROX(m2*x22, rhs2);
208 VERIFY_IS_APPROX(x21, x22);
209
210 // Now check with a rank deficient matrix
211 typedef Matrix<Scalar, RowsAtCompileTime3, ColsAtCompileTime> MatrixType3;
212 typedef Matrix<Scalar, RowsAtCompileTime3, 1> RhsType3;
213 Index rows3 = RowsAtCompileTime3==Dynamic ? internal::random<Index>(rank+1,2*cols) : Index(RowsAtCompileTime3);
214 Matrix<Scalar,RowsAtCompileTime3,Dynamic> C = Matrix<Scalar,RowsAtCompileTime3,Dynamic>::Random(rows3,rank);
215 MatrixType3 m3 = C * m2;
216 RhsType3 rhs3 = C * rhs2;
217 SVD_FOR_MIN_NORM(MatrixType3) svd3(m3, computationOptions);
218 SolutionType x3 = svd3.solve(rhs3);
219 VERIFY_IS_APPROX(m3*x3, rhs3);
220 VERIFY_IS_APPROX(m3*x21, rhs3);
221 VERIFY_IS_APPROX(m2*x3, rhs2);
222 VERIFY_IS_APPROX(x21, x3);
223 }
224
225 // Check full, compare_to_full, least_square, and min_norm for all possible compute-options
226 template<typename SvdType, typename MatrixType>
svd_test_all_computation_options(const MatrixType & m,bool full_only)227 void svd_test_all_computation_options(const MatrixType& m, bool full_only)
228 {
229 // if (QRPreconditioner == NoQRPreconditioner && m.rows() != m.cols())
230 // return;
231 SvdType fullSvd(m, ComputeFullU|ComputeFullV);
232 CALL_SUBTEST(( svd_check_full(m, fullSvd) ));
233 CALL_SUBTEST(( svd_least_square<SvdType>(m, ComputeFullU | ComputeFullV) ));
234 CALL_SUBTEST(( svd_min_norm(m, ComputeFullU | ComputeFullV) ));
235
236 #if defined __INTEL_COMPILER
237 // remark #111: statement is unreachable
238 #pragma warning disable 111
239 #endif
240 if(full_only)
241 return;
242
243 CALL_SUBTEST(( svd_compare_to_full(m, ComputeFullU, fullSvd) ));
244 CALL_SUBTEST(( svd_compare_to_full(m, ComputeFullV, fullSvd) ));
245 CALL_SUBTEST(( svd_compare_to_full(m, 0, fullSvd) ));
246
247 if (MatrixType::ColsAtCompileTime == Dynamic) {
248 // thin U/V are only available with dynamic number of columns
249 CALL_SUBTEST(( svd_compare_to_full(m, ComputeFullU|ComputeThinV, fullSvd) ));
250 CALL_SUBTEST(( svd_compare_to_full(m, ComputeThinV, fullSvd) ));
251 CALL_SUBTEST(( svd_compare_to_full(m, ComputeThinU|ComputeFullV, fullSvd) ));
252 CALL_SUBTEST(( svd_compare_to_full(m, ComputeThinU , fullSvd) ));
253 CALL_SUBTEST(( svd_compare_to_full(m, ComputeThinU|ComputeThinV, fullSvd) ));
254
255 CALL_SUBTEST(( svd_least_square<SvdType>(m, ComputeFullU | ComputeThinV) ));
256 CALL_SUBTEST(( svd_least_square<SvdType>(m, ComputeThinU | ComputeFullV) ));
257 CALL_SUBTEST(( svd_least_square<SvdType>(m, ComputeThinU | ComputeThinV) ));
258
259 CALL_SUBTEST(( svd_min_norm(m, ComputeFullU | ComputeThinV) ));
260 CALL_SUBTEST(( svd_min_norm(m, ComputeThinU | ComputeFullV) ));
261 CALL_SUBTEST(( svd_min_norm(m, ComputeThinU | ComputeThinV) ));
262
263 // test reconstruction
264 typedef typename MatrixType::Index Index;
265 Index diagSize = (std::min)(m.rows(), m.cols());
266 SvdType svd(m, ComputeThinU | ComputeThinV);
267 VERIFY_IS_APPROX(m, svd.matrixU().leftCols(diagSize) * svd.singularValues().asDiagonal() * svd.matrixV().leftCols(diagSize).adjoint());
268 }
269 }
270
271
272 // work around stupid msvc error when constructing at compile time an expression that involves
273 // a division by zero, even if the numeric type has floating point
274 template<typename Scalar>
zero()275 EIGEN_DONT_INLINE Scalar zero() { return Scalar(0); }
276
277 // workaround aggressive optimization in ICC
sub(T a,T b)278 template<typename T> EIGEN_DONT_INLINE T sub(T a, T b) { return a - b; }
279
280 // all this function does is verify we don't iterate infinitely on nan/inf values
281 template<typename SvdType, typename MatrixType>
svd_inf_nan()282 void svd_inf_nan()
283 {
284 SvdType svd;
285 typedef typename MatrixType::Scalar Scalar;
286 Scalar some_inf = Scalar(1) / zero<Scalar>();
287 VERIFY(sub(some_inf, some_inf) != sub(some_inf, some_inf));
288 svd.compute(MatrixType::Constant(10,10,some_inf), ComputeFullU | ComputeFullV);
289
290 Scalar nan = std::numeric_limits<Scalar>::quiet_NaN();
291 VERIFY(nan != nan);
292 svd.compute(MatrixType::Constant(10,10,nan), ComputeFullU | ComputeFullV);
293
294 MatrixType m = MatrixType::Zero(10,10);
295 m(internal::random<int>(0,9), internal::random<int>(0,9)) = some_inf;
296 svd.compute(m, ComputeFullU | ComputeFullV);
297
298 m = MatrixType::Zero(10,10);
299 m(internal::random<int>(0,9), internal::random<int>(0,9)) = nan;
300 svd.compute(m, ComputeFullU | ComputeFullV);
301
302 // regression test for bug 791
303 m.resize(3,3);
304 m << 0, 2*NumTraits<Scalar>::epsilon(), 0.5,
305 0, -0.5, 0,
306 nan, 0, 0;
307 svd.compute(m, ComputeFullU | ComputeFullV);
308
309 m.resize(4,4);
310 m << 1, 0, 0, 0,
311 0, 3, 1, 2e-308,
312 1, 0, 1, nan,
313 0, nan, nan, 0;
314 svd.compute(m, ComputeFullU | ComputeFullV);
315 }
316
317 // Regression test for bug 286: JacobiSVD loops indefinitely with some
318 // matrices containing denormal numbers.
319 template<typename>
svd_underoverflow()320 void svd_underoverflow()
321 {
322 #if defined __INTEL_COMPILER
323 // shut up warning #239: floating point underflow
324 #pragma warning push
325 #pragma warning disable 239
326 #endif
327 Matrix2d M;
328 M << -7.90884e-313, -4.94e-324,
329 0, 5.60844e-313;
330 SVD_DEFAULT(Matrix2d) svd;
331 svd.compute(M,ComputeFullU|ComputeFullV);
332 CALL_SUBTEST( svd_check_full(M,svd) );
333
334 // Check all 2x2 matrices made with the following coefficients:
335 VectorXd value_set(9);
336 value_set << 0, 1, -1, 5.60844e-313, -5.60844e-313, 4.94e-324, -4.94e-324, -4.94e-223, 4.94e-223;
337 Array4i id(0,0,0,0);
338 int k = 0;
339 do
340 {
341 M << value_set(id(0)), value_set(id(1)), value_set(id(2)), value_set(id(3));
342 svd.compute(M,ComputeFullU|ComputeFullV);
343 CALL_SUBTEST( svd_check_full(M,svd) );
344
345 id(k)++;
346 if(id(k)>=value_set.size())
347 {
348 while(k<3 && id(k)>=value_set.size()) id(++k)++;
349 id.head(k).setZero();
350 k=0;
351 }
352
353 } while((id<int(value_set.size())).all());
354
355 #if defined __INTEL_COMPILER
356 #pragma warning pop
357 #endif
358
359 // Check for overflow:
360 Matrix3d M3;
361 M3 << 4.4331978442502944e+307, -5.8585363752028680e+307, 6.4527017443412964e+307,
362 3.7841695601406358e+307, 2.4331702789740617e+306, -3.5235707140272905e+307,
363 -8.7190887618028355e+307, -7.3453213709232193e+307, -2.4367363684472105e+307;
364
365 SVD_DEFAULT(Matrix3d) svd3;
366 svd3.compute(M3,ComputeFullU|ComputeFullV); // just check we don't loop indefinitely
367 CALL_SUBTEST( svd_check_full(M3,svd3) );
368 }
369
370 // void jacobisvd(const MatrixType& a = MatrixType(), bool pickrandom = true)
371
372 template<typename MatrixType>
svd_all_trivial_2x2(void (* cb)(const MatrixType &,bool))373 void svd_all_trivial_2x2( void (*cb)(const MatrixType&,bool) )
374 {
375 MatrixType M;
376 VectorXd value_set(3);
377 value_set << 0, 1, -1;
378 Array4i id(0,0,0,0);
379 int k = 0;
380 do
381 {
382 M << value_set(id(0)), value_set(id(1)), value_set(id(2)), value_set(id(3));
383
384 cb(M,false);
385
386 id(k)++;
387 if(id(k)>=value_set.size())
388 {
389 while(k<3 && id(k)>=value_set.size()) id(++k)++;
390 id.head(k).setZero();
391 k=0;
392 }
393
394 } while((id<int(value_set.size())).all());
395 }
396
397 template<typename>
svd_preallocate()398 void svd_preallocate()
399 {
400 Vector3f v(3.f, 2.f, 1.f);
401 MatrixXf m = v.asDiagonal();
402
403 internal::set_is_malloc_allowed(false);
404 VERIFY_RAISES_ASSERT(VectorXf tmp(10);)
405 SVD_DEFAULT(MatrixXf) svd;
406 internal::set_is_malloc_allowed(true);
407 svd.compute(m);
408 VERIFY_IS_APPROX(svd.singularValues(), v);
409
410 SVD_DEFAULT(MatrixXf) svd2(3,3);
411 internal::set_is_malloc_allowed(false);
412 svd2.compute(m);
413 internal::set_is_malloc_allowed(true);
414 VERIFY_IS_APPROX(svd2.singularValues(), v);
415 VERIFY_RAISES_ASSERT(svd2.matrixU());
416 VERIFY_RAISES_ASSERT(svd2.matrixV());
417 svd2.compute(m, ComputeFullU | ComputeFullV);
418 VERIFY_IS_APPROX(svd2.matrixU(), Matrix3f::Identity());
419 VERIFY_IS_APPROX(svd2.matrixV(), Matrix3f::Identity());
420 internal::set_is_malloc_allowed(false);
421 svd2.compute(m);
422 internal::set_is_malloc_allowed(true);
423
424 SVD_DEFAULT(MatrixXf) svd3(3,3,ComputeFullU|ComputeFullV);
425 internal::set_is_malloc_allowed(false);
426 svd2.compute(m);
427 internal::set_is_malloc_allowed(true);
428 VERIFY_IS_APPROX(svd2.singularValues(), v);
429 VERIFY_IS_APPROX(svd2.matrixU(), Matrix3f::Identity());
430 VERIFY_IS_APPROX(svd2.matrixV(), Matrix3f::Identity());
431 internal::set_is_malloc_allowed(false);
432 svd2.compute(m, ComputeFullU|ComputeFullV);
433 internal::set_is_malloc_allowed(true);
434 }
435
436 template<typename SvdType,typename MatrixType>
svd_verify_assert(const MatrixType & m)437 void svd_verify_assert(const MatrixType& m)
438 {
439 typedef typename MatrixType::Scalar Scalar;
440 typedef typename MatrixType::Index Index;
441 Index rows = m.rows();
442 Index cols = m.cols();
443
444 enum {
445 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
446 ColsAtCompileTime = MatrixType::ColsAtCompileTime
447 };
448
449 typedef Matrix<Scalar, RowsAtCompileTime, 1> RhsType;
450 RhsType rhs(rows);
451 SvdType svd;
452 VERIFY_RAISES_ASSERT(svd.matrixU())
453 VERIFY_RAISES_ASSERT(svd.singularValues())
454 VERIFY_RAISES_ASSERT(svd.matrixV())
455 VERIFY_RAISES_ASSERT(svd.solve(rhs))
456 MatrixType a = MatrixType::Zero(rows, cols);
457 a.setZero();
458 svd.compute(a, 0);
459 VERIFY_RAISES_ASSERT(svd.matrixU())
460 VERIFY_RAISES_ASSERT(svd.matrixV())
461 svd.singularValues();
462 VERIFY_RAISES_ASSERT(svd.solve(rhs))
463
464 if (ColsAtCompileTime == Dynamic)
465 {
466 svd.compute(a, ComputeThinU);
467 svd.matrixU();
468 VERIFY_RAISES_ASSERT(svd.matrixV())
469 VERIFY_RAISES_ASSERT(svd.solve(rhs))
470 svd.compute(a, ComputeThinV);
471 svd.matrixV();
472 VERIFY_RAISES_ASSERT(svd.matrixU())
473 VERIFY_RAISES_ASSERT(svd.solve(rhs))
474 }
475 else
476 {
477 VERIFY_RAISES_ASSERT(svd.compute(a, ComputeThinU))
478 VERIFY_RAISES_ASSERT(svd.compute(a, ComputeThinV))
479 }
480 }
481
482 #undef SVD_DEFAULT
483 #undef SVD_FOR_MIN_NORM
484