Eigen  3.4.0 (git rev e3e74001f7c4bf95f0dde572e8a08c5b2918a3ab)
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 namespace Eigen {
14 
15 namespace internal {
16 template <typename _MatrixType>
17 struct traits<CompleteOrthogonalDecomposition<_MatrixType> >
18  : traits<_MatrixType> {
19  typedef MatrixXpr XprKind;
20  typedef SolverStorage StorageKind;
21  typedef int StorageIndex;
22  enum { Flags = 0 };
23 };
24 
25 } // end namespace internal
26 
50 template <typename _MatrixType> class CompleteOrthogonalDecomposition
51  : public SolverBase<CompleteOrthogonalDecomposition<_MatrixType> >
52 {
53  public:
54  typedef _MatrixType MatrixType;
56 
57  template<typename Derived>
58  friend struct internal::solve_assertion;
59 
60  EIGEN_GENERIC_PUBLIC_INTERFACE(CompleteOrthogonalDecomposition)
61  enum {
62  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
63  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
64  };
65  typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
68  typedef typename internal::plain_row_type<MatrixType, Index>::type
69  IntRowVectorType;
70  typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
71  typedef typename internal::plain_row_type<MatrixType, RealScalar>::type
72  RealRowVectorType;
73  typedef HouseholderSequence<
74  MatrixType, typename internal::remove_all<
75  typename HCoeffsType::ConjugateReturnType>::type>
77  typedef typename MatrixType::PlainObject PlainObject;
78 
79  private:
80  typedef typename PermutationType::Index PermIndexType;
81 
82  public:
90  CompleteOrthogonalDecomposition() : m_cpqr(), m_zCoeffs(), m_temp() {}
91 
99  : m_cpqr(rows, cols), m_zCoeffs((std::min)(rows, cols)), m_temp(cols) {}
100 
117  template <typename InputType>
119  : m_cpqr(matrix.rows(), matrix.cols()),
120  m_zCoeffs((std::min)(matrix.rows(), matrix.cols())),
121  m_temp(matrix.cols())
122  {
123  compute(matrix.derived());
124  }
125 
132  template<typename InputType>
134  : m_cpqr(matrix.derived()),
135  m_zCoeffs((std::min)(matrix.rows(), matrix.cols())),
136  m_temp(matrix.cols())
137  {
138  computeInPlace();
139  }
140 
141  #ifdef EIGEN_PARSED_BY_DOXYGEN
151  template <typename Rhs>
153  const MatrixBase<Rhs>& b) const;
154  #endif
155 
157  HouseholderSequenceType matrixQ(void) const { return m_cpqr.householderQ(); }
158 
161  MatrixType matrixZ() const {
162  MatrixType Z = MatrixType::Identity(m_cpqr.cols(), m_cpqr.cols());
163  applyZOnTheLeftInPlace<false>(Z);
164  return Z;
165  }
166 
170  const MatrixType& matrixQTZ() const { return m_cpqr.matrixQR(); }
171 
183  const MatrixType& matrixT() const { return m_cpqr.matrixQR(); }
184 
185  template <typename InputType>
186  CompleteOrthogonalDecomposition& compute(const EigenBase<InputType>& matrix) {
187  // Compute the column pivoted QR factorization A P = Q R.
188  m_cpqr.compute(matrix);
189  computeInPlace();
190  return *this;
191  }
192 
195  return m_cpqr.colsPermutation();
196  }
197 
211  typename MatrixType::RealScalar absDeterminant() const;
212 
226  typename MatrixType::RealScalar logAbsDeterminant() const;
227 
235  inline Index rank() const { return m_cpqr.rank(); }
236 
244  inline Index dimensionOfKernel() const { return m_cpqr.dimensionOfKernel(); }
245 
253  inline bool isInjective() const { return m_cpqr.isInjective(); }
254 
262  inline bool isSurjective() const { return m_cpqr.isSurjective(); }
263 
271  inline bool isInvertible() const { return m_cpqr.isInvertible(); }
272 
279  {
280  eigen_assert(m_cpqr.m_isInitialized && "CompleteOrthogonalDecomposition is not initialized.");
282  }
283 
284  inline Index rows() const { return m_cpqr.rows(); }
285  inline Index cols() const { return m_cpqr.cols(); }
286 
292  inline const HCoeffsType& hCoeffs() const { return m_cpqr.hCoeffs(); }
293 
299  const HCoeffsType& zCoeffs() const { return m_zCoeffs; }
300 
321  m_cpqr.setThreshold(threshold);
322  return *this;
323  }
324 
334  m_cpqr.setThreshold(Default);
335  return *this;
336  }
337 
342  RealScalar threshold() const { return m_cpqr.threshold(); }
343 
351  inline Index nonzeroPivots() const { return m_cpqr.nonzeroPivots(); }
352 
356  inline RealScalar maxPivot() const { return m_cpqr.maxPivot(); }
357 
367  eigen_assert(m_cpqr.m_isInitialized && "Decomposition is not initialized.");
368  return Success;
369  }
370 
371 #ifndef EIGEN_PARSED_BY_DOXYGEN
372  template <typename RhsType, typename DstType>
373  void _solve_impl(const RhsType& rhs, DstType& dst) const;
374 
375  template<bool Conjugate, typename RhsType, typename DstType>
376  void _solve_impl_transposed(const RhsType &rhs, DstType &dst) const;
377 #endif
378 
379  protected:
380  static void check_template_parameters() {
381  EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
382  }
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  check_template_parameters();
433 
434  // the column permutation is stored as int indices, so just to be sure:
435  eigen_assert(m_cpqr.cols() <= NumTraits<int>::highest());
436 
437  const Index rank = m_cpqr.rank();
438  const Index cols = m_cpqr.cols();
439  const Index rows = m_cpqr.rows();
440  m_zCoeffs.resize((std::min)(rows, cols));
441  m_temp.resize(cols);
442 
443  if (rank < cols) {
444  // We have reduced the (permuted) matrix to the form
445  // [R11 R12]
446  // [ 0 R22]
447  // where R11 is r-by-r (r = rank) upper triangular, R12 is
448  // r-by-(n-r), and R22 is empty or the norm of R22 is negligible.
449  // We now compute the complete orthogonal decomposition by applying
450  // Householder transformations from the right to the upper trapezoidal
451  // matrix X = [R11 R12] to zero out R12 and obtain the factorization
452  // [R11 R12] = [T11 0] * Z, where T11 is r-by-r upper triangular and
453  // Z = Z(0) * Z(1) ... Z(r-1) is an n-by-n orthogonal matrix.
454  // We store the data representing Z in R12 and m_zCoeffs.
455  for (Index k = rank - 1; k >= 0; --k) {
456  if (k != rank - 1) {
457  // Given the API for Householder reflectors, it is more convenient if
458  // we swap the leading parts of columns k and r-1 (zero-based) to form
459  // the matrix X_k = [X(0:k, k), X(0:k, r:n)]
460  m_cpqr.m_qr.col(k).head(k + 1).swap(
461  m_cpqr.m_qr.col(rank - 1).head(k + 1));
462  }
463  // Construct Householder reflector Z(k) to zero out the last row of X_k,
464  // i.e. choose Z(k) such that
465  // [X(k, k), X(k, r:n)] * Z(k) = [beta, 0, .., 0].
466  RealScalar beta;
467  m_cpqr.m_qr.row(k)
468  .tail(cols - rank + 1)
469  .makeHouseholderInPlace(m_zCoeffs(k), beta);
470  m_cpqr.m_qr(k, rank - 1) = beta;
471  if (k > 0) {
472  // Apply Z(k) to the first k rows of X_k
473  m_cpqr.m_qr.topRightCorner(k, cols - rank + 1)
474  .applyHouseholderOnTheRight(
475  m_cpqr.m_qr.row(k).tail(cols - rank).adjoint(), m_zCoeffs(k),
476  &m_temp(0));
477  }
478  if (k != rank - 1) {
479  // Swap X(0:k,k) back to its proper location.
480  m_cpqr.m_qr.col(k).head(k + 1).swap(
481  m_cpqr.m_qr.col(rank - 1).head(k + 1));
482  }
483  }
484  }
485 }
486 
487 template <typename MatrixType>
488 template <bool Conjugate, typename Rhs>
490  Rhs& rhs) const {
491  const Index cols = this->cols();
492  const Index nrhs = rhs.cols();
493  const Index rank = this->rank();
494  Matrix<typename Rhs::Scalar, Dynamic, 1> temp((std::max)(cols, nrhs));
495  for (Index k = rank-1; k >= 0; --k) {
496  if (k != rank - 1) {
497  rhs.row(k).swap(rhs.row(rank - 1));
498  }
499  rhs.middleRows(rank - 1, cols - rank + 1)
500  .applyHouseholderOnTheLeft(
501  matrixQTZ().row(k).tail(cols - rank).transpose().template conjugateIf<!Conjugate>(), zCoeffs().template conjugateIf<Conjugate>()(k),
502  &temp(0));
503  if (k != rank - 1) {
504  rhs.row(k).swap(rhs.row(rank - 1));
505  }
506  }
507 }
508 
509 template <typename MatrixType>
510 template <typename Rhs>
512  Rhs& rhs) const {
513  const Index cols = this->cols();
514  const Index nrhs = rhs.cols();
515  const Index rank = this->rank();
516  Matrix<typename Rhs::Scalar, Dynamic, 1> temp((std::max)(cols, nrhs));
517  for (Index k = 0; k < rank; ++k) {
518  if (k != rank - 1) {
519  rhs.row(k).swap(rhs.row(rank - 1));
520  }
521  rhs.middleRows(rank - 1, cols - rank + 1)
522  .applyHouseholderOnTheLeft(
523  matrixQTZ().row(k).tail(cols - rank).adjoint(), zCoeffs()(k),
524  &temp(0));
525  if (k != rank - 1) {
526  rhs.row(k).swap(rhs.row(rank - 1));
527  }
528  }
529 }
530 
531 #ifndef EIGEN_PARSED_BY_DOXYGEN
532 template <typename _MatrixType>
533 template <typename RhsType, typename DstType>
535  const RhsType& rhs, DstType& dst) const {
536  const Index rank = this->rank();
537  if (rank == 0) {
538  dst.setZero();
539  return;
540  }
541 
542  // Compute c = Q^* * rhs
543  typename RhsType::PlainObject c(rhs);
544  c.applyOnTheLeft(matrixQ().setLength(rank).adjoint());
545 
546  // Solve T z = c(1:rank, :)
547  dst.topRows(rank) = matrixT()
548  .topLeftCorner(rank, rank)
549  .template triangularView<Upper>()
550  .solve(c.topRows(rank));
551 
552  const Index cols = this->cols();
553  if (rank < cols) {
554  // Compute y = Z^* * [ z ]
555  // [ 0 ]
556  dst.bottomRows(cols - rank).setZero();
557  applyZAdjointOnTheLeftInPlace(dst);
558  }
559 
560  // Undo permutation to get x = P^{-1} * y.
561  dst = colsPermutation() * dst;
562 }
563 
564 template<typename _MatrixType>
565 template<bool Conjugate, typename RhsType, typename DstType>
566 void CompleteOrthogonalDecomposition<_MatrixType>::_solve_impl_transposed(const RhsType &rhs, DstType &dst) const
567 {
568  const Index rank = this->rank();
569 
570  if (rank == 0) {
571  dst.setZero();
572  return;
573  }
574 
575  typename RhsType::PlainObject c(colsPermutation().transpose()*rhs);
576 
577  if (rank < cols()) {
578  applyZOnTheLeftInPlace<!Conjugate>(c);
579  }
580 
581  matrixT().topLeftCorner(rank, rank)
582  .template triangularView<Upper>()
583  .transpose().template conjugateIf<Conjugate>()
584  .solveInPlace(c.topRows(rank));
585 
586  dst.topRows(rank) = c.topRows(rank);
587  dst.bottomRows(rows()-rank).setZero();
588 
589  dst.applyOnTheLeft(householderQ().setLength(rank).template conjugateIf<!Conjugate>() );
590 }
591 #endif
592 
593 namespace internal {
594 
595 template<typename MatrixType>
596 struct traits<Inverse<CompleteOrthogonalDecomposition<MatrixType> > >
597  : traits<typename Transpose<typename MatrixType::PlainObject>::PlainObject>
598 {
599  enum { Flags = 0 };
600 };
601 
602 template<typename DstXprType, typename MatrixType>
603 struct Assignment<DstXprType, Inverse<CompleteOrthogonalDecomposition<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename CompleteOrthogonalDecomposition<MatrixType>::Scalar>, Dense2Dense>
604 {
605  typedef CompleteOrthogonalDecomposition<MatrixType> CodType;
606  typedef Inverse<CodType> SrcXprType;
607  static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename CodType::Scalar> &)
608  {
609  typedef Matrix<typename CodType::Scalar, CodType::RowsAtCompileTime, CodType::RowsAtCompileTime, 0, CodType::MaxRowsAtCompileTime, CodType::MaxRowsAtCompileTime> IdentityMatrixType;
610  dst = src.nestedExpression().solve(IdentityMatrixType::Identity(src.cols(), src.cols()));
611  }
612 };
613 
614 } // end namespace internal
615 
617 template <typename MatrixType>
618 typename CompleteOrthogonalDecomposition<MatrixType>::HouseholderSequenceType
620  return m_cpqr.householderQ();
621 }
622 
627 template <typename Derived>
631 }
632 
633 } // end namespace Eigen
634 
635 #endif // EIGEN_COMPLETEORTHOGONALDECOMPOSITION_H
bool isInjective() const
Definition: ColPivHouseholderQR.h:285
HouseholderSequenceType householderQ() const
Definition: ColPivHouseholderQR.h:655
Index rank() const
Definition: ColPivHouseholderQR.h:255
RealScalar threshold() const
Definition: ColPivHouseholderQR.h:378
Index nonzeroPivots() const
Definition: ColPivHouseholderQR.h:394
Index dimensionOfKernel() const
Definition: ColPivHouseholderQR.h:272
bool isSurjective() const
Definition: ColPivHouseholderQR.h:298
bool isInvertible() const
Definition: ColPivHouseholderQR.h:310
const MatrixType & matrixQR() const
Definition: ColPivHouseholderQR.h:189
RealScalar maxPivot() const
Definition: ColPivHouseholderQR.h:403
const PermutationType & colsPermutation() const
Definition: ColPivHouseholderQR.h:214
const HCoeffsType & hCoeffs() const
Definition: ColPivHouseholderQR.h:334
ColPivHouseholderQR & setThreshold(const RealScalar &threshold)
Definition: ColPivHouseholderQR.h:353
Complete orthogonal decomposition (COD) of a matrix.
Definition: CompleteOrthogonalDecomposition.h:52
CompleteOrthogonalDecomposition(EigenBase< InputType > &matrix)
Constructs a complete orthogonal decomposition from a given matrix.
Definition: CompleteOrthogonalDecomposition.h:133
void applyZAdjointOnTheLeftInPlace(Rhs &rhs) const
Definition: CompleteOrthogonalDecomposition.h:511
const HCoeffsType & zCoeffs() const
Definition: CompleteOrthogonalDecomposition.h:299
ComputationInfo info() const
Reports whether the complete orthogonal decomposition was successful.
Definition: CompleteOrthogonalDecomposition.h:366
bool isInjective() const
Definition: CompleteOrthogonalDecomposition.h:253
CompleteOrthogonalDecomposition & setThreshold(Default_t)
Definition: CompleteOrthogonalDecomposition.h:333
RealScalar threshold() const
Definition: CompleteOrthogonalDecomposition.h:342
const Inverse< CompleteOrthogonalDecomposition > pseudoInverse() const
Definition: CompleteOrthogonalDecomposition.h:278
MatrixType matrixZ() const
Definition: CompleteOrthogonalDecomposition.h:161
bool isSurjective() const
Definition: CompleteOrthogonalDecomposition.h:262
RealScalar maxPivot() const
Definition: CompleteOrthogonalDecomposition.h:356
const PermutationType & colsPermutation() const
Definition: CompleteOrthogonalDecomposition.h:194
CompleteOrthogonalDecomposition()
Default Constructor.
Definition: CompleteOrthogonalDecomposition.h:90
bool isInvertible() const
Definition: CompleteOrthogonalDecomposition.h:271
const MatrixType & matrixT() const
Definition: CompleteOrthogonalDecomposition.h:183
CompleteOrthogonalDecomposition(Index rows, Index cols)
Default Constructor with memory preallocation.
Definition: CompleteOrthogonalDecomposition.h:98
CompleteOrthogonalDecomposition & setThreshold(const RealScalar &threshold)
Definition: CompleteOrthogonalDecomposition.h:320
const Solve< CompleteOrthogonalDecomposition, Rhs > solve(const MatrixBase< Rhs > &b) const
MatrixType::RealScalar absDeterminant() const
Definition: CompleteOrthogonalDecomposition.h:412
HouseholderSequenceType householderQ(void) const
Definition: CompleteOrthogonalDecomposition.h:619
Index dimensionOfKernel() const
Definition: CompleteOrthogonalDecomposition.h:244
const HCoeffsType & hCoeffs() const
Definition: CompleteOrthogonalDecomposition.h:292
MatrixType::RealScalar logAbsDeterminant() const
Definition: CompleteOrthogonalDecomposition.h:418
const MatrixType & matrixQTZ() const
Definition: CompleteOrthogonalDecomposition.h:170
void computeInPlace()
Definition: CompleteOrthogonalDecomposition.h:430
Index rank() const
Definition: CompleteOrthogonalDecomposition.h:235
Index nonzeroPivots() const
Definition: CompleteOrthogonalDecomposition.h:351
CompleteOrthogonalDecomposition(const EigenBase< InputType > &matrix)
Constructs a complete orthogonal decomposition from a given matrix.
Definition: CompleteOrthogonalDecomposition.h:118
void applyZOnTheLeftInPlace(Rhs &rhs) const
Definition: CompleteOrthogonalDecomposition.h:489
Sequence of Householder reflections acting on subspaces with decreasing size.
Definition: HouseholderSequence.h:121
Expression of the inverse of another expression.
Definition: Inverse.h:44
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:50
const CompleteOrthogonalDecomposition< PlainObject > completeOrthogonalDecomposition() const
Definition: CompleteOrthogonalDecomposition.h:629
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:180
Pseudo expression representing a solving operation.
Definition: Solve.h:63
A base class for matrix decomposition and solvers.
Definition: SolverBase.h:69
CompleteOrthogonalDecomposition< _MatrixType > & derived()
Definition: EigenBase.h:46
ComputationInfo
Definition: Constants.h:440
@ Success
Definition: Constants.h:442
Namespace containing all symbols from the Eigen library.
Definition: Core:141
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
Definition: EigenBase.h:30
Derived & derived()
Definition: EigenBase.h:46
Eigen::Index Index
The interface type of indices.
Definition: EigenBase.h:39
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:233