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