Eigen  3.3.3
FullPivHouseholderQR.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_FULLPIVOTINGHOUSEHOLDERQR_H
00012 #define EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H
00013 
00014 namespace Eigen { 
00015 
00016 namespace internal {
00017 
00018 template<typename _MatrixType> struct traits<FullPivHouseholderQR<_MatrixType> >
00019  : traits<_MatrixType>
00020 {
00021   enum { Flags = 0 };
00022 };
00023 
00024 template<typename MatrixType> struct FullPivHouseholderQRMatrixQReturnType;
00025 
00026 template<typename MatrixType>
00027 struct traits<FullPivHouseholderQRMatrixQReturnType<MatrixType> >
00028 {
00029   typedef typename MatrixType::PlainObject ReturnType;
00030 };
00031 
00032 } // end namespace internal
00033 
00057 template<typename _MatrixType> class FullPivHouseholderQR
00058 {
00059   public:
00060 
00061     typedef _MatrixType MatrixType;
00062     enum {
00063       RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00064       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00065       MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00066       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00067     };
00068     typedef typename MatrixType::Scalar Scalar;
00069     typedef typename MatrixType::RealScalar RealScalar;
00070     // FIXME should be int
00071     typedef typename MatrixType::StorageIndex StorageIndex;
00072     typedef internal::FullPivHouseholderQRMatrixQReturnType<MatrixType> MatrixQReturnType;
00073     typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
00074     typedef Matrix<StorageIndex, 1,
00075                    EIGEN_SIZE_MIN_PREFER_DYNAMIC(ColsAtCompileTime,RowsAtCompileTime), RowMajor, 1,
00076                    EIGEN_SIZE_MIN_PREFER_FIXED(MaxColsAtCompileTime,MaxRowsAtCompileTime)> IntDiagSizeVectorType;
00077     typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationType;
00078     typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
00079     typedef typename internal::plain_col_type<MatrixType>::type ColVectorType;
00080     typedef typename MatrixType::PlainObject PlainObject;
00081 
00087     FullPivHouseholderQR()
00088       : m_qr(),
00089         m_hCoeffs(),
00090         m_rows_transpositions(),
00091         m_cols_transpositions(),
00092         m_cols_permutation(),
00093         m_temp(),
00094         m_isInitialized(false),
00095         m_usePrescribedThreshold(false) {}
00096 
00103     FullPivHouseholderQR(Index rows, Index cols)
00104       : m_qr(rows, cols),
00105         m_hCoeffs((std::min)(rows,cols)),
00106         m_rows_transpositions((std::min)(rows,cols)),
00107         m_cols_transpositions((std::min)(rows,cols)),
00108         m_cols_permutation(cols),
00109         m_temp(cols),
00110         m_isInitialized(false),
00111         m_usePrescribedThreshold(false) {}
00112 
00125     template<typename InputType>
00126     explicit FullPivHouseholderQR(const EigenBase<InputType>& matrix)
00127       : m_qr(matrix.rows(), matrix.cols()),
00128         m_hCoeffs((std::min)(matrix.rows(), matrix.cols())),
00129         m_rows_transpositions((std::min)(matrix.rows(), matrix.cols())),
00130         m_cols_transpositions((std::min)(matrix.rows(), matrix.cols())),
00131         m_cols_permutation(matrix.cols()),
00132         m_temp(matrix.cols()),
00133         m_isInitialized(false),
00134         m_usePrescribedThreshold(false)
00135     {
00136       compute(matrix.derived());
00137     }
00138 
00145     template<typename InputType>
00146     explicit FullPivHouseholderQR(EigenBase<InputType>& matrix)
00147       : m_qr(matrix.derived()),
00148         m_hCoeffs((std::min)(matrix.rows(), matrix.cols())),
00149         m_rows_transpositions((std::min)(matrix.rows(), matrix.cols())),
00150         m_cols_transpositions((std::min)(matrix.rows(), matrix.cols())),
00151         m_cols_permutation(matrix.cols()),
00152         m_temp(matrix.cols()),
00153         m_isInitialized(false),
00154         m_usePrescribedThreshold(false)
00155     {
00156       computeInPlace();
00157     }
00158 
00174     template<typename Rhs>
00175     inline const Solve<FullPivHouseholderQR, Rhs>
00176     solve(const MatrixBase<Rhs>& b) const
00177     {
00178       eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00179       return Solve<FullPivHouseholderQR, Rhs>(*this, b.derived());
00180     }
00181 
00184     MatrixQReturnType matrixQ(void) const;
00185 
00188     const MatrixType& matrixQR() const
00189     {
00190       eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00191       return m_qr;
00192     }
00193 
00194     template<typename InputType>
00195     FullPivHouseholderQR& compute(const EigenBase<InputType>& matrix);
00196 
00198     const PermutationType& colsPermutation() const
00199     {
00200       eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00201       return m_cols_permutation;
00202     }
00203 
00205     const IntDiagSizeVectorType& rowsTranspositions() const
00206     {
00207       eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00208       return m_rows_transpositions;
00209     }
00210 
00224     typename MatrixType::RealScalar absDeterminant() const;
00225 
00238     typename MatrixType::RealScalar logAbsDeterminant() const;
00239 
00246     inline Index rank() const
00247     {
00248       using std::abs;
00249       eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00250       RealScalar premultiplied_threshold = abs(m_maxpivot) * threshold();
00251       Index result = 0;
00252       for(Index i = 0; i < m_nonzero_pivots; ++i)
00253         result += (abs(m_qr.coeff(i,i)) > premultiplied_threshold);
00254       return result;
00255     }
00256 
00263     inline Index dimensionOfKernel() const
00264     {
00265       eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00266       return cols() - rank();
00267     }
00268 
00276     inline bool isInjective() const
00277     {
00278       eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00279       return rank() == cols();
00280     }
00281 
00289     inline bool isSurjective() const
00290     {
00291       eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00292       return rank() == rows();
00293     }
00294 
00301     inline bool isInvertible() const
00302     {
00303       eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00304       return isInjective() && isSurjective();
00305     }
00306 
00312     inline const Inverse<FullPivHouseholderQR> inverse() const
00313     {
00314       eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00315       return Inverse<FullPivHouseholderQR>(*this);
00316     }
00317 
00318     inline Index rows() const { return m_qr.rows(); }
00319     inline Index cols() const { return m_qr.cols(); }
00320     
00325     const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
00326 
00344     FullPivHouseholderQR& setThreshold(const RealScalar& threshold)
00345     {
00346       m_usePrescribedThreshold = true;
00347       m_prescribedThreshold = threshold;
00348       return *this;
00349     }
00350 
00359     FullPivHouseholderQR& setThreshold(Default_t)
00360     {
00361       m_usePrescribedThreshold = false;
00362       return *this;
00363     }
00364 
00369     RealScalar threshold() const
00370     {
00371       eigen_assert(m_isInitialized || m_usePrescribedThreshold);
00372       return m_usePrescribedThreshold ? m_prescribedThreshold
00373       // this formula comes from experimenting (see "LU precision tuning" thread on the list)
00374       // and turns out to be identical to Higham's formula used already in LDLt.
00375                                       : NumTraits<Scalar>::epsilon() * RealScalar(m_qr.diagonalSize());
00376     }
00377 
00385     inline Index nonzeroPivots() const
00386     {
00387       eigen_assert(m_isInitialized && "LU is not initialized.");
00388       return m_nonzero_pivots;
00389     }
00390 
00394     RealScalar maxPivot() const { return m_maxpivot; }
00395     
00396     #ifndef EIGEN_PARSED_BY_DOXYGEN
00397     template<typename RhsType, typename DstType>
00398     EIGEN_DEVICE_FUNC
00399     void _solve_impl(const RhsType &rhs, DstType &dst) const;
00400     #endif
00401 
00402   protected:
00403     
00404     static void check_template_parameters()
00405     {
00406       EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
00407     }
00408     
00409     void computeInPlace();
00410     
00411     MatrixType m_qr;
00412     HCoeffsType m_hCoeffs;
00413     IntDiagSizeVectorType m_rows_transpositions;
00414     IntDiagSizeVectorType m_cols_transpositions;
00415     PermutationType m_cols_permutation;
00416     RowVectorType m_temp;
00417     bool m_isInitialized, m_usePrescribedThreshold;
00418     RealScalar m_prescribedThreshold, m_maxpivot;
00419     Index m_nonzero_pivots;
00420     RealScalar m_precision;
00421     Index m_det_pq;
00422 };
00423 
00424 template<typename MatrixType>
00425 typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType>::absDeterminant() const
00426 {
00427   using std::abs;
00428   eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00429   eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
00430   return abs(m_qr.diagonal().prod());
00431 }
00432 
00433 template<typename MatrixType>
00434 typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType>::logAbsDeterminant() const
00435 {
00436   eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00437   eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
00438   return m_qr.diagonal().cwiseAbs().array().log().sum();
00439 }
00440 
00447 template<typename MatrixType>
00448 template<typename InputType>
00449 FullPivHouseholderQR<MatrixType>& FullPivHouseholderQR<MatrixType>::compute(const EigenBase<InputType>& matrix)
00450 {
00451   m_qr = matrix.derived();
00452   computeInPlace();
00453   return *this;
00454 }
00455 
00456 template<typename MatrixType>
00457 void FullPivHouseholderQR<MatrixType>::computeInPlace()
00458 {
00459   check_template_parameters();
00460 
00461   using std::abs;
00462   Index rows = m_qr.rows();
00463   Index cols = m_qr.cols();
00464   Index size = (std::min)(rows,cols);
00465 
00466   
00467   m_hCoeffs.resize(size);
00468 
00469   m_temp.resize(cols);
00470 
00471   m_precision = NumTraits<Scalar>::epsilon() * RealScalar(size);
00472 
00473   m_rows_transpositions.resize(size);
00474   m_cols_transpositions.resize(size);
00475   Index number_of_transpositions = 0;
00476 
00477   RealScalar biggest(0);
00478 
00479   m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
00480   m_maxpivot = RealScalar(0);
00481 
00482   for (Index k = 0; k < size; ++k)
00483   {
00484     Index row_of_biggest_in_corner, col_of_biggest_in_corner;
00485     typedef internal::scalar_score_coeff_op<Scalar> Scoring;
00486     typedef typename Scoring::result_type Score;
00487 
00488     Score score = m_qr.bottomRightCorner(rows-k, cols-k)
00489                       .unaryExpr(Scoring())
00490                       .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner);
00491     row_of_biggest_in_corner += k;
00492     col_of_biggest_in_corner += k;
00493     RealScalar biggest_in_corner = internal::abs_knowing_score<Scalar>()(m_qr(row_of_biggest_in_corner, col_of_biggest_in_corner), score);
00494     if(k==0) biggest = biggest_in_corner;
00495 
00496     // if the corner is negligible, then we have less than full rank, and we can finish early
00497     if(internal::isMuchSmallerThan(biggest_in_corner, biggest, m_precision))
00498     {
00499       m_nonzero_pivots = k;
00500       for(Index i = k; i < size; i++)
00501       {
00502         m_rows_transpositions.coeffRef(i) = i;
00503         m_cols_transpositions.coeffRef(i) = i;
00504         m_hCoeffs.coeffRef(i) = Scalar(0);
00505       }
00506       break;
00507     }
00508 
00509     m_rows_transpositions.coeffRef(k) = row_of_biggest_in_corner;
00510     m_cols_transpositions.coeffRef(k) = col_of_biggest_in_corner;
00511     if(k != row_of_biggest_in_corner) {
00512       m_qr.row(k).tail(cols-k).swap(m_qr.row(row_of_biggest_in_corner).tail(cols-k));
00513       ++number_of_transpositions;
00514     }
00515     if(k != col_of_biggest_in_corner) {
00516       m_qr.col(k).swap(m_qr.col(col_of_biggest_in_corner));
00517       ++number_of_transpositions;
00518     }
00519 
00520     RealScalar beta;
00521     m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta);
00522     m_qr.coeffRef(k,k) = beta;
00523 
00524     // remember the maximum absolute value of diagonal coefficients
00525     if(abs(beta) > m_maxpivot) m_maxpivot = abs(beta);
00526 
00527     m_qr.bottomRightCorner(rows-k, cols-k-1)
00528         .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1));
00529   }
00530 
00531   m_cols_permutation.setIdentity(cols);
00532   for(Index k = 0; k < size; ++k)
00533     m_cols_permutation.applyTranspositionOnTheRight(k, m_cols_transpositions.coeff(k));
00534 
00535   m_det_pq = (number_of_transpositions%2) ? -1 : 1;
00536   m_isInitialized = true;
00537 }
00538 
00539 #ifndef EIGEN_PARSED_BY_DOXYGEN
00540 template<typename _MatrixType>
00541 template<typename RhsType, typename DstType>
00542 void FullPivHouseholderQR<_MatrixType>::_solve_impl(const RhsType &rhs, DstType &dst) const
00543 {
00544   eigen_assert(rhs.rows() == rows());
00545   const Index l_rank = rank();
00546 
00547   // FIXME introduce nonzeroPivots() and use it here. and more generally,
00548   // make the same improvements in this dec as in FullPivLU.
00549   if(l_rank==0)
00550   {
00551     dst.setZero();
00552     return;
00553   }
00554 
00555   typename RhsType::PlainObject c(rhs);
00556 
00557   Matrix<Scalar,1,RhsType::ColsAtCompileTime> temp(rhs.cols());
00558   for (Index k = 0; k < l_rank; ++k)
00559   {
00560     Index remainingSize = rows()-k;
00561     c.row(k).swap(c.row(m_rows_transpositions.coeff(k)));
00562     c.bottomRightCorner(remainingSize, rhs.cols())
00563       .applyHouseholderOnTheLeft(m_qr.col(k).tail(remainingSize-1),
00564                                m_hCoeffs.coeff(k), &temp.coeffRef(0));
00565   }
00566 
00567   m_qr.topLeftCorner(l_rank, l_rank)
00568       .template triangularView<Upper>()
00569       .solveInPlace(c.topRows(l_rank));
00570 
00571   for(Index i = 0; i < l_rank; ++i) dst.row(m_cols_permutation.indices().coeff(i)) = c.row(i);
00572   for(Index i = l_rank; i < cols(); ++i) dst.row(m_cols_permutation.indices().coeff(i)).setZero();
00573 }
00574 #endif
00575 
00576 namespace internal {
00577   
00578 template<typename DstXprType, typename MatrixType>
00579 struct Assignment<DstXprType, Inverse<FullPivHouseholderQR<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename FullPivHouseholderQR<MatrixType>::Scalar>, Dense2Dense>
00580 {
00581   typedef FullPivHouseholderQR<MatrixType> QrType;
00582   typedef Inverse<QrType> SrcXprType;
00583   static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename QrType::Scalar> &)
00584   {    
00585     dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.cols()));
00586   }
00587 };
00588 
00595 template<typename MatrixType> struct FullPivHouseholderQRMatrixQReturnType
00596   : public ReturnByValue<FullPivHouseholderQRMatrixQReturnType<MatrixType> >
00597 {
00598 public:
00599   typedef typename FullPivHouseholderQR<MatrixType>::IntDiagSizeVectorType IntDiagSizeVectorType;
00600   typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
00601   typedef Matrix<typename MatrixType::Scalar, 1, MatrixType::RowsAtCompileTime, RowMajor, 1,
00602                  MatrixType::MaxRowsAtCompileTime> WorkVectorType;
00603 
00604   FullPivHouseholderQRMatrixQReturnType(const MatrixType&       qr,
00605                                         const HCoeffsType&      hCoeffs,
00606                                         const IntDiagSizeVectorType& rowsTranspositions)
00607     : m_qr(qr),
00608       m_hCoeffs(hCoeffs),
00609       m_rowsTranspositions(rowsTranspositions)
00610   {}
00611 
00612   template <typename ResultType>
00613   void evalTo(ResultType& result) const
00614   {
00615     const Index rows = m_qr.rows();
00616     WorkVectorType workspace(rows);
00617     evalTo(result, workspace);
00618   }
00619 
00620   template <typename ResultType>
00621   void evalTo(ResultType& result, WorkVectorType& workspace) const
00622   {
00623     using numext::conj;
00624     // compute the product H'_0 H'_1 ... H'_n-1,
00625     // where H_k is the k-th Householder transformation I - h_k v_k v_k'
00626     // and v_k is the k-th Householder vector [1,m_qr(k+1,k), m_qr(k+2,k), ...]
00627     const Index rows = m_qr.rows();
00628     const Index cols = m_qr.cols();
00629     const Index size = (std::min)(rows, cols);
00630     workspace.resize(rows);
00631     result.setIdentity(rows, rows);
00632     for (Index k = size-1; k >= 0; k--)
00633     {
00634       result.block(k, k, rows-k, rows-k)
00635             .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), conj(m_hCoeffs.coeff(k)), &workspace.coeffRef(k));
00636       result.row(k).swap(result.row(m_rowsTranspositions.coeff(k)));
00637     }
00638   }
00639 
00640   Index rows() const { return m_qr.rows(); }
00641   Index cols() const { return m_qr.rows(); }
00642 
00643 protected:
00644   typename MatrixType::Nested m_qr;
00645   typename HCoeffsType::Nested m_hCoeffs;
00646   typename IntDiagSizeVectorType::Nested m_rowsTranspositions;
00647 };
00648 
00649 // template<typename MatrixType>
00650 // struct evaluator<FullPivHouseholderQRMatrixQReturnType<MatrixType> >
00651 //  : public evaluator<ReturnByValue<FullPivHouseholderQRMatrixQReturnType<MatrixType> > >
00652 // {};
00653 
00654 } // end namespace internal
00655 
00656 template<typename MatrixType>
00657 inline typename FullPivHouseholderQR<MatrixType>::MatrixQReturnType FullPivHouseholderQR<MatrixType>::matrixQ() const
00658 {
00659   eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00660   return MatrixQReturnType(m_qr, m_hCoeffs, m_rows_transpositions);
00661 }
00662 
00667 template<typename Derived>
00668 const FullPivHouseholderQR<typename MatrixBase<Derived>::PlainObject>
00669 MatrixBase<Derived>::fullPivHouseholderQr() const
00670 {
00671   return FullPivHouseholderQR<PlainObject>(eval());
00672 }
00673 
00674 } // end namespace Eigen
00675 
00676 #endif // EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H
 All Classes Functions Variables Typedefs Enumerations Enumerator Friends