00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 #ifndef EIGEN_COLPIVOTINGHOUSEHOLDERQR_H
00027 #define EIGEN_COLPIVOTINGHOUSEHOLDERQR_H
00028
00029 namespace Eigen {
00030
00052 template<typename _MatrixType> class ColPivHouseholderQR
00053 {
00054 public:
00055
00056 typedef _MatrixType MatrixType;
00057 enum {
00058 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00059 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00060 Options = MatrixType::Options,
00061 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00062 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00063 };
00064 typedef typename MatrixType::Scalar Scalar;
00065 typedef typename MatrixType::RealScalar RealScalar;
00066 typedef typename MatrixType::Index Index;
00067 typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, Options, MaxRowsAtCompileTime, MaxRowsAtCompileTime> MatrixQType;
00068 typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
00069 typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationType;
00070 typedef typename internal::plain_row_type<MatrixType, Index>::type IntRowVectorType;
00071 typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
00072 typedef typename internal::plain_row_type<MatrixType, RealScalar>::type RealRowVectorType;
00073 typedef typename HouseholderSequence<MatrixType,HCoeffsType>::ConjugateReturnType HouseholderSequenceType;
00074
00081 ColPivHouseholderQR()
00082 : m_qr(),
00083 m_hCoeffs(),
00084 m_colsPermutation(),
00085 m_colsTranspositions(),
00086 m_temp(),
00087 m_colSqNorms(),
00088 m_isInitialized(false) {}
00089
00096 ColPivHouseholderQR(Index rows, Index cols)
00097 : m_qr(rows, cols),
00098 m_hCoeffs((std::min)(rows,cols)),
00099 m_colsPermutation(cols),
00100 m_colsTranspositions(cols),
00101 m_temp(cols),
00102 m_colSqNorms(cols),
00103 m_isInitialized(false),
00104 m_usePrescribedThreshold(false) {}
00105
00106 ColPivHouseholderQR(const MatrixType& matrix)
00107 : m_qr(matrix.rows(), matrix.cols()),
00108 m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
00109 m_colsPermutation(matrix.cols()),
00110 m_colsTranspositions(matrix.cols()),
00111 m_temp(matrix.cols()),
00112 m_colSqNorms(matrix.cols()),
00113 m_isInitialized(false),
00114 m_usePrescribedThreshold(false)
00115 {
00116 compute(matrix);
00117 }
00118
00136 template<typename Rhs>
00137 inline const internal::solve_retval<ColPivHouseholderQR, Rhs>
00138 solve(const MatrixBase<Rhs>& b) const
00139 {
00140 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00141 return internal::solve_retval<ColPivHouseholderQR, Rhs>(*this, b.derived());
00142 }
00143
00144 HouseholderSequenceType householderQ(void) const;
00145
00148 const MatrixType& matrixQR() const
00149 {
00150 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00151 return m_qr;
00152 }
00153
00154 ColPivHouseholderQR& compute(const MatrixType& matrix);
00155
00156 const PermutationType& colsPermutation() const
00157 {
00158 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00159 return m_colsPermutation;
00160 }
00161
00175 typename MatrixType::RealScalar absDeterminant() const;
00176
00189 typename MatrixType::RealScalar logAbsDeterminant() const;
00190
00197 inline Index rank() const
00198 {
00199 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00200 RealScalar premultiplied_threshold = internal::abs(m_maxpivot) * threshold();
00201 Index result = 0;
00202 for(Index i = 0; i < m_nonzero_pivots; ++i)
00203 result += (internal::abs(m_qr.coeff(i,i)) > premultiplied_threshold);
00204 return result;
00205 }
00206
00213 inline Index dimensionOfKernel() const
00214 {
00215 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00216 return cols() - rank();
00217 }
00218
00226 inline bool isInjective() const
00227 {
00228 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00229 return rank() == cols();
00230 }
00231
00239 inline bool isSurjective() const
00240 {
00241 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00242 return rank() == rows();
00243 }
00244
00251 inline bool isInvertible() const
00252 {
00253 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00254 return isInjective() && isSurjective();
00255 }
00256
00262 inline const
00263 internal::solve_retval<ColPivHouseholderQR, typename MatrixType::IdentityReturnType>
00264 inverse() const
00265 {
00266 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00267 return internal::solve_retval<ColPivHouseholderQR,typename MatrixType::IdentityReturnType>
00268 (*this, MatrixType::Identity(m_qr.rows(), m_qr.cols()));
00269 }
00270
00271 inline Index rows() const { return m_qr.rows(); }
00272 inline Index cols() const { return m_qr.cols(); }
00273 const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
00274
00292 ColPivHouseholderQR& setThreshold(const RealScalar& threshold)
00293 {
00294 m_usePrescribedThreshold = true;
00295 m_prescribedThreshold = threshold;
00296 return *this;
00297 }
00298
00307 ColPivHouseholderQR& setThreshold(Default_t)
00308 {
00309 m_usePrescribedThreshold = false;
00310 return *this;
00311 }
00312
00317 RealScalar threshold() const
00318 {
00319 eigen_assert(m_isInitialized || m_usePrescribedThreshold);
00320 return m_usePrescribedThreshold ? m_prescribedThreshold
00321
00322
00323 : NumTraits<Scalar>::epsilon() * m_qr.diagonalSize();
00324 }
00325
00333 inline Index nonzeroPivots() const
00334 {
00335 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00336 return m_nonzero_pivots;
00337 }
00338
00342 RealScalar maxPivot() const { return m_maxpivot; }
00343
00344 protected:
00345 MatrixType m_qr;
00346 HCoeffsType m_hCoeffs;
00347 PermutationType m_colsPermutation;
00348 IntRowVectorType m_colsTranspositions;
00349 RowVectorType m_temp;
00350 RealRowVectorType m_colSqNorms;
00351 bool m_isInitialized, m_usePrescribedThreshold;
00352 RealScalar m_prescribedThreshold, m_maxpivot;
00353 Index m_nonzero_pivots;
00354 Index m_det_pq;
00355 };
00356
00357 template<typename MatrixType>
00358 typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::absDeterminant() const
00359 {
00360 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00361 eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
00362 return internal::abs(m_qr.diagonal().prod());
00363 }
00364
00365 template<typename MatrixType>
00366 typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::logAbsDeterminant() const
00367 {
00368 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00369 eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
00370 return m_qr.diagonal().cwiseAbs().array().log().sum();
00371 }
00372
00373 template<typename MatrixType>
00374 ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const MatrixType& matrix)
00375 {
00376 Index rows = matrix.rows();
00377 Index cols = matrix.cols();
00378 Index size = matrix.diagonalSize();
00379
00380 m_qr = matrix;
00381 m_hCoeffs.resize(size);
00382
00383 m_temp.resize(cols);
00384
00385 m_colsTranspositions.resize(matrix.cols());
00386 Index number_of_transpositions = 0;
00387
00388 m_colSqNorms.resize(cols);
00389 for(Index k = 0; k < cols; ++k)
00390 m_colSqNorms.coeffRef(k) = m_qr.col(k).squaredNorm();
00391
00392 RealScalar threshold_helper = m_colSqNorms.maxCoeff() * internal::abs2(NumTraits<Scalar>::epsilon()) / RealScalar(rows);
00393
00394 m_nonzero_pivots = size;
00395 m_maxpivot = RealScalar(0);
00396
00397 for(Index k = 0; k < size; ++k)
00398 {
00399
00400 Index biggest_col_index;
00401 RealScalar biggest_col_sq_norm = m_colSqNorms.tail(cols-k).maxCoeff(&biggest_col_index);
00402 biggest_col_index += k;
00403
00404
00405
00406
00407
00408 biggest_col_sq_norm = m_qr.col(biggest_col_index).tail(rows-k).squaredNorm();
00409
00410
00411 m_colSqNorms.coeffRef(biggest_col_index) = biggest_col_sq_norm;
00412
00413
00414
00415
00416
00417
00418 if(biggest_col_sq_norm < threshold_helper * RealScalar(rows-k))
00419 {
00420 m_nonzero_pivots = k;
00421 m_hCoeffs.tail(size-k).setZero();
00422 m_qr.bottomRightCorner(rows-k,cols-k)
00423 .template triangularView<StrictlyLower>()
00424 .setZero();
00425 break;
00426 }
00427
00428
00429 m_colsTranspositions.coeffRef(k) = biggest_col_index;
00430 if(k != biggest_col_index) {
00431 m_qr.col(k).swap(m_qr.col(biggest_col_index));
00432 std::swap(m_colSqNorms.coeffRef(k), m_colSqNorms.coeffRef(biggest_col_index));
00433 ++number_of_transpositions;
00434 }
00435
00436
00437 RealScalar beta;
00438 m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta);
00439
00440
00441 m_qr.coeffRef(k,k) = beta;
00442
00443
00444 if(internal::abs(beta) > m_maxpivot) m_maxpivot = internal::abs(beta);
00445
00446
00447 m_qr.bottomRightCorner(rows-k, cols-k-1)
00448 .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1));
00449
00450
00451 m_colSqNorms.tail(cols-k-1) -= m_qr.row(k).tail(cols-k-1).cwiseAbs2();
00452 }
00453
00454 m_colsPermutation.setIdentity(cols);
00455 for(Index k = 0; k < m_nonzero_pivots; ++k)
00456 m_colsPermutation.applyTranspositionOnTheRight(k, m_colsTranspositions.coeff(k));
00457
00458 m_det_pq = (number_of_transpositions%2) ? -1 : 1;
00459 m_isInitialized = true;
00460
00461 return *this;
00462 }
00463
00464 namespace internal {
00465
00466 template<typename _MatrixType, typename Rhs>
00467 struct solve_retval<ColPivHouseholderQR<_MatrixType>, Rhs>
00468 : solve_retval_base<ColPivHouseholderQR<_MatrixType>, Rhs>
00469 {
00470 EIGEN_MAKE_SOLVE_HELPERS(ColPivHouseholderQR<_MatrixType>,Rhs)
00471
00472 template<typename Dest> void evalTo(Dest& dst) const
00473 {
00474 eigen_assert(rhs().rows() == dec().rows());
00475
00476 const int cols = dec().cols(),
00477 nonzero_pivots = dec().nonzeroPivots();
00478
00479 if(nonzero_pivots == 0)
00480 {
00481 dst.setZero();
00482 return;
00483 }
00484
00485 typename Rhs::PlainObject c(rhs());
00486
00487
00488 c.applyOnTheLeft(householderSequence(dec().matrixQR(), dec().hCoeffs())
00489 .setLength(dec().nonzeroPivots())
00490 .transpose()
00491 );
00492
00493 dec().matrixQR()
00494 .topLeftCorner(nonzero_pivots, nonzero_pivots)
00495 .template triangularView<Upper>()
00496 .solveInPlace(c.topRows(nonzero_pivots));
00497
00498
00499 typename Rhs::PlainObject d(c);
00500 d.topRows(nonzero_pivots)
00501 = dec().matrixQR()
00502 .topLeftCorner(nonzero_pivots, nonzero_pivots)
00503 .template triangularView<Upper>()
00504 * c.topRows(nonzero_pivots);
00505
00506 for(Index i = 0; i < nonzero_pivots; ++i) dst.row(dec().colsPermutation().indices().coeff(i)) = c.row(i);
00507 for(Index i = nonzero_pivots; i < cols; ++i) dst.row(dec().colsPermutation().indices().coeff(i)).setZero();
00508 }
00509 };
00510
00511 }
00512
00514 template<typename MatrixType>
00515 typename ColPivHouseholderQR<MatrixType>::HouseholderSequenceType ColPivHouseholderQR<MatrixType>
00516 ::householderQ() const
00517 {
00518 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00519 return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate()).setLength(m_nonzero_pivots);
00520 }
00521
00526 template<typename Derived>
00527 const ColPivHouseholderQR<typename MatrixBase<Derived>::PlainObject>
00528 MatrixBase<Derived>::colPivHouseholderQr() const
00529 {
00530 return ColPivHouseholderQR<PlainObject>(eval());
00531 }
00532
00533 }
00534
00535 #endif // EIGEN_COLPIVOTINGHOUSEHOLDERQR_H