Eigen  3.3.90 (mercurial changeset d63bd3cb7378)
 All Classes Namespaces Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
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  enum { Flags = 0 };
20 };
21 
22 } // end namespace internal
23 
47 template <typename _MatrixType>
48 class CompleteOrthogonalDecomposition {
49  public:
50  typedef _MatrixType MatrixType;
51  enum {
52  RowsAtCompileTime = MatrixType::RowsAtCompileTime,
53  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
54  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
55  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
56  };
57  typedef typename MatrixType::Scalar Scalar;
58  typedef typename MatrixType::RealScalar RealScalar;
59  typedef typename MatrixType::StorageIndex StorageIndex;
60  typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
61  typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime>
62  PermutationType;
63  typedef typename internal::plain_row_type<MatrixType, Index>::type
64  IntRowVectorType;
65  typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
66  typedef typename internal::plain_row_type<MatrixType, RealScalar>::type
67  RealRowVectorType;
68  typedef HouseholderSequence<
69  MatrixType, typename internal::remove_all<
70  typename HCoeffsType::ConjugateReturnType>::type>
71  HouseholderSequenceType;
72  typedef typename MatrixType::PlainObject PlainObject;
73 
74  private:
75  typedef typename PermutationType::Index PermIndexType;
76 
77  public:
85  CompleteOrthogonalDecomposition() : m_cpqr(), m_zCoeffs(), m_temp() {}
86 
94  : m_cpqr(rows, cols), m_zCoeffs((std::min)(rows, cols)), m_temp(cols) {}
95 
112  template <typename InputType>
114  : m_cpqr(matrix.rows(), matrix.cols()),
115  m_zCoeffs((std::min)(matrix.rows(), matrix.cols())),
116  m_temp(matrix.cols())
117  {
118  compute(matrix.derived());
119  }
120 
127  template<typename InputType>
129  : m_cpqr(matrix.derived()),
130  m_zCoeffs((std::min)(matrix.rows(), matrix.cols())),
131  m_temp(matrix.cols())
132  {
133  computeInPlace();
134  }
135 
136 
146  template <typename Rhs>
148  const MatrixBase<Rhs>& b) const {
149  eigen_assert(m_cpqr.m_isInitialized &&
150  "CompleteOrthogonalDecomposition is not initialized.");
152  }
153 
154  HouseholderSequenceType householderQ(void) const;
155  HouseholderSequenceType matrixQ(void) const { return m_cpqr.householderQ(); }
156 
159  MatrixType matrixZ() const {
160  MatrixType Z = MatrixType::Identity(m_cpqr.cols(), m_cpqr.cols());
162  return Z.adjoint();
163  }
164 
168  const MatrixType& matrixQTZ() const { return m_cpqr.matrixQR(); }
169 
181  const MatrixType& matrixT() const { return m_cpqr.matrixQR(); }
182 
183  template <typename InputType>
184  CompleteOrthogonalDecomposition& compute(const EigenBase<InputType>& matrix) {
185  // Compute the column pivoted QR factorization A P = Q R.
186  m_cpqr.compute(matrix);
187  computeInPlace();
188  return *this;
189  }
190 
193  return m_cpqr.colsPermutation();
194  }
195 
209  typename MatrixType::RealScalar absDeterminant() const;
210 
224  typename MatrixType::RealScalar logAbsDeterminant() const;
225 
233  inline Index rank() const { return m_cpqr.rank(); }
234 
242  inline Index dimensionOfKernel() const { return m_cpqr.dimensionOfKernel(); }
243 
251  inline bool isInjective() const { return m_cpqr.isInjective(); }
252 
260  inline bool isSurjective() const { return m_cpqr.isSurjective(); }
261 
269  inline bool isInvertible() const { return m_cpqr.isInvertible(); }
270 
277  {
279  }
280 
281  inline Index rows() const { return m_cpqr.rows(); }
282  inline Index cols() const { return m_cpqr.cols(); }
283 
289  inline const HCoeffsType& hCoeffs() const { return m_cpqr.hCoeffs(); }
290 
296  const HCoeffsType& zCoeffs() const { return m_zCoeffs; }
297 
318  m_cpqr.setThreshold(threshold);
319  return *this;
320  }
321 
331  m_cpqr.setThreshold(Default);
332  return *this;
333  }
334 
339  RealScalar threshold() const { return m_cpqr.threshold(); }
340 
348  inline Index nonzeroPivots() const { return m_cpqr.nonzeroPivots(); }
349 
353  inline RealScalar maxPivot() const { return m_cpqr.maxPivot(); }
354 
364  eigen_assert(m_cpqr.m_isInitialized && "Decomposition is not initialized.");
365  return Success;
366  }
367 
368 #ifndef EIGEN_PARSED_BY_DOXYGEN
369  template <typename RhsType, typename DstType>
370  void _solve_impl(const RhsType& rhs, DstType& dst) const;
371 #endif
372 
373  protected:
374  static void check_template_parameters() {
375  EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
376  }
377 
378  void computeInPlace();
379 
382  template <typename Rhs>
383  void applyZAdjointOnTheLeftInPlace(Rhs& rhs) const;
384 
385  ColPivHouseholderQR<MatrixType> m_cpqr;
386  HCoeffsType m_zCoeffs;
387  RowVectorType m_temp;
388 };
389 
390 template <typename MatrixType>
391 typename MatrixType::RealScalar
393  return m_cpqr.absDeterminant();
394 }
395 
396 template <typename MatrixType>
397 typename MatrixType::RealScalar
399  return m_cpqr.logAbsDeterminant();
400 }
401 
409 template <typename MatrixType>
411 {
412  check_template_parameters();
413 
414  // the column permutation is stored as int indices, so just to be sure:
415  eigen_assert(m_cpqr.cols() <= NumTraits<int>::highest());
416 
417  const Index rank = m_cpqr.rank();
418  const Index cols = m_cpqr.cols();
419  const Index rows = m_cpqr.rows();
420  m_zCoeffs.resize((std::min)(rows, cols));
421  m_temp.resize(cols);
422 
423  if (rank < cols) {
424  // We have reduced the (permuted) matrix to the form
425  // [R11 R12]
426  // [ 0 R22]
427  // where R11 is r-by-r (r = rank) upper triangular, R12 is
428  // r-by-(n-r), and R22 is empty or the norm of R22 is negligible.
429  // We now compute the complete orthogonal decomposition by applying
430  // Householder transformations from the right to the upper trapezoidal
431  // matrix X = [R11 R12] to zero out R12 and obtain the factorization
432  // [R11 R12] = [T11 0] * Z, where T11 is r-by-r upper triangular and
433  // Z = Z(0) * Z(1) ... Z(r-1) is an n-by-n orthogonal matrix.
434  // We store the data representing Z in R12 and m_zCoeffs.
435  for (Index k = rank - 1; k >= 0; --k) {
436  if (k != rank - 1) {
437  // Given the API for Householder reflectors, it is more convenient if
438  // we swap the leading parts of columns k and r-1 (zero-based) to form
439  // the matrix X_k = [X(0:k, k), X(0:k, r:n)]
440  m_cpqr.m_qr.col(k).head(k + 1).swap(
441  m_cpqr.m_qr.col(rank - 1).head(k + 1));
442  }
443  // Construct Householder reflector Z(k) to zero out the last row of X_k,
444  // i.e. choose Z(k) such that
445  // [X(k, k), X(k, r:n)] * Z(k) = [beta, 0, .., 0].
446  RealScalar beta;
447  m_cpqr.m_qr.row(k)
448  .tail(cols - rank + 1)
449  .makeHouseholderInPlace(m_zCoeffs(k), beta);
450  m_cpqr.m_qr(k, rank - 1) = beta;
451  if (k > 0) {
452  // Apply Z(k) to the first k rows of X_k
453  m_cpqr.m_qr.topRightCorner(k, cols - rank + 1)
454  .applyHouseholderOnTheRight(
455  m_cpqr.m_qr.row(k).tail(cols - rank).transpose(), m_zCoeffs(k),
456  &m_temp(0));
457  }
458  if (k != rank - 1) {
459  // Swap X(0:k,k) back to its proper location.
460  m_cpqr.m_qr.col(k).head(k + 1).swap(
461  m_cpqr.m_qr.col(rank - 1).head(k + 1));
462  }
463  }
464  }
465 }
466 
467 template <typename MatrixType>
468 template <typename Rhs>
470  Rhs& rhs) const {
471  const Index cols = this->cols();
472  const Index nrhs = rhs.cols();
473  const Index rank = this->rank();
474  Matrix<typename MatrixType::Scalar, Dynamic, 1> temp((std::max)(cols, nrhs));
475  for (Index k = 0; k < rank; ++k) {
476  if (k != rank - 1) {
477  rhs.row(k).swap(rhs.row(rank - 1));
478  }
479  rhs.middleRows(rank - 1, cols - rank + 1)
480  .applyHouseholderOnTheLeft(
481  matrixQTZ().row(k).tail(cols - rank).adjoint(), zCoeffs()(k),
482  &temp(0));
483  if (k != rank - 1) {
484  rhs.row(k).swap(rhs.row(rank - 1));
485  }
486  }
487 }
488 
489 #ifndef EIGEN_PARSED_BY_DOXYGEN
490 template <typename _MatrixType>
491 template <typename RhsType, typename DstType>
493  const RhsType& rhs, DstType& dst) const {
494  eigen_assert(rhs.rows() == this->rows());
495 
496  const Index rank = this->rank();
497  if (rank == 0) {
498  dst.setZero();
499  return;
500  }
501 
502  // Compute c = Q^* * rhs
503  // Note that the matrix Q = H_0^* H_1^*... so its inverse is
504  // Q^* = (H_0 H_1 ...)^T
505  typename RhsType::PlainObject c(rhs);
506  c.applyOnTheLeft(
507  householderSequence(matrixQTZ(), hCoeffs()).setLength(rank).transpose());
508 
509  // Solve T z = c(1:rank, :)
510  dst.topRows(rank) = matrixT()
511  .topLeftCorner(rank, rank)
512  .template triangularView<Upper>()
513  .solve(c.topRows(rank));
514 
515  const Index cols = this->cols();
516  if (rank < cols) {
517  // Compute y = Z^* * [ z ]
518  // [ 0 ]
519  dst.bottomRows(cols - rank).setZero();
520  applyZAdjointOnTheLeftInPlace(dst);
521  }
522 
523  // Undo permutation to get x = P^{-1} * y.
524  dst = colsPermutation() * dst;
525 }
526 #endif
527 
528 namespace internal {
529 
530 template<typename DstXprType, typename MatrixType>
531 struct Assignment<DstXprType, Inverse<CompleteOrthogonalDecomposition<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename CompleteOrthogonalDecomposition<MatrixType>::Scalar>, Dense2Dense>
532 {
533  typedef CompleteOrthogonalDecomposition<MatrixType> CodType;
534  typedef Inverse<CodType> SrcXprType;
535  static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename CodType::Scalar> &)
536  {
537  dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.rows()));
538  }
539 };
540 
541 } // end namespace internal
542 
544 template <typename MatrixType>
545 typename CompleteOrthogonalDecomposition<MatrixType>::HouseholderSequenceType
547  return m_cpqr.householderQ();
548 }
549 
554 template <typename Derived>
558 }
559 
560 } // end namespace Eigen
561 
562 #endif // EIGEN_COMPLETEORTHOGONALDECOMPOSITION_H
const PermutationType & colsPermutation() const
Definition: CompleteOrthogonalDecomposition.h:192
CompleteOrthogonalDecomposition(const EigenBase< InputType > &matrix)
Constructs a complete orthogonal decomposition from a given matrix.
Definition: CompleteOrthogonalDecomposition.h:113
CompleteOrthogonalDecomposition & setThreshold(const RealScalar &threshold)
Definition: CompleteOrthogonalDecomposition.h:317
CompleteOrthogonalDecomposition()
Default Constructor.
Definition: CompleteOrthogonalDecomposition.h:85
bool isSurjective() const
Definition: CompleteOrthogonalDecomposition.h:260
MatrixType::RealScalar absDeterminant() const
Definition: CompleteOrthogonalDecomposition.h:392
CompleteOrthogonalDecomposition & setThreshold(Default_t)
Definition: CompleteOrthogonalDecomposition.h:330
bool isInvertible() const
Definition: CompleteOrthogonalDecomposition.h:269
MatrixType::RealScalar logAbsDeterminant() const
Definition: CompleteOrthogonalDecomposition.h:398
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:150
void applyZAdjointOnTheLeftInPlace(Rhs &rhs) const
Definition: CompleteOrthogonalDecomposition.h:469
Derived & derived()
Definition: EigenBase.h:45
HouseholderSequenceType householderQ(void) const
Definition: CompleteOrthogonalDecomposition.h:546
Complete orthogonal decomposition (COD) of a matrix.
Definition: ForwardDeclarations.h:258
Eigen::Index Index
The interface type of indices.
Definition: EigenBase.h:38
HouseholderSequence< VectorsType, CoeffsType > householderSequence(const VectorsType &v, const CoeffsType &h)
Convenience function for constructing a Householder sequence.
Definition: HouseholderSequence.h:451
Definition: EigenBase.h:29
Expression of the inverse of another expression.
Definition: Inverse.h:43
void computeInPlace()
Definition: CompleteOrthogonalDecomposition.h:410
const MatrixType & matrixT() const
Definition: CompleteOrthogonalDecomposition.h:181
MatrixType matrixZ() const
Definition: CompleteOrthogonalDecomposition.h:159
const MatrixType & matrixQTZ() const
Definition: CompleteOrthogonalDecomposition.h:168
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:33
Index nonzeroPivots() const
Definition: CompleteOrthogonalDecomposition.h:348
CompleteOrthogonalDecomposition(Index rows, Index cols)
Default Constructor with memory preallocation.
Definition: CompleteOrthogonalDecomposition.h:93
const CompleteOrthogonalDecomposition< PlainObject > completeOrthogonalDecomposition() const
Definition: CompleteOrthogonalDecomposition.h:556
const HCoeffsType & hCoeffs() const
Definition: CompleteOrthogonalDecomposition.h:289
Definition: Constants.h:436
Index dimensionOfKernel() const
Definition: CompleteOrthogonalDecomposition.h:242
RealScalar maxPivot() const
Definition: CompleteOrthogonalDecomposition.h:353
const Solve< CompleteOrthogonalDecomposition, Rhs > solve(const MatrixBase< Rhs > &b) const
Definition: CompleteOrthogonalDecomposition.h:147
const HCoeffsType & zCoeffs() const
Definition: CompleteOrthogonalDecomposition.h:296
CompleteOrthogonalDecomposition(EigenBase< InputType > &matrix)
Constructs a complete orthogonal decomposition from a given matrix.
Definition: CompleteOrthogonalDecomposition.h:128
Index rank() const
Definition: CompleteOrthogonalDecomposition.h:233
Pseudo expression representing a solving operation.
Definition: Solve.h:62
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:178
ComputationInfo
Definition: Constants.h:434
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
const Inverse< CompleteOrthogonalDecomposition > pseudoInverse() const
Definition: CompleteOrthogonalDecomposition.h:276
RealScalar threshold() const
Definition: CompleteOrthogonalDecomposition.h:339
ComputationInfo info() const
Reports whether the complete orthogonal decomposition was succesful.
Definition: CompleteOrthogonalDecomposition.h:363
bool isInjective() const
Definition: CompleteOrthogonalDecomposition.h:251