qr_colpivoting.cpp (13867B)
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 // 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 #include "main.h" 12 #include <Eigen/QR> 13 #include <Eigen/SVD> 14 #include "solverbase.h" 15 16 template <typename MatrixType> 17 void cod() { 18 STATIC_CHECK(( internal::is_same<typename CompleteOrthogonalDecomposition<MatrixType>::StorageIndex,int>::value )); 19 20 Index rows = internal::random<Index>(2, EIGEN_TEST_MAX_SIZE); 21 Index cols = internal::random<Index>(2, EIGEN_TEST_MAX_SIZE); 22 Index cols2 = internal::random<Index>(2, EIGEN_TEST_MAX_SIZE); 23 Index rank = internal::random<Index>(1, (std::min)(rows, cols) - 1); 24 25 typedef typename MatrixType::Scalar Scalar; 26 typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 27 MatrixType::RowsAtCompileTime> 28 MatrixQType; 29 MatrixType matrix; 30 createRandomPIMatrixOfRank(rank, rows, cols, matrix); 31 CompleteOrthogonalDecomposition<MatrixType> cod(matrix); 32 VERIFY(rank == cod.rank()); 33 VERIFY(cols - cod.rank() == cod.dimensionOfKernel()); 34 VERIFY(!cod.isInjective()); 35 VERIFY(!cod.isInvertible()); 36 VERIFY(!cod.isSurjective()); 37 38 MatrixQType q = cod.householderQ(); 39 VERIFY_IS_UNITARY(q); 40 41 MatrixType z = cod.matrixZ(); 42 VERIFY_IS_UNITARY(z); 43 44 MatrixType t; 45 t.setZero(rows, cols); 46 t.topLeftCorner(rank, rank) = 47 cod.matrixT().topLeftCorner(rank, rank).template triangularView<Upper>(); 48 49 MatrixType c = q * t * z * cod.colsPermutation().inverse(); 50 VERIFY_IS_APPROX(matrix, c); 51 52 check_solverbase<MatrixType, MatrixType>(matrix, cod, rows, cols, cols2); 53 54 // Verify that we get the same minimum-norm solution as the SVD. 55 MatrixType exact_solution = MatrixType::Random(cols, cols2); 56 MatrixType rhs = matrix * exact_solution; 57 MatrixType cod_solution = cod.solve(rhs); 58 JacobiSVD<MatrixType> svd(matrix, ComputeThinU | ComputeThinV); 59 MatrixType svd_solution = svd.solve(rhs); 60 VERIFY_IS_APPROX(cod_solution, svd_solution); 61 62 MatrixType pinv = cod.pseudoInverse(); 63 VERIFY_IS_APPROX(cod_solution, pinv * rhs); 64 } 65 66 template <typename MatrixType, int Cols2> 67 void cod_fixedsize() { 68 enum { 69 Rows = MatrixType::RowsAtCompileTime, 70 Cols = MatrixType::ColsAtCompileTime 71 }; 72 typedef typename MatrixType::Scalar Scalar; 73 typedef CompleteOrthogonalDecomposition<Matrix<Scalar, Rows, Cols> > COD; 74 int rank = internal::random<int>(1, (std::min)(int(Rows), int(Cols)) - 1); 75 Matrix<Scalar, Rows, Cols> matrix; 76 createRandomPIMatrixOfRank(rank, Rows, Cols, matrix); 77 COD cod(matrix); 78 VERIFY(rank == cod.rank()); 79 VERIFY(Cols - cod.rank() == cod.dimensionOfKernel()); 80 VERIFY(cod.isInjective() == (rank == Rows)); 81 VERIFY(cod.isSurjective() == (rank == Cols)); 82 VERIFY(cod.isInvertible() == (cod.isInjective() && cod.isSurjective())); 83 84 check_solverbase<Matrix<Scalar, Cols, Cols2>, Matrix<Scalar, Rows, Cols2> >(matrix, cod, Rows, Cols, Cols2); 85 86 // Verify that we get the same minimum-norm solution as the SVD. 87 Matrix<Scalar, Cols, Cols2> exact_solution; 88 exact_solution.setRandom(Cols, Cols2); 89 Matrix<Scalar, Rows, Cols2> rhs = matrix * exact_solution; 90 Matrix<Scalar, Cols, Cols2> cod_solution = cod.solve(rhs); 91 JacobiSVD<MatrixType> svd(matrix, ComputeFullU | ComputeFullV); 92 Matrix<Scalar, Cols, Cols2> svd_solution = svd.solve(rhs); 93 VERIFY_IS_APPROX(cod_solution, svd_solution); 94 95 typename Inverse<COD>::PlainObject pinv = cod.pseudoInverse(); 96 VERIFY_IS_APPROX(cod_solution, pinv * rhs); 97 } 98 99 template<typename MatrixType> void qr() 100 { 101 using std::sqrt; 102 103 STATIC_CHECK(( internal::is_same<typename ColPivHouseholderQR<MatrixType>::StorageIndex,int>::value )); 104 105 Index rows = internal::random<Index>(2,EIGEN_TEST_MAX_SIZE), cols = internal::random<Index>(2,EIGEN_TEST_MAX_SIZE), cols2 = internal::random<Index>(2,EIGEN_TEST_MAX_SIZE); 106 Index rank = internal::random<Index>(1, (std::min)(rows, cols)-1); 107 108 typedef typename MatrixType::Scalar Scalar; 109 typedef typename MatrixType::RealScalar RealScalar; 110 typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> MatrixQType; 111 MatrixType m1; 112 createRandomPIMatrixOfRank(rank,rows,cols,m1); 113 ColPivHouseholderQR<MatrixType> qr(m1); 114 VERIFY_IS_EQUAL(rank, qr.rank()); 115 VERIFY_IS_EQUAL(cols - qr.rank(), qr.dimensionOfKernel()); 116 VERIFY(!qr.isInjective()); 117 VERIFY(!qr.isInvertible()); 118 VERIFY(!qr.isSurjective()); 119 120 MatrixQType q = qr.householderQ(); 121 VERIFY_IS_UNITARY(q); 122 123 MatrixType r = qr.matrixQR().template triangularView<Upper>(); 124 MatrixType c = q * r * qr.colsPermutation().inverse(); 125 VERIFY_IS_APPROX(m1, c); 126 127 // Verify that the absolute value of the diagonal elements in R are 128 // non-increasing until they reach the singularity threshold. 129 RealScalar threshold = 130 sqrt(RealScalar(rows)) * numext::abs(r(0, 0)) * NumTraits<Scalar>::epsilon(); 131 for (Index i = 0; i < (std::min)(rows, cols) - 1; ++i) { 132 RealScalar x = numext::abs(r(i, i)); 133 RealScalar y = numext::abs(r(i + 1, i + 1)); 134 if (x < threshold && y < threshold) continue; 135 if (!test_isApproxOrLessThan(y, x)) { 136 for (Index j = 0; j < (std::min)(rows, cols); ++j) { 137 std::cout << "i = " << j << ", |r_ii| = " << numext::abs(r(j, j)) << std::endl; 138 } 139 std::cout << "Failure at i=" << i << ", rank=" << rank 140 << ", threshold=" << threshold << std::endl; 141 } 142 VERIFY_IS_APPROX_OR_LESS_THAN(y, x); 143 } 144 145 check_solverbase<MatrixType, MatrixType>(m1, qr, rows, cols, cols2); 146 147 { 148 MatrixType m2, m3; 149 Index size = rows; 150 do { 151 m1 = MatrixType::Random(size,size); 152 qr.compute(m1); 153 } while(!qr.isInvertible()); 154 MatrixType m1_inv = qr.inverse(); 155 m3 = m1 * MatrixType::Random(size,cols2); 156 m2 = qr.solve(m3); 157 VERIFY_IS_APPROX(m2, m1_inv*m3); 158 } 159 } 160 161 template<typename MatrixType, int Cols2> void qr_fixedsize() 162 { 163 using std::sqrt; 164 using std::abs; 165 enum { Rows = MatrixType::RowsAtCompileTime, Cols = MatrixType::ColsAtCompileTime }; 166 typedef typename MatrixType::Scalar Scalar; 167 typedef typename MatrixType::RealScalar RealScalar; 168 int rank = internal::random<int>(1, (std::min)(int(Rows), int(Cols))-1); 169 Matrix<Scalar,Rows,Cols> m1; 170 createRandomPIMatrixOfRank(rank,Rows,Cols,m1); 171 ColPivHouseholderQR<Matrix<Scalar,Rows,Cols> > qr(m1); 172 VERIFY_IS_EQUAL(rank, qr.rank()); 173 VERIFY_IS_EQUAL(Cols - qr.rank(), qr.dimensionOfKernel()); 174 VERIFY_IS_EQUAL(qr.isInjective(), (rank == Rows)); 175 VERIFY_IS_EQUAL(qr.isSurjective(), (rank == Cols)); 176 VERIFY_IS_EQUAL(qr.isInvertible(), (qr.isInjective() && qr.isSurjective())); 177 178 Matrix<Scalar,Rows,Cols> r = qr.matrixQR().template triangularView<Upper>(); 179 Matrix<Scalar,Rows,Cols> c = qr.householderQ() * r * qr.colsPermutation().inverse(); 180 VERIFY_IS_APPROX(m1, c); 181 182 check_solverbase<Matrix<Scalar,Cols,Cols2>, Matrix<Scalar,Rows,Cols2> >(m1, qr, Rows, Cols, Cols2); 183 184 // Verify that the absolute value of the diagonal elements in R are 185 // non-increasing until they reache the singularity threshold. 186 RealScalar threshold = 187 sqrt(RealScalar(Rows)) * (std::abs)(r(0, 0)) * NumTraits<Scalar>::epsilon(); 188 for (Index i = 0; i < (std::min)(int(Rows), int(Cols)) - 1; ++i) { 189 RealScalar x = numext::abs(r(i, i)); 190 RealScalar y = numext::abs(r(i + 1, i + 1)); 191 if (x < threshold && y < threshold) continue; 192 if (!test_isApproxOrLessThan(y, x)) { 193 for (Index j = 0; j < (std::min)(int(Rows), int(Cols)); ++j) { 194 std::cout << "i = " << j << ", |r_ii| = " << numext::abs(r(j, j)) << std::endl; 195 } 196 std::cout << "Failure at i=" << i << ", rank=" << rank 197 << ", threshold=" << threshold << std::endl; 198 } 199 VERIFY_IS_APPROX_OR_LESS_THAN(y, x); 200 } 201 } 202 203 // This test is meant to verify that pivots are chosen such that 204 // even for a graded matrix, the diagonal of R falls of roughly 205 // monotonically until it reaches the threshold for singularity. 206 // We use the so-called Kahan matrix, which is a famous counter-example 207 // for rank-revealing QR. See 208 // http://www.netlib.org/lapack/lawnspdf/lawn176.pdf 209 // page 3 for more detail. 210 template<typename MatrixType> void qr_kahan_matrix() 211 { 212 using std::sqrt; 213 using std::abs; 214 typedef typename MatrixType::Scalar Scalar; 215 typedef typename MatrixType::RealScalar RealScalar; 216 217 Index rows = 300, cols = rows; 218 219 MatrixType m1; 220 m1.setZero(rows,cols); 221 RealScalar s = std::pow(NumTraits<RealScalar>::epsilon(), 1.0 / rows); 222 RealScalar c = std::sqrt(1 - s*s); 223 RealScalar pow_s_i(1.0); // pow(s,i) 224 for (Index i = 0; i < rows; ++i) { 225 m1(i, i) = pow_s_i; 226 m1.row(i).tail(rows - i - 1) = -pow_s_i * c * MatrixType::Ones(1, rows - i - 1); 227 pow_s_i *= s; 228 } 229 m1 = (m1 + m1.transpose()).eval(); 230 ColPivHouseholderQR<MatrixType> qr(m1); 231 MatrixType r = qr.matrixQR().template triangularView<Upper>(); 232 233 RealScalar threshold = 234 std::sqrt(RealScalar(rows)) * numext::abs(r(0, 0)) * NumTraits<Scalar>::epsilon(); 235 for (Index i = 0; i < (std::min)(rows, cols) - 1; ++i) { 236 RealScalar x = numext::abs(r(i, i)); 237 RealScalar y = numext::abs(r(i + 1, i + 1)); 238 if (x < threshold && y < threshold) continue; 239 if (!test_isApproxOrLessThan(y, x)) { 240 for (Index j = 0; j < (std::min)(rows, cols); ++j) { 241 std::cout << "i = " << j << ", |r_ii| = " << numext::abs(r(j, j)) << std::endl; 242 } 243 std::cout << "Failure at i=" << i << ", rank=" << qr.rank() 244 << ", threshold=" << threshold << std::endl; 245 } 246 VERIFY_IS_APPROX_OR_LESS_THAN(y, x); 247 } 248 } 249 250 template<typename MatrixType> void qr_invertible() 251 { 252 using std::log; 253 using std::abs; 254 typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; 255 typedef typename MatrixType::Scalar Scalar; 256 257 int size = internal::random<int>(10,50); 258 259 MatrixType m1(size, size), m2(size, size), m3(size, size); 260 m1 = MatrixType::Random(size,size); 261 262 if (internal::is_same<RealScalar,float>::value) 263 { 264 // let's build a matrix more stable to inverse 265 MatrixType a = MatrixType::Random(size,size*2); 266 m1 += a * a.adjoint(); 267 } 268 269 ColPivHouseholderQR<MatrixType> qr(m1); 270 271 check_solverbase<MatrixType, MatrixType>(m1, qr, size, size, size); 272 273 // now construct a matrix with prescribed determinant 274 m1.setZero(); 275 for(int i = 0; i < size; i++) m1(i,i) = internal::random<Scalar>(); 276 RealScalar absdet = abs(m1.diagonal().prod()); 277 m3 = qr.householderQ(); // get a unitary 278 m1 = m3 * m1 * m3; 279 qr.compute(m1); 280 VERIFY_IS_APPROX(absdet, qr.absDeterminant()); 281 VERIFY_IS_APPROX(log(absdet), qr.logAbsDeterminant()); 282 } 283 284 template<typename MatrixType> void qr_verify_assert() 285 { 286 MatrixType tmp; 287 288 ColPivHouseholderQR<MatrixType> qr; 289 VERIFY_RAISES_ASSERT(qr.matrixQR()) 290 VERIFY_RAISES_ASSERT(qr.solve(tmp)) 291 VERIFY_RAISES_ASSERT(qr.transpose().solve(tmp)) 292 VERIFY_RAISES_ASSERT(qr.adjoint().solve(tmp)) 293 VERIFY_RAISES_ASSERT(qr.householderQ()) 294 VERIFY_RAISES_ASSERT(qr.dimensionOfKernel()) 295 VERIFY_RAISES_ASSERT(qr.isInjective()) 296 VERIFY_RAISES_ASSERT(qr.isSurjective()) 297 VERIFY_RAISES_ASSERT(qr.isInvertible()) 298 VERIFY_RAISES_ASSERT(qr.inverse()) 299 VERIFY_RAISES_ASSERT(qr.absDeterminant()) 300 VERIFY_RAISES_ASSERT(qr.logAbsDeterminant()) 301 } 302 303 template<typename MatrixType> void cod_verify_assert() 304 { 305 MatrixType tmp; 306 307 CompleteOrthogonalDecomposition<MatrixType> cod; 308 VERIFY_RAISES_ASSERT(cod.matrixQTZ()) 309 VERIFY_RAISES_ASSERT(cod.solve(tmp)) 310 VERIFY_RAISES_ASSERT(cod.transpose().solve(tmp)) 311 VERIFY_RAISES_ASSERT(cod.adjoint().solve(tmp)) 312 VERIFY_RAISES_ASSERT(cod.householderQ()) 313 VERIFY_RAISES_ASSERT(cod.dimensionOfKernel()) 314 VERIFY_RAISES_ASSERT(cod.isInjective()) 315 VERIFY_RAISES_ASSERT(cod.isSurjective()) 316 VERIFY_RAISES_ASSERT(cod.isInvertible()) 317 VERIFY_RAISES_ASSERT(cod.pseudoInverse()) 318 VERIFY_RAISES_ASSERT(cod.absDeterminant()) 319 VERIFY_RAISES_ASSERT(cod.logAbsDeterminant()) 320 } 321 322 EIGEN_DECLARE_TEST(qr_colpivoting) 323 { 324 for(int i = 0; i < g_repeat; i++) { 325 CALL_SUBTEST_1( qr<MatrixXf>() ); 326 CALL_SUBTEST_2( qr<MatrixXd>() ); 327 CALL_SUBTEST_3( qr<MatrixXcd>() ); 328 CALL_SUBTEST_4(( qr_fixedsize<Matrix<float,3,5>, 4 >() )); 329 CALL_SUBTEST_5(( qr_fixedsize<Matrix<double,6,2>, 3 >() )); 330 CALL_SUBTEST_5(( qr_fixedsize<Matrix<double,1,1>, 1 >() )); 331 } 332 333 for(int i = 0; i < g_repeat; i++) { 334 CALL_SUBTEST_1( cod<MatrixXf>() ); 335 CALL_SUBTEST_2( cod<MatrixXd>() ); 336 CALL_SUBTEST_3( cod<MatrixXcd>() ); 337 CALL_SUBTEST_4(( cod_fixedsize<Matrix<float,3,5>, 4 >() )); 338 CALL_SUBTEST_5(( cod_fixedsize<Matrix<double,6,2>, 3 >() )); 339 CALL_SUBTEST_5(( cod_fixedsize<Matrix<double,1,1>, 1 >() )); 340 } 341 342 for(int i = 0; i < g_repeat; i++) { 343 CALL_SUBTEST_1( qr_invertible<MatrixXf>() ); 344 CALL_SUBTEST_2( qr_invertible<MatrixXd>() ); 345 CALL_SUBTEST_6( qr_invertible<MatrixXcf>() ); 346 CALL_SUBTEST_3( qr_invertible<MatrixXcd>() ); 347 } 348 349 CALL_SUBTEST_7(qr_verify_assert<Matrix3f>()); 350 CALL_SUBTEST_8(qr_verify_assert<Matrix3d>()); 351 CALL_SUBTEST_1(qr_verify_assert<MatrixXf>()); 352 CALL_SUBTEST_2(qr_verify_assert<MatrixXd>()); 353 CALL_SUBTEST_6(qr_verify_assert<MatrixXcf>()); 354 CALL_SUBTEST_3(qr_verify_assert<MatrixXcd>()); 355 356 CALL_SUBTEST_7(cod_verify_assert<Matrix3f>()); 357 CALL_SUBTEST_8(cod_verify_assert<Matrix3d>()); 358 CALL_SUBTEST_1(cod_verify_assert<MatrixXf>()); 359 CALL_SUBTEST_2(cod_verify_assert<MatrixXd>()); 360 CALL_SUBTEST_6(cod_verify_assert<MatrixXcf>()); 361 CALL_SUBTEST_3(cod_verify_assert<MatrixXcd>()); 362 363 // Test problem size constructors 364 CALL_SUBTEST_9(ColPivHouseholderQR<MatrixXf>(10, 20)); 365 366 CALL_SUBTEST_1( qr_kahan_matrix<MatrixXf>() ); 367 CALL_SUBTEST_2( qr_kahan_matrix<MatrixXd>() ); 368 }