$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_FULLPIVOTINGHOUSEHOLDERQR_H 00012 #define EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H 00013 00014 namespace Eigen { 00015 00016 namespace internal { 00017 00018 template<typename MatrixType> struct FullPivHouseholderQRMatrixQReturnType; 00019 00020 template<typename MatrixType> 00021 struct traits<FullPivHouseholderQRMatrixQReturnType<MatrixType> > 00022 { 00023 typedef typename MatrixType::PlainObject ReturnType; 00024 }; 00025 00026 } 00027 00049 template<typename _MatrixType> class FullPivHouseholderQR 00050 { 00051 public: 00052 00053 typedef _MatrixType MatrixType; 00054 enum { 00055 RowsAtCompileTime = MatrixType::RowsAtCompileTime, 00056 ColsAtCompileTime = MatrixType::ColsAtCompileTime, 00057 Options = MatrixType::Options, 00058 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 00059 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime 00060 }; 00061 typedef typename MatrixType::Scalar Scalar; 00062 typedef typename MatrixType::RealScalar RealScalar; 00063 typedef typename MatrixType::Index Index; 00064 typedef internal::FullPivHouseholderQRMatrixQReturnType<MatrixType> MatrixQReturnType; 00065 typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType; 00066 typedef Matrix<Index, 1, 00067 EIGEN_SIZE_MIN_PREFER_DYNAMIC(ColsAtCompileTime,RowsAtCompileTime), RowMajor, 1, 00068 EIGEN_SIZE_MIN_PREFER_FIXED(MaxColsAtCompileTime,MaxRowsAtCompileTime)> IntDiagSizeVectorType; 00069 typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationType; 00070 typedef typename internal::plain_row_type<MatrixType>::type RowVectorType; 00071 typedef typename internal::plain_col_type<MatrixType>::type ColVectorType; 00072 00078 FullPivHouseholderQR() 00079 : m_qr(), 00080 m_hCoeffs(), 00081 m_rows_transpositions(), 00082 m_cols_transpositions(), 00083 m_cols_permutation(), 00084 m_temp(), 00085 m_isInitialized(false), 00086 m_usePrescribedThreshold(false) {} 00087 00094 FullPivHouseholderQR(Index rows, Index cols) 00095 : m_qr(rows, cols), 00096 m_hCoeffs((std::min)(rows,cols)), 00097 m_rows_transpositions((std::min)(rows,cols)), 00098 m_cols_transpositions((std::min)(rows,cols)), 00099 m_cols_permutation(cols), 00100 m_temp(cols), 00101 m_isInitialized(false), 00102 m_usePrescribedThreshold(false) {} 00103 00116 FullPivHouseholderQR(const MatrixType& matrix) 00117 : m_qr(matrix.rows(), matrix.cols()), 00118 m_hCoeffs((std::min)(matrix.rows(), matrix.cols())), 00119 m_rows_transpositions((std::min)(matrix.rows(), matrix.cols())), 00120 m_cols_transpositions((std::min)(matrix.rows(), matrix.cols())), 00121 m_cols_permutation(matrix.cols()), 00122 m_temp(matrix.cols()), 00123 m_isInitialized(false), 00124 m_usePrescribedThreshold(false) 00125 { 00126 compute(matrix); 00127 } 00128 00147 template<typename Rhs> 00148 inline const internal::solve_retval<FullPivHouseholderQR, Rhs> 00149 solve(const MatrixBase<Rhs>& b) const 00150 { 00151 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); 00152 return internal::solve_retval<FullPivHouseholderQR, Rhs>(*this, b.derived()); 00153 } 00154 00157 MatrixQReturnType matrixQ(void) const; 00158 00161 const MatrixType& matrixQR() const 00162 { 00163 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); 00164 return m_qr; 00165 } 00166 00167 FullPivHouseholderQR& compute(const MatrixType& matrix); 00168 00170 const PermutationType& colsPermutation() const 00171 { 00172 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); 00173 return m_cols_permutation; 00174 } 00175 00177 const IntDiagSizeVectorType& rowsTranspositions() const 00178 { 00179 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); 00180 return m_rows_transpositions; 00181 } 00182 00196 typename MatrixType::RealScalar absDeterminant() const; 00197 00210 typename MatrixType::RealScalar logAbsDeterminant() const; 00211 00218 inline Index rank() const 00219 { 00220 using std::abs; 00221 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); 00222 RealScalar premultiplied_threshold = abs(m_maxpivot) * threshold(); 00223 Index result = 0; 00224 for(Index i = 0; i < m_nonzero_pivots; ++i) 00225 result += (abs(m_qr.coeff(i,i)) > premultiplied_threshold); 00226 return result; 00227 } 00228 00235 inline Index dimensionOfKernel() const 00236 { 00237 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); 00238 return cols() - rank(); 00239 } 00240 00248 inline bool isInjective() const 00249 { 00250 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); 00251 return rank() == cols(); 00252 } 00253 00261 inline bool isSurjective() const 00262 { 00263 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); 00264 return rank() == rows(); 00265 } 00266 00273 inline bool isInvertible() const 00274 { 00275 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); 00276 return isInjective() && isSurjective(); 00277 } 00278 inline const 00284 internal::solve_retval<FullPivHouseholderQR, typename MatrixType::IdentityReturnType> 00285 inverse() const 00286 { 00287 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); 00288 return internal::solve_retval<FullPivHouseholderQR,typename MatrixType::IdentityReturnType> 00289 (*this, MatrixType::Identity(m_qr.rows(), m_qr.cols())); 00290 } 00291 00292 inline Index rows() const { return m_qr.rows(); } 00293 inline Index cols() const { return m_qr.cols(); } 00294 00299 const HCoeffsType& hCoeffs() const { return m_hCoeffs; } 00300 00318 FullPivHouseholderQR& setThreshold(const RealScalar& threshold) 00319 { 00320 m_usePrescribedThreshold = true; 00321 m_prescribedThreshold = threshold; 00322 return *this; 00323 } 00324 00333 FullPivHouseholderQR& setThreshold(Default_t) 00334 { 00335 m_usePrescribedThreshold = false; 00336 return *this; 00337 } 00338 00343 RealScalar threshold() const 00344 { 00345 eigen_assert(m_isInitialized || m_usePrescribedThreshold); 00346 return m_usePrescribedThreshold ? m_prescribedThreshold 00347 // this formula comes from experimenting (see "LU precision tuning" thread on the list) 00348 // and turns out to be identical to Higham's formula used already in LDLt. 00349 : NumTraits<Scalar>::epsilon() * RealScalar(m_qr.diagonalSize()); 00350 } 00351 00359 inline Index nonzeroPivots() const 00360 { 00361 eigen_assert(m_isInitialized && "LU is not initialized."); 00362 return m_nonzero_pivots; 00363 } 00364 00368 RealScalar maxPivot() const { return m_maxpivot; } 00369 00370 protected: 00371 00372 static void check_template_parameters() 00373 { 00374 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar); 00375 } 00376 00377 MatrixType m_qr; 00378 HCoeffsType m_hCoeffs; 00379 IntDiagSizeVectorType m_rows_transpositions; 00380 IntDiagSizeVectorType m_cols_transpositions; 00381 PermutationType m_cols_permutation; 00382 RowVectorType m_temp; 00383 bool m_isInitialized, m_usePrescribedThreshold; 00384 RealScalar m_prescribedThreshold, m_maxpivot; 00385 Index m_nonzero_pivots; 00386 RealScalar m_precision; 00387 Index m_det_pq; 00388 }; 00389 00390 template<typename MatrixType> 00391 typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType>::absDeterminant() const 00392 { 00393 using std::abs; 00394 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); 00395 eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!"); 00396 return abs(m_qr.diagonal().prod()); 00397 } 00398 00399 template<typename MatrixType> 00400 typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType>::logAbsDeterminant() const 00401 { 00402 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); 00403 eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!"); 00404 return m_qr.diagonal().cwiseAbs().array().log().sum(); 00405 } 00406 00413 template<typename MatrixType> 00414 FullPivHouseholderQR<MatrixType>& FullPivHouseholderQR<MatrixType>::compute(const MatrixType& matrix) 00415 { 00416 check_template_parameters(); 00417 00418 using std::abs; 00419 Index rows = matrix.rows(); 00420 Index cols = matrix.cols(); 00421 Index size = (std::min)(rows,cols); 00422 00423 m_qr = matrix; 00424 m_hCoeffs.resize(size); 00425 00426 m_temp.resize(cols); 00427 00428 m_precision = NumTraits<Scalar>::epsilon() * RealScalar(size); 00429 00430 m_rows_transpositions.resize(size); 00431 m_cols_transpositions.resize(size); 00432 Index number_of_transpositions = 0; 00433 00434 RealScalar biggest(0); 00435 00436 m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case) 00437 m_maxpivot = RealScalar(0); 00438 00439 for (Index k = 0; k < size; ++k) 00440 { 00441 Index row_of_biggest_in_corner, col_of_biggest_in_corner; 00442 RealScalar biggest_in_corner; 00443 00444 biggest_in_corner = m_qr.bottomRightCorner(rows-k, cols-k) 00445 .cwiseAbs() 00446 .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner); 00447 row_of_biggest_in_corner += k; 00448 col_of_biggest_in_corner += k; 00449 if(k==0) biggest = biggest_in_corner; 00450 00451 // if the corner is negligible, then we have less than full rank, and we can finish early 00452 if(internal::isMuchSmallerThan(biggest_in_corner, biggest, m_precision)) 00453 { 00454 m_nonzero_pivots = k; 00455 for(Index i = k; i < size; i++) 00456 { 00457 m_rows_transpositions.coeffRef(i) = i; 00458 m_cols_transpositions.coeffRef(i) = i; 00459 m_hCoeffs.coeffRef(i) = Scalar(0); 00460 } 00461 break; 00462 } 00463 00464 m_rows_transpositions.coeffRef(k) = row_of_biggest_in_corner; 00465 m_cols_transpositions.coeffRef(k) = col_of_biggest_in_corner; 00466 if(k != row_of_biggest_in_corner) { 00467 m_qr.row(k).tail(cols-k).swap(m_qr.row(row_of_biggest_in_corner).tail(cols-k)); 00468 ++number_of_transpositions; 00469 } 00470 if(k != col_of_biggest_in_corner) { 00471 m_qr.col(k).swap(m_qr.col(col_of_biggest_in_corner)); 00472 ++number_of_transpositions; 00473 } 00474 00475 RealScalar beta; 00476 m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta); 00477 m_qr.coeffRef(k,k) = beta; 00478 00479 // remember the maximum absolute value of diagonal coefficients 00480 if(abs(beta) > m_maxpivot) m_maxpivot = abs(beta); 00481 00482 m_qr.bottomRightCorner(rows-k, cols-k-1) 00483 .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1)); 00484 } 00485 00486 m_cols_permutation.setIdentity(cols); 00487 for(Index k = 0; k < size; ++k) 00488 m_cols_permutation.applyTranspositionOnTheRight(k, m_cols_transpositions.coeff(k)); 00489 00490 m_det_pq = (number_of_transpositions%2) ? -1 : 1; 00491 m_isInitialized = true; 00492 00493 return *this; 00494 } 00495 00496 namespace internal { 00497 00498 template<typename _MatrixType, typename Rhs> 00499 struct solve_retval<FullPivHouseholderQR<_MatrixType>, Rhs> 00500 : solve_retval_base<FullPivHouseholderQR<_MatrixType>, Rhs> 00501 { 00502 EIGEN_MAKE_SOLVE_HELPERS(FullPivHouseholderQR<_MatrixType>,Rhs) 00503 00504 template<typename Dest> void evalTo(Dest& dst) const 00505 { 00506 const Index rows = dec().rows(), cols = dec().cols(); 00507 eigen_assert(rhs().rows() == rows); 00508 00509 // FIXME introduce nonzeroPivots() and use it here. and more generally, 00510 // make the same improvements in this dec as in FullPivLU. 00511 if(dec().rank()==0) 00512 { 00513 dst.setZero(); 00514 return; 00515 } 00516 00517 typename Rhs::PlainObject c(rhs()); 00518 00519 Matrix<Scalar,1,Rhs::ColsAtCompileTime> temp(rhs().cols()); 00520 for (Index k = 0; k < dec().rank(); ++k) 00521 { 00522 Index remainingSize = rows-k; 00523 c.row(k).swap(c.row(dec().rowsTranspositions().coeff(k))); 00524 c.bottomRightCorner(remainingSize, rhs().cols()) 00525 .applyHouseholderOnTheLeft(dec().matrixQR().col(k).tail(remainingSize-1), 00526 dec().hCoeffs().coeff(k), &temp.coeffRef(0)); 00527 } 00528 00529 dec().matrixQR() 00530 .topLeftCorner(dec().rank(), dec().rank()) 00531 .template triangularView<Upper>() 00532 .solveInPlace(c.topRows(dec().rank())); 00533 00534 for(Index i = 0; i < dec().rank(); ++i) dst.row(dec().colsPermutation().indices().coeff(i)) = c.row(i); 00535 for(Index i = dec().rank(); i < cols; ++i) dst.row(dec().colsPermutation().indices().coeff(i)).setZero(); 00536 } 00537 }; 00538 00545 template<typename MatrixType> struct FullPivHouseholderQRMatrixQReturnType 00546 : public ReturnByValue<FullPivHouseholderQRMatrixQReturnType<MatrixType> > 00547 { 00548 public: 00549 typedef typename MatrixType::Index Index; 00550 typedef typename FullPivHouseholderQR<MatrixType>::IntDiagSizeVectorType IntDiagSizeVectorType; 00551 typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType; 00552 typedef Matrix<typename MatrixType::Scalar, 1, MatrixType::RowsAtCompileTime, RowMajor, 1, 00553 MatrixType::MaxRowsAtCompileTime> WorkVectorType; 00554 00555 FullPivHouseholderQRMatrixQReturnType(const MatrixType& qr, 00556 const HCoeffsType& hCoeffs, 00557 const IntDiagSizeVectorType& rowsTranspositions) 00558 : m_qr(qr), 00559 m_hCoeffs(hCoeffs), 00560 m_rowsTranspositions(rowsTranspositions) 00561 {} 00562 00563 template <typename ResultType> 00564 void evalTo(ResultType& result) const 00565 { 00566 const Index rows = m_qr.rows(); 00567 WorkVectorType workspace(rows); 00568 evalTo(result, workspace); 00569 } 00570 00571 template <typename ResultType> 00572 void evalTo(ResultType& result, WorkVectorType& workspace) const 00573 { 00574 using numext::conj; 00575 // compute the product H'_0 H'_1 ... H'_n-1, 00576 // where H_k is the k-th Householder transformation I - h_k v_k v_k' 00577 // and v_k is the k-th Householder vector [1,m_qr(k+1,k), m_qr(k+2,k), ...] 00578 const Index rows = m_qr.rows(); 00579 const Index cols = m_qr.cols(); 00580 const Index size = (std::min)(rows, cols); 00581 workspace.resize(rows); 00582 result.setIdentity(rows, rows); 00583 for (Index k = size-1; k >= 0; k--) 00584 { 00585 result.block(k, k, rows-k, rows-k) 00586 .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), conj(m_hCoeffs.coeff(k)), &workspace.coeffRef(k)); 00587 result.row(k).swap(result.row(m_rowsTranspositions.coeff(k))); 00588 } 00589 } 00590 00591 Index rows() const { return m_qr.rows(); } 00592 Index cols() const { return m_qr.rows(); } 00593 00594 protected: 00595 typename MatrixType::Nested m_qr; 00596 typename HCoeffsType::Nested m_hCoeffs; 00597 typename IntDiagSizeVectorType::Nested m_rowsTranspositions; 00598 }; 00599 00600 } // end namespace internal 00601 00602 template<typename MatrixType> 00603 inline typename FullPivHouseholderQR<MatrixType>::MatrixQReturnType FullPivHouseholderQR<MatrixType>::matrixQ() const 00604 { 00605 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); 00606 return MatrixQReturnType(m_qr, m_hCoeffs, m_rows_transpositions); 00607 } 00608 00613 template<typename Derived> 00614 const FullPivHouseholderQR<typename MatrixBase<Derived>::PlainObject> 00615 MatrixBase<Derived>::fullPivHouseholderQr() const 00616 { 00617 return FullPivHouseholderQR<PlainObject>(eval()); 00618 } 00619 00620 } // end namespace Eigen 00621 00622 #endif // EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H