Home | History | Annotate | Download | only in QR
      1 // This file is part of Eigen, a lightweight C++ template library
      2 // for linear algebra.
      3 //
      4 // Copyright (C) 2008-2009 Gael Guennebaud <gael.guennebaud (at) inria.fr>
      5 // Copyright (C) 2009 Benoit Jacob <jacob.benoit.1 (at) gmail.com>
      6 //
      7 // This Source Code Form is subject to the terms of the Mozilla
      8 // Public License v. 2.0. If a copy of the MPL was not distributed
      9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
     10 
     11 #ifndef EIGEN_COLPIVOTINGHOUSEHOLDERQR_H
     12 #define EIGEN_COLPIVOTINGHOUSEHOLDERQR_H
     13 
     14 namespace Eigen {
     15 
     16 namespace internal {
     17 template<typename _MatrixType> struct traits<ColPivHouseholderQR<_MatrixType> >
     18  : traits<_MatrixType>
     19 {
     20   enum { Flags = 0 };
     21 };
     22 
     23 } // end namespace internal
     24 
     25 /** \ingroup QR_Module
     26   *
     27   * \class ColPivHouseholderQR
     28   *
     29   * \brief Householder rank-revealing QR decomposition of a matrix with column-pivoting
     30   *
     31   * \tparam _MatrixType the type of the matrix of which we are computing the QR decomposition
     32   *
     33   * This class performs a rank-revealing QR decomposition of a matrix \b A into matrices \b P, \b Q and \b R
     34   * such that
     35   * \f[
     36   *  \mathbf{A} \, \mathbf{P} = \mathbf{Q} \, \mathbf{R}
     37   * \f]
     38   * by using Householder transformations. Here, \b P is a permutation matrix, \b Q a unitary matrix and \b R an
     39   * upper triangular matrix.
     40   *
     41   * This decomposition performs column pivoting in order to be rank-revealing and improve
     42   * numerical stability. It is slower than HouseholderQR, and faster than FullPivHouseholderQR.
     43   *
     44   * This class supports the \link InplaceDecomposition inplace decomposition \endlink mechanism.
     45   *
     46   * \sa MatrixBase::colPivHouseholderQr()
     47   */
     48 template<typename _MatrixType> class ColPivHouseholderQR
     49 {
     50   public:
     51 
     52     typedef _MatrixType MatrixType;
     53     enum {
     54       RowsAtCompileTime = MatrixType::RowsAtCompileTime,
     55       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
     56       MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
     57       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
     58     };
     59     typedef typename MatrixType::Scalar Scalar;
     60     typedef typename MatrixType::RealScalar RealScalar;
     61     // FIXME should be int
     62     typedef typename MatrixType::StorageIndex StorageIndex;
     63     typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
     64     typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationType;
     65     typedef typename internal::plain_row_type<MatrixType, Index>::type IntRowVectorType;
     66     typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
     67     typedef typename internal::plain_row_type<MatrixType, RealScalar>::type RealRowVectorType;
     68     typedef HouseholderSequence<MatrixType,typename internal::remove_all<typename HCoeffsType::ConjugateReturnType>::type> HouseholderSequenceType;
     69     typedef typename MatrixType::PlainObject PlainObject;
     70 
     71   private:
     72 
     73     typedef typename PermutationType::StorageIndex PermIndexType;
     74 
     75   public:
     76 
     77     /**
     78     * \brief Default Constructor.
     79     *
     80     * The default constructor is useful in cases in which the user intends to
     81     * perform decompositions via ColPivHouseholderQR::compute(const MatrixType&).
     82     */
     83     ColPivHouseholderQR()
     84       : m_qr(),
     85         m_hCoeffs(),
     86         m_colsPermutation(),
     87         m_colsTranspositions(),
     88         m_temp(),
     89         m_colNormsUpdated(),
     90         m_colNormsDirect(),
     91         m_isInitialized(false),
     92         m_usePrescribedThreshold(false) {}
     93 
     94     /** \brief Default Constructor with memory preallocation
     95       *
     96       * Like the default constructor but with preallocation of the internal data
     97       * according to the specified problem \a size.
     98       * \sa ColPivHouseholderQR()
     99       */
    100     ColPivHouseholderQR(Index rows, Index cols)
    101       : m_qr(rows, cols),
    102         m_hCoeffs((std::min)(rows,cols)),
    103         m_colsPermutation(PermIndexType(cols)),
    104         m_colsTranspositions(cols),
    105         m_temp(cols),
    106         m_colNormsUpdated(cols),
    107         m_colNormsDirect(cols),
    108         m_isInitialized(false),
    109         m_usePrescribedThreshold(false) {}
    110 
    111     /** \brief Constructs a QR factorization from a given matrix
    112       *
    113       * This constructor computes the QR factorization of the matrix \a matrix by calling
    114       * the method compute(). It is a short cut for:
    115       *
    116       * \code
    117       * ColPivHouseholderQR<MatrixType> qr(matrix.rows(), matrix.cols());
    118       * qr.compute(matrix);
    119       * \endcode
    120       *
    121       * \sa compute()
    122       */
    123     template<typename InputType>
    124     explicit ColPivHouseholderQR(const EigenBase<InputType>& matrix)
    125       : m_qr(matrix.rows(), matrix.cols()),
    126         m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
    127         m_colsPermutation(PermIndexType(matrix.cols())),
    128         m_colsTranspositions(matrix.cols()),
    129         m_temp(matrix.cols()),
    130         m_colNormsUpdated(matrix.cols()),
    131         m_colNormsDirect(matrix.cols()),
    132         m_isInitialized(false),
    133         m_usePrescribedThreshold(false)
    134     {
    135       compute(matrix.derived());
    136     }
    137 
    138     /** \brief Constructs a QR factorization from a given matrix
    139       *
    140       * This overloaded constructor is provided for \link InplaceDecomposition inplace decomposition \endlink when \c MatrixType is a Eigen::Ref.
    141       *
    142       * \sa ColPivHouseholderQR(const EigenBase&)
    143       */
    144     template<typename InputType>
    145     explicit ColPivHouseholderQR(EigenBase<InputType>& matrix)
    146       : m_qr(matrix.derived()),
    147         m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
    148         m_colsPermutation(PermIndexType(matrix.cols())),
    149         m_colsTranspositions(matrix.cols()),
    150         m_temp(matrix.cols()),
    151         m_colNormsUpdated(matrix.cols()),
    152         m_colNormsDirect(matrix.cols()),
    153         m_isInitialized(false),
    154         m_usePrescribedThreshold(false)
    155     {
    156       computeInPlace();
    157     }
    158 
    159     /** This method finds a solution x to the equation Ax=b, where A is the matrix of which
    160       * *this is the QR decomposition, if any exists.
    161       *
    162       * \param b the right-hand-side of the equation to solve.
    163       *
    164       * \returns a solution.
    165       *
    166       * \note_about_checking_solutions
    167       *
    168       * \note_about_arbitrary_choice_of_solution
    169       *
    170       * Example: \include ColPivHouseholderQR_solve.cpp
    171       * Output: \verbinclude ColPivHouseholderQR_solve.out
    172       */
    173     template<typename Rhs>
    174     inline const Solve<ColPivHouseholderQR, Rhs>
    175     solve(const MatrixBase<Rhs>& b) const
    176     {
    177       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
    178       return Solve<ColPivHouseholderQR, Rhs>(*this, b.derived());
    179     }
    180 
    181     HouseholderSequenceType householderQ() const;
    182     HouseholderSequenceType matrixQ() const
    183     {
    184       return householderQ();
    185     }
    186 
    187     /** \returns a reference to the matrix where the Householder QR decomposition is stored
    188       */
    189     const MatrixType& matrixQR() const
    190     {
    191       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
    192       return m_qr;
    193     }
    194 
    195     /** \returns a reference to the matrix where the result Householder QR is stored
    196      * \warning The strict lower part of this matrix contains internal values.
    197      * Only the upper triangular part should be referenced. To get it, use
    198      * \code matrixR().template triangularView<Upper>() \endcode
    199      * For rank-deficient matrices, use
    200      * \code
    201      * matrixR().topLeftCorner(rank(), rank()).template triangularView<Upper>()
    202      * \endcode
    203      */
    204     const MatrixType& matrixR() const
    205     {
    206       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
    207       return m_qr;
    208     }
    209 
    210     template<typename InputType>
    211     ColPivHouseholderQR& compute(const EigenBase<InputType>& matrix);
    212 
    213     /** \returns a const reference to the column permutation matrix */
    214     const PermutationType& colsPermutation() const
    215     {
    216       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
    217       return m_colsPermutation;
    218     }
    219 
    220     /** \returns the absolute value of the determinant of the matrix of which
    221       * *this is the QR decomposition. It has only linear complexity
    222       * (that is, O(n) where n is the dimension of the square matrix)
    223       * as the QR decomposition has already been computed.
    224       *
    225       * \note This is only for square matrices.
    226       *
    227       * \warning a determinant can be very big or small, so for matrices
    228       * of large enough dimension, there is a risk of overflow/underflow.
    229       * One way to work around that is to use logAbsDeterminant() instead.
    230       *
    231       * \sa logAbsDeterminant(), MatrixBase::determinant()
    232       */
    233     typename MatrixType::RealScalar absDeterminant() const;
    234 
    235     /** \returns the natural log of the absolute value of the determinant of the matrix of which
    236       * *this is the QR decomposition. It has only linear complexity
    237       * (that is, O(n) where n is the dimension of the square matrix)
    238       * as the QR decomposition has already been computed.
    239       *
    240       * \note This is only for square matrices.
    241       *
    242       * \note This method is useful to work around the risk of overflow/underflow that's inherent
    243       * to determinant computation.
    244       *
    245       * \sa absDeterminant(), MatrixBase::determinant()
    246       */
    247     typename MatrixType::RealScalar logAbsDeterminant() const;
    248 
    249     /** \returns the rank of the matrix of which *this is the QR decomposition.
    250       *
    251       * \note This method has to determine which pivots should be considered nonzero.
    252       *       For that, it uses the threshold value that you can control by calling
    253       *       setThreshold(const RealScalar&).
    254       */
    255     inline Index rank() const
    256     {
    257       using std::abs;
    258       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
    259       RealScalar premultiplied_threshold = abs(m_maxpivot) * threshold();
    260       Index result = 0;
    261       for(Index i = 0; i < m_nonzero_pivots; ++i)
    262         result += (abs(m_qr.coeff(i,i)) > premultiplied_threshold);
    263       return result;
    264     }
    265 
    266     /** \returns the dimension of the kernel of the matrix of which *this is the QR decomposition.
    267       *
    268       * \note This method has to determine which pivots should be considered nonzero.
    269       *       For that, it uses the threshold value that you can control by calling
    270       *       setThreshold(const RealScalar&).
    271       */
    272     inline Index dimensionOfKernel() const
    273     {
    274       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
    275       return cols() - rank();
    276     }
    277 
    278     /** \returns true if the matrix of which *this is the QR decomposition represents an injective
    279       *          linear map, i.e. has trivial kernel; false otherwise.
    280       *
    281       * \note This method has to determine which pivots should be considered nonzero.
    282       *       For that, it uses the threshold value that you can control by calling
    283       *       setThreshold(const RealScalar&).
    284       */
    285     inline bool isInjective() const
    286     {
    287       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
    288       return rank() == cols();
    289     }
    290 
    291     /** \returns true if the matrix of which *this is the QR decomposition represents a surjective
    292       *          linear map; false otherwise.
    293       *
    294       * \note This method has to determine which pivots should be considered nonzero.
    295       *       For that, it uses the threshold value that you can control by calling
    296       *       setThreshold(const RealScalar&).
    297       */
    298     inline bool isSurjective() const
    299     {
    300       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
    301       return rank() == rows();
    302     }
    303 
    304     /** \returns true if the matrix of which *this is the QR decomposition is invertible.
    305       *
    306       * \note This method has to determine which pivots should be considered nonzero.
    307       *       For that, it uses the threshold value that you can control by calling
    308       *       setThreshold(const RealScalar&).
    309       */
    310     inline bool isInvertible() const
    311     {
    312       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
    313       return isInjective() && isSurjective();
    314     }
    315 
    316     /** \returns the inverse of the matrix of which *this is the QR decomposition.
    317       *
    318       * \note If this matrix is not invertible, the returned matrix has undefined coefficients.
    319       *       Use isInvertible() to first determine whether this matrix is invertible.
    320       */
    321     inline const Inverse<ColPivHouseholderQR> inverse() const
    322     {
    323       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
    324       return Inverse<ColPivHouseholderQR>(*this);
    325     }
    326 
    327     inline Index rows() const { return m_qr.rows(); }
    328     inline Index cols() const { return m_qr.cols(); }
    329 
    330     /** \returns a const reference to the vector of Householder coefficients used to represent the factor \c Q.
    331       *
    332       * For advanced uses only.
    333       */
    334     const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
    335 
    336     /** Allows to prescribe a threshold to be used by certain methods, such as rank(),
    337       * who need to determine when pivots are to be considered nonzero. This is not used for the
    338       * QR decomposition itself.
    339       *
    340       * When it needs to get the threshold value, Eigen calls threshold(). By default, this
    341       * uses a formula to automatically determine a reasonable threshold.
    342       * Once you have called the present method setThreshold(const RealScalar&),
    343       * your value is used instead.
    344       *
    345       * \param threshold The new value to use as the threshold.
    346       *
    347       * A pivot will be considered nonzero if its absolute value is strictly greater than
    348       *  \f$ \vert pivot \vert \leqslant threshold \times \vert maxpivot \vert \f$
    349       * where maxpivot is the biggest pivot.
    350       *
    351       * If you want to come back to the default behavior, call setThreshold(Default_t)
    352       */
    353     ColPivHouseholderQR& setThreshold(const RealScalar& threshold)
    354     {
    355       m_usePrescribedThreshold = true;
    356       m_prescribedThreshold = threshold;
    357       return *this;
    358     }
    359 
    360     /** Allows to come back to the default behavior, letting Eigen use its default formula for
    361       * determining the threshold.
    362       *
    363       * You should pass the special object Eigen::Default as parameter here.
    364       * \code qr.setThreshold(Eigen::Default); \endcode
    365       *
    366       * See the documentation of setThreshold(const RealScalar&).
    367       */
    368     ColPivHouseholderQR& setThreshold(Default_t)
    369     {
    370       m_usePrescribedThreshold = false;
    371       return *this;
    372     }
    373 
    374     /** Returns the threshold that will be used by certain methods such as rank().
    375       *
    376       * See the documentation of setThreshold(const RealScalar&).
    377       */
    378     RealScalar threshold() const
    379     {
    380       eigen_assert(m_isInitialized || m_usePrescribedThreshold);
    381       return m_usePrescribedThreshold ? m_prescribedThreshold
    382       // this formula comes from experimenting (see "LU precision tuning" thread on the list)
    383       // and turns out to be identical to Higham's formula used already in LDLt.
    384                                       : NumTraits<Scalar>::epsilon() * RealScalar(m_qr.diagonalSize());
    385     }
    386 
    387     /** \returns the number of nonzero pivots in the QR decomposition.
    388       * Here nonzero is meant in the exact sense, not in a fuzzy sense.
    389       * So that notion isn't really intrinsically interesting, but it is
    390       * still useful when implementing algorithms.
    391       *
    392       * \sa rank()
    393       */
    394     inline Index nonzeroPivots() const
    395     {
    396       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
    397       return m_nonzero_pivots;
    398     }
    399 
    400     /** \returns the absolute value of the biggest pivot, i.e. the biggest
    401       *          diagonal coefficient of R.
    402       */
    403     RealScalar maxPivot() const { return m_maxpivot; }
    404 
    405     /** \brief Reports whether the QR factorization was succesful.
    406       *
    407       * \note This function always returns \c Success. It is provided for compatibility
    408       * with other factorization routines.
    409       * \returns \c Success
    410       */
    411     ComputationInfo info() const
    412     {
    413       eigen_assert(m_isInitialized && "Decomposition is not initialized.");
    414       return Success;
    415     }
    416 
    417     #ifndef EIGEN_PARSED_BY_DOXYGEN
    418     template<typename RhsType, typename DstType>
    419     EIGEN_DEVICE_FUNC
    420     void _solve_impl(const RhsType &rhs, DstType &dst) const;
    421     #endif
    422 
    423   protected:
    424 
    425     friend class CompleteOrthogonalDecomposition<MatrixType>;
    426 
    427     static void check_template_parameters()
    428     {
    429       EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
    430     }
    431 
    432     void computeInPlace();
    433 
    434     MatrixType m_qr;
    435     HCoeffsType m_hCoeffs;
    436     PermutationType m_colsPermutation;
    437     IntRowVectorType m_colsTranspositions;
    438     RowVectorType m_temp;
    439     RealRowVectorType m_colNormsUpdated;
    440     RealRowVectorType m_colNormsDirect;
    441     bool m_isInitialized, m_usePrescribedThreshold;
    442     RealScalar m_prescribedThreshold, m_maxpivot;
    443     Index m_nonzero_pivots;
    444     Index m_det_pq;
    445 };
    446 
    447 template<typename MatrixType>
    448 typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::absDeterminant() const
    449 {
    450   using std::abs;
    451   eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
    452   eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
    453   return abs(m_qr.diagonal().prod());
    454 }
    455 
    456 template<typename MatrixType>
    457 typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::logAbsDeterminant() const
    458 {
    459   eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
    460   eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
    461   return m_qr.diagonal().cwiseAbs().array().log().sum();
    462 }
    463 
    464 /** Performs the QR factorization of the given matrix \a matrix. The result of
    465   * the factorization is stored into \c *this, and a reference to \c *this
    466   * is returned.
    467   *
    468   * \sa class ColPivHouseholderQR, ColPivHouseholderQR(const MatrixType&)
    469   */
    470 template<typename MatrixType>
    471 template<typename InputType>
    472 ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const EigenBase<InputType>& matrix)
    473 {
    474   m_qr = matrix.derived();
    475   computeInPlace();
    476   return *this;
    477 }
    478 
    479 template<typename MatrixType>
    480 void ColPivHouseholderQR<MatrixType>::computeInPlace()
    481 {
    482   check_template_parameters();
    483 
    484   // the column permutation is stored as int indices, so just to be sure:
    485   eigen_assert(m_qr.cols()<=NumTraits<int>::highest());
    486 
    487   using std::abs;
    488 
    489   Index rows = m_qr.rows();
    490   Index cols = m_qr.cols();
    491   Index size = m_qr.diagonalSize();
    492 
    493   m_hCoeffs.resize(size);
    494 
    495   m_temp.resize(cols);
    496 
    497   m_colsTranspositions.resize(m_qr.cols());
    498   Index number_of_transpositions = 0;
    499 
    500   m_colNormsUpdated.resize(cols);
    501   m_colNormsDirect.resize(cols);
    502   for (Index k = 0; k < cols; ++k) {
    503     // colNormsDirect(k) caches the most recent directly computed norm of
    504     // column k.
    505     m_colNormsDirect.coeffRef(k) = m_qr.col(k).norm();
    506     m_colNormsUpdated.coeffRef(k) = m_colNormsDirect.coeffRef(k);
    507   }
    508 
    509   RealScalar threshold_helper =  numext::abs2<Scalar>(m_colNormsUpdated.maxCoeff() * NumTraits<Scalar>::epsilon()) / RealScalar(rows);
    510   RealScalar norm_downdate_threshold = numext::sqrt(NumTraits<Scalar>::epsilon());
    511 
    512   m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
    513   m_maxpivot = RealScalar(0);
    514 
    515   for(Index k = 0; k < size; ++k)
    516   {
    517     // first, we look up in our table m_colNormsUpdated which column has the biggest norm
    518     Index biggest_col_index;
    519     RealScalar biggest_col_sq_norm = numext::abs2(m_colNormsUpdated.tail(cols-k).maxCoeff(&biggest_col_index));
    520     biggest_col_index += k;
    521 
    522     // Track the number of meaningful pivots but do not stop the decomposition to make
    523     // sure that the initial matrix is properly reproduced. See bug 941.
    524     if(m_nonzero_pivots==size && biggest_col_sq_norm < threshold_helper * RealScalar(rows-k))
    525       m_nonzero_pivots = k;
    526 
    527     // apply the transposition to the columns
    528     m_colsTranspositions.coeffRef(k) = biggest_col_index;
    529     if(k != biggest_col_index) {
    530       m_qr.col(k).swap(m_qr.col(biggest_col_index));
    531       std::swap(m_colNormsUpdated.coeffRef(k), m_colNormsUpdated.coeffRef(biggest_col_index));
    532       std::swap(m_colNormsDirect.coeffRef(k), m_colNormsDirect.coeffRef(biggest_col_index));
    533       ++number_of_transpositions;
    534     }
    535 
    536     // generate the householder vector, store it below the diagonal
    537     RealScalar beta;
    538     m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta);
    539 
    540     // apply the householder transformation to the diagonal coefficient
    541     m_qr.coeffRef(k,k) = beta;
    542 
    543     // remember the maximum absolute value of diagonal coefficients
    544     if(abs(beta) > m_maxpivot) m_maxpivot = abs(beta);
    545 
    546     // apply the householder transformation
    547     m_qr.bottomRightCorner(rows-k, cols-k-1)
    548         .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1));
    549 
    550     // update our table of norms of the columns
    551     for (Index j = k + 1; j < cols; ++j) {
    552       // The following implements the stable norm downgrade step discussed in
    553       // http://www.netlib.org/lapack/lawnspdf/lawn176.pdf
    554       // and used in LAPACK routines xGEQPF and xGEQP3.
    555       // See lines 278-297 in http://www.netlib.org/lapack/explore-html/dc/df4/sgeqpf_8f_source.html
    556       if (m_colNormsUpdated.coeffRef(j) != 0) {
    557         RealScalar temp = abs(m_qr.coeffRef(k, j)) / m_colNormsUpdated.coeffRef(j);
    558         temp = (RealScalar(1) + temp) * (RealScalar(1) - temp);
    559         temp = temp < 0 ? 0 : temp;
    560         RealScalar temp2 = temp * numext::abs2<Scalar>(m_colNormsUpdated.coeffRef(j) /
    561                                                        m_colNormsDirect.coeffRef(j));
    562         if (temp2 <= norm_downdate_threshold) {
    563           // The updated norm has become too inaccurate so re-compute the column
    564           // norm directly.
    565           m_colNormsDirect.coeffRef(j) = m_qr.col(j).tail(rows - k - 1).norm();
    566           m_colNormsUpdated.coeffRef(j) = m_colNormsDirect.coeffRef(j);
    567         } else {
    568           m_colNormsUpdated.coeffRef(j) *= numext::sqrt(temp);
    569         }
    570       }
    571     }
    572   }
    573 
    574   m_colsPermutation.setIdentity(PermIndexType(cols));
    575   for(PermIndexType k = 0; k < size/*m_nonzero_pivots*/; ++k)
    576     m_colsPermutation.applyTranspositionOnTheRight(k, PermIndexType(m_colsTranspositions.coeff(k)));
    577 
    578   m_det_pq = (number_of_transpositions%2) ? -1 : 1;
    579   m_isInitialized = true;
    580 }
    581 
    582 #ifndef EIGEN_PARSED_BY_DOXYGEN
    583 template<typename _MatrixType>
    584 template<typename RhsType, typename DstType>
    585 void ColPivHouseholderQR<_MatrixType>::_solve_impl(const RhsType &rhs, DstType &dst) const
    586 {
    587   eigen_assert(rhs.rows() == rows());
    588 
    589   const Index nonzero_pivots = nonzeroPivots();
    590 
    591   if(nonzero_pivots == 0)
    592   {
    593     dst.setZero();
    594     return;
    595   }
    596 
    597   typename RhsType::PlainObject c(rhs);
    598 
    599   // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
    600   c.applyOnTheLeft(householderSequence(m_qr, m_hCoeffs)
    601                     .setLength(nonzero_pivots)
    602                     .transpose()
    603     );
    604 
    605   m_qr.topLeftCorner(nonzero_pivots, nonzero_pivots)
    606       .template triangularView<Upper>()
    607       .solveInPlace(c.topRows(nonzero_pivots));
    608 
    609   for(Index i = 0; i < nonzero_pivots; ++i) dst.row(m_colsPermutation.indices().coeff(i)) = c.row(i);
    610   for(Index i = nonzero_pivots; i < cols(); ++i) dst.row(m_colsPermutation.indices().coeff(i)).setZero();
    611 }
    612 #endif
    613 
    614 namespace internal {
    615 
    616 template<typename DstXprType, typename MatrixType>
    617 struct Assignment<DstXprType, Inverse<ColPivHouseholderQR<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename ColPivHouseholderQR<MatrixType>::Scalar>, Dense2Dense>
    618 {
    619   typedef ColPivHouseholderQR<MatrixType> QrType;
    620   typedef Inverse<QrType> SrcXprType;
    621   static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename QrType::Scalar> &)
    622   {
    623     dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.cols()));
    624   }
    625 };
    626 
    627 } // end namespace internal
    628 
    629 /** \returns the matrix Q as a sequence of householder transformations.
    630   * You can extract the meaningful part only by using:
    631   * \code qr.householderQ().setLength(qr.nonzeroPivots()) \endcode*/
    632 template<typename MatrixType>
    633 typename ColPivHouseholderQR<MatrixType>::HouseholderSequenceType ColPivHouseholderQR<MatrixType>
    634   ::householderQ() const
    635 {
    636   eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
    637   return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate());
    638 }
    639 
    640 /** \return the column-pivoting Householder QR decomposition of \c *this.
    641   *
    642   * \sa class ColPivHouseholderQR
    643   */
    644 template<typename Derived>
    645 const ColPivHouseholderQR<typename MatrixBase<Derived>::PlainObject>
    646 MatrixBase<Derived>::colPivHouseholderQr() const
    647 {
    648   return ColPivHouseholderQR<PlainObject>(eval());
    649 }
    650 
    651 } // end namespace Eigen
    652 
    653 #endif // EIGEN_COLPIVOTINGHOUSEHOLDERQR_H
    654