Eigen  3.3.3
CompleteOrthogonalDecomposition.h
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2016 Rasmus Munk Larsen <rmlarsen@google.com>
00005 //
00006 // This Source Code Form is subject to the terms of the Mozilla
00007 // Public License v. 2.0. If a copy of the MPL was not distributed
00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00009 
00010 #ifndef EIGEN_COMPLETEORTHOGONALDECOMPOSITION_H
00011 #define EIGEN_COMPLETEORTHOGONALDECOMPOSITION_H
00012 
00013 namespace Eigen {
00014 
00015 namespace internal {
00016 template <typename _MatrixType>
00017 struct traits<CompleteOrthogonalDecomposition<_MatrixType> >
00018     : traits<_MatrixType> {
00019   enum { Flags = 0 };
00020 };
00021 
00022 }  // end namespace internal
00023 
00047 template <typename _MatrixType>
00048 class CompleteOrthogonalDecomposition {
00049  public:
00050   typedef _MatrixType MatrixType;
00051   enum {
00052     RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00053     ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00054     MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00055     MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00056   };
00057   typedef typename MatrixType::Scalar Scalar;
00058   typedef typename MatrixType::RealScalar RealScalar;
00059   typedef typename MatrixType::StorageIndex StorageIndex;
00060   typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
00061   typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime>
00062       PermutationType;
00063   typedef typename internal::plain_row_type<MatrixType, Index>::type
00064       IntRowVectorType;
00065   typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
00066   typedef typename internal::plain_row_type<MatrixType, RealScalar>::type
00067       RealRowVectorType;
00068   typedef HouseholderSequence<
00069       MatrixType, typename internal::remove_all<
00070                       typename HCoeffsType::ConjugateReturnType>::type>
00071       HouseholderSequenceType;
00072   typedef typename MatrixType::PlainObject PlainObject;
00073 
00074  private:
00075   typedef typename PermutationType::Index PermIndexType;
00076 
00077  public:
00085   CompleteOrthogonalDecomposition() : m_cpqr(), m_zCoeffs(), m_temp() {}
00086 
00093   CompleteOrthogonalDecomposition(Index rows, Index cols)
00094       : m_cpqr(rows, cols), m_zCoeffs((std::min)(rows, cols)), m_temp(cols) {}
00095 
00112   template <typename InputType>
00113   explicit CompleteOrthogonalDecomposition(const EigenBase<InputType>& matrix)
00114       : m_cpqr(matrix.rows(), matrix.cols()),
00115         m_zCoeffs((std::min)(matrix.rows(), matrix.cols())),
00116         m_temp(matrix.cols())
00117   {
00118     compute(matrix.derived());
00119   }
00120 
00127   template<typename InputType>
00128   explicit CompleteOrthogonalDecomposition(EigenBase<InputType>& matrix)
00129     : m_cpqr(matrix.derived()),
00130       m_zCoeffs((std::min)(matrix.rows(), matrix.cols())),
00131       m_temp(matrix.cols())
00132   {
00133     computeInPlace();
00134   }
00135 
00136 
00146   template <typename Rhs>
00147   inline const Solve<CompleteOrthogonalDecomposition, Rhs> solve(
00148       const MatrixBase<Rhs>& b) const {
00149     eigen_assert(m_cpqr.m_isInitialized &&
00150                  "CompleteOrthogonalDecomposition is not initialized.");
00151     return Solve<CompleteOrthogonalDecomposition, Rhs>(*this, b.derived());
00152   }
00153 
00154   HouseholderSequenceType householderQ(void) const;
00155   HouseholderSequenceType matrixQ(void) const { return m_cpqr.householderQ(); }
00156 
00159   MatrixType matrixZ() const {
00160     MatrixType Z = MatrixType::Identity(m_cpqr.cols(), m_cpqr.cols());
00161     applyZAdjointOnTheLeftInPlace(Z);
00162     return Z.adjoint();
00163   }
00164 
00168   const MatrixType& matrixQTZ() const { return m_cpqr.matrixQR(); }
00169 
00181   const MatrixType& matrixT() const { return m_cpqr.matrixQR(); }
00182 
00183   template <typename InputType>
00184   CompleteOrthogonalDecomposition& compute(const EigenBase<InputType>& matrix) {
00185     // Compute the column pivoted QR factorization A P = Q R.
00186     m_cpqr.compute(matrix);
00187     computeInPlace();
00188     return *this;
00189   }
00190 
00192   const PermutationType& colsPermutation() const {
00193     return m_cpqr.colsPermutation();
00194   }
00195 
00209   typename MatrixType::RealScalar absDeterminant() const;
00210 
00224   typename MatrixType::RealScalar logAbsDeterminant() const;
00225 
00233   inline Index rank() const { return m_cpqr.rank(); }
00234 
00242   inline Index dimensionOfKernel() const { return m_cpqr.dimensionOfKernel(); }
00243 
00251   inline bool isInjective() const { return m_cpqr.isInjective(); }
00252 
00260   inline bool isSurjective() const { return m_cpqr.isSurjective(); }
00261 
00269   inline bool isInvertible() const { return m_cpqr.isInvertible(); }
00270 
00276   inline const Inverse<CompleteOrthogonalDecomposition> pseudoInverse() const
00277   {
00278     return Inverse<CompleteOrthogonalDecomposition>(*this);
00279   }
00280 
00281   inline Index rows() const { return m_cpqr.rows(); }
00282   inline Index cols() const { return m_cpqr.cols(); }
00283 
00289   inline const HCoeffsType& hCoeffs() const { return m_cpqr.hCoeffs(); }
00290 
00296   const HCoeffsType& zCoeffs() const { return m_zCoeffs; }
00297 
00317   CompleteOrthogonalDecomposition& setThreshold(const RealScalar& threshold) {
00318     m_cpqr.setThreshold(threshold);
00319     return *this;
00320   }
00321 
00330   CompleteOrthogonalDecomposition& setThreshold(Default_t) {
00331     m_cpqr.setThreshold(Default);
00332     return *this;
00333   }
00334 
00339   RealScalar threshold() const { return m_cpqr.threshold(); }
00340 
00348   inline Index nonzeroPivots() const { return m_cpqr.nonzeroPivots(); }
00349 
00353   inline RealScalar maxPivot() const { return m_cpqr.maxPivot(); }
00354 
00363   ComputationInfo info() const {
00364     eigen_assert(m_cpqr.m_isInitialized && "Decomposition is not initialized.");
00365     return Success;
00366   }
00367 
00368 #ifndef EIGEN_PARSED_BY_DOXYGEN
00369   template <typename RhsType, typename DstType>
00370   EIGEN_DEVICE_FUNC void _solve_impl(const RhsType& rhs, DstType& dst) const;
00371 #endif
00372 
00373  protected:
00374   static void check_template_parameters() {
00375     EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
00376   }
00377 
00378   void computeInPlace();
00379 
00382   template <typename Rhs>
00383   void applyZAdjointOnTheLeftInPlace(Rhs& rhs) const;
00384 
00385   ColPivHouseholderQR<MatrixType> m_cpqr;
00386   HCoeffsType m_zCoeffs;
00387   RowVectorType m_temp;
00388 };
00389 
00390 template <typename MatrixType>
00391 typename MatrixType::RealScalar
00392 CompleteOrthogonalDecomposition<MatrixType>::absDeterminant() const {
00393   return m_cpqr.absDeterminant();
00394 }
00395 
00396 template <typename MatrixType>
00397 typename MatrixType::RealScalar
00398 CompleteOrthogonalDecomposition<MatrixType>::logAbsDeterminant() const {
00399   return m_cpqr.logAbsDeterminant();
00400 }
00401 
00409 template <typename MatrixType>
00410 void CompleteOrthogonalDecomposition<MatrixType>::computeInPlace()
00411 {
00412   check_template_parameters();
00413 
00414   // the column permutation is stored as int indices, so just to be sure:
00415   eigen_assert(m_cpqr.cols() <= NumTraits<int>::highest());
00416 
00417   const Index rank = m_cpqr.rank();
00418   const Index cols = m_cpqr.cols();
00419   const Index rows = m_cpqr.rows();
00420   m_zCoeffs.resize((std::min)(rows, cols));
00421   m_temp.resize(cols);
00422 
00423   if (rank < cols) {
00424     // We have reduced the (permuted) matrix to the form
00425     //   [R11 R12]
00426     //   [ 0  R22]
00427     // where R11 is r-by-r (r = rank) upper triangular, R12 is
00428     // r-by-(n-r), and R22 is empty or the norm of R22 is negligible.
00429     // We now compute the complete orthogonal decomposition by applying
00430     // Householder transformations from the right to the upper trapezoidal
00431     // matrix X = [R11 R12] to zero out R12 and obtain the factorization
00432     // [R11 R12] = [T11 0] * Z, where T11 is r-by-r upper triangular and
00433     // Z = Z(0) * Z(1) ... Z(r-1) is an n-by-n orthogonal matrix.
00434     // We store the data representing Z in R12 and m_zCoeffs.
00435     for (Index k = rank - 1; k >= 0; --k) {
00436       if (k != rank - 1) {
00437         // Given the API for Householder reflectors, it is more convenient if
00438         // we swap the leading parts of columns k and r-1 (zero-based) to form
00439         // the matrix X_k = [X(0:k, k), X(0:k, r:n)]
00440         m_cpqr.m_qr.col(k).head(k + 1).swap(
00441             m_cpqr.m_qr.col(rank - 1).head(k + 1));
00442       }
00443       // Construct Householder reflector Z(k) to zero out the last row of X_k,
00444       // i.e. choose Z(k) such that
00445       // [X(k, k), X(k, r:n)] * Z(k) = [beta, 0, .., 0].
00446       RealScalar beta;
00447       m_cpqr.m_qr.row(k)
00448           .tail(cols - rank + 1)
00449           .makeHouseholderInPlace(m_zCoeffs(k), beta);
00450       m_cpqr.m_qr(k, rank - 1) = beta;
00451       if (k > 0) {
00452         // Apply Z(k) to the first k rows of X_k
00453         m_cpqr.m_qr.topRightCorner(k, cols - rank + 1)
00454             .applyHouseholderOnTheRight(
00455                 m_cpqr.m_qr.row(k).tail(cols - rank).transpose(), m_zCoeffs(k),
00456                 &m_temp(0));
00457       }
00458       if (k != rank - 1) {
00459         // Swap X(0:k,k) back to its proper location.
00460         m_cpqr.m_qr.col(k).head(k + 1).swap(
00461             m_cpqr.m_qr.col(rank - 1).head(k + 1));
00462       }
00463     }
00464   }
00465 }
00466 
00467 template <typename MatrixType>
00468 template <typename Rhs>
00469 void CompleteOrthogonalDecomposition<MatrixType>::applyZAdjointOnTheLeftInPlace(
00470     Rhs& rhs) const {
00471   const Index cols = this->cols();
00472   const Index nrhs = rhs.cols();
00473   const Index rank = this->rank();
00474   Matrix<typename MatrixType::Scalar, Dynamic, 1> temp((std::max)(cols, nrhs));
00475   for (Index k = 0; k < rank; ++k) {
00476     if (k != rank - 1) {
00477       rhs.row(k).swap(rhs.row(rank - 1));
00478     }
00479     rhs.middleRows(rank - 1, cols - rank + 1)
00480         .applyHouseholderOnTheLeft(
00481             matrixQTZ().row(k).tail(cols - rank).adjoint(), zCoeffs()(k),
00482             &temp(0));
00483     if (k != rank - 1) {
00484       rhs.row(k).swap(rhs.row(rank - 1));
00485     }
00486   }
00487 }
00488 
00489 #ifndef EIGEN_PARSED_BY_DOXYGEN
00490 template <typename _MatrixType>
00491 template <typename RhsType, typename DstType>
00492 void CompleteOrthogonalDecomposition<_MatrixType>::_solve_impl(
00493     const RhsType& rhs, DstType& dst) const {
00494   eigen_assert(rhs.rows() == this->rows());
00495 
00496   const Index rank = this->rank();
00497   if (rank == 0) {
00498     dst.setZero();
00499     return;
00500   }
00501 
00502   // Compute c = Q^* * rhs
00503   // Note that the matrix Q = H_0^* H_1^*... so its inverse is
00504   // Q^* = (H_0 H_1 ...)^T
00505   typename RhsType::PlainObject c(rhs);
00506   c.applyOnTheLeft(
00507       householderSequence(matrixQTZ(), hCoeffs()).setLength(rank).transpose());
00508 
00509   // Solve T z = c(1:rank, :)
00510   dst.topRows(rank) = matrixT()
00511                           .topLeftCorner(rank, rank)
00512                           .template triangularView<Upper>()
00513                           .solve(c.topRows(rank));
00514 
00515   const Index cols = this->cols();
00516   if (rank < cols) {
00517     // Compute y = Z^* * [ z ]
00518     //                   [ 0 ]
00519     dst.bottomRows(cols - rank).setZero();
00520     applyZAdjointOnTheLeftInPlace(dst);
00521   }
00522 
00523   // Undo permutation to get x = P^{-1} * y.
00524   dst = colsPermutation() * dst;
00525 }
00526 #endif
00527 
00528 namespace internal {
00529 
00530 template<typename DstXprType, typename MatrixType>
00531 struct Assignment<DstXprType, Inverse<CompleteOrthogonalDecomposition<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename CompleteOrthogonalDecomposition<MatrixType>::Scalar>, Dense2Dense>
00532 {
00533   typedef CompleteOrthogonalDecomposition<MatrixType> CodType;
00534   typedef Inverse<CodType> SrcXprType;
00535   static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename CodType::Scalar> &)
00536   {
00537     dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.rows()));
00538   }
00539 };
00540 
00541 } // end namespace internal
00542 
00544 template <typename MatrixType>
00545 typename CompleteOrthogonalDecomposition<MatrixType>::HouseholderSequenceType
00546 CompleteOrthogonalDecomposition<MatrixType>::householderQ() const {
00547   return m_cpqr.householderQ();
00548 }
00549 
00554 template <typename Derived>
00555 const CompleteOrthogonalDecomposition<typename MatrixBase<Derived>::PlainObject>
00556 MatrixBase<Derived>::completeOrthogonalDecomposition() const {
00557   return CompleteOrthogonalDecomposition<PlainObject>(eval());
00558 }
00559 
00560 }  // end namespace Eigen
00561 
00562 #endif  // EIGEN_COMPLETEORTHOGONALDECOMPOSITION_H
 All Classes Functions Variables Typedefs Enumerations Enumerator Friends