diff options
Diffstat (limited to 'Eigen/src/QR/ColPivHouseholderQR.h')
| -rw-r--r-- | Eigen/src/QR/ColPivHouseholderQR.h | 110 |
1 files changed, 85 insertions, 25 deletions
diff --git a/Eigen/src/QR/ColPivHouseholderQR.h b/Eigen/src/QR/ColPivHouseholderQR.h index 2daa23cc3..773d1df8f 100644 --- a/Eigen/src/QR/ColPivHouseholderQR.h +++ b/Eigen/src/QR/ColPivHouseholderQR.h @@ -55,7 +55,13 @@ template<typename _MatrixType> class ColPivHouseholderQR typedef typename internal::plain_row_type<MatrixType, Index>::type IntRowVectorType; typedef typename internal::plain_row_type<MatrixType>::type RowVectorType; typedef typename internal::plain_row_type<MatrixType, RealScalar>::type RealRowVectorType; - typedef typename HouseholderSequence<MatrixType,HCoeffsType>::ConjugateReturnType HouseholderSequenceType; + typedef HouseholderSequence<MatrixType,typename internal::remove_all<typename HCoeffsType::ConjugateReturnType>::type> HouseholderSequenceType; + + private: + + typedef typename PermutationType::Index PermIndexType; + + public: /** * \brief Default Constructor. @@ -70,7 +76,8 @@ template<typename _MatrixType> class ColPivHouseholderQR m_colsTranspositions(), m_temp(), m_colSqNorms(), - m_isInitialized(false) {} + m_isInitialized(false), + m_usePrescribedThreshold(false) {} /** \brief Default Constructor with memory preallocation * @@ -81,17 +88,29 @@ template<typename _MatrixType> class ColPivHouseholderQR ColPivHouseholderQR(Index rows, Index cols) : m_qr(rows, cols), m_hCoeffs((std::min)(rows,cols)), - m_colsPermutation(cols), + m_colsPermutation(PermIndexType(cols)), m_colsTranspositions(cols), m_temp(cols), m_colSqNorms(cols), m_isInitialized(false), m_usePrescribedThreshold(false) {} + /** \brief Constructs a QR factorization from a given matrix + * + * This constructor computes the QR factorization of the matrix \a matrix by calling + * the method compute(). It is a short cut for: + * + * \code + * ColPivHouseholderQR<MatrixType> qr(matrix.rows(), matrix.cols()); + * qr.compute(matrix); + * \endcode + * + * \sa compute() + */ ColPivHouseholderQR(const MatrixType& matrix) : m_qr(matrix.rows(), matrix.cols()), m_hCoeffs((std::min)(matrix.rows(),matrix.cols())), - m_colsPermutation(matrix.cols()), + m_colsPermutation(PermIndexType(matrix.cols())), m_colsTranspositions(matrix.cols()), m_temp(matrix.cols()), m_colSqNorms(matrix.cols()), @@ -127,6 +146,10 @@ template<typename _MatrixType> class ColPivHouseholderQR } HouseholderSequenceType householderQ(void) const; + HouseholderSequenceType matrixQ(void) const + { + return householderQ(); + } /** \returns a reference to the matrix where the Householder QR decomposition is stored */ @@ -135,9 +158,25 @@ template<typename _MatrixType> class ColPivHouseholderQR eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); return m_qr; } - + + /** \returns a reference to the matrix where the result Householder QR is stored + * \warning The strict lower part of this matrix contains internal values. + * Only the upper triangular part should be referenced. To get it, use + * \code matrixR().template triangularView<Upper>() \endcode + * For rank-deficient matrices, use + * \code + * matrixR().topLeftCorner(rank(), rank()).template triangularView<Upper>() + * \endcode + */ + const MatrixType& matrixR() const + { + eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); + return m_qr; + } + ColPivHouseholderQR& compute(const MatrixType& matrix); + /** \returns a const reference to the column permutation matrix */ const PermutationType& colsPermutation() const { eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); @@ -181,11 +220,12 @@ template<typename _MatrixType> class ColPivHouseholderQR */ inline Index rank() const { + using std::abs; eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); - RealScalar premultiplied_threshold = internal::abs(m_maxpivot) * threshold(); + RealScalar premultiplied_threshold = abs(m_maxpivot) * threshold(); Index result = 0; for(Index i = 0; i < m_nonzero_pivots; ++i) - result += (internal::abs(m_qr.coeff(i,i)) > premultiplied_threshold); + result += (abs(m_qr.coeff(i,i)) > premultiplied_threshold); return result; } @@ -255,6 +295,11 @@ template<typename _MatrixType> class ColPivHouseholderQR inline Index rows() const { return m_qr.rows(); } inline Index cols() const { return m_qr.cols(); } + + /** \returns a const reference to the vector of Householder coefficients used to represent the factor \c Q. + * + * For advanced uses only. + */ const HCoeffsType& hCoeffs() const { return m_hCoeffs; } /** Allows to prescribe a threshold to be used by certain methods, such as rank(), @@ -305,7 +350,7 @@ template<typename _MatrixType> class ColPivHouseholderQR return m_usePrescribedThreshold ? m_prescribedThreshold // this formula comes from experimenting (see "LU precision tuning" thread on the list) // and turns out to be identical to Higham's formula used already in LDLt. - : NumTraits<Scalar>::epsilon() * m_qr.diagonalSize(); + : NumTraits<Scalar>::epsilon() * RealScalar(m_qr.diagonalSize()); } /** \returns the number of nonzero pivots in the QR decomposition. @@ -325,6 +370,18 @@ template<typename _MatrixType> class ColPivHouseholderQR * diagonal coefficient of R. */ RealScalar maxPivot() const { return m_maxpivot; } + + /** \brief Reports whether the QR factorization was succesful. + * + * \note This function always returns \c Success. It is provided for compatibility + * with other factorization routines. + * \returns \c Success + */ + ComputationInfo info() const + { + eigen_assert(m_isInitialized && "Decomposition is not initialized."); + return Success; + } protected: MatrixType m_qr; @@ -342,9 +399,10 @@ template<typename _MatrixType> class ColPivHouseholderQR template<typename MatrixType> typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::absDeterminant() const { + using std::abs; eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!"); - return internal::abs(m_qr.diagonal().prod()); + return abs(m_qr.diagonal().prod()); } template<typename MatrixType> @@ -355,12 +413,22 @@ typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::logAbsDetermina return m_qr.diagonal().cwiseAbs().array().log().sum(); } +/** Performs the QR factorization of the given matrix \a matrix. The result of + * the factorization is stored into \c *this, and a reference to \c *this + * is returned. + * + * \sa class ColPivHouseholderQR, ColPivHouseholderQR(const MatrixType&) + */ template<typename MatrixType> ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const MatrixType& matrix) { + using std::abs; Index rows = matrix.rows(); Index cols = matrix.cols(); Index size = matrix.diagonalSize(); + + // the column permutation is stored as int indices, so just to be sure: + eigen_assert(cols<=NumTraits<int>::highest()); m_qr = matrix; m_hCoeffs.resize(size); @@ -374,7 +442,7 @@ ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const for(Index k = 0; k < cols; ++k) m_colSqNorms.coeffRef(k) = m_qr.col(k).squaredNorm(); - RealScalar threshold_helper = m_colSqNorms.maxCoeff() * internal::abs2(NumTraits<Scalar>::epsilon()) / RealScalar(rows); + RealScalar threshold_helper = m_colSqNorms.maxCoeff() * numext::abs2(NumTraits<Scalar>::epsilon()) / RealScalar(rows); m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case) m_maxpivot = RealScalar(0); @@ -426,7 +494,7 @@ ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const m_qr.coeffRef(k,k) = beta; // remember the maximum absolute value of diagonal coefficients - if(internal::abs(beta) > m_maxpivot) m_maxpivot = internal::abs(beta); + if(abs(beta) > m_maxpivot) m_maxpivot = abs(beta); // apply the householder transformation m_qr.bottomRightCorner(rows-k, cols-k-1) @@ -436,9 +504,9 @@ ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const m_colSqNorms.tail(cols-k-1) -= m_qr.row(k).tail(cols-k-1).cwiseAbs2(); } - m_colsPermutation.setIdentity(cols); - for(Index k = 0; k < m_nonzero_pivots; ++k) - m_colsPermutation.applyTranspositionOnTheRight(k, m_colsTranspositions.coeff(k)); + m_colsPermutation.setIdentity(PermIndexType(cols)); + for(PermIndexType k = 0; k < m_nonzero_pivots; ++k) + m_colsPermutation.applyTranspositionOnTheRight(k, PermIndexType(m_colsTranspositions.coeff(k))); m_det_pq = (number_of_transpositions%2) ? -1 : 1; m_isInitialized = true; @@ -458,8 +526,8 @@ struct solve_retval<ColPivHouseholderQR<_MatrixType>, Rhs> { eigen_assert(rhs().rows() == dec().rows()); - const int cols = dec().cols(), - nonzero_pivots = dec().nonzeroPivots(); + const Index cols = dec().cols(), + nonzero_pivots = dec().nonzeroPivots(); if(nonzero_pivots == 0) { @@ -475,19 +543,11 @@ struct solve_retval<ColPivHouseholderQR<_MatrixType>, Rhs> .transpose() ); - dec().matrixQR() + dec().matrixR() .topLeftCorner(nonzero_pivots, nonzero_pivots) .template triangularView<Upper>() .solveInPlace(c.topRows(nonzero_pivots)); - - typename Rhs::PlainObject d(c); - d.topRows(nonzero_pivots) - = dec().matrixQR() - .topLeftCorner(nonzero_pivots, nonzero_pivots) - .template triangularView<Upper>() - * c.topRows(nonzero_pivots); - for(Index i = 0; i < nonzero_pivots; ++i) dst.row(dec().colsPermutation().indices().coeff(i)) = c.row(i); for(Index i = nonzero_pivots; i < cols; ++i) dst.row(dec().colsPermutation().indices().coeff(i)).setZero(); } |