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 
00050 template<typename _MatrixType> class FullPivHouseholderQR
00051 {
00052   public:
00053 
00054     typedef _MatrixType MatrixType;
00055     enum {
00056       RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00057       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00058       Options = MatrixType::Options,
00059       MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00060       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00061     };
00062     typedef typename MatrixType::Scalar Scalar;
00063     typedef typename MatrixType::RealScalar RealScalar;
00064     typedef typename MatrixType::Index Index;
00065     typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, Options, MaxRowsAtCompileTime, MaxRowsAtCompileTime> MatrixQType;
00066     typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
00067     typedef Matrix<Index, 1, ColsAtCompileTime, RowMajor, 1, MaxColsAtCompileTime> IntRowVectorType;
00068     typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationType;
00069     typedef typename internal::plain_col_type<MatrixType, Index>::type IntColVectorType;
00070     typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
00071     typedef typename internal::plain_col_type<MatrixType>::type ColVectorType;
00072 
00078     FullPivHouseholderQR()
00079       : m_qr(),
00080         m_hCoeffs(),
00081         m_rows_transpositions(),
00082         m_cols_transpositions(),
00083         m_cols_permutation(),
00084         m_temp(),
00085         m_isInitialized(false),
00086         m_usePrescribedThreshold(false) {}
00087 
00094     FullPivHouseholderQR(Index rows, Index cols)
00095       : m_qr(rows, cols),
00096         m_hCoeffs((std::min)(rows,cols)),
00097         m_rows_transpositions(rows),
00098         m_cols_transpositions(cols),
00099         m_cols_permutation(cols),
00100         m_temp((std::min)(rows,cols)),
00101         m_isInitialized(false),
00102         m_usePrescribedThreshold(false) {}
00103 
00104     FullPivHouseholderQR(const MatrixType& matrix)
00105       : m_qr(matrix.rows(), matrix.cols()),
00106         m_hCoeffs((std::min)(matrix.rows(), matrix.cols())),
00107         m_rows_transpositions(matrix.rows()),
00108         m_cols_transpositions(matrix.cols()),
00109         m_cols_permutation(matrix.cols()),
00110         m_temp((std::min)(matrix.rows(), matrix.cols())),
00111         m_isInitialized(false),
00112         m_usePrescribedThreshold(false)
00113     {
00114       compute(matrix);
00115     }
00116 
00134     template<typename Rhs>
00135     inline const internal::solve_retval<FullPivHouseholderQR, Rhs>
00136     solve(const MatrixBase<Rhs>& b) const
00137     {
00138       eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00139       return internal::solve_retval<FullPivHouseholderQR, Rhs>(*this, b.derived());
00140     }
00141 
00142     MatrixQType matrixQ(void) const;
00143 
00146     const MatrixType& matrixQR() const
00147     {
00148       eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00149       return m_qr;
00150     }
00151 
00152     FullPivHouseholderQR& compute(const MatrixType& matrix);
00153 
00154     const PermutationType& colsPermutation() const
00155     {
00156       eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00157       return m_cols_permutation;
00158     }
00159 
00160     const IntColVectorType& rowsTranspositions() const
00161     {
00162       eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00163       return m_rows_transpositions;
00164     }
00165 
00179     typename MatrixType::RealScalar absDeterminant() const;
00180 
00193     typename MatrixType::RealScalar logAbsDeterminant() const;
00194 
00201     inline Index rank() const
00202     {
00203       eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00204       RealScalar premultiplied_threshold = internal::abs(m_maxpivot) * threshold();
00205       Index result = 0;
00206       for(Index i = 0; i < m_nonzero_pivots; ++i)
00207         result += (internal::abs(m_qr.coeff(i,i)) > premultiplied_threshold);
00208       return result;
00209     }
00210 
00217     inline Index dimensionOfKernel() const
00218     {
00219       eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00220       return cols() - rank();
00221     }
00222 
00230     inline bool isInjective() const
00231     {
00232       eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00233       return rank() == cols();
00234     }
00235 
00243     inline bool isSurjective() const
00244     {
00245       eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00246       return rank() == rows();
00247     }
00248 
00255     inline bool isInvertible() const
00256     {
00257       eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00258       return isInjective() && isSurjective();
00259     }
00260     inline const
00266     internal::solve_retval<FullPivHouseholderQR, typename MatrixType::IdentityReturnType>
00267     inverse() const
00268     {
00269       eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00270       return internal::solve_retval<FullPivHouseholderQR,typename MatrixType::IdentityReturnType>
00271                (*this, MatrixType::Identity(m_qr.rows(), m_qr.cols()));
00272     }
00273 
00274     inline Index rows() const { return m_qr.rows(); }
00275     inline Index cols() const { return m_qr.cols(); }
00276     const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
00277 
00295     FullPivHouseholderQR& setThreshold(const RealScalar& threshold)
00296     {
00297       m_usePrescribedThreshold = true;
00298       m_prescribedThreshold = threshold;
00299       return *this;
00300     }
00301 
00310     FullPivHouseholderQR& setThreshold(Default_t)
00311     {
00312       m_usePrescribedThreshold = false;
00313       return *this;
00314     }
00315 
00320     RealScalar threshold() const
00321     {
00322       eigen_assert(m_isInitialized || m_usePrescribedThreshold);
00323       return m_usePrescribedThreshold ? m_prescribedThreshold
00324       // this formula comes from experimenting (see "LU precision tuning" thread on the list)
00325       // and turns out to be identical to Higham's formula used already in LDLt.
00326                                       : NumTraits<Scalar>::epsilon() * m_qr.diagonalSize();
00327     }
00328 
00336     inline Index nonzeroPivots() const
00337     {
00338       eigen_assert(m_isInitialized && "LU is not initialized.");
00339       return m_nonzero_pivots;
00340     }
00341 
00345     RealScalar maxPivot() const { return m_maxpivot; }
00346 
00347   protected:
00348     MatrixType m_qr;
00349     HCoeffsType m_hCoeffs;
00350     IntColVectorType m_rows_transpositions;
00351     IntRowVectorType m_cols_transpositions;
00352     PermutationType m_cols_permutation;
00353     RowVectorType m_temp;
00354     bool m_isInitialized, m_usePrescribedThreshold;
00355     RealScalar m_prescribedThreshold, m_maxpivot;
00356     Index m_nonzero_pivots;
00357     RealScalar m_precision;
00358     Index m_det_pq;
00359 };
00360 
00361 template<typename MatrixType>
00362 typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType>::absDeterminant() const
00363 {
00364   eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00365   eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
00366   return internal::abs(m_qr.diagonal().prod());
00367 }
00368 
00369 template<typename MatrixType>
00370 typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType>::logAbsDeterminant() const
00371 {
00372   eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00373   eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
00374   return m_qr.diagonal().cwiseAbs().array().log().sum();
00375 }
00376 
00377 template<typename MatrixType>
00378 FullPivHouseholderQR<MatrixType>& FullPivHouseholderQR<MatrixType>::compute(const MatrixType& matrix)
00379 {
00380   Index rows = matrix.rows();
00381   Index cols = matrix.cols();
00382   Index size = (std::min)(rows,cols);
00383 
00384   m_qr = matrix;
00385   m_hCoeffs.resize(size);
00386 
00387   m_temp.resize(cols);
00388 
00389   m_precision = NumTraits<Scalar>::epsilon() * size;
00390 
00391   m_rows_transpositions.resize(matrix.rows());
00392   m_cols_transpositions.resize(matrix.cols());
00393   Index number_of_transpositions = 0;
00394 
00395   RealScalar biggest(0);
00396 
00397   m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
00398   m_maxpivot = RealScalar(0);
00399 
00400   for (Index k = 0; k < size; ++k)
00401   {
00402     Index row_of_biggest_in_corner, col_of_biggest_in_corner;
00403     RealScalar biggest_in_corner;
00404 
00405     biggest_in_corner = m_qr.bottomRightCorner(rows-k, cols-k)
00406                             .cwiseAbs()
00407                             .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner);
00408     row_of_biggest_in_corner += k;
00409     col_of_biggest_in_corner += k;
00410     if(k==0) biggest = biggest_in_corner;
00411 
00412     // if the corner is negligible, then we have less than full rank, and we can finish early
00413     if(internal::isMuchSmallerThan(biggest_in_corner, biggest, m_precision))
00414     {
00415       m_nonzero_pivots = k;
00416       for(Index i = k; i < size; i++)
00417       {
00418         m_rows_transpositions.coeffRef(i) = i;
00419         m_cols_transpositions.coeffRef(i) = i;
00420         m_hCoeffs.coeffRef(i) = Scalar(0);
00421       }
00422       break;
00423     }
00424 
00425     m_rows_transpositions.coeffRef(k) = row_of_biggest_in_corner;
00426     m_cols_transpositions.coeffRef(k) = col_of_biggest_in_corner;
00427     if(k != row_of_biggest_in_corner) {
00428       m_qr.row(k).tail(cols-k).swap(m_qr.row(row_of_biggest_in_corner).tail(cols-k));
00429       ++number_of_transpositions;
00430     }
00431     if(k != col_of_biggest_in_corner) {
00432       m_qr.col(k).swap(m_qr.col(col_of_biggest_in_corner));
00433       ++number_of_transpositions;
00434     }
00435 
00436     RealScalar beta;
00437     m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta);
00438     m_qr.coeffRef(k,k) = beta;
00439 
00440     // remember the maximum absolute value of diagonal coefficients
00441     if(internal::abs(beta) > m_maxpivot) m_maxpivot = internal::abs(beta);
00442 
00443     m_qr.bottomRightCorner(rows-k, cols-k-1)
00444         .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1));
00445   }
00446 
00447   m_cols_permutation.setIdentity(cols);
00448   for(Index k = 0; k < size; ++k)
00449     m_cols_permutation.applyTranspositionOnTheRight(k, m_cols_transpositions.coeff(k));
00450 
00451   m_det_pq = (number_of_transpositions%2) ? -1 : 1;
00452   m_isInitialized = true;
00453 
00454   return *this;
00455 }
00456 
00457 namespace internal {
00458 
00459 template<typename _MatrixType, typename Rhs>
00460 struct solve_retval<FullPivHouseholderQR<_MatrixType>, Rhs>
00461   : solve_retval_base<FullPivHouseholderQR<_MatrixType>, Rhs>
00462 {
00463   EIGEN_MAKE_SOLVE_HELPERS(FullPivHouseholderQR<_MatrixType>,Rhs)
00464 
00465   template<typename Dest> void evalTo(Dest& dst) const
00466   {
00467     const Index rows = dec().rows(), cols = dec().cols();
00468     eigen_assert(rhs().rows() == rows);
00469 
00470     // FIXME introduce nonzeroPivots() and use it here. and more generally,
00471     // make the same improvements in this dec as in FullPivLU.
00472     if(dec().rank()==0)
00473     {
00474       dst.setZero();
00475       return;
00476     }
00477 
00478     typename Rhs::PlainObject c(rhs());
00479 
00480     Matrix<Scalar,1,Rhs::ColsAtCompileTime> temp(rhs().cols());
00481     for (Index k = 0; k < dec().rank(); ++k)
00482     {
00483       Index remainingSize = rows-k;
00484       c.row(k).swap(c.row(dec().rowsTranspositions().coeff(k)));
00485       c.bottomRightCorner(remainingSize, rhs().cols())
00486        .applyHouseholderOnTheLeft(dec().matrixQR().col(k).tail(remainingSize-1),
00487                                   dec().hCoeffs().coeff(k), &temp.coeffRef(0));
00488     }
00489 
00490     if(!dec().isSurjective())
00491     {
00492       // is c is in the image of R ?
00493       RealScalar biggest_in_upper_part_of_c = c.topRows(   dec().rank()     ).cwiseAbs().maxCoeff();
00494       RealScalar biggest_in_lower_part_of_c = c.bottomRows(rows-dec().rank()).cwiseAbs().maxCoeff();
00495       // FIXME brain dead
00496       const RealScalar m_precision = NumTraits<Scalar>::epsilon() * (std::min)(rows,cols);
00497       // this internal:: prefix is needed by at least gcc 3.4 and ICC
00498       if(!internal::isMuchSmallerThan(biggest_in_lower_part_of_c, biggest_in_upper_part_of_c, m_precision))
00499         return;
00500     }
00501     dec().matrixQR()
00502        .topLeftCorner(dec().rank(), dec().rank())
00503        .template triangularView<Upper>()
00504        .solveInPlace(c.topRows(dec().rank()));
00505 
00506     for(Index i = 0; i < dec().rank(); ++i) dst.row(dec().colsPermutation().indices().coeff(i)) = c.row(i);
00507     for(Index i = dec().rank(); i < cols; ++i) dst.row(dec().colsPermutation().indices().coeff(i)).setZero();
00508   }
00509 };
00510 
00511 } // end namespace internal
00512 
00514 template<typename MatrixType>
00515 typename FullPivHouseholderQR<MatrixType>::MatrixQType FullPivHouseholderQR<MatrixType>::matrixQ() const
00516 {
00517   eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00518   // compute the product H'_0 H'_1 ... H'_n-1,
00519   // where H_k is the k-th Householder transformation I - h_k v_k v_k'
00520   // and v_k is the k-th Householder vector [1,m_qr(k+1,k), m_qr(k+2,k), ...]
00521   Index rows = m_qr.rows();
00522   Index cols = m_qr.cols();
00523   Index size = (std::min)(rows,cols);
00524   MatrixQType res = MatrixQType::Identity(rows, rows);
00525   Matrix<Scalar,1,MatrixType::RowsAtCompileTime> temp(rows);
00526   for (Index k = size-1; k >= 0; k--)
00527   {
00528     res.block(k, k, rows-k, rows-k)
00529        .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), internal::conj(m_hCoeffs.coeff(k)), &temp.coeffRef(k));
00530     res.row(k).swap(res.row(m_rows_transpositions.coeff(k)));
00531   }
00532   return res;
00533 }
00534 
00539 template<typename Derived>
00540 const FullPivHouseholderQR<typename MatrixBase<Derived>::PlainObject>
00541 MatrixBase<Derived>::fullPivHouseholderQr() const
00542 {
00543   return FullPivHouseholderQR<PlainObject>(eval());
00544 }
00545 
00546 #endif // EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H


libicr
Author(s): Robert Krug
autogenerated on Mon Jan 6 2014 11:32:42