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