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 (at) inria.fr> 5 // Copyright (C) 2009 Benoit Jacob <jacob.benoit.1 (at) 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> 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> 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> 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> 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> 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> 275 EIGEN_DONT_INLINE Scalar zero() { return Scalar(0); } 276 277 // workaround aggressive optimization in ICC 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> 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> 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> 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> 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> 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
