Eigen  3.3.3
ColPivHouseholderQR.h
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 
00016 namespace internal {
00017 template<typename _MatrixType> struct traits<ColPivHouseholderQR<_MatrixType> >
00018  : traits<_MatrixType>
00019 {
00020   enum { Flags = 0 };
00021 };
00022 
00023 } // end namespace internal
00024 
00048 template<typename _MatrixType> class ColPivHouseholderQR
00049 {
00050   public:
00051 
00052     typedef _MatrixType MatrixType;
00053     enum {
00054       RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00055       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00056       MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00057       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00058     };
00059     typedef typename MatrixType::Scalar Scalar;
00060     typedef typename MatrixType::RealScalar RealScalar;
00061     // FIXME should be int
00062     typedef typename MatrixType::StorageIndex StorageIndex;
00063     typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
00064     typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationType;
00065     typedef typename internal::plain_row_type<MatrixType, Index>::type IntRowVectorType;
00066     typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
00067     typedef typename internal::plain_row_type<MatrixType, RealScalar>::type RealRowVectorType;
00068     typedef HouseholderSequence<MatrixType,typename internal::remove_all<typename HCoeffsType::ConjugateReturnType>::type> HouseholderSequenceType;
00069     typedef typename MatrixType::PlainObject PlainObject;
00070 
00071   private:
00072 
00073     typedef typename PermutationType::StorageIndex PermIndexType;
00074 
00075   public:
00076 
00083     ColPivHouseholderQR()
00084       : m_qr(),
00085         m_hCoeffs(),
00086         m_colsPermutation(),
00087         m_colsTranspositions(),
00088         m_temp(),
00089         m_colNormsUpdated(),
00090         m_colNormsDirect(),
00091         m_isInitialized(false),
00092         m_usePrescribedThreshold(false) {}
00093 
00100     ColPivHouseholderQR(Index rows, Index cols)
00101       : m_qr(rows, cols),
00102         m_hCoeffs((std::min)(rows,cols)),
00103         m_colsPermutation(PermIndexType(cols)),
00104         m_colsTranspositions(cols),
00105         m_temp(cols),
00106         m_colNormsUpdated(cols),
00107         m_colNormsDirect(cols),
00108         m_isInitialized(false),
00109         m_usePrescribedThreshold(false) {}
00110 
00123     template<typename InputType>
00124     explicit ColPivHouseholderQR(const EigenBase<InputType>& matrix)
00125       : m_qr(matrix.rows(), matrix.cols()),
00126         m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
00127         m_colsPermutation(PermIndexType(matrix.cols())),
00128         m_colsTranspositions(matrix.cols()),
00129         m_temp(matrix.cols()),
00130         m_colNormsUpdated(matrix.cols()),
00131         m_colNormsDirect(matrix.cols()),
00132         m_isInitialized(false),
00133         m_usePrescribedThreshold(false)
00134     {
00135       compute(matrix.derived());
00136     }
00137 
00144     template<typename InputType>
00145     explicit ColPivHouseholderQR(EigenBase<InputType>& matrix)
00146       : m_qr(matrix.derived()),
00147         m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
00148         m_colsPermutation(PermIndexType(matrix.cols())),
00149         m_colsTranspositions(matrix.cols()),
00150         m_temp(matrix.cols()),
00151         m_colNormsUpdated(matrix.cols()),
00152         m_colNormsDirect(matrix.cols()),
00153         m_isInitialized(false),
00154         m_usePrescribedThreshold(false)
00155     {
00156       computeInPlace();
00157     }
00158 
00173     template<typename Rhs>
00174     inline const Solve<ColPivHouseholderQR, Rhs>
00175     solve(const MatrixBase<Rhs>& b) const
00176     {
00177       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00178       return Solve<ColPivHouseholderQR, Rhs>(*this, b.derived());
00179     }
00180 
00181     HouseholderSequenceType householderQ() const;
00182     HouseholderSequenceType matrixQ() const
00183     {
00184       return householderQ();
00185     }
00186 
00189     const MatrixType& matrixQR() const
00190     {
00191       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00192       return m_qr;
00193     }
00194 
00204     const MatrixType& matrixR() const
00205     {
00206       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00207       return m_qr;
00208     }
00209 
00210     template<typename InputType>
00211     ColPivHouseholderQR& compute(const EigenBase<InputType>& matrix);
00212 
00214     const PermutationType& colsPermutation() const
00215     {
00216       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00217       return m_colsPermutation;
00218     }
00219 
00233     typename MatrixType::RealScalar absDeterminant() const;
00234 
00247     typename MatrixType::RealScalar logAbsDeterminant() const;
00248 
00255     inline Index rank() const
00256     {
00257       using std::abs;
00258       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00259       RealScalar premultiplied_threshold = abs(m_maxpivot) * threshold();
00260       Index result = 0;
00261       for(Index i = 0; i < m_nonzero_pivots; ++i)
00262         result += (abs(m_qr.coeff(i,i)) > premultiplied_threshold);
00263       return result;
00264     }
00265 
00272     inline Index dimensionOfKernel() const
00273     {
00274       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00275       return cols() - rank();
00276     }
00277 
00285     inline bool isInjective() const
00286     {
00287       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00288       return rank() == cols();
00289     }
00290 
00298     inline bool isSurjective() const
00299     {
00300       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00301       return rank() == rows();
00302     }
00303 
00310     inline bool isInvertible() const
00311     {
00312       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00313       return isInjective() && isSurjective();
00314     }
00315 
00321     inline const Inverse<ColPivHouseholderQR> inverse() const
00322     {
00323       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00324       return Inverse<ColPivHouseholderQR>(*this);
00325     }
00326 
00327     inline Index rows() const { return m_qr.rows(); }
00328     inline Index cols() const { return m_qr.cols(); }
00329 
00334     const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
00335 
00353     ColPivHouseholderQR& setThreshold(const RealScalar& threshold)
00354     {
00355       m_usePrescribedThreshold = true;
00356       m_prescribedThreshold = threshold;
00357       return *this;
00358     }
00359 
00368     ColPivHouseholderQR& setThreshold(Default_t)
00369     {
00370       m_usePrescribedThreshold = false;
00371       return *this;
00372     }
00373 
00378     RealScalar threshold() const
00379     {
00380       eigen_assert(m_isInitialized || m_usePrescribedThreshold);
00381       return m_usePrescribedThreshold ? m_prescribedThreshold
00382       // this formula comes from experimenting (see "LU precision tuning" thread on the list)
00383       // and turns out to be identical to Higham's formula used already in LDLt.
00384                                       : NumTraits<Scalar>::epsilon() * RealScalar(m_qr.diagonalSize());
00385     }
00386 
00394     inline Index nonzeroPivots() const
00395     {
00396       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00397       return m_nonzero_pivots;
00398     }
00399 
00403     RealScalar maxPivot() const { return m_maxpivot; }
00404 
00411     ComputationInfo info() const
00412     {
00413       eigen_assert(m_isInitialized && "Decomposition is not initialized.");
00414       return Success;
00415     }
00416 
00417     #ifndef EIGEN_PARSED_BY_DOXYGEN
00418     template<typename RhsType, typename DstType>
00419     EIGEN_DEVICE_FUNC
00420     void _solve_impl(const RhsType &rhs, DstType &dst) const;
00421     #endif
00422 
00423   protected:
00424 
00425     friend class CompleteOrthogonalDecomposition<MatrixType>;
00426 
00427     static void check_template_parameters()
00428     {
00429       EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
00430     }
00431 
00432     void computeInPlace();
00433 
00434     MatrixType m_qr;
00435     HCoeffsType m_hCoeffs;
00436     PermutationType m_colsPermutation;
00437     IntRowVectorType m_colsTranspositions;
00438     RowVectorType m_temp;
00439     RealRowVectorType m_colNormsUpdated;
00440     RealRowVectorType m_colNormsDirect;
00441     bool m_isInitialized, m_usePrescribedThreshold;
00442     RealScalar m_prescribedThreshold, m_maxpivot;
00443     Index m_nonzero_pivots;
00444     Index m_det_pq;
00445 };
00446 
00447 template<typename MatrixType>
00448 typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::absDeterminant() const
00449 {
00450   using std::abs;
00451   eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00452   eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
00453   return abs(m_qr.diagonal().prod());
00454 }
00455 
00456 template<typename MatrixType>
00457 typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::logAbsDeterminant() const
00458 {
00459   eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00460   eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
00461   return m_qr.diagonal().cwiseAbs().array().log().sum();
00462 }
00463 
00470 template<typename MatrixType>
00471 template<typename InputType>
00472 ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const EigenBase<InputType>& matrix)
00473 {
00474   m_qr = matrix.derived();
00475   computeInPlace();
00476   return *this;
00477 }
00478 
00479 template<typename MatrixType>
00480 void ColPivHouseholderQR<MatrixType>::computeInPlace()
00481 {
00482   check_template_parameters();
00483 
00484   // the column permutation is stored as int indices, so just to be sure:
00485   eigen_assert(m_qr.cols()<=NumTraits<int>::highest());
00486 
00487   using std::abs;
00488 
00489   Index rows = m_qr.rows();
00490   Index cols = m_qr.cols();
00491   Index size = m_qr.diagonalSize();
00492 
00493   m_hCoeffs.resize(size);
00494 
00495   m_temp.resize(cols);
00496 
00497   m_colsTranspositions.resize(m_qr.cols());
00498   Index number_of_transpositions = 0;
00499 
00500   m_colNormsUpdated.resize(cols);
00501   m_colNormsDirect.resize(cols);
00502   for (Index k = 0; k < cols; ++k) {
00503     // colNormsDirect(k) caches the most recent directly computed norm of
00504     // column k.
00505     m_colNormsDirect.coeffRef(k) = m_qr.col(k).norm();
00506     m_colNormsUpdated.coeffRef(k) = m_colNormsDirect.coeffRef(k);
00507   }
00508 
00509   RealScalar threshold_helper =  numext::abs2<Scalar>(m_colNormsUpdated.maxCoeff() * NumTraits<Scalar>::epsilon()) / RealScalar(rows);
00510   RealScalar norm_downdate_threshold = numext::sqrt(NumTraits<Scalar>::epsilon());
00511 
00512   m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
00513   m_maxpivot = RealScalar(0);
00514 
00515   for(Index k = 0; k < size; ++k)
00516   {
00517     // first, we look up in our table m_colNormsUpdated which column has the biggest norm
00518     Index biggest_col_index;
00519     RealScalar biggest_col_sq_norm = numext::abs2(m_colNormsUpdated.tail(cols-k).maxCoeff(&biggest_col_index));
00520     biggest_col_index += k;
00521 
00522     // Track the number of meaningful pivots but do not stop the decomposition to make
00523     // sure that the initial matrix is properly reproduced. See bug 941.
00524     if(m_nonzero_pivots==size && biggest_col_sq_norm < threshold_helper * RealScalar(rows-k))
00525       m_nonzero_pivots = k;
00526 
00527     // apply the transposition to the columns
00528     m_colsTranspositions.coeffRef(k) = biggest_col_index;
00529     if(k != biggest_col_index) {
00530       m_qr.col(k).swap(m_qr.col(biggest_col_index));
00531       std::swap(m_colNormsUpdated.coeffRef(k), m_colNormsUpdated.coeffRef(biggest_col_index));
00532       std::swap(m_colNormsDirect.coeffRef(k), m_colNormsDirect.coeffRef(biggest_col_index));
00533       ++number_of_transpositions;
00534     }
00535 
00536     // generate the householder vector, store it below the diagonal
00537     RealScalar beta;
00538     m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta);
00539 
00540     // apply the householder transformation to the diagonal coefficient
00541     m_qr.coeffRef(k,k) = beta;
00542 
00543     // remember the maximum absolute value of diagonal coefficients
00544     if(abs(beta) > m_maxpivot) m_maxpivot = abs(beta);
00545 
00546     // apply the householder transformation
00547     m_qr.bottomRightCorner(rows-k, cols-k-1)
00548         .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1));
00549 
00550     // update our table of norms of the columns
00551     for (Index j = k + 1; j < cols; ++j) {
00552       // The following implements the stable norm downgrade step discussed in
00553       // http://www.netlib.org/lapack/lawnspdf/lawn176.pdf
00554       // and used in LAPACK routines xGEQPF and xGEQP3.
00555       // See lines 278-297 in http://www.netlib.org/lapack/explore-html/dc/df4/sgeqpf_8f_source.html
00556       if (m_colNormsUpdated.coeffRef(j) != 0) {
00557         RealScalar temp = abs(m_qr.coeffRef(k, j)) / m_colNormsUpdated.coeffRef(j);
00558         temp = (RealScalar(1) + temp) * (RealScalar(1) - temp);
00559         temp = temp < 0 ? 0 : temp;
00560         RealScalar temp2 = temp * numext::abs2<Scalar>(m_colNormsUpdated.coeffRef(j) /
00561                                                        m_colNormsDirect.coeffRef(j));
00562         if (temp2 <= norm_downdate_threshold) {
00563           // The updated norm has become too inaccurate so re-compute the column
00564           // norm directly.
00565           m_colNormsDirect.coeffRef(j) = m_qr.col(j).tail(rows - k - 1).norm();
00566           m_colNormsUpdated.coeffRef(j) = m_colNormsDirect.coeffRef(j);
00567         } else {
00568           m_colNormsUpdated.coeffRef(j) *= numext::sqrt(temp);
00569         }
00570       }
00571     }
00572   }
00573 
00574   m_colsPermutation.setIdentity(PermIndexType(cols));
00575   for(PermIndexType k = 0; k < size/*m_nonzero_pivots*/; ++k)
00576     m_colsPermutation.applyTranspositionOnTheRight(k, PermIndexType(m_colsTranspositions.coeff(k)));
00577 
00578   m_det_pq = (number_of_transpositions%2) ? -1 : 1;
00579   m_isInitialized = true;
00580 }
00581 
00582 #ifndef EIGEN_PARSED_BY_DOXYGEN
00583 template<typename _MatrixType>
00584 template<typename RhsType, typename DstType>
00585 void ColPivHouseholderQR<_MatrixType>::_solve_impl(const RhsType &rhs, DstType &dst) const
00586 {
00587   eigen_assert(rhs.rows() == rows());
00588 
00589   const Index nonzero_pivots = nonzeroPivots();
00590 
00591   if(nonzero_pivots == 0)
00592   {
00593     dst.setZero();
00594     return;
00595   }
00596 
00597   typename RhsType::PlainObject c(rhs);
00598 
00599   // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
00600   c.applyOnTheLeft(householderSequence(m_qr, m_hCoeffs)
00601                     .setLength(nonzero_pivots)
00602                     .transpose()
00603     );
00604 
00605   m_qr.topLeftCorner(nonzero_pivots, nonzero_pivots)
00606       .template triangularView<Upper>()
00607       .solveInPlace(c.topRows(nonzero_pivots));
00608 
00609   for(Index i = 0; i < nonzero_pivots; ++i) dst.row(m_colsPermutation.indices().coeff(i)) = c.row(i);
00610   for(Index i = nonzero_pivots; i < cols(); ++i) dst.row(m_colsPermutation.indices().coeff(i)).setZero();
00611 }
00612 #endif
00613 
00614 namespace internal {
00615 
00616 template<typename DstXprType, typename MatrixType>
00617 struct Assignment<DstXprType, Inverse<ColPivHouseholderQR<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename ColPivHouseholderQR<MatrixType>::Scalar>, Dense2Dense>
00618 {
00619   typedef ColPivHouseholderQR<MatrixType> QrType;
00620   typedef Inverse<QrType> SrcXprType;
00621   static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename QrType::Scalar> &)
00622   {
00623     dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.cols()));
00624   }
00625 };
00626 
00627 } // end namespace internal
00628 
00632 template<typename MatrixType>
00633 typename ColPivHouseholderQR<MatrixType>::HouseholderSequenceType ColPivHouseholderQR<MatrixType>
00634   ::householderQ() const
00635 {
00636   eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00637   return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate());
00638 }
00639 
00644 template<typename Derived>
00645 const ColPivHouseholderQR<typename MatrixBase<Derived>::PlainObject>
00646 MatrixBase<Derived>::colPivHouseholderQr() const
00647 {
00648   return ColPivHouseholderQR<PlainObject>(eval());
00649 }
00650 
00651 } // end namespace Eigen
00652 
00653 #endif // EIGEN_COLPIVOTINGHOUSEHOLDERQR_H
 All Classes Functions Variables Typedefs Enumerations Enumerator Friends