Eigen  3.4.90 (git rev 67eeba6e720c5745abc77ae6c92ce0a44aa7b7ae)
CompleteOrthogonalDecomposition.h
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2016 Rasmus Munk Larsen <rmlarsen@google.com>
5 //
6 // This Source Code Form is subject to the terms of the Mozilla
7 // Public License v. 2.0. If a copy of the MPL was not distributed
8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9 
10 #ifndef EIGEN_COMPLETEORTHOGONALDECOMPOSITION_H
11 #define EIGEN_COMPLETEORTHOGONALDECOMPOSITION_H
12 
13 #include "./InternalHeaderCheck.h"
14 
15 namespace Eigen {
16 
17 namespace internal {
18 template <typename MatrixType_>
19 struct traits<CompleteOrthogonalDecomposition<MatrixType_> >
20  : traits<MatrixType_> {
21  typedef MatrixXpr XprKind;
22  typedef SolverStorage StorageKind;
23  typedef int StorageIndex;
24  enum { Flags = 0 };
25 };
26 
27 } // end namespace internal
28 
52 template <typename MatrixType_> class CompleteOrthogonalDecomposition
53  : public SolverBase<CompleteOrthogonalDecomposition<MatrixType_> >
54 {
55  public:
56  typedef MatrixType_ MatrixType;
58 
59  template<typename Derived>
60  friend struct internal::solve_assertion;
61 
62  EIGEN_GENERIC_PUBLIC_INTERFACE(CompleteOrthogonalDecomposition)
63  enum {
64  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
65  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
66  };
67  typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
70  typedef typename internal::plain_row_type<MatrixType, Index>::type
71  IntRowVectorType;
72  typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
73  typedef typename internal::plain_row_type<MatrixType, RealScalar>::type
74  RealRowVectorType;
75  typedef HouseholderSequence<
76  MatrixType, internal::remove_all_t<
77  typename HCoeffsType::ConjugateReturnType>>
79  typedef typename MatrixType::PlainObject PlainObject;
80 
81  private:
82  typedef typename PermutationType::Index PermIndexType;
83 
84  public:
92  CompleteOrthogonalDecomposition() : m_cpqr(), m_zCoeffs(), m_temp() {}
93 
101  : m_cpqr(rows, cols), m_zCoeffs((std::min)(rows, cols)), m_temp(cols) {}
102 
119  template <typename InputType>
121  : m_cpqr(matrix.rows(), matrix.cols()),
122  m_zCoeffs((std::min)(matrix.rows(), matrix.cols())),
123  m_temp(matrix.cols())
124  {
125  compute(matrix.derived());
126  }
127 
134  template<typename InputType>
136  : m_cpqr(matrix.derived()),
137  m_zCoeffs((std::min)(matrix.rows(), matrix.cols())),
138  m_temp(matrix.cols())
139  {
140  computeInPlace();
141  }
142 
143  #ifdef EIGEN_PARSED_BY_DOXYGEN
153  template <typename Rhs>
155  const MatrixBase<Rhs>& b) const;
156  #endif
157 
159  HouseholderSequenceType matrixQ(void) const { return m_cpqr.householderQ(); }
160 
163  MatrixType matrixZ() const {
164  MatrixType Z = MatrixType::Identity(m_cpqr.cols(), m_cpqr.cols());
165  applyZOnTheLeftInPlace<false>(Z);
166  return Z;
167  }
168 
172  const MatrixType& matrixQTZ() const { return m_cpqr.matrixQR(); }
173 
185  const MatrixType& matrixT() const { return m_cpqr.matrixQR(); }
186 
187  template <typename InputType>
188  CompleteOrthogonalDecomposition& compute(const EigenBase<InputType>& matrix) {
189  // Compute the column pivoted QR factorization A P = Q R.
190  m_cpqr.compute(matrix);
191  computeInPlace();
192  return *this;
193  }
194 
197  return m_cpqr.colsPermutation();
198  }
199 
213  typename MatrixType::RealScalar absDeterminant() const;
214 
228  typename MatrixType::RealScalar logAbsDeterminant() const;
229 
237  inline Index rank() const { return m_cpqr.rank(); }
238 
246  inline Index dimensionOfKernel() const { return m_cpqr.dimensionOfKernel(); }
247 
255  inline bool isInjective() const { return m_cpqr.isInjective(); }
256 
264  inline bool isSurjective() const { return m_cpqr.isSurjective(); }
265 
273  inline bool isInvertible() const { return m_cpqr.isInvertible(); }
274 
281  {
282  eigen_assert(m_cpqr.m_isInitialized && "CompleteOrthogonalDecomposition is not initialized.");
284  }
285 
286  inline Index rows() const { return m_cpqr.rows(); }
287  inline Index cols() const { return m_cpqr.cols(); }
288 
294  inline const HCoeffsType& hCoeffs() const { return m_cpqr.hCoeffs(); }
295 
301  const HCoeffsType& zCoeffs() const { return m_zCoeffs; }
302 
323  m_cpqr.setThreshold(threshold);
324  return *this;
325  }
326 
336  m_cpqr.setThreshold(Default);
337  return *this;
338  }
339 
344  RealScalar threshold() const { return m_cpqr.threshold(); }
345 
353  inline Index nonzeroPivots() const { return m_cpqr.nonzeroPivots(); }
354 
358  inline RealScalar maxPivot() const { return m_cpqr.maxPivot(); }
359 
369  eigen_assert(m_cpqr.m_isInitialized && "Decomposition is not initialized.");
370  return Success;
371  }
372 
373 #ifndef EIGEN_PARSED_BY_DOXYGEN
374  template <typename RhsType, typename DstType>
375  void _solve_impl(const RhsType& rhs, DstType& dst) const;
376 
377  template<bool Conjugate, typename RhsType, typename DstType>
378  void _solve_impl_transposed(const RhsType &rhs, DstType &dst) const;
379 #endif
380 
381  protected:
382  EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar)
383 
384  template<bool Transpose_, typename Rhs>
385  void _check_solve_assertion(const Rhs& b) const {
386  EIGEN_ONLY_USED_FOR_DEBUG(b);
387  eigen_assert(m_cpqr.m_isInitialized && "CompleteOrthogonalDecomposition is not initialized.");
388  eigen_assert((Transpose_?derived().cols():derived().rows())==b.rows() && "CompleteOrthogonalDecomposition::solve(): invalid number of rows of the right hand side matrix b");
389  }
390 
391  void computeInPlace();
392 
397  template <bool Conjugate, typename Rhs>
398  void applyZOnTheLeftInPlace(Rhs& rhs) const;
399 
402  template <typename Rhs>
403  void applyZAdjointOnTheLeftInPlace(Rhs& rhs) const;
404 
405  ColPivHouseholderQR<MatrixType> m_cpqr;
406  HCoeffsType m_zCoeffs;
407  RowVectorType m_temp;
408 };
409 
410 template <typename MatrixType>
411 typename MatrixType::RealScalar
413  return m_cpqr.absDeterminant();
414 }
415 
416 template <typename MatrixType>
417 typename MatrixType::RealScalar
419  return m_cpqr.logAbsDeterminant();
420 }
421 
429 template <typename MatrixType>
431 {
432  // the column permutation is stored as int indices, so just to be sure:
433  eigen_assert(m_cpqr.cols() <= NumTraits<int>::highest());
434 
435  const Index rank = m_cpqr.rank();
436  const Index cols = m_cpqr.cols();
437  const Index rows = m_cpqr.rows();
438  m_zCoeffs.resize((std::min)(rows, cols));
439  m_temp.resize(cols);
440 
441  if (rank < cols) {
442  // We have reduced the (permuted) matrix to the form
443  // [R11 R12]
444  // [ 0 R22]
445  // where R11 is r-by-r (r = rank) upper triangular, R12 is
446  // r-by-(n-r), and R22 is empty or the norm of R22 is negligible.
447  // We now compute the complete orthogonal decomposition by applying
448  // Householder transformations from the right to the upper trapezoidal
449  // matrix X = [R11 R12] to zero out R12 and obtain the factorization
450  // [R11 R12] = [T11 0] * Z, where T11 is r-by-r upper triangular and
451  // Z = Z(0) * Z(1) ... Z(r-1) is an n-by-n orthogonal matrix.
452  // We store the data representing Z in R12 and m_zCoeffs.
453  for (Index k = rank - 1; k >= 0; --k) {
454  if (k != rank - 1) {
455  // Given the API for Householder reflectors, it is more convenient if
456  // we swap the leading parts of columns k and r-1 (zero-based) to form
457  // the matrix X_k = [X(0:k, k), X(0:k, r:n)]
458  m_cpqr.m_qr.col(k).head(k + 1).swap(
459  m_cpqr.m_qr.col(rank - 1).head(k + 1));
460  }
461  // Construct Householder reflector Z(k) to zero out the last row of X_k,
462  // i.e. choose Z(k) such that
463  // [X(k, k), X(k, r:n)] * Z(k) = [beta, 0, .., 0].
464  RealScalar beta;
465  m_cpqr.m_qr.row(k)
466  .tail(cols - rank + 1)
467  .makeHouseholderInPlace(m_zCoeffs(k), beta);
468  m_cpqr.m_qr(k, rank - 1) = beta;
469  if (k > 0) {
470  // Apply Z(k) to the first k rows of X_k
471  m_cpqr.m_qr.topRightCorner(k, cols - rank + 1)
472  .applyHouseholderOnTheRight(
473  m_cpqr.m_qr.row(k).tail(cols - rank).adjoint(), m_zCoeffs(k),
474  &m_temp(0));
475  }
476  if (k != rank - 1) {
477  // Swap X(0:k,k) back to its proper location.
478  m_cpqr.m_qr.col(k).head(k + 1).swap(
479  m_cpqr.m_qr.col(rank - 1).head(k + 1));
480  }
481  }
482  }
483 }
484 
485 template <typename MatrixType>
486 template <bool Conjugate, typename Rhs>
488  Rhs& rhs) const {
489  const Index cols = this->cols();
490  const Index nrhs = rhs.cols();
491  const Index rank = this->rank();
492  Matrix<typename Rhs::Scalar, Dynamic, 1> temp((std::max)(cols, nrhs));
493  for (Index k = rank-1; k >= 0; --k) {
494  if (k != rank - 1) {
495  rhs.row(k).swap(rhs.row(rank - 1));
496  }
497  rhs.middleRows(rank - 1, cols - rank + 1)
498  .applyHouseholderOnTheLeft(
499  matrixQTZ().row(k).tail(cols - rank).transpose().template conjugateIf<!Conjugate>(), zCoeffs().template conjugateIf<Conjugate>()(k),
500  &temp(0));
501  if (k != rank - 1) {
502  rhs.row(k).swap(rhs.row(rank - 1));
503  }
504  }
505 }
506 
507 template <typename MatrixType>
508 template <typename Rhs>
510  Rhs& rhs) const {
511  const Index cols = this->cols();
512  const Index nrhs = rhs.cols();
513  const Index rank = this->rank();
514  Matrix<typename Rhs::Scalar, Dynamic, 1> temp((std::max)(cols, nrhs));
515  for (Index k = 0; k < rank; ++k) {
516  if (k != rank - 1) {
517  rhs.row(k).swap(rhs.row(rank - 1));
518  }
519  rhs.middleRows(rank - 1, cols - rank + 1)
520  .applyHouseholderOnTheLeft(
521  matrixQTZ().row(k).tail(cols - rank).adjoint(), zCoeffs()(k),
522  &temp(0));
523  if (k != rank - 1) {
524  rhs.row(k).swap(rhs.row(rank - 1));
525  }
526  }
527 }
528 
529 #ifndef EIGEN_PARSED_BY_DOXYGEN
530 template <typename MatrixType_>
531 template <typename RhsType, typename DstType>
533  const RhsType& rhs, DstType& dst) const {
534  const Index rank = this->rank();
535  if (rank == 0) {
536  dst.setZero();
537  return;
538  }
539 
540  // Compute c = Q^* * rhs
541  typename RhsType::PlainObject c(rhs);
542  c.applyOnTheLeft(matrixQ().setLength(rank).adjoint());
543 
544  // Solve T z = c(1:rank, :)
545  dst.topRows(rank) = matrixT()
546  .topLeftCorner(rank, rank)
547  .template triangularView<Upper>()
548  .solve(c.topRows(rank));
549 
550  const Index cols = this->cols();
551  if (rank < cols) {
552  // Compute y = Z^* * [ z ]
553  // [ 0 ]
554  dst.bottomRows(cols - rank).setZero();
555  applyZAdjointOnTheLeftInPlace(dst);
556  }
557 
558  // Undo permutation to get x = P^{-1} * y.
559  dst = colsPermutation() * dst;
560 }
561 
562 template<typename MatrixType_>
563 template<bool Conjugate, typename RhsType, typename DstType>
564 void CompleteOrthogonalDecomposition<MatrixType_>::_solve_impl_transposed(const RhsType &rhs, DstType &dst) const
565 {
566  const Index rank = this->rank();
567 
568  if (rank == 0) {
569  dst.setZero();
570  return;
571  }
572 
573  typename RhsType::PlainObject c(colsPermutation().transpose()*rhs);
574 
575  if (rank < cols()) {
576  applyZOnTheLeftInPlace<!Conjugate>(c);
577  }
578 
579  matrixT().topLeftCorner(rank, rank)
580  .template triangularView<Upper>()
581  .transpose().template conjugateIf<Conjugate>()
582  .solveInPlace(c.topRows(rank));
583 
584  dst.topRows(rank) = c.topRows(rank);
585  dst.bottomRows(rows()-rank).setZero();
586 
587  dst.applyOnTheLeft(householderQ().setLength(rank).template conjugateIf<!Conjugate>() );
588 }
589 #endif
590 
591 namespace internal {
592 
593 template<typename MatrixType>
594 struct traits<Inverse<CompleteOrthogonalDecomposition<MatrixType> > >
595  : traits<typename Transpose<typename MatrixType::PlainObject>::PlainObject>
596 {
597  enum { Flags = 0 };
598 };
599 
600 template<typename DstXprType, typename MatrixType>
601 struct Assignment<DstXprType, Inverse<CompleteOrthogonalDecomposition<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename CompleteOrthogonalDecomposition<MatrixType>::Scalar>, Dense2Dense>
602 {
603  typedef CompleteOrthogonalDecomposition<MatrixType> CodType;
604  typedef Inverse<CodType> SrcXprType;
605  static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename CodType::Scalar> &)
606  {
607  typedef Matrix<typename CodType::Scalar, CodType::RowsAtCompileTime, CodType::RowsAtCompileTime, 0, CodType::MaxRowsAtCompileTime, CodType::MaxRowsAtCompileTime> IdentityMatrixType;
608  dst = src.nestedExpression().solve(IdentityMatrixType::Identity(src.cols(), src.cols()));
609  }
610 };
611 
612 } // end namespace internal
613 
615 template <typename MatrixType>
616 typename CompleteOrthogonalDecomposition<MatrixType>::HouseholderSequenceType
618  return m_cpqr.householderQ();
619 }
620 
625 template <typename Derived>
629 }
630 
631 } // end namespace Eigen
632 
633 #endif // EIGEN_COMPLETEORTHOGONALDECOMPOSITION_H
const HCoeffsType & hCoeffs() const
Definition: ColPivHouseholderQR.h:336
Index nonzeroPivots() const
Definition: ColPivHouseholderQR.h:396
bool isInvertible() const
Definition: ColPivHouseholderQR.h:312
HouseholderSequenceType householderQ() const
Definition: ColPivHouseholderQR.h:652
const MatrixType & matrixQR() const
Definition: ColPivHouseholderQR.h:191
Index rank() const
Definition: ColPivHouseholderQR.h:257
RealScalar threshold() const
Definition: ColPivHouseholderQR.h:380
const PermutationType & colsPermutation() const
Definition: ColPivHouseholderQR.h:216
RealScalar maxPivot() const
Definition: ColPivHouseholderQR.h:405
bool isSurjective() const
Definition: ColPivHouseholderQR.h:300
Index dimensionOfKernel() const
Definition: ColPivHouseholderQR.h:274
ColPivHouseholderQR & setThreshold(const RealScalar &threshold)
Definition: ColPivHouseholderQR.h:355
bool isInjective() const
Definition: ColPivHouseholderQR.h:287
Complete orthogonal decomposition (COD) of a matrix.
Definition: CompleteOrthogonalDecomposition.h:54
void applyZAdjointOnTheLeftInPlace(Rhs &rhs) const
Definition: CompleteOrthogonalDecomposition.h:509
Index rank() const
Definition: CompleteOrthogonalDecomposition.h:237
RealScalar threshold() const
Definition: CompleteOrthogonalDecomposition.h:344
const HCoeffsType & zCoeffs() const
Definition: CompleteOrthogonalDecomposition.h:301
Index nonzeroPivots() const
Definition: CompleteOrthogonalDecomposition.h:353
const PermutationType & colsPermutation() const
Definition: CompleteOrthogonalDecomposition.h:196
bool isInvertible() const
Definition: CompleteOrthogonalDecomposition.h:273
ComputationInfo info() const
Reports whether the complete orthogonal decomposition was successful.
Definition: CompleteOrthogonalDecomposition.h:368
bool isInjective() const
Definition: CompleteOrthogonalDecomposition.h:255
CompleteOrthogonalDecomposition & setThreshold(Default_t)
Definition: CompleteOrthogonalDecomposition.h:335
MatrixType matrixZ() const
Definition: CompleteOrthogonalDecomposition.h:163
bool isSurjective() const
Definition: CompleteOrthogonalDecomposition.h:264
const Solve< CompleteOrthogonalDecomposition, Rhs > solve(const MatrixBase< Rhs > &b) const
CompleteOrthogonalDecomposition()
Default Constructor.
Definition: CompleteOrthogonalDecomposition.h:92
const MatrixType & matrixT() const
Definition: CompleteOrthogonalDecomposition.h:185
CompleteOrthogonalDecomposition & setThreshold(const RealScalar &threshold)
Definition: CompleteOrthogonalDecomposition.h:322
CompleteOrthogonalDecomposition(const EigenBase< InputType > &matrix)
Constructs a complete orthogonal decomposition from a given matrix.
Definition: CompleteOrthogonalDecomposition.h:120
Index dimensionOfKernel() const
Definition: CompleteOrthogonalDecomposition.h:246
const Inverse< CompleteOrthogonalDecomposition > pseudoInverse() const
Definition: CompleteOrthogonalDecomposition.h:280
MatrixType::RealScalar absDeterminant() const
Definition: CompleteOrthogonalDecomposition.h:412
RealScalar maxPivot() const
Definition: CompleteOrthogonalDecomposition.h:358
HouseholderSequenceType householderQ(void) const
Definition: CompleteOrthogonalDecomposition.h:617
MatrixType::RealScalar logAbsDeterminant() const
Definition: CompleteOrthogonalDecomposition.h:418
CompleteOrthogonalDecomposition(Index rows, Index cols)
Default Constructor with memory preallocation.
Definition: CompleteOrthogonalDecomposition.h:100
void computeInPlace()
Definition: CompleteOrthogonalDecomposition.h:430
const MatrixType & matrixQTZ() const
Definition: CompleteOrthogonalDecomposition.h:172
CompleteOrthogonalDecomposition(EigenBase< InputType > &matrix)
Constructs a complete orthogonal decomposition from a given matrix.
Definition: CompleteOrthogonalDecomposition.h:135
const HCoeffsType & hCoeffs() const
Definition: CompleteOrthogonalDecomposition.h:294
void applyZOnTheLeftInPlace(Rhs &rhs) const
Definition: CompleteOrthogonalDecomposition.h:487
Sequence of Householder reflections acting on subspaces with decreasing size.
Definition: HouseholderSequence.h:123
Expression of the inverse of another expression.
Definition: Inverse.h:46
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:52
const CompleteOrthogonalDecomposition< PlainObject > completeOrthogonalDecomposition() const
Definition: CompleteOrthogonalDecomposition.h:627
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:182
Pseudo expression representing a solving operation.
Definition: Solve.h:65
A base class for matrix decomposition and solvers.
Definition: SolverBase.h:71
CompleteOrthogonalDecomposition< MatrixType_ > & derived()
Definition: EigenBase.h:48
ComputationInfo
Definition: Constants.h:442
@ Success
Definition: Constants.h:444
Namespace containing all symbols from the Eigen library.
Definition: Core:139
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:59
Definition: EigenBase.h:32
Derived & derived()
Definition: EigenBase.h:48
Eigen::Index Index
The interface type of indices.
Definition: EigenBase.h:41
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:231