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 // 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 FullPivHouseholderQRMatrixQReturnType;
00019 
00020 template<typename MatrixType>
00021 struct traits<FullPivHouseholderQRMatrixQReturnType<MatrixType> >
00022 {
00023   typedef typename MatrixType::PlainObject ReturnType;
00024 };
00025 
00026 }
00027 
00049 template<typename _MatrixType> class FullPivHouseholderQR
00050 {
00051   public:
00052 
00053     typedef _MatrixType MatrixType;
00054     enum {
00055       RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00056       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00057       Options = MatrixType::Options,
00058       MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00059       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00060     };
00061     typedef typename MatrixType::Scalar Scalar;
00062     typedef typename MatrixType::RealScalar RealScalar;
00063     typedef typename MatrixType::Index Index;
00064     typedef internal::FullPivHouseholderQRMatrixQReturnType<MatrixType> MatrixQReturnType;
00065     typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
00066     typedef Matrix<Index, 1, ColsAtCompileTime, RowMajor, 1, MaxColsAtCompileTime> IntRowVectorType;
00067     typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationType;
00068     typedef typename internal::plain_col_type<MatrixType, Index>::type IntColVectorType;
00069     typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
00070     typedef typename internal::plain_col_type<MatrixType>::type ColVectorType;
00071 
00077     FullPivHouseholderQR()
00078       : m_qr(),
00079         m_hCoeffs(),
00080         m_rows_transpositions(),
00081         m_cols_transpositions(),
00082         m_cols_permutation(),
00083         m_temp(),
00084         m_isInitialized(false),
00085         m_usePrescribedThreshold(false) {}
00086 
00093     FullPivHouseholderQR(Index rows, Index cols)
00094       : m_qr(rows, cols),
00095         m_hCoeffs((std::min)(rows,cols)),
00096         m_rows_transpositions(rows),
00097         m_cols_transpositions(cols),
00098         m_cols_permutation(cols),
00099         m_temp((std::min)(rows,cols)),
00100         m_isInitialized(false),
00101         m_usePrescribedThreshold(false) {}
00102 
00115     FullPivHouseholderQR(const MatrixType& matrix)
00116       : m_qr(matrix.rows(), matrix.cols()),
00117         m_hCoeffs((std::min)(matrix.rows(), matrix.cols())),
00118         m_rows_transpositions(matrix.rows()),
00119         m_cols_transpositions(matrix.cols()),
00120         m_cols_permutation(matrix.cols()),
00121         m_temp((std::min)(matrix.rows(), matrix.cols())),
00122         m_isInitialized(false),
00123         m_usePrescribedThreshold(false)
00124     {
00125       compute(matrix);
00126     }
00127 
00145     template<typename Rhs>
00146     inline const internal::solve_retval<FullPivHouseholderQR, Rhs>
00147     solve(const MatrixBase<Rhs>& b) const
00148     {
00149       eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00150       return internal::solve_retval<FullPivHouseholderQR, Rhs>(*this, b.derived());
00151     }
00152 
00155     MatrixQReturnType matrixQ(void) const;
00156 
00159     const MatrixType& matrixQR() const
00160     {
00161       eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00162       return m_qr;
00163     }
00164 
00165     FullPivHouseholderQR& compute(const MatrixType& matrix);
00166 
00168     const PermutationType& colsPermutation() const
00169     {
00170       eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00171       return m_cols_permutation;
00172     }
00173 
00175     const IntColVectorType& rowsTranspositions() const
00176     {
00177       eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00178       return m_rows_transpositions;
00179     }
00180 
00194     typename MatrixType::RealScalar absDeterminant() const;
00195 
00208     typename MatrixType::RealScalar logAbsDeterminant() const;
00209 
00216     inline Index rank() const
00217     {
00218       using std::abs;
00219       eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00220       RealScalar premultiplied_threshold = abs(m_maxpivot) * threshold();
00221       Index result = 0;
00222       for(Index i = 0; i < m_nonzero_pivots; ++i)
00223         result += (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     
00297     const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
00298 
00316     FullPivHouseholderQR& setThreshold(const RealScalar& threshold)
00317     {
00318       m_usePrescribedThreshold = true;
00319       m_prescribedThreshold = threshold;
00320       return *this;
00321     }
00322 
00331     FullPivHouseholderQR& setThreshold(Default_t)
00332     {
00333       m_usePrescribedThreshold = false;
00334       return *this;
00335     }
00336 
00341     RealScalar threshold() const
00342     {
00343       eigen_assert(m_isInitialized || m_usePrescribedThreshold);
00344       return m_usePrescribedThreshold ? m_prescribedThreshold
00345       // this formula comes from experimenting (see "LU precision tuning" thread on the list)
00346       // and turns out to be identical to Higham's formula used already in LDLt.
00347                                       : NumTraits<Scalar>::epsilon() * m_qr.diagonalSize();
00348     }
00349 
00357     inline Index nonzeroPivots() const
00358     {
00359       eigen_assert(m_isInitialized && "LU is not initialized.");
00360       return m_nonzero_pivots;
00361     }
00362 
00366     RealScalar maxPivot() const { return m_maxpivot; }
00367 
00368   protected:
00369     MatrixType m_qr;
00370     HCoeffsType m_hCoeffs;
00371     IntColVectorType m_rows_transpositions;
00372     IntRowVectorType m_cols_transpositions;
00373     PermutationType m_cols_permutation;
00374     RowVectorType m_temp;
00375     bool m_isInitialized, m_usePrescribedThreshold;
00376     RealScalar m_prescribedThreshold, m_maxpivot;
00377     Index m_nonzero_pivots;
00378     RealScalar m_precision;
00379     Index m_det_pq;
00380 };
00381 
00382 template<typename MatrixType>
00383 typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType>::absDeterminant() const
00384 {
00385   using std::abs;
00386   eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00387   eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
00388   return abs(m_qr.diagonal().prod());
00389 }
00390 
00391 template<typename MatrixType>
00392 typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType>::logAbsDeterminant() const
00393 {
00394   eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00395   eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
00396   return m_qr.diagonal().cwiseAbs().array().log().sum();
00397 }
00398 
00405 template<typename MatrixType>
00406 FullPivHouseholderQR<MatrixType>& FullPivHouseholderQR<MatrixType>::compute(const MatrixType& matrix)
00407 {
00408   using std::abs;
00409   Index rows = matrix.rows();
00410   Index cols = matrix.cols();
00411   Index size = (std::min)(rows,cols);
00412 
00413   m_qr = matrix;
00414   m_hCoeffs.resize(size);
00415 
00416   m_temp.resize(cols);
00417 
00418   m_precision = NumTraits<Scalar>::epsilon() * size;
00419 
00420   m_rows_transpositions.resize(matrix.rows());
00421   m_cols_transpositions.resize(matrix.cols());
00422   Index number_of_transpositions = 0;
00423 
00424   RealScalar biggest(0);
00425 
00426   m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
00427   m_maxpivot = RealScalar(0);
00428 
00429   for (Index k = 0; k < size; ++k)
00430   {
00431     Index row_of_biggest_in_corner, col_of_biggest_in_corner;
00432     RealScalar biggest_in_corner;
00433 
00434     biggest_in_corner = m_qr.bottomRightCorner(rows-k, cols-k)
00435                             .cwiseAbs()
00436                             .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner);
00437     row_of_biggest_in_corner += k;
00438     col_of_biggest_in_corner += k;
00439     if(k==0) biggest = biggest_in_corner;
00440 
00441     // if the corner is negligible, then we have less than full rank, and we can finish early
00442     if(internal::isMuchSmallerThan(biggest_in_corner, biggest, m_precision))
00443     {
00444       m_nonzero_pivots = k;
00445       for(Index i = k; i < size; i++)
00446       {
00447         m_rows_transpositions.coeffRef(i) = i;
00448         m_cols_transpositions.coeffRef(i) = i;
00449         m_hCoeffs.coeffRef(i) = Scalar(0);
00450       }
00451       break;
00452     }
00453 
00454     m_rows_transpositions.coeffRef(k) = row_of_biggest_in_corner;
00455     m_cols_transpositions.coeffRef(k) = col_of_biggest_in_corner;
00456     if(k != row_of_biggest_in_corner) {
00457       m_qr.row(k).tail(cols-k).swap(m_qr.row(row_of_biggest_in_corner).tail(cols-k));
00458       ++number_of_transpositions;
00459     }
00460     if(k != col_of_biggest_in_corner) {
00461       m_qr.col(k).swap(m_qr.col(col_of_biggest_in_corner));
00462       ++number_of_transpositions;
00463     }
00464 
00465     RealScalar beta;
00466     m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta);
00467     m_qr.coeffRef(k,k) = beta;
00468 
00469     // remember the maximum absolute value of diagonal coefficients
00470     if(abs(beta) > m_maxpivot) m_maxpivot = abs(beta);
00471 
00472     m_qr.bottomRightCorner(rows-k, cols-k-1)
00473         .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1));
00474   }
00475 
00476   m_cols_permutation.setIdentity(cols);
00477   for(Index k = 0; k < size; ++k)
00478     m_cols_permutation.applyTranspositionOnTheRight(k, m_cols_transpositions.coeff(k));
00479 
00480   m_det_pq = (number_of_transpositions%2) ? -1 : 1;
00481   m_isInitialized = true;
00482 
00483   return *this;
00484 }
00485 
00486 namespace internal {
00487 
00488 template<typename _MatrixType, typename Rhs>
00489 struct solve_retval<FullPivHouseholderQR<_MatrixType>, Rhs>
00490   : solve_retval_base<FullPivHouseholderQR<_MatrixType>, Rhs>
00491 {
00492   EIGEN_MAKE_SOLVE_HELPERS(FullPivHouseholderQR<_MatrixType>,Rhs)
00493 
00494   template<typename Dest> void evalTo(Dest& dst) const
00495   {
00496     const Index rows = dec().rows(), cols = dec().cols();
00497     eigen_assert(rhs().rows() == rows);
00498 
00499     // FIXME introduce nonzeroPivots() and use it here. and more generally,
00500     // make the same improvements in this dec as in FullPivLU.
00501     if(dec().rank()==0)
00502     {
00503       dst.setZero();
00504       return;
00505     }
00506 
00507     typename Rhs::PlainObject c(rhs());
00508 
00509     Matrix<Scalar,1,Rhs::ColsAtCompileTime> temp(rhs().cols());
00510     for (Index k = 0; k < dec().rank(); ++k)
00511     {
00512       Index remainingSize = rows-k;
00513       c.row(k).swap(c.row(dec().rowsTranspositions().coeff(k)));
00514       c.bottomRightCorner(remainingSize, rhs().cols())
00515        .applyHouseholderOnTheLeft(dec().matrixQR().col(k).tail(remainingSize-1),
00516                                   dec().hCoeffs().coeff(k), &temp.coeffRef(0));
00517     }
00518 
00519     if(!dec().isSurjective())
00520     {
00521       // is c is in the image of R ?
00522       RealScalar biggest_in_upper_part_of_c = c.topRows(   dec().rank()     ).cwiseAbs().maxCoeff();
00523       RealScalar biggest_in_lower_part_of_c = c.bottomRows(rows-dec().rank()).cwiseAbs().maxCoeff();
00524       // FIXME brain dead
00525       const RealScalar m_precision = NumTraits<Scalar>::epsilon() * (std::min)(rows,cols);
00526       // this internal:: prefix is needed by at least gcc 3.4 and ICC
00527       if(!internal::isMuchSmallerThan(biggest_in_lower_part_of_c, biggest_in_upper_part_of_c, m_precision))
00528         return;
00529     }
00530     dec().matrixQR()
00531        .topLeftCorner(dec().rank(), dec().rank())
00532        .template triangularView<Upper>()
00533        .solveInPlace(c.topRows(dec().rank()));
00534 
00535     for(Index i = 0; i < dec().rank(); ++i) dst.row(dec().colsPermutation().indices().coeff(i)) = c.row(i);
00536     for(Index i = dec().rank(); i < cols; ++i) dst.row(dec().colsPermutation().indices().coeff(i)).setZero();
00537   }
00538 };
00539 
00546 template<typename MatrixType> struct FullPivHouseholderQRMatrixQReturnType
00547   : public ReturnByValue<FullPivHouseholderQRMatrixQReturnType<MatrixType> >
00548 {
00549 public:
00550   typedef typename MatrixType::Index Index;
00551   typedef typename internal::plain_col_type<MatrixType, Index>::type IntColVectorType;
00552   typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
00553   typedef Matrix<typename MatrixType::Scalar, 1, MatrixType::RowsAtCompileTime, RowMajor, 1,
00554                  MatrixType::MaxRowsAtCompileTime> WorkVectorType;
00555 
00556   FullPivHouseholderQRMatrixQReturnType(const MatrixType&       qr,
00557                                         const HCoeffsType&      hCoeffs,
00558                                         const IntColVectorType& rowsTranspositions)
00559     : m_qr(qr),
00560       m_hCoeffs(hCoeffs),
00561       m_rowsTranspositions(rowsTranspositions)
00562       {}
00563 
00564   template <typename ResultType>
00565   void evalTo(ResultType& result) const
00566   {
00567     const Index rows = m_qr.rows();
00568     WorkVectorType workspace(rows);
00569     evalTo(result, workspace);
00570   }
00571 
00572   template <typename ResultType>
00573   void evalTo(ResultType& result, WorkVectorType& workspace) const
00574   {
00575     using numext::conj;
00576     // compute the product H'_0 H'_1 ... H'_n-1,
00577     // where H_k is the k-th Householder transformation I - h_k v_k v_k'
00578     // and v_k is the k-th Householder vector [1,m_qr(k+1,k), m_qr(k+2,k), ...]
00579     const Index rows = m_qr.rows();
00580     const Index cols = m_qr.cols();
00581     const Index size = (std::min)(rows, cols);
00582     workspace.resize(rows);
00583     result.setIdentity(rows, rows);
00584     for (Index k = size-1; k >= 0; k--)
00585     {
00586       result.block(k, k, rows-k, rows-k)
00587             .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), conj(m_hCoeffs.coeff(k)), &workspace.coeffRef(k));
00588       result.row(k).swap(result.row(m_rowsTranspositions.coeff(k)));
00589     }
00590   }
00591 
00592     Index rows() const { return m_qr.rows(); }
00593     Index cols() const { return m_qr.rows(); }
00594 
00595 protected:
00596   typename MatrixType::Nested m_qr;
00597   typename HCoeffsType::Nested m_hCoeffs;
00598   typename IntColVectorType::Nested m_rowsTranspositions;
00599 };
00600 
00601 } // end namespace internal
00602 
00603 template<typename MatrixType>
00604 inline typename FullPivHouseholderQR<MatrixType>::MatrixQReturnType FullPivHouseholderQR<MatrixType>::matrixQ() const
00605 {
00606   eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00607   return MatrixQReturnType(m_qr, m_hCoeffs, m_rows_transpositions);
00608 }
00609 
00614 template<typename Derived>
00615 const FullPivHouseholderQR<typename MatrixBase<Derived>::PlainObject>
00616 MatrixBase<Derived>::fullPivHouseholderQr() const
00617 {
00618   return FullPivHouseholderQR<PlainObject>(eval());
00619 }
00620 
00621 } // end namespace Eigen
00622 
00623 #endif // EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Sat Jun 8 2019 19:37:05