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