Eigen  3.4.90 (git rev a4098ac676528a83cfb73d4d26ce1b42ec05f47c)
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
15namespace Eigen {
16
17namespace internal {
18template <typename MatrixType_>
19struct 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
52template <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, typename internal::remove_all<
77 typename HCoeffsType::ConjugateReturnType>::type>
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 {
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>
189 // Compute the column pivoted QR factorization A P = Q R.
190 m_cpqr.compute(matrix);
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
410template <typename MatrixType>
411typename MatrixType::RealScalar
413 return m_cpqr.absDeterminant();
414}
415
416template <typename MatrixType>
417typename MatrixType::RealScalar
419 return m_cpqr.logAbsDeterminant();
420}
421
429template <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
485template <typename MatrixType>
486template <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
507template <typename MatrixType>
508template <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
530template <typename MatrixType_>
531template <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
562template<typename MatrixType_>
563template<bool Conjugate, typename RhsType, typename DstType>
564void 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
591namespace internal {
592
593template<typename MatrixType>
594struct traits<Inverse<CompleteOrthogonalDecomposition<MatrixType> > >
595 : traits<typename Transpose<typename MatrixType::PlainObject>::PlainObject>
596{
597 enum { Flags = 0 };
598};
599
600template<typename DstXprType, typename MatrixType>
601struct 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
615template <typename MatrixType>
616typename CompleteOrthogonalDecomposition<MatrixType>::HouseholderSequenceType
618 return m_cpqr.householderQ();
619}
620
625template <typename Derived>
629}
630
631} // end namespace Eigen
632
633#endif // EIGEN_COMPLETEORTHOGONALDECOMPOSITION_H
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
Index nonzeroPivots() const
Definition: CompleteOrthogonalDecomposition.h:353
const MatrixType & matrixT() const
Definition: CompleteOrthogonalDecomposition.h:185
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
MatrixType matrixZ() const
Definition: CompleteOrthogonalDecomposition.h:163
const PermutationType & colsPermutation() const
Definition: CompleteOrthogonalDecomposition.h:196
bool isSurjective() const
Definition: CompleteOrthogonalDecomposition.h:264
const HCoeffsType & hCoeffs() const
Definition: CompleteOrthogonalDecomposition.h:294
CompleteOrthogonalDecomposition & setThreshold(const RealScalar &threshold)
Definition: CompleteOrthogonalDecomposition.h:322
CompleteOrthogonalDecomposition()
Default Constructor.
Definition: CompleteOrthogonalDecomposition.h:92
const Solve< CompleteOrthogonalDecomposition, Rhs > solve(const MatrixBase< Rhs > &b) const
CompleteOrthogonalDecomposition(const EigenBase< InputType > &matrix)
Constructs a complete orthogonal decomposition from a given matrix.
Definition: CompleteOrthogonalDecomposition.h:120
const MatrixType & matrixQTZ() const
Definition: CompleteOrthogonalDecomposition.h:172
const Inverse< CompleteOrthogonalDecomposition > pseudoInverse() const
Definition: CompleteOrthogonalDecomposition.h:280
Index dimensionOfKernel() const
Definition: CompleteOrthogonalDecomposition.h:246
MatrixType::RealScalar absDeterminant() const
Definition: CompleteOrthogonalDecomposition.h:412
RealScalar maxPivot() const
Definition: CompleteOrthogonalDecomposition.h:358
HouseholderSequenceType householderQ(void) const
Definition: CompleteOrthogonalDecomposition.h:617
const HCoeffsType & zCoeffs() const
Definition: CompleteOrthogonalDecomposition.h:301
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
CompleteOrthogonalDecomposition(EigenBase< InputType > &matrix)
Constructs a complete orthogonal decomposition from a given matrix.
Definition: CompleteOrthogonalDecomposition.h:135
CompleteOrthogonalDecomposition & setThreshold(Default_t)
Definition: CompleteOrthogonalDecomposition.h:335
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: B01_Experimental.dox:1
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:59
Definition: EigenBase.h:32
Eigen::Index Index
The interface type of indices.
Definition: EigenBase.h:41
Derived & derived()
Definition: EigenBase.h:48
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:235