![]() |
Eigen
3.3.3
|
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 00016 namespace internal { 00017 template<typename _MatrixType> struct traits<ColPivHouseholderQR<_MatrixType> > 00018 : traits<_MatrixType> 00019 { 00020 enum { Flags = 0 }; 00021 }; 00022 00023 } // end namespace internal 00024 00048 template<typename _MatrixType> class ColPivHouseholderQR 00049 { 00050 public: 00051 00052 typedef _MatrixType MatrixType; 00053 enum { 00054 RowsAtCompileTime = MatrixType::RowsAtCompileTime, 00055 ColsAtCompileTime = MatrixType::ColsAtCompileTime, 00056 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 00057 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime 00058 }; 00059 typedef typename MatrixType::Scalar Scalar; 00060 typedef typename MatrixType::RealScalar RealScalar; 00061 // FIXME should be int 00062 typedef typename MatrixType::StorageIndex StorageIndex; 00063 typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType; 00064 typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationType; 00065 typedef typename internal::plain_row_type<MatrixType, Index>::type IntRowVectorType; 00066 typedef typename internal::plain_row_type<MatrixType>::type RowVectorType; 00067 typedef typename internal::plain_row_type<MatrixType, RealScalar>::type RealRowVectorType; 00068 typedef HouseholderSequence<MatrixType,typename internal::remove_all<typename HCoeffsType::ConjugateReturnType>::type> HouseholderSequenceType; 00069 typedef typename MatrixType::PlainObject PlainObject; 00070 00071 private: 00072 00073 typedef typename PermutationType::StorageIndex PermIndexType; 00074 00075 public: 00076 00083 ColPivHouseholderQR() 00084 : m_qr(), 00085 m_hCoeffs(), 00086 m_colsPermutation(), 00087 m_colsTranspositions(), 00088 m_temp(), 00089 m_colNormsUpdated(), 00090 m_colNormsDirect(), 00091 m_isInitialized(false), 00092 m_usePrescribedThreshold(false) {} 00093 00100 ColPivHouseholderQR(Index rows, Index cols) 00101 : m_qr(rows, cols), 00102 m_hCoeffs((std::min)(rows,cols)), 00103 m_colsPermutation(PermIndexType(cols)), 00104 m_colsTranspositions(cols), 00105 m_temp(cols), 00106 m_colNormsUpdated(cols), 00107 m_colNormsDirect(cols), 00108 m_isInitialized(false), 00109 m_usePrescribedThreshold(false) {} 00110 00123 template<typename InputType> 00124 explicit ColPivHouseholderQR(const EigenBase<InputType>& matrix) 00125 : m_qr(matrix.rows(), matrix.cols()), 00126 m_hCoeffs((std::min)(matrix.rows(),matrix.cols())), 00127 m_colsPermutation(PermIndexType(matrix.cols())), 00128 m_colsTranspositions(matrix.cols()), 00129 m_temp(matrix.cols()), 00130 m_colNormsUpdated(matrix.cols()), 00131 m_colNormsDirect(matrix.cols()), 00132 m_isInitialized(false), 00133 m_usePrescribedThreshold(false) 00134 { 00135 compute(matrix.derived()); 00136 } 00137 00144 template<typename InputType> 00145 explicit ColPivHouseholderQR(EigenBase<InputType>& matrix) 00146 : m_qr(matrix.derived()), 00147 m_hCoeffs((std::min)(matrix.rows(),matrix.cols())), 00148 m_colsPermutation(PermIndexType(matrix.cols())), 00149 m_colsTranspositions(matrix.cols()), 00150 m_temp(matrix.cols()), 00151 m_colNormsUpdated(matrix.cols()), 00152 m_colNormsDirect(matrix.cols()), 00153 m_isInitialized(false), 00154 m_usePrescribedThreshold(false) 00155 { 00156 computeInPlace(); 00157 } 00158 00173 template<typename Rhs> 00174 inline const Solve<ColPivHouseholderQR, Rhs> 00175 solve(const MatrixBase<Rhs>& b) const 00176 { 00177 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); 00178 return Solve<ColPivHouseholderQR, Rhs>(*this, b.derived()); 00179 } 00180 00181 HouseholderSequenceType householderQ() const; 00182 HouseholderSequenceType matrixQ() const 00183 { 00184 return householderQ(); 00185 } 00186 00189 const MatrixType& matrixQR() const 00190 { 00191 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); 00192 return m_qr; 00193 } 00194 00204 const MatrixType& matrixR() const 00205 { 00206 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); 00207 return m_qr; 00208 } 00209 00210 template<typename InputType> 00211 ColPivHouseholderQR& compute(const EigenBase<InputType>& matrix); 00212 00214 const PermutationType& colsPermutation() const 00215 { 00216 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); 00217 return m_colsPermutation; 00218 } 00219 00233 typename MatrixType::RealScalar absDeterminant() const; 00234 00247 typename MatrixType::RealScalar logAbsDeterminant() const; 00248 00255 inline Index rank() const 00256 { 00257 using std::abs; 00258 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); 00259 RealScalar premultiplied_threshold = abs(m_maxpivot) * threshold(); 00260 Index result = 0; 00261 for(Index i = 0; i < m_nonzero_pivots; ++i) 00262 result += (abs(m_qr.coeff(i,i)) > premultiplied_threshold); 00263 return result; 00264 } 00265 00272 inline Index dimensionOfKernel() const 00273 { 00274 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); 00275 return cols() - rank(); 00276 } 00277 00285 inline bool isInjective() const 00286 { 00287 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); 00288 return rank() == cols(); 00289 } 00290 00298 inline bool isSurjective() const 00299 { 00300 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); 00301 return rank() == rows(); 00302 } 00303 00310 inline bool isInvertible() const 00311 { 00312 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); 00313 return isInjective() && isSurjective(); 00314 } 00315 00321 inline const Inverse<ColPivHouseholderQR> inverse() const 00322 { 00323 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); 00324 return Inverse<ColPivHouseholderQR>(*this); 00325 } 00326 00327 inline Index rows() const { return m_qr.rows(); } 00328 inline Index cols() const { return m_qr.cols(); } 00329 00334 const HCoeffsType& hCoeffs() const { return m_hCoeffs; } 00335 00353 ColPivHouseholderQR& setThreshold(const RealScalar& threshold) 00354 { 00355 m_usePrescribedThreshold = true; 00356 m_prescribedThreshold = threshold; 00357 return *this; 00358 } 00359 00368 ColPivHouseholderQR& setThreshold(Default_t) 00369 { 00370 m_usePrescribedThreshold = false; 00371 return *this; 00372 } 00373 00378 RealScalar threshold() const 00379 { 00380 eigen_assert(m_isInitialized || m_usePrescribedThreshold); 00381 return m_usePrescribedThreshold ? m_prescribedThreshold 00382 // this formula comes from experimenting (see "LU precision tuning" thread on the list) 00383 // and turns out to be identical to Higham's formula used already in LDLt. 00384 : NumTraits<Scalar>::epsilon() * RealScalar(m_qr.diagonalSize()); 00385 } 00386 00394 inline Index nonzeroPivots() const 00395 { 00396 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); 00397 return m_nonzero_pivots; 00398 } 00399 00403 RealScalar maxPivot() const { return m_maxpivot; } 00404 00411 ComputationInfo info() const 00412 { 00413 eigen_assert(m_isInitialized && "Decomposition is not initialized."); 00414 return Success; 00415 } 00416 00417 #ifndef EIGEN_PARSED_BY_DOXYGEN 00418 template<typename RhsType, typename DstType> 00419 EIGEN_DEVICE_FUNC 00420 void _solve_impl(const RhsType &rhs, DstType &dst) const; 00421 #endif 00422 00423 protected: 00424 00425 friend class CompleteOrthogonalDecomposition<MatrixType>; 00426 00427 static void check_template_parameters() 00428 { 00429 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar); 00430 } 00431 00432 void computeInPlace(); 00433 00434 MatrixType m_qr; 00435 HCoeffsType m_hCoeffs; 00436 PermutationType m_colsPermutation; 00437 IntRowVectorType m_colsTranspositions; 00438 RowVectorType m_temp; 00439 RealRowVectorType m_colNormsUpdated; 00440 RealRowVectorType m_colNormsDirect; 00441 bool m_isInitialized, m_usePrescribedThreshold; 00442 RealScalar m_prescribedThreshold, m_maxpivot; 00443 Index m_nonzero_pivots; 00444 Index m_det_pq; 00445 }; 00446 00447 template<typename MatrixType> 00448 typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::absDeterminant() const 00449 { 00450 using std::abs; 00451 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); 00452 eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!"); 00453 return abs(m_qr.diagonal().prod()); 00454 } 00455 00456 template<typename MatrixType> 00457 typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::logAbsDeterminant() const 00458 { 00459 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); 00460 eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!"); 00461 return m_qr.diagonal().cwiseAbs().array().log().sum(); 00462 } 00463 00470 template<typename MatrixType> 00471 template<typename InputType> 00472 ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const EigenBase<InputType>& matrix) 00473 { 00474 m_qr = matrix.derived(); 00475 computeInPlace(); 00476 return *this; 00477 } 00478 00479 template<typename MatrixType> 00480 void ColPivHouseholderQR<MatrixType>::computeInPlace() 00481 { 00482 check_template_parameters(); 00483 00484 // the column permutation is stored as int indices, so just to be sure: 00485 eigen_assert(m_qr.cols()<=NumTraits<int>::highest()); 00486 00487 using std::abs; 00488 00489 Index rows = m_qr.rows(); 00490 Index cols = m_qr.cols(); 00491 Index size = m_qr.diagonalSize(); 00492 00493 m_hCoeffs.resize(size); 00494 00495 m_temp.resize(cols); 00496 00497 m_colsTranspositions.resize(m_qr.cols()); 00498 Index number_of_transpositions = 0; 00499 00500 m_colNormsUpdated.resize(cols); 00501 m_colNormsDirect.resize(cols); 00502 for (Index k = 0; k < cols; ++k) { 00503 // colNormsDirect(k) caches the most recent directly computed norm of 00504 // column k. 00505 m_colNormsDirect.coeffRef(k) = m_qr.col(k).norm(); 00506 m_colNormsUpdated.coeffRef(k) = m_colNormsDirect.coeffRef(k); 00507 } 00508 00509 RealScalar threshold_helper = numext::abs2<Scalar>(m_colNormsUpdated.maxCoeff() * NumTraits<Scalar>::epsilon()) / RealScalar(rows); 00510 RealScalar norm_downdate_threshold = numext::sqrt(NumTraits<Scalar>::epsilon()); 00511 00512 m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case) 00513 m_maxpivot = RealScalar(0); 00514 00515 for(Index k = 0; k < size; ++k) 00516 { 00517 // first, we look up in our table m_colNormsUpdated which column has the biggest norm 00518 Index biggest_col_index; 00519 RealScalar biggest_col_sq_norm = numext::abs2(m_colNormsUpdated.tail(cols-k).maxCoeff(&biggest_col_index)); 00520 biggest_col_index += k; 00521 00522 // Track the number of meaningful pivots but do not stop the decomposition to make 00523 // sure that the initial matrix is properly reproduced. See bug 941. 00524 if(m_nonzero_pivots==size && biggest_col_sq_norm < threshold_helper * RealScalar(rows-k)) 00525 m_nonzero_pivots = k; 00526 00527 // apply the transposition to the columns 00528 m_colsTranspositions.coeffRef(k) = biggest_col_index; 00529 if(k != biggest_col_index) { 00530 m_qr.col(k).swap(m_qr.col(biggest_col_index)); 00531 std::swap(m_colNormsUpdated.coeffRef(k), m_colNormsUpdated.coeffRef(biggest_col_index)); 00532 std::swap(m_colNormsDirect.coeffRef(k), m_colNormsDirect.coeffRef(biggest_col_index)); 00533 ++number_of_transpositions; 00534 } 00535 00536 // generate the householder vector, store it below the diagonal 00537 RealScalar beta; 00538 m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta); 00539 00540 // apply the householder transformation to the diagonal coefficient 00541 m_qr.coeffRef(k,k) = beta; 00542 00543 // remember the maximum absolute value of diagonal coefficients 00544 if(abs(beta) > m_maxpivot) m_maxpivot = abs(beta); 00545 00546 // apply the householder transformation 00547 m_qr.bottomRightCorner(rows-k, cols-k-1) 00548 .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1)); 00549 00550 // update our table of norms of the columns 00551 for (Index j = k + 1; j < cols; ++j) { 00552 // The following implements the stable norm downgrade step discussed in 00553 // http://www.netlib.org/lapack/lawnspdf/lawn176.pdf 00554 // and used in LAPACK routines xGEQPF and xGEQP3. 00555 // See lines 278-297 in http://www.netlib.org/lapack/explore-html/dc/df4/sgeqpf_8f_source.html 00556 if (m_colNormsUpdated.coeffRef(j) != 0) { 00557 RealScalar temp = abs(m_qr.coeffRef(k, j)) / m_colNormsUpdated.coeffRef(j); 00558 temp = (RealScalar(1) + temp) * (RealScalar(1) - temp); 00559 temp = temp < 0 ? 0 : temp; 00560 RealScalar temp2 = temp * numext::abs2<Scalar>(m_colNormsUpdated.coeffRef(j) / 00561 m_colNormsDirect.coeffRef(j)); 00562 if (temp2 <= norm_downdate_threshold) { 00563 // The updated norm has become too inaccurate so re-compute the column 00564 // norm directly. 00565 m_colNormsDirect.coeffRef(j) = m_qr.col(j).tail(rows - k - 1).norm(); 00566 m_colNormsUpdated.coeffRef(j) = m_colNormsDirect.coeffRef(j); 00567 } else { 00568 m_colNormsUpdated.coeffRef(j) *= numext::sqrt(temp); 00569 } 00570 } 00571 } 00572 } 00573 00574 m_colsPermutation.setIdentity(PermIndexType(cols)); 00575 for(PermIndexType k = 0; k < size/*m_nonzero_pivots*/; ++k) 00576 m_colsPermutation.applyTranspositionOnTheRight(k, PermIndexType(m_colsTranspositions.coeff(k))); 00577 00578 m_det_pq = (number_of_transpositions%2) ? -1 : 1; 00579 m_isInitialized = true; 00580 } 00581 00582 #ifndef EIGEN_PARSED_BY_DOXYGEN 00583 template<typename _MatrixType> 00584 template<typename RhsType, typename DstType> 00585 void ColPivHouseholderQR<_MatrixType>::_solve_impl(const RhsType &rhs, DstType &dst) const 00586 { 00587 eigen_assert(rhs.rows() == rows()); 00588 00589 const Index nonzero_pivots = nonzeroPivots(); 00590 00591 if(nonzero_pivots == 0) 00592 { 00593 dst.setZero(); 00594 return; 00595 } 00596 00597 typename RhsType::PlainObject c(rhs); 00598 00599 // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T 00600 c.applyOnTheLeft(householderSequence(m_qr, m_hCoeffs) 00601 .setLength(nonzero_pivots) 00602 .transpose() 00603 ); 00604 00605 m_qr.topLeftCorner(nonzero_pivots, nonzero_pivots) 00606 .template triangularView<Upper>() 00607 .solveInPlace(c.topRows(nonzero_pivots)); 00608 00609 for(Index i = 0; i < nonzero_pivots; ++i) dst.row(m_colsPermutation.indices().coeff(i)) = c.row(i); 00610 for(Index i = nonzero_pivots; i < cols(); ++i) dst.row(m_colsPermutation.indices().coeff(i)).setZero(); 00611 } 00612 #endif 00613 00614 namespace internal { 00615 00616 template<typename DstXprType, typename MatrixType> 00617 struct Assignment<DstXprType, Inverse<ColPivHouseholderQR<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename ColPivHouseholderQR<MatrixType>::Scalar>, Dense2Dense> 00618 { 00619 typedef ColPivHouseholderQR<MatrixType> QrType; 00620 typedef Inverse<QrType> SrcXprType; 00621 static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename QrType::Scalar> &) 00622 { 00623 dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.cols())); 00624 } 00625 }; 00626 00627 } // end namespace internal 00628 00632 template<typename MatrixType> 00633 typename ColPivHouseholderQR<MatrixType>::HouseholderSequenceType ColPivHouseholderQR<MatrixType> 00634 ::householderQ() const 00635 { 00636 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); 00637 return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate()); 00638 } 00639 00644 template<typename Derived> 00645 const ColPivHouseholderQR<typename MatrixBase<Derived>::PlainObject> 00646 MatrixBase<Derived>::colPivHouseholderQr() const 00647 { 00648 return ColPivHouseholderQR<PlainObject>(eval()); 00649 } 00650 00651 } // end namespace Eigen 00652 00653 #endif // EIGEN_COLPIVOTINGHOUSEHOLDERQR_H