![]() |
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_FULLPIVOTINGHOUSEHOLDERQR_H 00012 #define EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H 00013 00014 namespace Eigen { 00015 00016 namespace internal { 00017 00018 template<typename _MatrixType> struct traits<FullPivHouseholderQR<_MatrixType> > 00019 : traits<_MatrixType> 00020 { 00021 enum { Flags = 0 }; 00022 }; 00023 00024 template<typename MatrixType> struct FullPivHouseholderQRMatrixQReturnType; 00025 00026 template<typename MatrixType> 00027 struct traits<FullPivHouseholderQRMatrixQReturnType<MatrixType> > 00028 { 00029 typedef typename MatrixType::PlainObject ReturnType; 00030 }; 00031 00032 } // end namespace internal 00033 00057 template<typename _MatrixType> class FullPivHouseholderQR 00058 { 00059 public: 00060 00061 typedef _MatrixType MatrixType; 00062 enum { 00063 RowsAtCompileTime = MatrixType::RowsAtCompileTime, 00064 ColsAtCompileTime = MatrixType::ColsAtCompileTime, 00065 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 00066 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime 00067 }; 00068 typedef typename MatrixType::Scalar Scalar; 00069 typedef typename MatrixType::RealScalar RealScalar; 00070 // FIXME should be int 00071 typedef typename MatrixType::StorageIndex StorageIndex; 00072 typedef internal::FullPivHouseholderQRMatrixQReturnType<MatrixType> MatrixQReturnType; 00073 typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType; 00074 typedef Matrix<StorageIndex, 1, 00075 EIGEN_SIZE_MIN_PREFER_DYNAMIC(ColsAtCompileTime,RowsAtCompileTime), RowMajor, 1, 00076 EIGEN_SIZE_MIN_PREFER_FIXED(MaxColsAtCompileTime,MaxRowsAtCompileTime)> IntDiagSizeVectorType; 00077 typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationType; 00078 typedef typename internal::plain_row_type<MatrixType>::type RowVectorType; 00079 typedef typename internal::plain_col_type<MatrixType>::type ColVectorType; 00080 typedef typename MatrixType::PlainObject PlainObject; 00081 00087 FullPivHouseholderQR() 00088 : m_qr(), 00089 m_hCoeffs(), 00090 m_rows_transpositions(), 00091 m_cols_transpositions(), 00092 m_cols_permutation(), 00093 m_temp(), 00094 m_isInitialized(false), 00095 m_usePrescribedThreshold(false) {} 00096 00103 FullPivHouseholderQR(Index rows, Index cols) 00104 : m_qr(rows, cols), 00105 m_hCoeffs((std::min)(rows,cols)), 00106 m_rows_transpositions((std::min)(rows,cols)), 00107 m_cols_transpositions((std::min)(rows,cols)), 00108 m_cols_permutation(cols), 00109 m_temp(cols), 00110 m_isInitialized(false), 00111 m_usePrescribedThreshold(false) {} 00112 00125 template<typename InputType> 00126 explicit FullPivHouseholderQR(const EigenBase<InputType>& matrix) 00127 : m_qr(matrix.rows(), matrix.cols()), 00128 m_hCoeffs((std::min)(matrix.rows(), matrix.cols())), 00129 m_rows_transpositions((std::min)(matrix.rows(), matrix.cols())), 00130 m_cols_transpositions((std::min)(matrix.rows(), matrix.cols())), 00131 m_cols_permutation(matrix.cols()), 00132 m_temp(matrix.cols()), 00133 m_isInitialized(false), 00134 m_usePrescribedThreshold(false) 00135 { 00136 compute(matrix.derived()); 00137 } 00138 00145 template<typename InputType> 00146 explicit FullPivHouseholderQR(EigenBase<InputType>& matrix) 00147 : m_qr(matrix.derived()), 00148 m_hCoeffs((std::min)(matrix.rows(), matrix.cols())), 00149 m_rows_transpositions((std::min)(matrix.rows(), matrix.cols())), 00150 m_cols_transpositions((std::min)(matrix.rows(), matrix.cols())), 00151 m_cols_permutation(matrix.cols()), 00152 m_temp(matrix.cols()), 00153 m_isInitialized(false), 00154 m_usePrescribedThreshold(false) 00155 { 00156 computeInPlace(); 00157 } 00158 00174 template<typename Rhs> 00175 inline const Solve<FullPivHouseholderQR, Rhs> 00176 solve(const MatrixBase<Rhs>& b) const 00177 { 00178 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); 00179 return Solve<FullPivHouseholderQR, Rhs>(*this, b.derived()); 00180 } 00181 00184 MatrixQReturnType matrixQ(void) const; 00185 00188 const MatrixType& matrixQR() const 00189 { 00190 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); 00191 return m_qr; 00192 } 00193 00194 template<typename InputType> 00195 FullPivHouseholderQR& compute(const EigenBase<InputType>& matrix); 00196 00198 const PermutationType& colsPermutation() const 00199 { 00200 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); 00201 return m_cols_permutation; 00202 } 00203 00205 const IntDiagSizeVectorType& rowsTranspositions() const 00206 { 00207 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); 00208 return m_rows_transpositions; 00209 } 00210 00224 typename MatrixType::RealScalar absDeterminant() const; 00225 00238 typename MatrixType::RealScalar logAbsDeterminant() const; 00239 00246 inline Index rank() const 00247 { 00248 using std::abs; 00249 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); 00250 RealScalar premultiplied_threshold = abs(m_maxpivot) * threshold(); 00251 Index result = 0; 00252 for(Index i = 0; i < m_nonzero_pivots; ++i) 00253 result += (abs(m_qr.coeff(i,i)) > premultiplied_threshold); 00254 return result; 00255 } 00256 00263 inline Index dimensionOfKernel() const 00264 { 00265 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); 00266 return cols() - rank(); 00267 } 00268 00276 inline bool isInjective() const 00277 { 00278 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); 00279 return rank() == cols(); 00280 } 00281 00289 inline bool isSurjective() const 00290 { 00291 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); 00292 return rank() == rows(); 00293 } 00294 00301 inline bool isInvertible() const 00302 { 00303 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); 00304 return isInjective() && isSurjective(); 00305 } 00306 00312 inline const Inverse<FullPivHouseholderQR> inverse() const 00313 { 00314 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); 00315 return Inverse<FullPivHouseholderQR>(*this); 00316 } 00317 00318 inline Index rows() const { return m_qr.rows(); } 00319 inline Index cols() const { return m_qr.cols(); } 00320 00325 const HCoeffsType& hCoeffs() const { return m_hCoeffs; } 00326 00344 FullPivHouseholderQR& setThreshold(const RealScalar& threshold) 00345 { 00346 m_usePrescribedThreshold = true; 00347 m_prescribedThreshold = threshold; 00348 return *this; 00349 } 00350 00359 FullPivHouseholderQR& setThreshold(Default_t) 00360 { 00361 m_usePrescribedThreshold = false; 00362 return *this; 00363 } 00364 00369 RealScalar threshold() const 00370 { 00371 eigen_assert(m_isInitialized || m_usePrescribedThreshold); 00372 return m_usePrescribedThreshold ? m_prescribedThreshold 00373 // this formula comes from experimenting (see "LU precision tuning" thread on the list) 00374 // and turns out to be identical to Higham's formula used already in LDLt. 00375 : NumTraits<Scalar>::epsilon() * RealScalar(m_qr.diagonalSize()); 00376 } 00377 00385 inline Index nonzeroPivots() const 00386 { 00387 eigen_assert(m_isInitialized && "LU is not initialized."); 00388 return m_nonzero_pivots; 00389 } 00390 00394 RealScalar maxPivot() const { return m_maxpivot; } 00395 00396 #ifndef EIGEN_PARSED_BY_DOXYGEN 00397 template<typename RhsType, typename DstType> 00398 EIGEN_DEVICE_FUNC 00399 void _solve_impl(const RhsType &rhs, DstType &dst) const; 00400 #endif 00401 00402 protected: 00403 00404 static void check_template_parameters() 00405 { 00406 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar); 00407 } 00408 00409 void computeInPlace(); 00410 00411 MatrixType m_qr; 00412 HCoeffsType m_hCoeffs; 00413 IntDiagSizeVectorType m_rows_transpositions; 00414 IntDiagSizeVectorType m_cols_transpositions; 00415 PermutationType m_cols_permutation; 00416 RowVectorType m_temp; 00417 bool m_isInitialized, m_usePrescribedThreshold; 00418 RealScalar m_prescribedThreshold, m_maxpivot; 00419 Index m_nonzero_pivots; 00420 RealScalar m_precision; 00421 Index m_det_pq; 00422 }; 00423 00424 template<typename MatrixType> 00425 typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType>::absDeterminant() const 00426 { 00427 using std::abs; 00428 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); 00429 eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!"); 00430 return abs(m_qr.diagonal().prod()); 00431 } 00432 00433 template<typename MatrixType> 00434 typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType>::logAbsDeterminant() const 00435 { 00436 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); 00437 eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!"); 00438 return m_qr.diagonal().cwiseAbs().array().log().sum(); 00439 } 00440 00447 template<typename MatrixType> 00448 template<typename InputType> 00449 FullPivHouseholderQR<MatrixType>& FullPivHouseholderQR<MatrixType>::compute(const EigenBase<InputType>& matrix) 00450 { 00451 m_qr = matrix.derived(); 00452 computeInPlace(); 00453 return *this; 00454 } 00455 00456 template<typename MatrixType> 00457 void FullPivHouseholderQR<MatrixType>::computeInPlace() 00458 { 00459 check_template_parameters(); 00460 00461 using std::abs; 00462 Index rows = m_qr.rows(); 00463 Index cols = m_qr.cols(); 00464 Index size = (std::min)(rows,cols); 00465 00466 00467 m_hCoeffs.resize(size); 00468 00469 m_temp.resize(cols); 00470 00471 m_precision = NumTraits<Scalar>::epsilon() * RealScalar(size); 00472 00473 m_rows_transpositions.resize(size); 00474 m_cols_transpositions.resize(size); 00475 Index number_of_transpositions = 0; 00476 00477 RealScalar biggest(0); 00478 00479 m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case) 00480 m_maxpivot = RealScalar(0); 00481 00482 for (Index k = 0; k < size; ++k) 00483 { 00484 Index row_of_biggest_in_corner, col_of_biggest_in_corner; 00485 typedef internal::scalar_score_coeff_op<Scalar> Scoring; 00486 typedef typename Scoring::result_type Score; 00487 00488 Score score = m_qr.bottomRightCorner(rows-k, cols-k) 00489 .unaryExpr(Scoring()) 00490 .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner); 00491 row_of_biggest_in_corner += k; 00492 col_of_biggest_in_corner += k; 00493 RealScalar biggest_in_corner = internal::abs_knowing_score<Scalar>()(m_qr(row_of_biggest_in_corner, col_of_biggest_in_corner), score); 00494 if(k==0) biggest = biggest_in_corner; 00495 00496 // if the corner is negligible, then we have less than full rank, and we can finish early 00497 if(internal::isMuchSmallerThan(biggest_in_corner, biggest, m_precision)) 00498 { 00499 m_nonzero_pivots = k; 00500 for(Index i = k; i < size; i++) 00501 { 00502 m_rows_transpositions.coeffRef(i) = i; 00503 m_cols_transpositions.coeffRef(i) = i; 00504 m_hCoeffs.coeffRef(i) = Scalar(0); 00505 } 00506 break; 00507 } 00508 00509 m_rows_transpositions.coeffRef(k) = row_of_biggest_in_corner; 00510 m_cols_transpositions.coeffRef(k) = col_of_biggest_in_corner; 00511 if(k != row_of_biggest_in_corner) { 00512 m_qr.row(k).tail(cols-k).swap(m_qr.row(row_of_biggest_in_corner).tail(cols-k)); 00513 ++number_of_transpositions; 00514 } 00515 if(k != col_of_biggest_in_corner) { 00516 m_qr.col(k).swap(m_qr.col(col_of_biggest_in_corner)); 00517 ++number_of_transpositions; 00518 } 00519 00520 RealScalar beta; 00521 m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta); 00522 m_qr.coeffRef(k,k) = beta; 00523 00524 // remember the maximum absolute value of diagonal coefficients 00525 if(abs(beta) > m_maxpivot) m_maxpivot = abs(beta); 00526 00527 m_qr.bottomRightCorner(rows-k, cols-k-1) 00528 .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1)); 00529 } 00530 00531 m_cols_permutation.setIdentity(cols); 00532 for(Index k = 0; k < size; ++k) 00533 m_cols_permutation.applyTranspositionOnTheRight(k, m_cols_transpositions.coeff(k)); 00534 00535 m_det_pq = (number_of_transpositions%2) ? -1 : 1; 00536 m_isInitialized = true; 00537 } 00538 00539 #ifndef EIGEN_PARSED_BY_DOXYGEN 00540 template<typename _MatrixType> 00541 template<typename RhsType, typename DstType> 00542 void FullPivHouseholderQR<_MatrixType>::_solve_impl(const RhsType &rhs, DstType &dst) const 00543 { 00544 eigen_assert(rhs.rows() == rows()); 00545 const Index l_rank = rank(); 00546 00547 // FIXME introduce nonzeroPivots() and use it here. and more generally, 00548 // make the same improvements in this dec as in FullPivLU. 00549 if(l_rank==0) 00550 { 00551 dst.setZero(); 00552 return; 00553 } 00554 00555 typename RhsType::PlainObject c(rhs); 00556 00557 Matrix<Scalar,1,RhsType::ColsAtCompileTime> temp(rhs.cols()); 00558 for (Index k = 0; k < l_rank; ++k) 00559 { 00560 Index remainingSize = rows()-k; 00561 c.row(k).swap(c.row(m_rows_transpositions.coeff(k))); 00562 c.bottomRightCorner(remainingSize, rhs.cols()) 00563 .applyHouseholderOnTheLeft(m_qr.col(k).tail(remainingSize-1), 00564 m_hCoeffs.coeff(k), &temp.coeffRef(0)); 00565 } 00566 00567 m_qr.topLeftCorner(l_rank, l_rank) 00568 .template triangularView<Upper>() 00569 .solveInPlace(c.topRows(l_rank)); 00570 00571 for(Index i = 0; i < l_rank; ++i) dst.row(m_cols_permutation.indices().coeff(i)) = c.row(i); 00572 for(Index i = l_rank; i < cols(); ++i) dst.row(m_cols_permutation.indices().coeff(i)).setZero(); 00573 } 00574 #endif 00575 00576 namespace internal { 00577 00578 template<typename DstXprType, typename MatrixType> 00579 struct Assignment<DstXprType, Inverse<FullPivHouseholderQR<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename FullPivHouseholderQR<MatrixType>::Scalar>, Dense2Dense> 00580 { 00581 typedef FullPivHouseholderQR<MatrixType> QrType; 00582 typedef Inverse<QrType> SrcXprType; 00583 static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename QrType::Scalar> &) 00584 { 00585 dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.cols())); 00586 } 00587 }; 00588 00595 template<typename MatrixType> struct FullPivHouseholderQRMatrixQReturnType 00596 : public ReturnByValue<FullPivHouseholderQRMatrixQReturnType<MatrixType> > 00597 { 00598 public: 00599 typedef typename FullPivHouseholderQR<MatrixType>::IntDiagSizeVectorType IntDiagSizeVectorType; 00600 typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType; 00601 typedef Matrix<typename MatrixType::Scalar, 1, MatrixType::RowsAtCompileTime, RowMajor, 1, 00602 MatrixType::MaxRowsAtCompileTime> WorkVectorType; 00603 00604 FullPivHouseholderQRMatrixQReturnType(const MatrixType& qr, 00605 const HCoeffsType& hCoeffs, 00606 const IntDiagSizeVectorType& rowsTranspositions) 00607 : m_qr(qr), 00608 m_hCoeffs(hCoeffs), 00609 m_rowsTranspositions(rowsTranspositions) 00610 {} 00611 00612 template <typename ResultType> 00613 void evalTo(ResultType& result) const 00614 { 00615 const Index rows = m_qr.rows(); 00616 WorkVectorType workspace(rows); 00617 evalTo(result, workspace); 00618 } 00619 00620 template <typename ResultType> 00621 void evalTo(ResultType& result, WorkVectorType& workspace) const 00622 { 00623 using numext::conj; 00624 // compute the product H'_0 H'_1 ... H'_n-1, 00625 // where H_k is the k-th Householder transformation I - h_k v_k v_k' 00626 // and v_k is the k-th Householder vector [1,m_qr(k+1,k), m_qr(k+2,k), ...] 00627 const Index rows = m_qr.rows(); 00628 const Index cols = m_qr.cols(); 00629 const Index size = (std::min)(rows, cols); 00630 workspace.resize(rows); 00631 result.setIdentity(rows, rows); 00632 for (Index k = size-1; k >= 0; k--) 00633 { 00634 result.block(k, k, rows-k, rows-k) 00635 .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), conj(m_hCoeffs.coeff(k)), &workspace.coeffRef(k)); 00636 result.row(k).swap(result.row(m_rowsTranspositions.coeff(k))); 00637 } 00638 } 00639 00640 Index rows() const { return m_qr.rows(); } 00641 Index cols() const { return m_qr.rows(); } 00642 00643 protected: 00644 typename MatrixType::Nested m_qr; 00645 typename HCoeffsType::Nested m_hCoeffs; 00646 typename IntDiagSizeVectorType::Nested m_rowsTranspositions; 00647 }; 00648 00649 // template<typename MatrixType> 00650 // struct evaluator<FullPivHouseholderQRMatrixQReturnType<MatrixType> > 00651 // : public evaluator<ReturnByValue<FullPivHouseholderQRMatrixQReturnType<MatrixType> > > 00652 // {}; 00653 00654 } // end namespace internal 00655 00656 template<typename MatrixType> 00657 inline typename FullPivHouseholderQR<MatrixType>::MatrixQReturnType FullPivHouseholderQR<MatrixType>::matrixQ() const 00658 { 00659 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); 00660 return MatrixQReturnType(m_qr, m_hCoeffs, m_rows_transpositions); 00661 } 00662 00667 template<typename Derived> 00668 const FullPivHouseholderQR<typename MatrixBase<Derived>::PlainObject> 00669 MatrixBase<Derived>::fullPivHouseholderQr() const 00670 { 00671 return FullPivHouseholderQR<PlainObject>(eval()); 00672 } 00673 00674 } // end namespace Eigen 00675 00676 #endif // EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H