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