diff options
Diffstat (limited to 'Eigen/src/QR/ColPivHouseholderQR.h')
| -rw-r--r-- | Eigen/src/QR/ColPivHouseholderQR.h | 281 |
1 files changed, 177 insertions, 104 deletions
diff --git a/Eigen/src/QR/ColPivHouseholderQR.h b/Eigen/src/QR/ColPivHouseholderQR.h index 567eab7cd..0e47c8332 100644 --- a/Eigen/src/QR/ColPivHouseholderQR.h +++ b/Eigen/src/QR/ColPivHouseholderQR.h @@ -11,7 +11,16 @@ #ifndef EIGEN_COLPIVOTINGHOUSEHOLDERQR_H #define EIGEN_COLPIVOTINGHOUSEHOLDERQR_H -namespace Eigen { +namespace Eigen { + +namespace internal { +template<typename _MatrixType> struct traits<ColPivHouseholderQR<_MatrixType> > + : traits<_MatrixType> +{ + enum { Flags = 0 }; +}; + +} // end namespace internal /** \ingroup QR_Module * @@ -19,19 +28,21 @@ namespace Eigen { * * \brief Householder rank-revealing QR decomposition of a matrix with column-pivoting * - * \param MatrixType the type of the matrix of which we are computing the QR decomposition + * \tparam _MatrixType the type of the matrix of which we are computing the QR decomposition * * This class performs a rank-revealing QR decomposition of a matrix \b A into matrices \b P, \b Q and \b R - * such that + * such that * \f[ * \mathbf{A} \, \mathbf{P} = \mathbf{Q} \, \mathbf{R} * \f] - * by using Householder transformations. Here, \b P is a permutation matrix, \b Q a unitary matrix and \b R an + * by using Householder transformations. Here, \b P is a permutation matrix, \b Q a unitary matrix and \b R an * upper triangular matrix. * * This decomposition performs column pivoting in order to be rank-revealing and improve * numerical stability. It is slower than HouseholderQR, and faster than FullPivHouseholderQR. * + * This class supports the \link InplaceDecomposition inplace decomposition \endlink mechanism. + * * \sa MatrixBase::colPivHouseholderQr() */ template<typename _MatrixType> class ColPivHouseholderQR @@ -42,25 +53,25 @@ template<typename _MatrixType> class ColPivHouseholderQR enum { RowsAtCompileTime = MatrixType::RowsAtCompileTime, ColsAtCompileTime = MatrixType::ColsAtCompileTime, - Options = MatrixType::Options, MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime }; typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::RealScalar RealScalar; - typedef typename MatrixType::Index Index; - typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, Options, MaxRowsAtCompileTime, MaxRowsAtCompileTime> MatrixQType; + // FIXME should be int + typedef typename MatrixType::StorageIndex StorageIndex; typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType; typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationType; 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 HouseholderSequence<MatrixType,typename internal::remove_all<typename HCoeffsType::ConjugateReturnType>::type> HouseholderSequenceType; - + typedef typename MatrixType::PlainObject PlainObject; + private: - - typedef typename PermutationType::Index PermIndexType; - + + typedef typename PermutationType::StorageIndex PermIndexType; + public: /** @@ -75,7 +86,8 @@ template<typename _MatrixType> class ColPivHouseholderQR m_colsPermutation(), m_colsTranspositions(), m_temp(), - m_colSqNorms(), + m_colNormsUpdated(), + m_colNormsDirect(), m_isInitialized(false), m_usePrescribedThreshold(false) {} @@ -91,7 +103,8 @@ template<typename _MatrixType> class ColPivHouseholderQR m_colsPermutation(PermIndexType(cols)), m_colsTranspositions(cols), m_temp(cols), - m_colSqNorms(cols), + m_colNormsUpdated(cols), + m_colNormsDirect(cols), m_isInitialized(false), m_usePrescribedThreshold(false) {} @@ -99,25 +112,48 @@ template<typename _MatrixType> class ColPivHouseholderQR * * 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) + template<typename InputType> + explicit ColPivHouseholderQR(const EigenBase<InputType>& matrix) : m_qr(matrix.rows(), matrix.cols()), m_hCoeffs((std::min)(matrix.rows(),matrix.cols())), m_colsPermutation(PermIndexType(matrix.cols())), m_colsTranspositions(matrix.cols()), m_temp(matrix.cols()), - m_colSqNorms(matrix.cols()), + m_colNormsUpdated(matrix.cols()), + m_colNormsDirect(matrix.cols()), m_isInitialized(false), m_usePrescribedThreshold(false) { - compute(matrix); + compute(matrix.derived()); + } + + /** \brief Constructs a QR factorization from a given matrix + * + * This overloaded constructor is provided for \link InplaceDecomposition inplace decomposition \endlink when \c MatrixType is a Eigen::Ref. + * + * \sa ColPivHouseholderQR(const EigenBase&) + */ + template<typename InputType> + explicit ColPivHouseholderQR(EigenBase<InputType>& matrix) + : m_qr(matrix.derived()), + m_hCoeffs((std::min)(matrix.rows(),matrix.cols())), + m_colsPermutation(PermIndexType(matrix.cols())), + m_colsTranspositions(matrix.cols()), + m_temp(matrix.cols()), + m_colNormsUpdated(matrix.cols()), + m_colNormsDirect(matrix.cols()), + m_isInitialized(false), + m_usePrescribedThreshold(false) + { + computeInPlace(); } /** This method finds a solution x to the equation Ax=b, where A is the matrix of which @@ -127,9 +163,6 @@ template<typename _MatrixType> class ColPivHouseholderQR * * \returns a solution. * - * \note The case where b is a matrix is not yet implemented. Also, this - * code is space inefficient. - * * \note_about_checking_solutions * * \note_about_arbitrary_choice_of_solution @@ -138,17 +171,17 @@ template<typename _MatrixType> class ColPivHouseholderQR * Output: \verbinclude ColPivHouseholderQR_solve.out */ template<typename Rhs> - inline const internal::solve_retval<ColPivHouseholderQR, Rhs> + inline const Solve<ColPivHouseholderQR, Rhs> solve(const MatrixBase<Rhs>& b) const { eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); - return internal::solve_retval<ColPivHouseholderQR, Rhs>(*this, b.derived()); + return Solve<ColPivHouseholderQR, Rhs>(*this, b.derived()); } - HouseholderSequenceType householderQ(void) const; - HouseholderSequenceType matrixQ(void) const + HouseholderSequenceType householderQ() const; + HouseholderSequenceType matrixQ() const { - return householderQ(); + return householderQ(); } /** \returns a reference to the matrix where the Householder QR decomposition is stored @@ -158,14 +191,14 @@ 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. + + /** \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>() + * For rank-deficient matrices, use + * \code + * matrixR().topLeftCorner(rank(), rank()).template triangularView<Upper>() * \endcode */ const MatrixType& matrixR() const @@ -173,8 +206,9 @@ template<typename _MatrixType> class ColPivHouseholderQR eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); return m_qr; } - - ColPivHouseholderQR& compute(const MatrixType& matrix); + + template<typename InputType> + ColPivHouseholderQR& compute(const EigenBase<InputType>& matrix); /** \returns a const reference to the column permutation matrix */ const PermutationType& colsPermutation() const @@ -284,20 +318,17 @@ template<typename _MatrixType> class ColPivHouseholderQR * \note If this matrix is not invertible, the returned matrix has undefined coefficients. * Use isInvertible() to first determine whether this matrix is invertible. */ - inline const - internal::solve_retval<ColPivHouseholderQR, typename MatrixType::IdentityReturnType> - inverse() const + inline const Inverse<ColPivHouseholderQR> inverse() const { eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); - return internal::solve_retval<ColPivHouseholderQR,typename MatrixType::IdentityReturnType> - (*this, MatrixType::Identity(m_qr.rows(), m_qr.cols())); + return Inverse<ColPivHouseholderQR>(*this); } 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; } @@ -370,12 +401,12 @@ 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 + * \note This function always returns \c Success. It is provided for compatibility * with other factorization routines. - * \returns \c Success + * \returns \c Success */ ComputationInfo info() const { @@ -383,19 +414,30 @@ template<typename _MatrixType> class ColPivHouseholderQR return Success; } + #ifndef EIGEN_PARSED_BY_DOXYGEN + template<typename RhsType, typename DstType> + EIGEN_DEVICE_FUNC + void _solve_impl(const RhsType &rhs, DstType &dst) const; + #endif + protected: - + + friend class CompleteOrthogonalDecomposition<MatrixType>; + static void check_template_parameters() { EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar); } - + + void computeInPlace(); + MatrixType m_qr; HCoeffsType m_hCoeffs; PermutationType m_colsPermutation; IntRowVectorType m_colsTranspositions; RowVectorType m_temp; - RealRowVectorType m_colSqNorms; + RealRowVectorType m_colNormsUpdated; + RealRowVectorType m_colNormsDirect; bool m_isInitialized, m_usePrescribedThreshold; RealScalar m_prescribedThreshold, m_maxpivot; Index m_nonzero_pivots; @@ -426,51 +468,57 @@ typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::logAbsDetermina * \sa class ColPivHouseholderQR, ColPivHouseholderQR(const MatrixType&) */ template<typename MatrixType> -ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const MatrixType& matrix) +template<typename InputType> +ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const EigenBase<InputType>& matrix) +{ + m_qr = matrix.derived(); + computeInPlace(); + return *this; +} + +template<typename MatrixType> +void ColPivHouseholderQR<MatrixType>::computeInPlace() { check_template_parameters(); - - 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()); + eigen_assert(m_qr.cols()<=NumTraits<int>::highest()); + + using std::abs; + + Index rows = m_qr.rows(); + Index cols = m_qr.cols(); + Index size = m_qr.diagonalSize(); - m_qr = matrix; m_hCoeffs.resize(size); m_temp.resize(cols); - m_colsTranspositions.resize(matrix.cols()); + m_colsTranspositions.resize(m_qr.cols()); Index number_of_transpositions = 0; - m_colSqNorms.resize(cols); - for(Index k = 0; k < cols; ++k) - m_colSqNorms.coeffRef(k) = m_qr.col(k).squaredNorm(); + m_colNormsUpdated.resize(cols); + m_colNormsDirect.resize(cols); + for (Index k = 0; k < cols; ++k) { + // colNormsDirect(k) caches the most recent directly computed norm of + // column k. + m_colNormsDirect.coeffRef(k) = m_qr.col(k).norm(); + m_colNormsUpdated.coeffRef(k) = m_colNormsDirect.coeffRef(k); + } - RealScalar threshold_helper = m_colSqNorms.maxCoeff() * numext::abs2(NumTraits<Scalar>::epsilon()) / RealScalar(rows); + RealScalar threshold_helper = numext::abs2<Scalar>(m_colNormsUpdated.maxCoeff() * NumTraits<Scalar>::epsilon()) / RealScalar(rows); + RealScalar norm_downdate_threshold = numext::sqrt(NumTraits<Scalar>::epsilon()); m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case) m_maxpivot = RealScalar(0); for(Index k = 0; k < size; ++k) { - // first, we look up in our table m_colSqNorms which column has the biggest squared norm + // first, we look up in our table m_colNormsUpdated which column has the biggest norm Index biggest_col_index; - RealScalar biggest_col_sq_norm = m_colSqNorms.tail(cols-k).maxCoeff(&biggest_col_index); + RealScalar biggest_col_sq_norm = numext::abs2(m_colNormsUpdated.tail(cols-k).maxCoeff(&biggest_col_index)); biggest_col_index += k; - // since our table m_colSqNorms accumulates imprecision at every step, we must now recompute - // the actual squared norm of the selected column. - // Note that not doing so does result in solve() sometimes returning inf/nan values - // when running the unit test with 1000 repetitions. - biggest_col_sq_norm = m_qr.col(biggest_col_index).tail(rows-k).squaredNorm(); - - // we store that back into our table: it can't hurt to correct our table. - m_colSqNorms.coeffRef(biggest_col_index) = biggest_col_sq_norm; - // Track the number of meaningful pivots but do not stop the decomposition to make // sure that the initial matrix is properly reproduced. See bug 941. if(m_nonzero_pivots==size && biggest_col_sq_norm < threshold_helper * RealScalar(rows-k)) @@ -480,7 +528,8 @@ ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const m_colsTranspositions.coeffRef(k) = biggest_col_index; if(k != biggest_col_index) { m_qr.col(k).swap(m_qr.col(biggest_col_index)); - std::swap(m_colSqNorms.coeffRef(k), m_colSqNorms.coeffRef(biggest_col_index)); + std::swap(m_colNormsUpdated.coeffRef(k), m_colNormsUpdated.coeffRef(biggest_col_index)); + std::swap(m_colNormsDirect.coeffRef(k), m_colNormsDirect.coeffRef(biggest_col_index)); ++number_of_transpositions; } @@ -498,8 +547,28 @@ ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const m_qr.bottomRightCorner(rows-k, cols-k-1) .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1)); - // update our table of squared norms of the columns - m_colSqNorms.tail(cols-k-1) -= m_qr.row(k).tail(cols-k-1).cwiseAbs2(); + // update our table of norms of the columns + for (Index j = k + 1; j < cols; ++j) { + // The following implements the stable norm downgrade step discussed in + // http://www.netlib.org/lapack/lawnspdf/lawn176.pdf + // and used in LAPACK routines xGEQPF and xGEQP3. + // See lines 278-297 in http://www.netlib.org/lapack/explore-html/dc/df4/sgeqpf_8f_source.html + if (m_colNormsUpdated.coeffRef(j) != 0) { + RealScalar temp = abs(m_qr.coeffRef(k, j)) / m_colNormsUpdated.coeffRef(j); + temp = (RealScalar(1) + temp) * (RealScalar(1) - temp); + temp = temp < 0 ? 0 : temp; + RealScalar temp2 = temp * numext::abs2<Scalar>(m_colNormsUpdated.coeffRef(j) / + m_colNormsDirect.coeffRef(j)); + if (temp2 <= norm_downdate_threshold) { + // The updated norm has become too inaccurate so re-compute the column + // norm directly. + m_colNormsDirect.coeffRef(j) = m_qr.col(j).tail(rows - k - 1).norm(); + m_colNormsUpdated.coeffRef(j) = m_colNormsDirect.coeffRef(j); + } else { + m_colNormsUpdated.coeffRef(j) *= numext::sqrt(temp); + } + } + } } m_colsPermutation.setIdentity(PermIndexType(cols)); @@ -508,46 +577,50 @@ ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const m_det_pq = (number_of_transpositions%2) ? -1 : 1; m_isInitialized = true; - - return *this; } -namespace internal { - -template<typename _MatrixType, typename Rhs> -struct solve_retval<ColPivHouseholderQR<_MatrixType>, Rhs> - : solve_retval_base<ColPivHouseholderQR<_MatrixType>, Rhs> +#ifndef EIGEN_PARSED_BY_DOXYGEN +template<typename _MatrixType> +template<typename RhsType, typename DstType> +void ColPivHouseholderQR<_MatrixType>::_solve_impl(const RhsType &rhs, DstType &dst) const { - EIGEN_MAKE_SOLVE_HELPERS(ColPivHouseholderQR<_MatrixType>,Rhs) + eigen_assert(rhs.rows() == rows()); + + const Index nonzero_pivots = nonzeroPivots(); - template<typename Dest> void evalTo(Dest& dst) const + if(nonzero_pivots == 0) { - eigen_assert(rhs().rows() == dec().rows()); + dst.setZero(); + return; + } - const Index cols = dec().cols(), - nonzero_pivots = dec().nonzeroPivots(); + typename RhsType::PlainObject c(rhs); - if(nonzero_pivots == 0) - { - dst.setZero(); - return; - } + // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T + c.applyOnTheLeft(householderSequence(m_qr, m_hCoeffs) + .setLength(nonzero_pivots) + .transpose() + ); - typename Rhs::PlainObject c(rhs()); + m_qr.topLeftCorner(nonzero_pivots, nonzero_pivots) + .template triangularView<Upper>() + .solveInPlace(c.topRows(nonzero_pivots)); - // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T - c.applyOnTheLeft(householderSequence(dec().matrixQR(), dec().hCoeffs()) - .setLength(dec().nonzeroPivots()) - .transpose() - ); + for(Index i = 0; i < nonzero_pivots; ++i) dst.row(m_colsPermutation.indices().coeff(i)) = c.row(i); + for(Index i = nonzero_pivots; i < cols(); ++i) dst.row(m_colsPermutation.indices().coeff(i)).setZero(); +} +#endif - dec().matrixR() - .topLeftCorner(nonzero_pivots, nonzero_pivots) - .template triangularView<Upper>() - .solveInPlace(c.topRows(nonzero_pivots)); +namespace internal { - 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(); +template<typename DstXprType, typename MatrixType> +struct Assignment<DstXprType, Inverse<ColPivHouseholderQR<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename ColPivHouseholderQR<MatrixType>::Scalar>, Dense2Dense> +{ + typedef ColPivHouseholderQR<MatrixType> QrType; + typedef Inverse<QrType> SrcXprType; + static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename QrType::Scalar> &) + { + dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.cols())); } }; |