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
FullPivHouseholderQR.h
1// This file is part of Eigen, a lightweight C++ template library
2// for linear algebra.
3//
4// Copyright (C) 2008-2009 Gael Guennebaud <gael.guennebaud@inria.fr>
5// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
6//
7// This Source Code Form is subject to the terms of the Mozilla
8// Public License v. 2.0. If a copy of the MPL was not distributed
9// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
10
11#ifndef EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H
12#define EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H
13
14namespace Eigen {
15
16namespace internal {
17
18template<typename _MatrixType> struct traits<FullPivHouseholderQR<_MatrixType> >
19 : traits<_MatrixType>
20{
21 typedef MatrixXpr XprKind;
22 typedef SolverStorage StorageKind;
23 typedef int StorageIndex;
24 enum { Flags = 0 };
25};
26
27template<typename MatrixType> struct FullPivHouseholderQRMatrixQReturnType;
28
29template<typename MatrixType>
30struct traits<FullPivHouseholderQRMatrixQReturnType<MatrixType> >
31{
32 typedef typename MatrixType::PlainObject ReturnType;
33};
34
35} // end namespace internal
36
60template<typename _MatrixType> class FullPivHouseholderQR
61 : public SolverBase<FullPivHouseholderQR<_MatrixType> >
62{
63 public:
64
65 typedef _MatrixType MatrixType;
67 friend class SolverBase<FullPivHouseholderQR>;
68
69 EIGEN_GENERIC_PUBLIC_INTERFACE(FullPivHouseholderQR)
70 enum {
71 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
72 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
73 };
74 typedef internal::FullPivHouseholderQRMatrixQReturnType<MatrixType> MatrixQReturnType;
75 typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
76 typedef Matrix<StorageIndex, 1,
77 EIGEN_SIZE_MIN_PREFER_DYNAMIC(ColsAtCompileTime,RowsAtCompileTime), RowMajor, 1,
78 EIGEN_SIZE_MIN_PREFER_FIXED(MaxColsAtCompileTime,MaxRowsAtCompileTime)> IntDiagSizeVectorType;
80 typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
81 typedef typename internal::plain_col_type<MatrixType>::type ColVectorType;
82 typedef typename MatrixType::PlainObject PlainObject;
83
90 : m_qr(),
91 m_hCoeffs(),
92 m_rows_transpositions(),
93 m_cols_transpositions(),
94 m_cols_permutation(),
95 m_temp(),
96 m_isInitialized(false),
97 m_usePrescribedThreshold(false) {}
98
106 : m_qr(rows, cols),
107 m_hCoeffs((std::min)(rows,cols)),
108 m_rows_transpositions((std::min)(rows,cols)),
109 m_cols_transpositions((std::min)(rows,cols)),
110 m_cols_permutation(cols),
111 m_temp(cols),
112 m_isInitialized(false),
113 m_usePrescribedThreshold(false) {}
114
127 template<typename InputType>
129 : m_qr(matrix.rows(), matrix.cols()),
130 m_hCoeffs((std::min)(matrix.rows(), matrix.cols())),
131 m_rows_transpositions((std::min)(matrix.rows(), matrix.cols())),
132 m_cols_transpositions((std::min)(matrix.rows(), matrix.cols())),
133 m_cols_permutation(matrix.cols()),
134 m_temp(matrix.cols()),
135 m_isInitialized(false),
136 m_usePrescribedThreshold(false)
137 {
138 compute(matrix.derived());
139 }
140
147 template<typename InputType>
149 : m_qr(matrix.derived()),
150 m_hCoeffs((std::min)(matrix.rows(), matrix.cols())),
151 m_rows_transpositions((std::min)(matrix.rows(), matrix.cols())),
152 m_cols_transpositions((std::min)(matrix.rows(), matrix.cols())),
153 m_cols_permutation(matrix.cols()),
154 m_temp(matrix.cols()),
155 m_isInitialized(false),
156 m_usePrescribedThreshold(false)
157 {
158 computeInPlace();
159 }
160
161 #ifdef EIGEN_PARSED_BY_DOXYGEN
177 template<typename Rhs>
179 solve(const MatrixBase<Rhs>& b) const;
180 #endif
181
184 MatrixQReturnType matrixQ(void) const;
185
188 const MatrixType& matrixQR() const
189 {
190 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
191 return m_qr;
192 }
193
194 template<typename InputType>
195 FullPivHouseholderQR& compute(const EigenBase<InputType>& matrix);
196
199 {
200 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
201 return m_cols_permutation;
202 }
203
206 {
207 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
208 return m_rows_transpositions;
209 }
210
224 typename MatrixType::RealScalar absDeterminant() const;
225
238 typename MatrixType::RealScalar logAbsDeterminant() const;
239
246 inline Index rank() const
247 {
248 using std::abs;
249 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
250 RealScalar premultiplied_threshold = abs(m_maxpivot) * threshold();
251 Index result = 0;
252 for(Index i = 0; i < m_nonzero_pivots; ++i)
253 result += (abs(m_qr.coeff(i,i)) > premultiplied_threshold);
254 return result;
255 }
256
264 {
265 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
266 return cols() - rank();
267 }
268
276 inline bool isInjective() const
277 {
278 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
279 return rank() == cols();
280 }
281
289 inline bool isSurjective() const
290 {
291 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
292 return rank() == rows();
293 }
294
301 inline bool isInvertible() const
302 {
303 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
304 return isInjective() && isSurjective();
305 }
306
313 {
314 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
315 return Inverse<FullPivHouseholderQR>(*this);
316 }
317
318 inline Index rows() const { return m_qr.rows(); }
319 inline Index cols() const { return m_qr.cols(); }
320
325 const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
326
345 {
346 m_usePrescribedThreshold = true;
347 m_prescribedThreshold = threshold;
348 return *this;
349 }
350
360 {
361 m_usePrescribedThreshold = false;
362 return *this;
363 }
364
369 RealScalar threshold() const
370 {
371 eigen_assert(m_isInitialized || m_usePrescribedThreshold);
372 return m_usePrescribedThreshold ? m_prescribedThreshold
373 // this formula comes from experimenting (see "LU precision tuning" thread on the list)
374 // and turns out to be identical to Higham's formula used already in LDLt.
375 : NumTraits<Scalar>::epsilon() * RealScalar(m_qr.diagonalSize());
376 }
377
385 inline Index nonzeroPivots() const
386 {
387 eigen_assert(m_isInitialized && "LU is not initialized.");
388 return m_nonzero_pivots;
389 }
390
394 RealScalar maxPivot() const { return m_maxpivot; }
395
396 #ifndef EIGEN_PARSED_BY_DOXYGEN
397 template<typename RhsType, typename DstType>
398 void _solve_impl(const RhsType &rhs, DstType &dst) const;
399
400 template<bool Conjugate, typename RhsType, typename DstType>
401 void _solve_impl_transposed(const RhsType &rhs, DstType &dst) const;
402 #endif
403
404 protected:
405
406 static void check_template_parameters()
407 {
408 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
409 }
410
411 void computeInPlace();
412
413 MatrixType m_qr;
414 HCoeffsType m_hCoeffs;
415 IntDiagSizeVectorType m_rows_transpositions;
416 IntDiagSizeVectorType m_cols_transpositions;
417 PermutationType m_cols_permutation;
418 RowVectorType m_temp;
419 bool m_isInitialized, m_usePrescribedThreshold;
420 RealScalar m_prescribedThreshold, m_maxpivot;
421 Index m_nonzero_pivots;
422 RealScalar m_precision;
423 Index m_det_pq;
424};
425
426template<typename MatrixType>
427typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType>::absDeterminant() const
428{
429 using std::abs;
430 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
431 eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
432 return abs(m_qr.diagonal().prod());
433}
434
435template<typename MatrixType>
436typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType>::logAbsDeterminant() const
437{
438 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
439 eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
440 return m_qr.diagonal().cwiseAbs().array().log().sum();
441}
442
449template<typename MatrixType>
450template<typename InputType>
452{
453 m_qr = matrix.derived();
454 computeInPlace();
455 return *this;
456}
457
458template<typename MatrixType>
460{
461 check_template_parameters();
462
463 using std::abs;
464 Index rows = m_qr.rows();
465 Index cols = m_qr.cols();
466 Index size = (std::min)(rows,cols);
467
468
469 m_hCoeffs.resize(size);
470
471 m_temp.resize(cols);
472
473 m_precision = NumTraits<Scalar>::epsilon() * RealScalar(size);
474
475 m_rows_transpositions.resize(size);
476 m_cols_transpositions.resize(size);
477 Index number_of_transpositions = 0;
478
479 RealScalar biggest(0);
480
481 m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
482 m_maxpivot = RealScalar(0);
483
484 for (Index k = 0; k < size; ++k)
485 {
486 Index row_of_biggest_in_corner, col_of_biggest_in_corner;
487 typedef internal::scalar_score_coeff_op<Scalar> Scoring;
488 typedef typename Scoring::result_type Score;
489
490 Score score = m_qr.bottomRightCorner(rows-k, cols-k)
491 .unaryExpr(Scoring())
492 .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner);
493 row_of_biggest_in_corner += k;
494 col_of_biggest_in_corner += k;
495 RealScalar biggest_in_corner = internal::abs_knowing_score<Scalar>()(m_qr(row_of_biggest_in_corner, col_of_biggest_in_corner), score);
496 if(k==0) biggest = biggest_in_corner;
497
498 // if the corner is negligible, then we have less than full rank, and we can finish early
499 if(internal::isMuchSmallerThan(biggest_in_corner, biggest, m_precision))
500 {
501 m_nonzero_pivots = k;
502 for(Index i = k; i < size; i++)
503 {
504 m_rows_transpositions.coeffRef(i) = internal::convert_index<StorageIndex>(i);
505 m_cols_transpositions.coeffRef(i) = internal::convert_index<StorageIndex>(i);
506 m_hCoeffs.coeffRef(i) = Scalar(0);
507 }
508 break;
509 }
510
511 m_rows_transpositions.coeffRef(k) = internal::convert_index<StorageIndex>(row_of_biggest_in_corner);
512 m_cols_transpositions.coeffRef(k) = internal::convert_index<StorageIndex>(col_of_biggest_in_corner);
513 if(k != row_of_biggest_in_corner) {
514 m_qr.row(k).tail(cols-k).swap(m_qr.row(row_of_biggest_in_corner).tail(cols-k));
515 ++number_of_transpositions;
516 }
517 if(k != col_of_biggest_in_corner) {
518 m_qr.col(k).swap(m_qr.col(col_of_biggest_in_corner));
519 ++number_of_transpositions;
520 }
521
522 RealScalar beta;
523 m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta);
524 m_qr.coeffRef(k,k) = beta;
525
526 // remember the maximum absolute value of diagonal coefficients
527 if(abs(beta) > m_maxpivot) m_maxpivot = abs(beta);
528
529 m_qr.bottomRightCorner(rows-k, cols-k-1)
530 .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1));
531 }
532
533 m_cols_permutation.setIdentity(cols);
534 for(Index k = 0; k < size; ++k)
535 m_cols_permutation.applyTranspositionOnTheRight(k, m_cols_transpositions.coeff(k));
536
537 m_det_pq = (number_of_transpositions%2) ? -1 : 1;
538 m_isInitialized = true;
539}
540
541#ifndef EIGEN_PARSED_BY_DOXYGEN
542template<typename _MatrixType>
543template<typename RhsType, typename DstType>
544void FullPivHouseholderQR<_MatrixType>::_solve_impl(const RhsType &rhs, DstType &dst) const
545{
546 const Index l_rank = rank();
547
548 // FIXME introduce nonzeroPivots() and use it here. and more generally,
549 // make the same improvements in this dec as in FullPivLU.
550 if(l_rank==0)
551 {
552 dst.setZero();
553 return;
554 }
555
556 typename RhsType::PlainObject c(rhs);
557
558 Matrix<typename RhsType::Scalar,1,RhsType::ColsAtCompileTime> temp(rhs.cols());
559 for (Index k = 0; k < l_rank; ++k)
560 {
561 Index remainingSize = rows()-k;
562 c.row(k).swap(c.row(m_rows_transpositions.coeff(k)));
563 c.bottomRightCorner(remainingSize, rhs.cols())
564 .applyHouseholderOnTheLeft(m_qr.col(k).tail(remainingSize-1),
565 m_hCoeffs.coeff(k), &temp.coeffRef(0));
566 }
567
568 m_qr.topLeftCorner(l_rank, l_rank)
569 .template triangularView<Upper>()
570 .solveInPlace(c.topRows(l_rank));
571
572 for(Index i = 0; i < l_rank; ++i) dst.row(m_cols_permutation.indices().coeff(i)) = c.row(i);
573 for(Index i = l_rank; i < cols(); ++i) dst.row(m_cols_permutation.indices().coeff(i)).setZero();
574}
575
576template<typename _MatrixType>
577template<bool Conjugate, typename RhsType, typename DstType>
578void FullPivHouseholderQR<_MatrixType>::_solve_impl_transposed(const RhsType &rhs, DstType &dst) const
579{
580 const Index l_rank = rank();
581
582 if(l_rank == 0)
583 {
584 dst.setZero();
585 return;
586 }
587
588 typename RhsType::PlainObject c(m_cols_permutation.transpose()*rhs);
589
590 m_qr.topLeftCorner(l_rank, l_rank)
591 .template triangularView<Upper>()
592 .transpose().template conjugateIf<Conjugate>()
593 .solveInPlace(c.topRows(l_rank));
594
595 dst.topRows(l_rank) = c.topRows(l_rank);
596 dst.bottomRows(rows()-l_rank).setZero();
597
598 Matrix<Scalar, 1, DstType::ColsAtCompileTime> temp(dst.cols());
599 const Index size = (std::min)(rows(), cols());
600 for (Index k = size-1; k >= 0; --k)
601 {
602 Index remainingSize = rows()-k;
603
604 dst.bottomRightCorner(remainingSize, dst.cols())
605 .applyHouseholderOnTheLeft(m_qr.col(k).tail(remainingSize-1).template conjugateIf<!Conjugate>(),
606 m_hCoeffs.template conjugateIf<Conjugate>().coeff(k), &temp.coeffRef(0));
607
608 dst.row(k).swap(dst.row(m_rows_transpositions.coeff(k)));
609 }
610}
611#endif
612
613namespace internal {
614
615template<typename DstXprType, typename MatrixType>
616struct Assignment<DstXprType, Inverse<FullPivHouseholderQR<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename FullPivHouseholderQR<MatrixType>::Scalar>, Dense2Dense>
617{
618 typedef FullPivHouseholderQR<MatrixType> QrType;
619 typedef Inverse<QrType> SrcXprType;
620 static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename QrType::Scalar> &)
621 {
622 dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.cols()));
623 }
624};
625
632template<typename MatrixType> struct FullPivHouseholderQRMatrixQReturnType
633 : public ReturnByValue<FullPivHouseholderQRMatrixQReturnType<MatrixType> >
634{
635public:
636 typedef typename FullPivHouseholderQR<MatrixType>::IntDiagSizeVectorType IntDiagSizeVectorType;
637 typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
638 typedef Matrix<typename MatrixType::Scalar, 1, MatrixType::RowsAtCompileTime, RowMajor, 1,
639 MatrixType::MaxRowsAtCompileTime> WorkVectorType;
640
641 FullPivHouseholderQRMatrixQReturnType(const MatrixType& qr,
642 const HCoeffsType& hCoeffs,
643 const IntDiagSizeVectorType& rowsTranspositions)
644 : m_qr(qr),
645 m_hCoeffs(hCoeffs),
646 m_rowsTranspositions(rowsTranspositions)
647 {}
648
649 template <typename ResultType>
650 void evalTo(ResultType& result) const
651 {
652 const Index rows = m_qr.rows();
653 WorkVectorType workspace(rows);
654 evalTo(result, workspace);
655 }
656
657 template <typename ResultType>
658 void evalTo(ResultType& result, WorkVectorType& workspace) const
659 {
660 using numext::conj;
661 // compute the product H'_0 H'_1 ... H'_n-1,
662 // where H_k is the k-th Householder transformation I - h_k v_k v_k'
663 // and v_k is the k-th Householder vector [1,m_qr(k+1,k), m_qr(k+2,k), ...]
664 const Index rows = m_qr.rows();
665 const Index cols = m_qr.cols();
666 const Index size = (std::min)(rows, cols);
667 workspace.resize(rows);
668 result.setIdentity(rows, rows);
669 for (Index k = size-1; k >= 0; k--)
670 {
671 result.block(k, k, rows-k, rows-k)
672 .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), conj(m_hCoeffs.coeff(k)), &workspace.coeffRef(k));
673 result.row(k).swap(result.row(m_rowsTranspositions.coeff(k)));
674 }
675 }
676
677 Index rows() const { return m_qr.rows(); }
678 Index cols() const { return m_qr.rows(); }
679
680protected:
681 typename MatrixType::Nested m_qr;
682 typename HCoeffsType::Nested m_hCoeffs;
683 typename IntDiagSizeVectorType::Nested m_rowsTranspositions;
684};
685
686// template<typename MatrixType>
687// struct evaluator<FullPivHouseholderQRMatrixQReturnType<MatrixType> >
688// : public evaluator<ReturnByValue<FullPivHouseholderQRMatrixQReturnType<MatrixType> > >
689// {};
690
691} // end namespace internal
692
693template<typename MatrixType>
694inline typename FullPivHouseholderQR<MatrixType>::MatrixQReturnType FullPivHouseholderQR<MatrixType>::matrixQ() const
695{
696 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
697 return MatrixQReturnType(m_qr, m_hCoeffs, m_rows_transpositions);
698}
699
704template<typename Derived>
707{
709}
710
711} // end namespace Eigen
712
713#endif // EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H
Householder rank-revealing QR decomposition of a matrix with full pivoting.
Definition: FullPivHouseholderQR.h:62
FullPivHouseholderQR & setThreshold(const RealScalar &threshold)
Definition: FullPivHouseholderQR.h:344
MatrixType::RealScalar absDeterminant() const
Definition: FullPivHouseholderQR.h:427
const MatrixType & matrixQR() const
Definition: FullPivHouseholderQR.h:188
FullPivHouseholderQR & setThreshold(Default_t)
Definition: FullPivHouseholderQR.h:359
const PermutationType & colsPermutation() const
Definition: FullPivHouseholderQR.h:198
Index dimensionOfKernel() const
Definition: FullPivHouseholderQR.h:263
const IntDiagSizeVectorType & rowsTranspositions() const
Definition: FullPivHouseholderQR.h:205
bool isInjective() const
Definition: FullPivHouseholderQR.h:276
const HCoeffsType & hCoeffs() const
Definition: FullPivHouseholderQR.h:325
RealScalar maxPivot() const
Definition: FullPivHouseholderQR.h:394
bool isSurjective() const
Definition: FullPivHouseholderQR.h:289
MatrixType::RealScalar logAbsDeterminant() const
Definition: FullPivHouseholderQR.h:436
FullPivHouseholderQR(Index rows, Index cols)
Default Constructor with memory preallocation.
Definition: FullPivHouseholderQR.h:105
FullPivHouseholderQR(EigenBase< InputType > &matrix)
Constructs a QR factorization from a given matrix.
Definition: FullPivHouseholderQR.h:148
const Solve< FullPivHouseholderQR, Rhs > solve(const MatrixBase< Rhs > &b) const
MatrixQReturnType matrixQ(void) const
Definition: FullPivHouseholderQR.h:694
const Inverse< FullPivHouseholderQR > inverse() const
Definition: FullPivHouseholderQR.h:312
Index rank() const
Definition: FullPivHouseholderQR.h:246
FullPivHouseholderQR()
Default Constructor.
Definition: FullPivHouseholderQR.h:89
bool isInvertible() const
Definition: FullPivHouseholderQR.h:301
FullPivHouseholderQR(const EigenBase< InputType > &matrix)
Constructs a QR factorization from a given matrix.
Definition: FullPivHouseholderQR.h:128
Index nonzeroPivots() const
Definition: FullPivHouseholderQR.h:385
RealScalar threshold() const
Definition: FullPivHouseholderQR.h:369
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
FullPivHouseholderQR< _MatrixType > & derived()
Definition: EigenBase.h:46
const Solve< Derived, Rhs > solve(const MatrixBase< Rhs > &b) const
Definition: SolverBase.h:106
@ RowMajor
Definition: Constants.h:321
Namespace containing all symbols from the Eigen library.
Definition: Core:141
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_conjugate_op< typename Derived::Scalar >, const Derived > conj(const Eigen::ArrayBase< Derived > &x)
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