ColPivHouseholderQR.h
Go to the documentation of this file.
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 
00037 template<typename _MatrixType> class ColPivHouseholderQR
00038 {
00039   public:
00040 
00041     typedef _MatrixType MatrixType;
00042     enum {
00043       RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00044       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00045       Options = MatrixType::Options,
00046       MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00047       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00048     };
00049     typedef typename MatrixType::Scalar Scalar;
00050     typedef typename MatrixType::RealScalar RealScalar;
00051     typedef typename MatrixType::Index Index;
00052     typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, Options, MaxRowsAtCompileTime, MaxRowsAtCompileTime> MatrixQType;
00053     typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
00054     typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationType;
00055     typedef typename internal::plain_row_type<MatrixType, Index>::type IntRowVectorType;
00056     typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
00057     typedef typename internal::plain_row_type<MatrixType, RealScalar>::type RealRowVectorType;
00058     typedef HouseholderSequence<MatrixType,typename internal::remove_all<typename HCoeffsType::ConjugateReturnType>::type> HouseholderSequenceType;
00059     
00060   private:
00061     
00062     typedef typename PermutationType::Index PermIndexType;
00063     
00064   public:
00065 
00072     ColPivHouseholderQR()
00073       : m_qr(),
00074         m_hCoeffs(),
00075         m_colsPermutation(),
00076         m_colsTranspositions(),
00077         m_temp(),
00078         m_colSqNorms(),
00079         m_isInitialized(false) {}
00080 
00087     ColPivHouseholderQR(Index rows, Index cols)
00088       : m_qr(rows, cols),
00089         m_hCoeffs((std::min)(rows,cols)),
00090         m_colsPermutation(PermIndexType(cols)),
00091         m_colsTranspositions(cols),
00092         m_temp(cols),
00093         m_colSqNorms(cols),
00094         m_isInitialized(false),
00095         m_usePrescribedThreshold(false) {}
00096 
00109     ColPivHouseholderQR(const MatrixType& matrix)
00110       : m_qr(matrix.rows(), matrix.cols()),
00111         m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
00112         m_colsPermutation(PermIndexType(matrix.cols())),
00113         m_colsTranspositions(matrix.cols()),
00114         m_temp(matrix.cols()),
00115         m_colSqNorms(matrix.cols()),
00116         m_isInitialized(false),
00117         m_usePrescribedThreshold(false)
00118     {
00119       compute(matrix);
00120     }
00121 
00139     template<typename Rhs>
00140     inline const internal::solve_retval<ColPivHouseholderQR, Rhs>
00141     solve(const MatrixBase<Rhs>& b) const
00142     {
00143       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00144       return internal::solve_retval<ColPivHouseholderQR, Rhs>(*this, b.derived());
00145     }
00146 
00147     HouseholderSequenceType householderQ(void) const;
00148     HouseholderSequenceType matrixQ(void) const
00149     {
00150       return householderQ(); 
00151     }
00152 
00155     const MatrixType& matrixQR() const
00156     {
00157       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00158       return m_qr;
00159     }
00160     
00170     const MatrixType& matrixR() const
00171     {
00172       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00173       return m_qr;
00174     }
00175     
00176     ColPivHouseholderQR& compute(const MatrixType& matrix);
00177 
00179     const PermutationType& colsPermutation() const
00180     {
00181       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00182       return m_colsPermutation;
00183     }
00184 
00198     typename MatrixType::RealScalar absDeterminant() const;
00199 
00212     typename MatrixType::RealScalar logAbsDeterminant() const;
00213 
00220     inline Index rank() const
00221     {
00222       using std::abs;
00223       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00224       RealScalar premultiplied_threshold = abs(m_maxpivot) * threshold();
00225       Index result = 0;
00226       for(Index i = 0; i < m_nonzero_pivots; ++i)
00227         result += (abs(m_qr.coeff(i,i)) > premultiplied_threshold);
00228       return result;
00229     }
00230 
00237     inline Index dimensionOfKernel() const
00238     {
00239       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00240       return cols() - rank();
00241     }
00242 
00250     inline bool isInjective() const
00251     {
00252       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00253       return rank() == cols();
00254     }
00255 
00263     inline bool isSurjective() const
00264     {
00265       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00266       return rank() == rows();
00267     }
00268 
00275     inline bool isInvertible() const
00276     {
00277       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00278       return isInjective() && isSurjective();
00279     }
00280 
00286     inline const
00287     internal::solve_retval<ColPivHouseholderQR, typename MatrixType::IdentityReturnType>
00288     inverse() const
00289     {
00290       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00291       return internal::solve_retval<ColPivHouseholderQR,typename MatrixType::IdentityReturnType>
00292                (*this, MatrixType::Identity(m_qr.rows(), m_qr.cols()));
00293     }
00294 
00295     inline Index rows() const { return m_qr.rows(); }
00296     inline Index cols() const { return m_qr.cols(); }
00297     
00302     const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
00303 
00321     ColPivHouseholderQR& setThreshold(const RealScalar& threshold)
00322     {
00323       m_usePrescribedThreshold = true;
00324       m_prescribedThreshold = threshold;
00325       return *this;
00326     }
00327 
00336     ColPivHouseholderQR& setThreshold(Default_t)
00337     {
00338       m_usePrescribedThreshold = false;
00339       return *this;
00340     }
00341 
00346     RealScalar threshold() const
00347     {
00348       eigen_assert(m_isInitialized || m_usePrescribedThreshold);
00349       return m_usePrescribedThreshold ? m_prescribedThreshold
00350       // this formula comes from experimenting (see "LU precision tuning" thread on the list)
00351       // and turns out to be identical to Higham's formula used already in LDLt.
00352                                       : NumTraits<Scalar>::epsilon() * m_qr.diagonalSize();
00353     }
00354 
00362     inline Index nonzeroPivots() const
00363     {
00364       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00365       return m_nonzero_pivots;
00366     }
00367 
00371     RealScalar maxPivot() const { return m_maxpivot; }
00372     
00379     ComputationInfo info() const
00380     {
00381       eigen_assert(m_isInitialized && "Decomposition is not initialized.");
00382       return Success;
00383     }
00384 
00385   protected:
00386     MatrixType m_qr;
00387     HCoeffsType m_hCoeffs;
00388     PermutationType m_colsPermutation;
00389     IntRowVectorType m_colsTranspositions;
00390     RowVectorType m_temp;
00391     RealRowVectorType m_colSqNorms;
00392     bool m_isInitialized, m_usePrescribedThreshold;
00393     RealScalar m_prescribedThreshold, m_maxpivot;
00394     Index m_nonzero_pivots;
00395     Index m_det_pq;
00396 };
00397 
00398 template<typename MatrixType>
00399 typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::absDeterminant() const
00400 {
00401   using std::abs;
00402   eigen_assert(m_isInitialized && "ColPivHouseholderQR 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 abs(m_qr.diagonal().prod());
00405 }
00406 
00407 template<typename MatrixType>
00408 typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::logAbsDeterminant() const
00409 {
00410   eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00411   eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
00412   return m_qr.diagonal().cwiseAbs().array().log().sum();
00413 }
00414 
00421 template<typename MatrixType>
00422 ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const MatrixType& matrix)
00423 {
00424   using std::abs;
00425   Index rows = matrix.rows();
00426   Index cols = matrix.cols();
00427   Index size = matrix.diagonalSize();
00428   
00429   // the column permutation is stored as int indices, so just to be sure:
00430   eigen_assert(cols<=NumTraits<int>::highest());
00431 
00432   m_qr = matrix;
00433   m_hCoeffs.resize(size);
00434 
00435   m_temp.resize(cols);
00436 
00437   m_colsTranspositions.resize(matrix.cols());
00438   Index number_of_transpositions = 0;
00439 
00440   m_colSqNorms.resize(cols);
00441   for(Index k = 0; k < cols; ++k)
00442     m_colSqNorms.coeffRef(k) = m_qr.col(k).squaredNorm();
00443 
00444   RealScalar threshold_helper = m_colSqNorms.maxCoeff() * numext::abs2(NumTraits<Scalar>::epsilon()) / RealScalar(rows);
00445 
00446   m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
00447   m_maxpivot = RealScalar(0);
00448 
00449   for(Index k = 0; k < size; ++k)
00450   {
00451     // first, we look up in our table m_colSqNorms which column has the biggest squared norm
00452     Index biggest_col_index;
00453     RealScalar biggest_col_sq_norm = m_colSqNorms.tail(cols-k).maxCoeff(&biggest_col_index);
00454     biggest_col_index += k;
00455 
00456     // since our table m_colSqNorms accumulates imprecision at every step, we must now recompute
00457     // the actual squared norm of the selected column.
00458     // Note that not doing so does result in solve() sometimes returning inf/nan values
00459     // when running the unit test with 1000 repetitions.
00460     biggest_col_sq_norm = m_qr.col(biggest_col_index).tail(rows-k).squaredNorm();
00461 
00462     // we store that back into our table: it can't hurt to correct our table.
00463     m_colSqNorms.coeffRef(biggest_col_index) = biggest_col_sq_norm;
00464 
00465     // if the current biggest column is smaller than epsilon times the initial biggest column,
00466     // terminate to avoid generating nan/inf values.
00467     // Note that here, if we test instead for "biggest == 0", we get a failure every 1000 (or so)
00468     // repetitions of the unit test, with the result of solve() filled with large values of the order
00469     // of 1/(size*epsilon).
00470     if(biggest_col_sq_norm < threshold_helper * RealScalar(rows-k))
00471     {
00472       m_nonzero_pivots = k;
00473       m_hCoeffs.tail(size-k).setZero();
00474       m_qr.bottomRightCorner(rows-k,cols-k)
00475           .template triangularView<StrictlyLower>()
00476           .setZero();
00477       break;
00478     }
00479 
00480     // apply the transposition to the columns
00481     m_colsTranspositions.coeffRef(k) = biggest_col_index;
00482     if(k != biggest_col_index) {
00483       m_qr.col(k).swap(m_qr.col(biggest_col_index));
00484       std::swap(m_colSqNorms.coeffRef(k), m_colSqNorms.coeffRef(biggest_col_index));
00485       ++number_of_transpositions;
00486     }
00487 
00488     // generate the householder vector, store it below the diagonal
00489     RealScalar beta;
00490     m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta);
00491 
00492     // apply the householder transformation to the diagonal coefficient
00493     m_qr.coeffRef(k,k) = beta;
00494 
00495     // remember the maximum absolute value of diagonal coefficients
00496     if(abs(beta) > m_maxpivot) m_maxpivot = abs(beta);
00497 
00498     // apply the householder transformation
00499     m_qr.bottomRightCorner(rows-k, cols-k-1)
00500         .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1));
00501 
00502     // update our table of squared norms of the columns
00503     m_colSqNorms.tail(cols-k-1) -= m_qr.row(k).tail(cols-k-1).cwiseAbs2();
00504   }
00505 
00506   m_colsPermutation.setIdentity(PermIndexType(cols));
00507   for(PermIndexType k = 0; k < m_nonzero_pivots; ++k)
00508     m_colsPermutation.applyTranspositionOnTheRight(k, PermIndexType(m_colsTranspositions.coeff(k)));
00509 
00510   m_det_pq = (number_of_transpositions%2) ? -1 : 1;
00511   m_isInitialized = true;
00512 
00513   return *this;
00514 }
00515 
00516 namespace internal {
00517 
00518 template<typename _MatrixType, typename Rhs>
00519 struct solve_retval<ColPivHouseholderQR<_MatrixType>, Rhs>
00520   : solve_retval_base<ColPivHouseholderQR<_MatrixType>, Rhs>
00521 {
00522   EIGEN_MAKE_SOLVE_HELPERS(ColPivHouseholderQR<_MatrixType>,Rhs)
00523 
00524   template<typename Dest> void evalTo(Dest& dst) const
00525   {
00526     eigen_assert(rhs().rows() == dec().rows());
00527 
00528     const Index cols = dec().cols(),
00529                                 nonzero_pivots = dec().nonzeroPivots();
00530 
00531     if(nonzero_pivots == 0)
00532     {
00533       dst.setZero();
00534       return;
00535     }
00536 
00537     typename Rhs::PlainObject c(rhs());
00538 
00539     // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
00540     c.applyOnTheLeft(householderSequence(dec().matrixQR(), dec().hCoeffs())
00541                      .setLength(dec().nonzeroPivots())
00542                      .transpose()
00543       );
00544 
00545     dec().matrixR()
00546        .topLeftCorner(nonzero_pivots, nonzero_pivots)
00547        .template triangularView<Upper>()
00548        .solveInPlace(c.topRows(nonzero_pivots));
00549 
00550     for(Index i = 0; i < nonzero_pivots; ++i) dst.row(dec().colsPermutation().indices().coeff(i)) = c.row(i);
00551     for(Index i = nonzero_pivots; i < cols; ++i) dst.row(dec().colsPermutation().indices().coeff(i)).setZero();
00552   }
00553 };
00554 
00555 } // end namespace internal
00556 
00558 template<typename MatrixType>
00559 typename ColPivHouseholderQR<MatrixType>::HouseholderSequenceType ColPivHouseholderQR<MatrixType>
00560   ::householderQ() const
00561 {
00562   eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00563   return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate()).setLength(m_nonzero_pivots);
00564 }
00565 
00570 template<typename Derived>
00571 const ColPivHouseholderQR<typename MatrixBase<Derived>::PlainObject>
00572 MatrixBase<Derived>::colPivHouseholderQr() const
00573 {
00574   return ColPivHouseholderQR<PlainObject>(eval());
00575 }
00576 
00577 } // end namespace Eigen
00578 
00579 #endif // EIGEN_COLPIVOTINGHOUSEHOLDERQR_H


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Sat Jun 8 2019 19:36:53