![]() |
Eigen
3.3.3
|
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