$treeview $search $mathjax
Eigen
3.2.5
$projectbrief
|
$projectbrief
|
$searchbox |
00001 // This file is part of Eigen, a lightweight C++ template library 00002 // for linear algebra. 00003 // 00004 // Copyright (C) 2008-2009 Gael Guennebaud <gael.guennebaud@inria.fr> 00005 // Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com> 00006 // 00007 // This Source Code Form is subject to the terms of the Mozilla 00008 // Public License v. 2.0. If a copy of the MPL was not distributed 00009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 00010 00011 #ifndef EIGEN_COLPIVOTINGHOUSEHOLDERQR_H 00012 #define EIGEN_COLPIVOTINGHOUSEHOLDERQR_H 00013 00014 namespace Eigen { 00015 00037 template<typename _MatrixType> class ColPivHouseholderQR 00038 { 00039 public: 00040 00041 typedef _MatrixType MatrixType; 00042 enum { 00043 RowsAtCompileTime = MatrixType::RowsAtCompileTime, 00044 ColsAtCompileTime = MatrixType::ColsAtCompileTime, 00045 Options = MatrixType::Options, 00046 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 00047 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime 00048 }; 00049 typedef typename MatrixType::Scalar Scalar; 00050 typedef typename MatrixType::RealScalar RealScalar; 00051 typedef typename MatrixType::Index Index; 00052 typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, Options, MaxRowsAtCompileTime, MaxRowsAtCompileTime> MatrixQType; 00053 typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType; 00054 typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationType; 00055 typedef typename internal::plain_row_type<MatrixType, Index>::type IntRowVectorType; 00056 typedef typename internal::plain_row_type<MatrixType>::type RowVectorType; 00057 typedef typename internal::plain_row_type<MatrixType, RealScalar>::type RealRowVectorType; 00058 typedef HouseholderSequence<MatrixType,typename internal::remove_all<typename HCoeffsType::ConjugateReturnType>::type> HouseholderSequenceType; 00059 00060 private: 00061 00062 typedef typename PermutationType::Index PermIndexType; 00063 00064 public: 00065 00072 ColPivHouseholderQR() 00073 : m_qr(), 00074 m_hCoeffs(), 00075 m_colsPermutation(), 00076 m_colsTranspositions(), 00077 m_temp(), 00078 m_colSqNorms(), 00079 m_isInitialized(false), 00080 m_usePrescribedThreshold(false) {} 00081 00088 ColPivHouseholderQR(Index rows, Index cols) 00089 : m_qr(rows, cols), 00090 m_hCoeffs((std::min)(rows,cols)), 00091 m_colsPermutation(PermIndexType(cols)), 00092 m_colsTranspositions(cols), 00093 m_temp(cols), 00094 m_colSqNorms(cols), 00095 m_isInitialized(false), 00096 m_usePrescribedThreshold(false) {} 00097 00110 ColPivHouseholderQR(const MatrixType& matrix) 00111 : m_qr(matrix.rows(), matrix.cols()), 00112 m_hCoeffs((std::min)(matrix.rows(),matrix.cols())), 00113 m_colsPermutation(PermIndexType(matrix.cols())), 00114 m_colsTranspositions(matrix.cols()), 00115 m_temp(matrix.cols()), 00116 m_colSqNorms(matrix.cols()), 00117 m_isInitialized(false), 00118 m_usePrescribedThreshold(false) 00119 { 00120 compute(matrix); 00121 } 00122 00140 template<typename Rhs> 00141 inline const internal::solve_retval<ColPivHouseholderQR, Rhs> 00142 solve(const MatrixBase<Rhs>& b) const 00143 { 00144 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); 00145 return internal::solve_retval<ColPivHouseholderQR, Rhs>(*this, b.derived()); 00146 } 00147 00148 HouseholderSequenceType householderQ(void) const; 00149 HouseholderSequenceType matrixQ(void) const 00150 { 00151 return householderQ(); 00152 } 00153 00156 const MatrixType& matrixQR() const 00157 { 00158 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); 00159 return m_qr; 00160 } 00161 00171 const MatrixType& matrixR() const 00172 { 00173 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); 00174 return m_qr; 00175 } 00176 00177 ColPivHouseholderQR& compute(const MatrixType& matrix); 00178 00180 const PermutationType& colsPermutation() const 00181 { 00182 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); 00183 return m_colsPermutation; 00184 } 00185 00199 typename MatrixType::RealScalar absDeterminant() const; 00200 00213 typename MatrixType::RealScalar logAbsDeterminant() const; 00214 00221 inline Index rank() const 00222 { 00223 using std::abs; 00224 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); 00225 RealScalar premultiplied_threshold = abs(m_maxpivot) * threshold(); 00226 Index result = 0; 00227 for(Index i = 0; i < m_nonzero_pivots; ++i) 00228 result += (abs(m_qr.coeff(i,i)) > premultiplied_threshold); 00229 return result; 00230 } 00231 00238 inline Index dimensionOfKernel() const 00239 { 00240 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); 00241 return cols() - rank(); 00242 } 00243 00251 inline bool isInjective() const 00252 { 00253 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); 00254 return rank() == cols(); 00255 } 00256 00264 inline bool isSurjective() const 00265 { 00266 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); 00267 return rank() == rows(); 00268 } 00269 00276 inline bool isInvertible() const 00277 { 00278 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); 00279 return isInjective() && isSurjective(); 00280 } 00281 00287 inline const 00288 internal::solve_retval<ColPivHouseholderQR, typename MatrixType::IdentityReturnType> 00289 inverse() const 00290 { 00291 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); 00292 return internal::solve_retval<ColPivHouseholderQR,typename MatrixType::IdentityReturnType> 00293 (*this, MatrixType::Identity(m_qr.rows(), m_qr.cols())); 00294 } 00295 00296 inline Index rows() const { return m_qr.rows(); } 00297 inline Index cols() const { return m_qr.cols(); } 00298 00303 const HCoeffsType& hCoeffs() const { return m_hCoeffs; } 00304 00322 ColPivHouseholderQR& setThreshold(const RealScalar& threshold) 00323 { 00324 m_usePrescribedThreshold = true; 00325 m_prescribedThreshold = threshold; 00326 return *this; 00327 } 00328 00337 ColPivHouseholderQR& setThreshold(Default_t) 00338 { 00339 m_usePrescribedThreshold = false; 00340 return *this; 00341 } 00342 00347 RealScalar threshold() const 00348 { 00349 eigen_assert(m_isInitialized || m_usePrescribedThreshold); 00350 return m_usePrescribedThreshold ? m_prescribedThreshold 00351 // this formula comes from experimenting (see "LU precision tuning" thread on the list) 00352 // and turns out to be identical to Higham's formula used already in LDLt. 00353 : NumTraits<Scalar>::epsilon() * RealScalar(m_qr.diagonalSize()); 00354 } 00355 00363 inline Index nonzeroPivots() const 00364 { 00365 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); 00366 return m_nonzero_pivots; 00367 } 00368 00372 RealScalar maxPivot() const { return m_maxpivot; } 00373 00380 ComputationInfo info() const 00381 { 00382 eigen_assert(m_isInitialized && "Decomposition is not initialized."); 00383 return Success; 00384 } 00385 00386 protected: 00387 00388 static void check_template_parameters() 00389 { 00390 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar); 00391 } 00392 00393 MatrixType m_qr; 00394 HCoeffsType m_hCoeffs; 00395 PermutationType m_colsPermutation; 00396 IntRowVectorType m_colsTranspositions; 00397 RowVectorType m_temp; 00398 RealRowVectorType m_colSqNorms; 00399 bool m_isInitialized, m_usePrescribedThreshold; 00400 RealScalar m_prescribedThreshold, m_maxpivot; 00401 Index m_nonzero_pivots; 00402 Index m_det_pq; 00403 }; 00404 00405 template<typename MatrixType> 00406 typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::absDeterminant() const 00407 { 00408 using std::abs; 00409 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); 00410 eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!"); 00411 return abs(m_qr.diagonal().prod()); 00412 } 00413 00414 template<typename MatrixType> 00415 typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::logAbsDeterminant() const 00416 { 00417 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); 00418 eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!"); 00419 return m_qr.diagonal().cwiseAbs().array().log().sum(); 00420 } 00421 00428 template<typename MatrixType> 00429 ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const MatrixType& matrix) 00430 { 00431 check_template_parameters(); 00432 00433 using std::abs; 00434 Index rows = matrix.rows(); 00435 Index cols = matrix.cols(); 00436 Index size = matrix.diagonalSize(); 00437 00438 // the column permutation is stored as int indices, so just to be sure: 00439 eigen_assert(cols<=NumTraits<int>::highest()); 00440 00441 m_qr = matrix; 00442 m_hCoeffs.resize(size); 00443 00444 m_temp.resize(cols); 00445 00446 m_colsTranspositions.resize(matrix.cols()); 00447 Index number_of_transpositions = 0; 00448 00449 m_colSqNorms.resize(cols); 00450 for(Index k = 0; k < cols; ++k) 00451 m_colSqNorms.coeffRef(k) = m_qr.col(k).squaredNorm(); 00452 00453 RealScalar threshold_helper = m_colSqNorms.maxCoeff() * numext::abs2(NumTraits<Scalar>::epsilon()) / RealScalar(rows); 00454 00455 m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case) 00456 m_maxpivot = RealScalar(0); 00457 00458 for(Index k = 0; k < size; ++k) 00459 { 00460 // first, we look up in our table m_colSqNorms which column has the biggest squared norm 00461 Index biggest_col_index; 00462 RealScalar biggest_col_sq_norm = m_colSqNorms.tail(cols-k).maxCoeff(&biggest_col_index); 00463 biggest_col_index += k; 00464 00465 // since our table m_colSqNorms accumulates imprecision at every step, we must now recompute 00466 // the actual squared norm of the selected column. 00467 // Note that not doing so does result in solve() sometimes returning inf/nan values 00468 // when running the unit test with 1000 repetitions. 00469 biggest_col_sq_norm = m_qr.col(biggest_col_index).tail(rows-k).squaredNorm(); 00470 00471 // we store that back into our table: it can't hurt to correct our table. 00472 m_colSqNorms.coeffRef(biggest_col_index) = biggest_col_sq_norm; 00473 00474 // Track the number of meaningful pivots but do not stop the decomposition to make 00475 // sure that the initial matrix is properly reproduced. See bug 941. 00476 if(m_nonzero_pivots==size && biggest_col_sq_norm < threshold_helper * RealScalar(rows-k)) 00477 m_nonzero_pivots = k; 00478 00479 // apply the transposition to the columns 00480 m_colsTranspositions.coeffRef(k) = biggest_col_index; 00481 if(k != biggest_col_index) { 00482 m_qr.col(k).swap(m_qr.col(biggest_col_index)); 00483 std::swap(m_colSqNorms.coeffRef(k), m_colSqNorms.coeffRef(biggest_col_index)); 00484 ++number_of_transpositions; 00485 } 00486 00487 // generate the householder vector, store it below the diagonal 00488 RealScalar beta; 00489 m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta); 00490 00491 // apply the householder transformation to the diagonal coefficient 00492 m_qr.coeffRef(k,k) = beta; 00493 00494 // remember the maximum absolute value of diagonal coefficients 00495 if(abs(beta) > m_maxpivot) m_maxpivot = abs(beta); 00496 00497 // apply the householder transformation 00498 m_qr.bottomRightCorner(rows-k, cols-k-1) 00499 .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1)); 00500 00501 // update our table of squared norms of the columns 00502 m_colSqNorms.tail(cols-k-1) -= m_qr.row(k).tail(cols-k-1).cwiseAbs2(); 00503 } 00504 00505 m_colsPermutation.setIdentity(PermIndexType(cols)); 00506 for(PermIndexType k = 0; k < size/*m_nonzero_pivots*/; ++k) 00507 m_colsPermutation.applyTranspositionOnTheRight(k, PermIndexType(m_colsTranspositions.coeff(k))); 00508 00509 m_det_pq = (number_of_transpositions%2) ? -1 : 1; 00510 m_isInitialized = true; 00511 00512 return *this; 00513 } 00514 00515 namespace internal { 00516 00517 template<typename _MatrixType, typename Rhs> 00518 struct solve_retval<ColPivHouseholderQR<_MatrixType>, Rhs> 00519 : solve_retval_base<ColPivHouseholderQR<_MatrixType>, Rhs> 00520 { 00521 EIGEN_MAKE_SOLVE_HELPERS(ColPivHouseholderQR<_MatrixType>,Rhs) 00522 00523 template<typename Dest> void evalTo(Dest& dst) const 00524 { 00525 eigen_assert(rhs().rows() == dec().rows()); 00526 00527 const Index cols = dec().cols(), 00528 nonzero_pivots = dec().nonzeroPivots(); 00529 00530 if(nonzero_pivots == 0) 00531 { 00532 dst.setZero(); 00533 return; 00534 } 00535 00536 typename Rhs::PlainObject c(rhs()); 00537 00538 // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T 00539 c.applyOnTheLeft(householderSequence(dec().matrixQR(), dec().hCoeffs()) 00540 .setLength(dec().nonzeroPivots()) 00541 .transpose() 00542 ); 00543 00544 dec().matrixR() 00545 .topLeftCorner(nonzero_pivots, nonzero_pivots) 00546 .template triangularView<Upper>() 00547 .solveInPlace(c.topRows(nonzero_pivots)); 00548 00549 for(Index i = 0; i < nonzero_pivots; ++i) dst.row(dec().colsPermutation().indices().coeff(i)) = c.row(i); 00550 for(Index i = nonzero_pivots; i < cols; ++i) dst.row(dec().colsPermutation().indices().coeff(i)).setZero(); 00551 } 00552 }; 00553 00554 } // end namespace internal 00555 00559 template<typename MatrixType> 00560 typename ColPivHouseholderQR<MatrixType>::HouseholderSequenceType ColPivHouseholderQR<MatrixType> 00561 ::householderQ() const 00562 { 00563 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); 00564 return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate()); 00565 } 00566 00571 template<typename Derived> 00572 const ColPivHouseholderQR<typename MatrixBase<Derived>::PlainObject> 00573 MatrixBase<Derived>::colPivHouseholderQr() const 00574 { 00575 return ColPivHouseholderQR<PlainObject>(eval()); 00576 } 00577 00578 } // end namespace Eigen 00579 00580 #endif // EIGEN_COLPIVOTINGHOUSEHOLDERQR_H