Please, help us to better know about our user community by answering the following short survey: https://forms.gle/wpyrxWi18ox9Z5ae9
Eigen  3.4.0
 
Loading...
Searching...
No Matches
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
13namespace Eigen {
14
15namespace internal {
16template <typename _MatrixType>
17struct 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
50template <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 {
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>
187 // Compute the column pivoted QR factorization A P = Q R.
188 m_cpqr.compute(matrix);
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
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 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
487template <typename MatrixType>
488template <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
509template <typename MatrixType>
510template <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
532template <typename _MatrixType>
533template <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
564template<typename _MatrixType>
565template<bool Conjugate, typename RhsType, typename DstType>
566void 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
593namespace internal {
594
595template<typename MatrixType>
596struct traits<Inverse<CompleteOrthogonalDecomposition<MatrixType> > >
597 : traits<typename Transpose<typename MatrixType::PlainObject>::PlainObject>
598{
599 enum { Flags = 0 };
600};
601
602template<typename DstXprType, typename MatrixType>
603struct 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
617template <typename MatrixType>
618typename CompleteOrthogonalDecomposition<MatrixType>::HouseholderSequenceType
620 return m_cpqr.householderQ();
621}
622
627template <typename Derived>
631}
632
633} // end namespace Eigen
634
635#endif // EIGEN_COMPLETEORTHOGONALDECOMPOSITION_H
bool isInjective() const
Definition: ColPivHouseholderQR.h:285
const HCoeffsType & hCoeffs() const
Definition: ColPivHouseholderQR.h:334
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
const PermutationType & colsPermutation() const
Definition: ColPivHouseholderQR.h:214
Index dimensionOfKernel() const
Definition: ColPivHouseholderQR.h:272
bool isSurjective() const
Definition: ColPivHouseholderQR.h:298
bool isInvertible() const
Definition: ColPivHouseholderQR.h:310
RealScalar maxPivot() const
Definition: ColPivHouseholderQR.h:403
ColPivHouseholderQR & setThreshold(const RealScalar &threshold)
Definition: ColPivHouseholderQR.h:353
const MatrixType & matrixQR() const
Definition: ColPivHouseholderQR.h:189
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 Inverse< CompleteOrthogonalDecomposition > pseudoInverse() const
Definition: CompleteOrthogonalDecomposition.h:278
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
RealScalar threshold() const
Definition: CompleteOrthogonalDecomposition.h:342
CompleteOrthogonalDecomposition & setThreshold(Default_t)
Definition: CompleteOrthogonalDecomposition.h:333
MatrixType matrixZ() const
Definition: CompleteOrthogonalDecomposition.h:161
bool isSurjective() const
Definition: CompleteOrthogonalDecomposition.h:262
RealScalar maxPivot() const
Definition: CompleteOrthogonalDecomposition.h:356
CompleteOrthogonalDecomposition()
Default Constructor.
Definition: CompleteOrthogonalDecomposition.h:90
bool isInvertible() const
Definition: CompleteOrthogonalDecomposition.h:271
const MatrixType & matrixQTZ() const
Definition: CompleteOrthogonalDecomposition.h:170
CompleteOrthogonalDecomposition(Index rows, Index cols)
Default Constructor with memory preallocation.
Definition: CompleteOrthogonalDecomposition.h:98
const PermutationType & colsPermutation() const
Definition: CompleteOrthogonalDecomposition.h:194
const MatrixType & matrixT() const
Definition: CompleteOrthogonalDecomposition.h:183
MatrixType::RealScalar absDeterminant() const
Definition: CompleteOrthogonalDecomposition.h:412
const HCoeffsType & hCoeffs() const
Definition: CompleteOrthogonalDecomposition.h:292
HouseholderSequenceType householderQ(void) const
Definition: CompleteOrthogonalDecomposition.h:619
Index dimensionOfKernel() const
Definition: CompleteOrthogonalDecomposition.h:244
MatrixType::RealScalar logAbsDeterminant() const
Definition: CompleteOrthogonalDecomposition.h:418
CompleteOrthogonalDecomposition & setThreshold(const RealScalar &threshold)
Definition: CompleteOrthogonalDecomposition.h:320
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
const Solve< CompleteOrthogonalDecomposition, Rhs > solve(const MatrixBase< Rhs > &b) const
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
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:180
Permutation matrix.
Definition: PermutationMatrix.h:298
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
Eigen::Index Index
The interface type of indices.
Definition: EigenBase.h:39
Derived & derived()
Definition: EigenBase.h:46
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:233