FullPivHouseholderQR.h
Go to the documentation of this file.
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 // Eigen is free software; you can redistribute it and/or
00008 // modify it under the terms of the GNU Lesser General Public
00009 // License as published by the Free Software Foundation; either
00010 // version 3 of the License, or (at your option) any later version.
00011 //
00012 // Alternatively, you can redistribute it and/or
00013 // modify it under the terms of the GNU General Public License as
00014 // published by the Free Software Foundation; either version 2 of
00015 // the License, or (at your option) any later version.
00016 //
00017 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00018 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00019 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00020 // GNU General Public License for more details.
00021 //
00022 // You should have received a copy of the GNU Lesser General Public
00023 // License and a copy of the GNU General Public License along with
00024 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00025 
00026 #ifndef EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H
00027 #define EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H
00028 
00029 namespace Eigen { 
00030 
00031 namespace internal {
00032 
00033 template<typename MatrixType> struct FullPivHouseholderQRMatrixQReturnType;
00034 
00035 template<typename MatrixType>
00036 struct traits<FullPivHouseholderQRMatrixQReturnType<MatrixType> >
00037 {
00038   typedef typename MatrixType::PlainObject ReturnType;
00039 };
00040 
00041 }
00042 
00064 template<typename _MatrixType> class FullPivHouseholderQR
00065 {
00066   public:
00067 
00068     typedef _MatrixType MatrixType;
00069     enum {
00070       RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00071       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00072       Options = MatrixType::Options,
00073       MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00074       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00075     };
00076     typedef typename MatrixType::Scalar Scalar;
00077     typedef typename MatrixType::RealScalar RealScalar;
00078     typedef typename MatrixType::Index Index;
00079     typedef internal::FullPivHouseholderQRMatrixQReturnType<MatrixType> MatrixQReturnType;
00080     typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
00081     typedef Matrix<Index, 1, ColsAtCompileTime, RowMajor, 1, MaxColsAtCompileTime> IntRowVectorType;
00082     typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationType;
00083     typedef typename internal::plain_col_type<MatrixType, Index>::type IntColVectorType;
00084     typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
00085     typedef typename internal::plain_col_type<MatrixType>::type ColVectorType;
00086 
00092     FullPivHouseholderQR()
00093       : m_qr(),
00094         m_hCoeffs(),
00095         m_rows_transpositions(),
00096         m_cols_transpositions(),
00097         m_cols_permutation(),
00098         m_temp(),
00099         m_isInitialized(false),
00100         m_usePrescribedThreshold(false) {}
00101 
00108     FullPivHouseholderQR(Index rows, Index cols)
00109       : m_qr(rows, cols),
00110         m_hCoeffs((std::min)(rows,cols)),
00111         m_rows_transpositions(rows),
00112         m_cols_transpositions(cols),
00113         m_cols_permutation(cols),
00114         m_temp((std::min)(rows,cols)),
00115         m_isInitialized(false),
00116         m_usePrescribedThreshold(false) {}
00117 
00118     FullPivHouseholderQR(const MatrixType& matrix)
00119       : m_qr(matrix.rows(), matrix.cols()),
00120         m_hCoeffs((std::min)(matrix.rows(), matrix.cols())),
00121         m_rows_transpositions(matrix.rows()),
00122         m_cols_transpositions(matrix.cols()),
00123         m_cols_permutation(matrix.cols()),
00124         m_temp((std::min)(matrix.rows(), matrix.cols())),
00125         m_isInitialized(false),
00126         m_usePrescribedThreshold(false)
00127     {
00128       compute(matrix);
00129     }
00130 
00148     template<typename Rhs>
00149     inline const internal::solve_retval<FullPivHouseholderQR, Rhs>
00150     solve(const MatrixBase<Rhs>& b) const
00151     {
00152       eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00153       return internal::solve_retval<FullPivHouseholderQR, Rhs>(*this, b.derived());
00154     }
00155 
00158     MatrixQReturnType matrixQ(void) const;
00159 
00162     const MatrixType& matrixQR() const
00163     {
00164       eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00165       return m_qr;
00166     }
00167 
00168     FullPivHouseholderQR& compute(const MatrixType& matrix);
00169 
00170     const PermutationType& colsPermutation() const
00171     {
00172       eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00173       return m_cols_permutation;
00174     }
00175 
00176     const IntColVectorType& rowsTranspositions() const
00177     {
00178       eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00179       return m_rows_transpositions;
00180     }
00181 
00195     typename MatrixType::RealScalar absDeterminant() const;
00196 
00209     typename MatrixType::RealScalar logAbsDeterminant() const;
00210 
00217     inline Index rank() const
00218     {
00219       eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00220       RealScalar premultiplied_threshold = internal::abs(m_maxpivot) * threshold();
00221       Index result = 0;
00222       for(Index i = 0; i < m_nonzero_pivots; ++i)
00223         result += (internal::abs(m_qr.coeff(i,i)) > premultiplied_threshold);
00224       return result;
00225     }
00226 
00233     inline Index dimensionOfKernel() const
00234     {
00235       eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00236       return cols() - rank();
00237     }
00238 
00246     inline bool isInjective() const
00247     {
00248       eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00249       return rank() == cols();
00250     }
00251 
00259     inline bool isSurjective() const
00260     {
00261       eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00262       return rank() == rows();
00263     }
00264 
00271     inline bool isInvertible() const
00272     {
00273       eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00274       return isInjective() && isSurjective();
00275     }
00276     inline const
00282     internal::solve_retval<FullPivHouseholderQR, typename MatrixType::IdentityReturnType>
00283     inverse() const
00284     {
00285       eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00286       return internal::solve_retval<FullPivHouseholderQR,typename MatrixType::IdentityReturnType>
00287                (*this, MatrixType::Identity(m_qr.rows(), m_qr.cols()));
00288     }
00289 
00290     inline Index rows() const { return m_qr.rows(); }
00291     inline Index cols() const { return m_qr.cols(); }
00292     const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
00293 
00311     FullPivHouseholderQR& setThreshold(const RealScalar& threshold)
00312     {
00313       m_usePrescribedThreshold = true;
00314       m_prescribedThreshold = threshold;
00315       return *this;
00316     }
00317 
00326     FullPivHouseholderQR& setThreshold(Default_t)
00327     {
00328       m_usePrescribedThreshold = false;
00329       return *this;
00330     }
00331 
00336     RealScalar threshold() const
00337     {
00338       eigen_assert(m_isInitialized || m_usePrescribedThreshold);
00339       return m_usePrescribedThreshold ? m_prescribedThreshold
00340       // this formula comes from experimenting (see "LU precision tuning" thread on the list)
00341       // and turns out to be identical to Higham's formula used already in LDLt.
00342                                       : NumTraits<Scalar>::epsilon() * m_qr.diagonalSize();
00343     }
00344 
00352     inline Index nonzeroPivots() const
00353     {
00354       eigen_assert(m_isInitialized && "LU is not initialized.");
00355       return m_nonzero_pivots;
00356     }
00357 
00361     RealScalar maxPivot() const { return m_maxpivot; }
00362 
00363   protected:
00364     MatrixType m_qr;
00365     HCoeffsType m_hCoeffs;
00366     IntColVectorType m_rows_transpositions;
00367     IntRowVectorType m_cols_transpositions;
00368     PermutationType m_cols_permutation;
00369     RowVectorType m_temp;
00370     bool m_isInitialized, m_usePrescribedThreshold;
00371     RealScalar m_prescribedThreshold, m_maxpivot;
00372     Index m_nonzero_pivots;
00373     RealScalar m_precision;
00374     Index m_det_pq;
00375 };
00376 
00377 template<typename MatrixType>
00378 typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType>::absDeterminant() const
00379 {
00380   eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00381   eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
00382   return internal::abs(m_qr.diagonal().prod());
00383 }
00384 
00385 template<typename MatrixType>
00386 typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType>::logAbsDeterminant() const
00387 {
00388   eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00389   eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
00390   return m_qr.diagonal().cwiseAbs().array().log().sum();
00391 }
00392 
00393 template<typename MatrixType>
00394 FullPivHouseholderQR<MatrixType>& FullPivHouseholderQR<MatrixType>::compute(const MatrixType& matrix)
00395 {
00396   Index rows = matrix.rows();
00397   Index cols = matrix.cols();
00398   Index size = (std::min)(rows,cols);
00399 
00400   m_qr = matrix;
00401   m_hCoeffs.resize(size);
00402 
00403   m_temp.resize(cols);
00404 
00405   m_precision = NumTraits<Scalar>::epsilon() * size;
00406 
00407   m_rows_transpositions.resize(matrix.rows());
00408   m_cols_transpositions.resize(matrix.cols());
00409   Index number_of_transpositions = 0;
00410 
00411   RealScalar biggest(0);
00412 
00413   m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
00414   m_maxpivot = RealScalar(0);
00415 
00416   for (Index k = 0; k < size; ++k)
00417   {
00418     Index row_of_biggest_in_corner, col_of_biggest_in_corner;
00419     RealScalar biggest_in_corner;
00420 
00421     biggest_in_corner = m_qr.bottomRightCorner(rows-k, cols-k)
00422                             .cwiseAbs()
00423                             .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner);
00424     row_of_biggest_in_corner += k;
00425     col_of_biggest_in_corner += k;
00426     if(k==0) biggest = biggest_in_corner;
00427 
00428     // if the corner is negligible, then we have less than full rank, and we can finish early
00429     if(internal::isMuchSmallerThan(biggest_in_corner, biggest, m_precision))
00430     {
00431       m_nonzero_pivots = k;
00432       for(Index i = k; i < size; i++)
00433       {
00434         m_rows_transpositions.coeffRef(i) = i;
00435         m_cols_transpositions.coeffRef(i) = i;
00436         m_hCoeffs.coeffRef(i) = Scalar(0);
00437       }
00438       break;
00439     }
00440 
00441     m_rows_transpositions.coeffRef(k) = row_of_biggest_in_corner;
00442     m_cols_transpositions.coeffRef(k) = col_of_biggest_in_corner;
00443     if(k != row_of_biggest_in_corner) {
00444       m_qr.row(k).tail(cols-k).swap(m_qr.row(row_of_biggest_in_corner).tail(cols-k));
00445       ++number_of_transpositions;
00446     }
00447     if(k != col_of_biggest_in_corner) {
00448       m_qr.col(k).swap(m_qr.col(col_of_biggest_in_corner));
00449       ++number_of_transpositions;
00450     }
00451 
00452     RealScalar beta;
00453     m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta);
00454     m_qr.coeffRef(k,k) = beta;
00455 
00456     // remember the maximum absolute value of diagonal coefficients
00457     if(internal::abs(beta) > m_maxpivot) m_maxpivot = internal::abs(beta);
00458 
00459     m_qr.bottomRightCorner(rows-k, cols-k-1)
00460         .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1));
00461   }
00462 
00463   m_cols_permutation.setIdentity(cols);
00464   for(Index k = 0; k < size; ++k)
00465     m_cols_permutation.applyTranspositionOnTheRight(k, m_cols_transpositions.coeff(k));
00466 
00467   m_det_pq = (number_of_transpositions%2) ? -1 : 1;
00468   m_isInitialized = true;
00469 
00470   return *this;
00471 }
00472 
00473 namespace internal {
00474 
00475 template<typename _MatrixType, typename Rhs>
00476 struct solve_retval<FullPivHouseholderQR<_MatrixType>, Rhs>
00477   : solve_retval_base<FullPivHouseholderQR<_MatrixType>, Rhs>
00478 {
00479   EIGEN_MAKE_SOLVE_HELPERS(FullPivHouseholderQR<_MatrixType>,Rhs)
00480 
00481   template<typename Dest> void evalTo(Dest& dst) const
00482   {
00483     const Index rows = dec().rows(), cols = dec().cols();
00484     eigen_assert(rhs().rows() == rows);
00485 
00486     // FIXME introduce nonzeroPivots() and use it here. and more generally,
00487     // make the same improvements in this dec as in FullPivLU.
00488     if(dec().rank()==0)
00489     {
00490       dst.setZero();
00491       return;
00492     }
00493 
00494     typename Rhs::PlainObject c(rhs());
00495 
00496     Matrix<Scalar,1,Rhs::ColsAtCompileTime> temp(rhs().cols());
00497     for (Index k = 0; k < dec().rank(); ++k)
00498     {
00499       Index remainingSize = rows-k;
00500       c.row(k).swap(c.row(dec().rowsTranspositions().coeff(k)));
00501       c.bottomRightCorner(remainingSize, rhs().cols())
00502        .applyHouseholderOnTheLeft(dec().matrixQR().col(k).tail(remainingSize-1),
00503                                   dec().hCoeffs().coeff(k), &temp.coeffRef(0));
00504     }
00505 
00506     if(!dec().isSurjective())
00507     {
00508       // is c is in the image of R ?
00509       RealScalar biggest_in_upper_part_of_c = c.topRows(   dec().rank()     ).cwiseAbs().maxCoeff();
00510       RealScalar biggest_in_lower_part_of_c = c.bottomRows(rows-dec().rank()).cwiseAbs().maxCoeff();
00511       // FIXME brain dead
00512       const RealScalar m_precision = NumTraits<Scalar>::epsilon() * (std::min)(rows,cols);
00513       // this internal:: prefix is needed by at least gcc 3.4 and ICC
00514       if(!internal::isMuchSmallerThan(biggest_in_lower_part_of_c, biggest_in_upper_part_of_c, m_precision))
00515         return;
00516     }
00517     dec().matrixQR()
00518        .topLeftCorner(dec().rank(), dec().rank())
00519        .template triangularView<Upper>()
00520        .solveInPlace(c.topRows(dec().rank()));
00521 
00522     for(Index i = 0; i < dec().rank(); ++i) dst.row(dec().colsPermutation().indices().coeff(i)) = c.row(i);
00523     for(Index i = dec().rank(); i < cols; ++i) dst.row(dec().colsPermutation().indices().coeff(i)).setZero();
00524   }
00525 };
00526 
00533 template<typename MatrixType> struct FullPivHouseholderQRMatrixQReturnType
00534   : public ReturnByValue<FullPivHouseholderQRMatrixQReturnType<MatrixType> >
00535 {
00536 public:
00537   typedef typename MatrixType::Index Index;
00538   typedef typename internal::plain_col_type<MatrixType, Index>::type IntColVectorType;
00539   typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
00540   typedef Matrix<typename MatrixType::Scalar, 1, MatrixType::RowsAtCompileTime, RowMajor, 1,
00541                  MatrixType::MaxRowsAtCompileTime> WorkVectorType;
00542 
00543   FullPivHouseholderQRMatrixQReturnType(const MatrixType&       qr,
00544                                         const HCoeffsType&      hCoeffs,
00545                                         const IntColVectorType& rowsTranspositions)
00546     : m_qr(qr),
00547       m_hCoeffs(hCoeffs),
00548       m_rowsTranspositions(rowsTranspositions)
00549       {}
00550 
00551   template <typename ResultType>
00552   void evalTo(ResultType& result) const
00553   {
00554     const Index rows = m_qr.rows();
00555     WorkVectorType workspace(rows);
00556     evalTo(result, workspace);
00557   }
00558 
00559   template <typename ResultType>
00560   void evalTo(ResultType& result, WorkVectorType& workspace) const
00561   {
00562     // compute the product H'_0 H'_1 ... H'_n-1,
00563     // where H_k is the k-th Householder transformation I - h_k v_k v_k'
00564     // and v_k is the k-th Householder vector [1,m_qr(k+1,k), m_qr(k+2,k), ...]
00565     const Index rows = m_qr.rows();
00566     const Index cols = m_qr.cols();
00567     const Index size = (std::min)(rows, cols);
00568     workspace.resize(rows);
00569     result.setIdentity(rows, rows);
00570     for (Index k = size-1; k >= 0; k--)
00571     {
00572       result.block(k, k, rows-k, rows-k)
00573             .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), internal::conj(m_hCoeffs.coeff(k)), &workspace.coeffRef(k));
00574       result.row(k).swap(result.row(m_rowsTranspositions.coeff(k)));
00575     }
00576   }
00577 
00578     Index rows() const { return m_qr.rows(); }
00579     Index cols() const { return m_qr.rows(); }
00580 
00581 protected:
00582   typename MatrixType::Nested m_qr;
00583   typename HCoeffsType::Nested m_hCoeffs;
00584   typename IntColVectorType::Nested m_rowsTranspositions;
00585 };
00586 
00587 } // end namespace internal
00588 
00589 template<typename MatrixType>
00590 inline typename FullPivHouseholderQR<MatrixType>::MatrixQReturnType FullPivHouseholderQR<MatrixType>::matrixQ() const
00591 {
00592   eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00593   return MatrixQReturnType(m_qr, m_hCoeffs, m_rows_transpositions);
00594 }
00595 
00600 template<typename Derived>
00601 const FullPivHouseholderQR<typename MatrixBase<Derived>::PlainObject>
00602 MatrixBase<Derived>::fullPivHouseholderQr() const
00603 {
00604   return FullPivHouseholderQR<PlainObject>(eval());
00605 }
00606 
00607 } // end namespace Eigen
00608 
00609 #endif // EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H