1// This file is part of Eigen, a lightweight C++ template library 2// for linear algebra. 3// 4// Copyright (C) 2008-2009 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 EIGEN_COLPIVOTINGHOUSEHOLDERQR_H 12#define EIGEN_COLPIVOTINGHOUSEHOLDERQR_H 13 14namespace Eigen { 15 16/** \ingroup QR_Module 17 * 18 * \class ColPivHouseholderQR 19 * 20 * \brief Householder rank-revealing QR decomposition of a matrix with column-pivoting 21 * 22 * \param MatrixType the type of the matrix of which we are computing the QR decomposition 23 * 24 * This class performs a rank-revealing QR decomposition of a matrix \b A into matrices \b P, \b Q and \b R 25 * such that 26 * \f[ 27 * \mathbf{A} \, \mathbf{P} = \mathbf{Q} \, \mathbf{R} 28 * \f] 29 * by using Householder transformations. Here, \b P is a permutation matrix, \b Q a unitary matrix and \b R an 30 * upper triangular matrix. 31 * 32 * This decomposition performs column pivoting in order to be rank-revealing and improve 33 * numerical stability. It is slower than HouseholderQR, and faster than FullPivHouseholderQR. 34 * 35 * \sa MatrixBase::colPivHouseholderQr() 36 */ 37template<typename _MatrixType> class ColPivHouseholderQR 38{ 39 public: 40 41 typedef _MatrixType MatrixType; 42 enum { 43 RowsAtCompileTime = MatrixType::RowsAtCompileTime, 44 ColsAtCompileTime = MatrixType::ColsAtCompileTime, 45 Options = MatrixType::Options, 46 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 47 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime 48 }; 49 typedef typename MatrixType::Scalar Scalar; 50 typedef typename MatrixType::RealScalar RealScalar; 51 typedef typename MatrixType::Index Index; 52 typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, Options, MaxRowsAtCompileTime, MaxRowsAtCompileTime> MatrixQType; 53 typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType; 54 typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationType; 55 typedef typename internal::plain_row_type<MatrixType, Index>::type IntRowVectorType; 56 typedef typename internal::plain_row_type<MatrixType>::type RowVectorType; 57 typedef typename internal::plain_row_type<MatrixType, RealScalar>::type RealRowVectorType; 58 typedef HouseholderSequence<MatrixType,typename internal::remove_all<typename HCoeffsType::ConjugateReturnType>::type> HouseholderSequenceType; 59 60 private: 61 62 typedef typename PermutationType::Index PermIndexType; 63 64 public: 65 66 /** 67 * \brief Default Constructor. 68 * 69 * The default constructor is useful in cases in which the user intends to 70 * perform decompositions via ColPivHouseholderQR::compute(const MatrixType&). 71 */ 72 ColPivHouseholderQR() 73 : m_qr(), 74 m_hCoeffs(), 75 m_colsPermutation(), 76 m_colsTranspositions(), 77 m_temp(), 78 m_colSqNorms(), 79 m_isInitialized(false), 80 m_usePrescribedThreshold(false) {} 81 82 /** \brief Default Constructor with memory preallocation 83 * 84 * Like the default constructor but with preallocation of the internal data 85 * according to the specified problem \a size. 86 * \sa ColPivHouseholderQR() 87 */ 88 ColPivHouseholderQR(Index rows, Index cols) 89 : m_qr(rows, cols), 90 m_hCoeffs((std::min)(rows,cols)), 91 m_colsPermutation(PermIndexType(cols)), 92 m_colsTranspositions(cols), 93 m_temp(cols), 94 m_colSqNorms(cols), 95 m_isInitialized(false), 96 m_usePrescribedThreshold(false) {} 97 98 /** \brief Constructs a QR factorization from a given matrix 99 * 100 * This constructor computes the QR factorization of the matrix \a matrix by calling 101 * the method compute(). It is a short cut for: 102 * 103 * \code 104 * ColPivHouseholderQR<MatrixType> qr(matrix.rows(), matrix.cols()); 105 * qr.compute(matrix); 106 * \endcode 107 * 108 * \sa compute() 109 */ 110 ColPivHouseholderQR(const MatrixType& matrix) 111 : m_qr(matrix.rows(), matrix.cols()), 112 m_hCoeffs((std::min)(matrix.rows(),matrix.cols())), 113 m_colsPermutation(PermIndexType(matrix.cols())), 114 m_colsTranspositions(matrix.cols()), 115 m_temp(matrix.cols()), 116 m_colSqNorms(matrix.cols()), 117 m_isInitialized(false), 118 m_usePrescribedThreshold(false) 119 { 120 compute(matrix); 121 } 122 123 /** This method finds a solution x to the equation Ax=b, where A is the matrix of which 124 * *this is the QR decomposition, if any exists. 125 * 126 * \param b the right-hand-side of the equation to solve. 127 * 128 * \returns a solution. 129 * 130 * \note The case where b is a matrix is not yet implemented. Also, this 131 * code is space inefficient. 132 * 133 * \note_about_checking_solutions 134 * 135 * \note_about_arbitrary_choice_of_solution 136 * 137 * Example: \include ColPivHouseholderQR_solve.cpp 138 * Output: \verbinclude ColPivHouseholderQR_solve.out 139 */ 140 template<typename Rhs> 141 inline const internal::solve_retval<ColPivHouseholderQR, Rhs> 142 solve(const MatrixBase<Rhs>& b) const 143 { 144 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); 145 return internal::solve_retval<ColPivHouseholderQR, Rhs>(*this, b.derived()); 146 } 147 148 HouseholderSequenceType householderQ(void) const; 149 HouseholderSequenceType matrixQ(void) const 150 { 151 return householderQ(); 152 } 153 154 /** \returns a reference to the matrix where the Householder QR decomposition is stored 155 */ 156 const MatrixType& matrixQR() const 157 { 158 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); 159 return m_qr; 160 } 161 162 /** \returns a reference to the matrix where the result Householder QR is stored 163 * \warning The strict lower part of this matrix contains internal values. 164 * Only the upper triangular part should be referenced. To get it, use 165 * \code matrixR().template triangularView<Upper>() \endcode 166 * For rank-deficient matrices, use 167 * \code 168 * matrixR().topLeftCorner(rank(), rank()).template triangularView<Upper>() 169 * \endcode 170 */ 171 const MatrixType& matrixR() const 172 { 173 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); 174 return m_qr; 175 } 176 177 ColPivHouseholderQR& compute(const MatrixType& matrix); 178 179 /** \returns a const reference to the column permutation matrix */ 180 const PermutationType& colsPermutation() const 181 { 182 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); 183 return m_colsPermutation; 184 } 185 186 /** \returns the absolute value of the determinant of the matrix of which 187 * *this is the QR decomposition. It has only linear complexity 188 * (that is, O(n) where n is the dimension of the square matrix) 189 * as the QR decomposition has already been computed. 190 * 191 * \note This is only for square matrices. 192 * 193 * \warning a determinant can be very big or small, so for matrices 194 * of large enough dimension, there is a risk of overflow/underflow. 195 * One way to work around that is to use logAbsDeterminant() instead. 196 * 197 * \sa logAbsDeterminant(), MatrixBase::determinant() 198 */ 199 typename MatrixType::RealScalar absDeterminant() const; 200 201 /** \returns the natural log of the absolute value of the determinant of the matrix of which 202 * *this is the QR decomposition. It has only linear complexity 203 * (that is, O(n) where n is the dimension of the square matrix) 204 * as the QR decomposition has already been computed. 205 * 206 * \note This is only for square matrices. 207 * 208 * \note This method is useful to work around the risk of overflow/underflow that's inherent 209 * to determinant computation. 210 * 211 * \sa absDeterminant(), MatrixBase::determinant() 212 */ 213 typename MatrixType::RealScalar logAbsDeterminant() const; 214 215 /** \returns the rank of the matrix of which *this is the QR decomposition. 216 * 217 * \note This method has to determine which pivots should be considered nonzero. 218 * For that, it uses the threshold value that you can control by calling 219 * setThreshold(const RealScalar&). 220 */ 221 inline Index rank() const 222 { 223 using std::abs; 224 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); 225 RealScalar premultiplied_threshold = abs(m_maxpivot) * threshold(); 226 Index result = 0; 227 for(Index i = 0; i < m_nonzero_pivots; ++i) 228 result += (abs(m_qr.coeff(i,i)) > premultiplied_threshold); 229 return result; 230 } 231 232 /** \returns the dimension of the kernel of the matrix of which *this is the QR decomposition. 233 * 234 * \note This method has to determine which pivots should be considered nonzero. 235 * For that, it uses the threshold value that you can control by calling 236 * setThreshold(const RealScalar&). 237 */ 238 inline Index dimensionOfKernel() const 239 { 240 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); 241 return cols() - rank(); 242 } 243 244 /** \returns true if the matrix of which *this is the QR decomposition represents an injective 245 * linear map, i.e. has trivial kernel; false otherwise. 246 * 247 * \note This method has to determine which pivots should be considered nonzero. 248 * For that, it uses the threshold value that you can control by calling 249 * setThreshold(const RealScalar&). 250 */ 251 inline bool isInjective() const 252 { 253 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); 254 return rank() == cols(); 255 } 256 257 /** \returns true if the matrix of which *this is the QR decomposition represents a surjective 258 * linear map; false otherwise. 259 * 260 * \note This method has to determine which pivots should be considered nonzero. 261 * For that, it uses the threshold value that you can control by calling 262 * setThreshold(const RealScalar&). 263 */ 264 inline bool isSurjective() const 265 { 266 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); 267 return rank() == rows(); 268 } 269 270 /** \returns true if the matrix of which *this is the QR decomposition is invertible. 271 * 272 * \note This method has to determine which pivots should be considered nonzero. 273 * For that, it uses the threshold value that you can control by calling 274 * setThreshold(const RealScalar&). 275 */ 276 inline bool isInvertible() const 277 { 278 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); 279 return isInjective() && isSurjective(); 280 } 281 282 /** \returns the inverse of the matrix of which *this is the QR decomposition. 283 * 284 * \note If this matrix is not invertible, the returned matrix has undefined coefficients. 285 * Use isInvertible() to first determine whether this matrix is invertible. 286 */ 287 inline const 288 internal::solve_retval<ColPivHouseholderQR, typename MatrixType::IdentityReturnType> 289 inverse() const 290 { 291 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); 292 return internal::solve_retval<ColPivHouseholderQR,typename MatrixType::IdentityReturnType> 293 (*this, MatrixType::Identity(m_qr.rows(), m_qr.cols())); 294 } 295 296 inline Index rows() const { return m_qr.rows(); } 297 inline Index cols() const { return m_qr.cols(); } 298 299 /** \returns a const reference to the vector of Householder coefficients used to represent the factor \c Q. 300 * 301 * For advanced uses only. 302 */ 303 const HCoeffsType& hCoeffs() const { return m_hCoeffs; } 304 305 /** Allows to prescribe a threshold to be used by certain methods, such as rank(), 306 * who need to determine when pivots are to be considered nonzero. This is not used for the 307 * QR decomposition itself. 308 * 309 * When it needs to get the threshold value, Eigen calls threshold(). By default, this 310 * uses a formula to automatically determine a reasonable threshold. 311 * Once you have called the present method setThreshold(const RealScalar&), 312 * your value is used instead. 313 * 314 * \param threshold The new value to use as the threshold. 315 * 316 * A pivot will be considered nonzero if its absolute value is strictly greater than 317 * \f$ \vert pivot \vert \leqslant threshold \times \vert maxpivot \vert \f$ 318 * where maxpivot is the biggest pivot. 319 * 320 * If you want to come back to the default behavior, call setThreshold(Default_t) 321 */ 322 ColPivHouseholderQR& setThreshold(const RealScalar& threshold) 323 { 324 m_usePrescribedThreshold = true; 325 m_prescribedThreshold = threshold; 326 return *this; 327 } 328 329 /** Allows to come back to the default behavior, letting Eigen use its default formula for 330 * determining the threshold. 331 * 332 * You should pass the special object Eigen::Default as parameter here. 333 * \code qr.setThreshold(Eigen::Default); \endcode 334 * 335 * See the documentation of setThreshold(const RealScalar&). 336 */ 337 ColPivHouseholderQR& setThreshold(Default_t) 338 { 339 m_usePrescribedThreshold = false; 340 return *this; 341 } 342 343 /** Returns the threshold that will be used by certain methods such as rank(). 344 * 345 * See the documentation of setThreshold(const RealScalar&). 346 */ 347 RealScalar threshold() const 348 { 349 eigen_assert(m_isInitialized || m_usePrescribedThreshold); 350 return m_usePrescribedThreshold ? m_prescribedThreshold 351 // this formula comes from experimenting (see "LU precision tuning" thread on the list) 352 // and turns out to be identical to Higham's formula used already in LDLt. 353 : NumTraits<Scalar>::epsilon() * RealScalar(m_qr.diagonalSize()); 354 } 355 356 /** \returns the number of nonzero pivots in the QR decomposition. 357 * Here nonzero is meant in the exact sense, not in a fuzzy sense. 358 * So that notion isn't really intrinsically interesting, but it is 359 * still useful when implementing algorithms. 360 * 361 * \sa rank() 362 */ 363 inline Index nonzeroPivots() const 364 { 365 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); 366 return m_nonzero_pivots; 367 } 368 369 /** \returns the absolute value of the biggest pivot, i.e. the biggest 370 * diagonal coefficient of R. 371 */ 372 RealScalar maxPivot() const { return m_maxpivot; } 373 374 /** \brief Reports whether the QR factorization was succesful. 375 * 376 * \note This function always returns \c Success. It is provided for compatibility 377 * with other factorization routines. 378 * \returns \c Success 379 */ 380 ComputationInfo info() const 381 { 382 eigen_assert(m_isInitialized && "Decomposition is not initialized."); 383 return Success; 384 } 385 386 protected: 387 MatrixType m_qr; 388 HCoeffsType m_hCoeffs; 389 PermutationType m_colsPermutation; 390 IntRowVectorType m_colsTranspositions; 391 RowVectorType m_temp; 392 RealRowVectorType m_colSqNorms; 393 bool m_isInitialized, m_usePrescribedThreshold; 394 RealScalar m_prescribedThreshold, m_maxpivot; 395 Index m_nonzero_pivots; 396 Index m_det_pq; 397}; 398 399template<typename MatrixType> 400typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::absDeterminant() const 401{ 402 using std::abs; 403 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); 404 eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!"); 405 return abs(m_qr.diagonal().prod()); 406} 407 408template<typename MatrixType> 409typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::logAbsDeterminant() const 410{ 411 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); 412 eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!"); 413 return m_qr.diagonal().cwiseAbs().array().log().sum(); 414} 415 416/** Performs the QR factorization of the given matrix \a matrix. The result of 417 * the factorization is stored into \c *this, and a reference to \c *this 418 * is returned. 419 * 420 * \sa class ColPivHouseholderQR, ColPivHouseholderQR(const MatrixType&) 421 */ 422template<typename MatrixType> 423ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const MatrixType& matrix) 424{ 425 using std::abs; 426 Index rows = matrix.rows(); 427 Index cols = matrix.cols(); 428 Index size = matrix.diagonalSize(); 429 430 // the column permutation is stored as int indices, so just to be sure: 431 eigen_assert(cols<=NumTraits<int>::highest()); 432 433 m_qr = matrix; 434 m_hCoeffs.resize(size); 435 436 m_temp.resize(cols); 437 438 m_colsTranspositions.resize(matrix.cols()); 439 Index number_of_transpositions = 0; 440 441 m_colSqNorms.resize(cols); 442 for(Index k = 0; k < cols; ++k) 443 m_colSqNorms.coeffRef(k) = m_qr.col(k).squaredNorm(); 444 445 RealScalar threshold_helper = m_colSqNorms.maxCoeff() * numext::abs2(NumTraits<Scalar>::epsilon()) / RealScalar(rows); 446 447 m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case) 448 m_maxpivot = RealScalar(0); 449 450 for(Index k = 0; k < size; ++k) 451 { 452 // first, we look up in our table m_colSqNorms which column has the biggest squared norm 453 Index biggest_col_index; 454 RealScalar biggest_col_sq_norm = m_colSqNorms.tail(cols-k).maxCoeff(&biggest_col_index); 455 biggest_col_index += k; 456 457 // since our table m_colSqNorms accumulates imprecision at every step, we must now recompute 458 // the actual squared norm of the selected column. 459 // Note that not doing so does result in solve() sometimes returning inf/nan values 460 // when running the unit test with 1000 repetitions. 461 biggest_col_sq_norm = m_qr.col(biggest_col_index).tail(rows-k).squaredNorm(); 462 463 // we store that back into our table: it can't hurt to correct our table. 464 m_colSqNorms.coeffRef(biggest_col_index) = biggest_col_sq_norm; 465 466 // if the current biggest column is smaller than epsilon times the initial biggest column, 467 // terminate to avoid generating nan/inf values. 468 // Note that here, if we test instead for "biggest == 0", we get a failure every 1000 (or so) 469 // repetitions of the unit test, with the result of solve() filled with large values of the order 470 // of 1/(size*epsilon). 471 if(biggest_col_sq_norm < threshold_helper * RealScalar(rows-k)) 472 { 473 m_nonzero_pivots = k; 474 m_hCoeffs.tail(size-k).setZero(); 475 m_qr.bottomRightCorner(rows-k,cols-k) 476 .template triangularView<StrictlyLower>() 477 .setZero(); 478 break; 479 } 480 481 // apply the transposition to the columns 482 m_colsTranspositions.coeffRef(k) = biggest_col_index; 483 if(k != biggest_col_index) { 484 m_qr.col(k).swap(m_qr.col(biggest_col_index)); 485 std::swap(m_colSqNorms.coeffRef(k), m_colSqNorms.coeffRef(biggest_col_index)); 486 ++number_of_transpositions; 487 } 488 489 // generate the householder vector, store it below the diagonal 490 RealScalar beta; 491 m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta); 492 493 // apply the householder transformation to the diagonal coefficient 494 m_qr.coeffRef(k,k) = beta; 495 496 // remember the maximum absolute value of diagonal coefficients 497 if(abs(beta) > m_maxpivot) m_maxpivot = abs(beta); 498 499 // apply the householder transformation 500 m_qr.bottomRightCorner(rows-k, cols-k-1) 501 .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1)); 502 503 // update our table of squared norms of the columns 504 m_colSqNorms.tail(cols-k-1) -= m_qr.row(k).tail(cols-k-1).cwiseAbs2(); 505 } 506 507 m_colsPermutation.setIdentity(PermIndexType(cols)); 508 for(PermIndexType k = 0; k < m_nonzero_pivots; ++k) 509 m_colsPermutation.applyTranspositionOnTheRight(k, PermIndexType(m_colsTranspositions.coeff(k))); 510 511 m_det_pq = (number_of_transpositions%2) ? -1 : 1; 512 m_isInitialized = true; 513 514 return *this; 515} 516 517namespace internal { 518 519template<typename _MatrixType, typename Rhs> 520struct solve_retval<ColPivHouseholderQR<_MatrixType>, Rhs> 521 : solve_retval_base<ColPivHouseholderQR<_MatrixType>, Rhs> 522{ 523 EIGEN_MAKE_SOLVE_HELPERS(ColPivHouseholderQR<_MatrixType>,Rhs) 524 525 template<typename Dest> void evalTo(Dest& dst) const 526 { 527 eigen_assert(rhs().rows() == dec().rows()); 528 529 const Index cols = dec().cols(), 530 nonzero_pivots = dec().nonzeroPivots(); 531 532 if(nonzero_pivots == 0) 533 { 534 dst.setZero(); 535 return; 536 } 537 538 typename Rhs::PlainObject c(rhs()); 539 540 // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T 541 c.applyOnTheLeft(householderSequence(dec().matrixQR(), dec().hCoeffs()) 542 .setLength(dec().nonzeroPivots()) 543 .transpose() 544 ); 545 546 dec().matrixR() 547 .topLeftCorner(nonzero_pivots, nonzero_pivots) 548 .template triangularView<Upper>() 549 .solveInPlace(c.topRows(nonzero_pivots)); 550 551 for(Index i = 0; i < nonzero_pivots; ++i) dst.row(dec().colsPermutation().indices().coeff(i)) = c.row(i); 552 for(Index i = nonzero_pivots; i < cols; ++i) dst.row(dec().colsPermutation().indices().coeff(i)).setZero(); 553 } 554}; 555 556} // end namespace internal 557 558/** \returns the matrix Q as a sequence of householder transformations */ 559template<typename MatrixType> 560typename ColPivHouseholderQR<MatrixType>::HouseholderSequenceType ColPivHouseholderQR<MatrixType> 561 ::householderQ() const 562{ 563 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); 564 return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate()).setLength(m_nonzero_pivots); 565} 566 567/** \return the column-pivoting Householder QR decomposition of \c *this. 568 * 569 * \sa class ColPivHouseholderQR 570 */ 571template<typename Derived> 572const ColPivHouseholderQR<typename MatrixBase<Derived>::PlainObject> 573MatrixBase<Derived>::colPivHouseholderQr() const 574{ 575 return ColPivHouseholderQR<PlainObject>(eval()); 576} 577 578} // end namespace Eigen 579 580#endif // EIGEN_COLPIVOTINGHOUSEHOLDERQR_H 581