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