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