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