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_COLPIVOTINGHOUSEHOLDERQR_H
12 #define EIGEN_COLPIVOTINGHOUSEHOLDERQR_H
13 
14 namespace Eigen {
15 
16 /** \ingroup QR_Module
17   *
18   * \class ColPivHouseholderQR
19   *
20   * \brief Householder rank-revealing QR decomposition of a matrix with column-pivoting
21   *
22   * \param MatrixType the type of the matrix of which we are computing the QR decomposition
23   *
24   * This class performs a rank-revealing QR decomposition of a matrix \b A into matrices \b P, \b Q and \b R
25   * such that
26   * \f[
27   *  \mathbf{A} \, \mathbf{P} = \mathbf{Q} \, \mathbf{R}
28   * \f]
29   * by using Householder transformations. Here, \b P is a permutation matrix, \b Q a unitary matrix and \b R an
30   * upper triangular matrix.
31   *
32   * This decomposition performs column pivoting in order to be rank-revealing and improve
33   * numerical stability. It is slower than HouseholderQR, and faster than FullPivHouseholderQR.
34   *
35   * \sa MatrixBase::colPivHouseholderQr()
36   */
37 template<typename _MatrixType> class ColPivHouseholderQR
38 {
39   public:
40 
41     typedef _MatrixType MatrixType;
42     enum {
43       RowsAtCompileTime = MatrixType::RowsAtCompileTime,
44       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
45       Options = MatrixType::Options,
46       MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
47       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
48     };
49     typedef typename MatrixType::Scalar Scalar;
50     typedef typename MatrixType::RealScalar RealScalar;
51     typedef typename MatrixType::Index Index;
52     typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, Options, MaxRowsAtCompileTime, MaxRowsAtCompileTime> MatrixQType;
53     typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
54     typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationType;
55     typedef typename internal::plain_row_type<MatrixType, Index>::type IntRowVectorType;
56     typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
57     typedef typename internal::plain_row_type<MatrixType, RealScalar>::type RealRowVectorType;
58     typedef HouseholderSequence<MatrixType,typename internal::remove_all<typename HCoeffsType::ConjugateReturnType>::type> HouseholderSequenceType;
59 
60   private:
61 
62     typedef typename PermutationType::Index PermIndexType;
63 
64   public:
65 
66     /**
67     * \brief Default Constructor.
68     *
69     * The default constructor is useful in cases in which the user intends to
70     * perform decompositions via ColPivHouseholderQR::compute(const MatrixType&).
71     */
ColPivHouseholderQR()72     ColPivHouseholderQR()
73       : m_qr(),
74         m_hCoeffs(),
75         m_colsPermutation(),
76         m_colsTranspositions(),
77         m_temp(),
78         m_colSqNorms(),
79         m_isInitialized(false),
80         m_usePrescribedThreshold(false) {}
81 
82     /** \brief Default Constructor with memory preallocation
83       *
84       * Like the default constructor but with preallocation of the internal data
85       * according to the specified problem \a size.
86       * \sa ColPivHouseholderQR()
87       */
ColPivHouseholderQR(Index rows,Index cols)88     ColPivHouseholderQR(Index rows, Index cols)
89       : m_qr(rows, cols),
90         m_hCoeffs((std::min)(rows,cols)),
91         m_colsPermutation(PermIndexType(cols)),
92         m_colsTranspositions(cols),
93         m_temp(cols),
94         m_colSqNorms(cols),
95         m_isInitialized(false),
96         m_usePrescribedThreshold(false) {}
97 
98     /** \brief Constructs a QR factorization from a given matrix
99       *
100       * This constructor computes the QR factorization of the matrix \a matrix by calling
101       * the method compute(). It is a short cut for:
102       *
103       * \code
104       * ColPivHouseholderQR<MatrixType> qr(matrix.rows(), matrix.cols());
105       * qr.compute(matrix);
106       * \endcode
107       *
108       * \sa compute()
109       */
ColPivHouseholderQR(const MatrixType & matrix)110     ColPivHouseholderQR(const MatrixType& matrix)
111       : m_qr(matrix.rows(), matrix.cols()),
112         m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
113         m_colsPermutation(PermIndexType(matrix.cols())),
114         m_colsTranspositions(matrix.cols()),
115         m_temp(matrix.cols()),
116         m_colSqNorms(matrix.cols()),
117         m_isInitialized(false),
118         m_usePrescribedThreshold(false)
119     {
120       compute(matrix);
121     }
122 
123     /** This method finds a solution x to the equation Ax=b, where A is the matrix of which
124       * *this is the QR decomposition, if any exists.
125       *
126       * \param b the right-hand-side of the equation to solve.
127       *
128       * \returns a solution.
129       *
130       * \note The case where b is a matrix is not yet implemented. Also, this
131       *       code is space inefficient.
132       *
133       * \note_about_checking_solutions
134       *
135       * \note_about_arbitrary_choice_of_solution
136       *
137       * Example: \include ColPivHouseholderQR_solve.cpp
138       * Output: \verbinclude ColPivHouseholderQR_solve.out
139       */
140     template<typename Rhs>
141     inline const internal::solve_retval<ColPivHouseholderQR, Rhs>
solve(const MatrixBase<Rhs> & b)142     solve(const MatrixBase<Rhs>& b) const
143     {
144       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
145       return internal::solve_retval<ColPivHouseholderQR, Rhs>(*this, b.derived());
146     }
147 
148     HouseholderSequenceType householderQ(void) const;
matrixQ(void)149     HouseholderSequenceType matrixQ(void) const
150     {
151       return householderQ();
152     }
153 
154     /** \returns a reference to the matrix where the Householder QR decomposition is stored
155       */
matrixQR()156     const MatrixType& matrixQR() const
157     {
158       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
159       return m_qr;
160     }
161 
162     /** \returns a reference to the matrix where the result Householder QR is stored
163      * \warning The strict lower part of this matrix contains internal values.
164      * Only the upper triangular part should be referenced. To get it, use
165      * \code matrixR().template triangularView<Upper>() \endcode
166      * For rank-deficient matrices, use
167      * \code
168      * matrixR().topLeftCorner(rank(), rank()).template triangularView<Upper>()
169      * \endcode
170      */
matrixR()171     const MatrixType& matrixR() const
172     {
173       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
174       return m_qr;
175     }
176 
177     ColPivHouseholderQR& compute(const MatrixType& matrix);
178 
179     /** \returns a const reference to the column permutation matrix */
colsPermutation()180     const PermutationType& colsPermutation() const
181     {
182       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
183       return m_colsPermutation;
184     }
185 
186     /** \returns the absolute value of the determinant of the matrix of which
187       * *this is the QR decomposition. It has only linear complexity
188       * (that is, O(n) where n is the dimension of the square matrix)
189       * as the QR decomposition has already been computed.
190       *
191       * \note This is only for square matrices.
192       *
193       * \warning a determinant can be very big or small, so for matrices
194       * of large enough dimension, there is a risk of overflow/underflow.
195       * One way to work around that is to use logAbsDeterminant() instead.
196       *
197       * \sa logAbsDeterminant(), MatrixBase::determinant()
198       */
199     typename MatrixType::RealScalar absDeterminant() const;
200 
201     /** \returns the natural log of the absolute value of the determinant of the matrix of which
202       * *this is the QR decomposition. It has only linear complexity
203       * (that is, O(n) where n is the dimension of the square matrix)
204       * as the QR decomposition has already been computed.
205       *
206       * \note This is only for square matrices.
207       *
208       * \note This method is useful to work around the risk of overflow/underflow that's inherent
209       * to determinant computation.
210       *
211       * \sa absDeterminant(), MatrixBase::determinant()
212       */
213     typename MatrixType::RealScalar logAbsDeterminant() const;
214 
215     /** \returns the rank of the matrix of which *this is the QR decomposition.
216       *
217       * \note This method has to determine which pivots should be considered nonzero.
218       *       For that, it uses the threshold value that you can control by calling
219       *       setThreshold(const RealScalar&).
220       */
rank()221     inline Index rank() const
222     {
223       using std::abs;
224       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
225       RealScalar premultiplied_threshold = abs(m_maxpivot) * threshold();
226       Index result = 0;
227       for(Index i = 0; i < m_nonzero_pivots; ++i)
228         result += (abs(m_qr.coeff(i,i)) > premultiplied_threshold);
229       return result;
230     }
231 
232     /** \returns the dimension of the kernel of the matrix of which *this is the QR decomposition.
233       *
234       * \note This method has to determine which pivots should be considered nonzero.
235       *       For that, it uses the threshold value that you can control by calling
236       *       setThreshold(const RealScalar&).
237       */
dimensionOfKernel()238     inline Index dimensionOfKernel() const
239     {
240       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
241       return cols() - rank();
242     }
243 
244     /** \returns true if the matrix of which *this is the QR decomposition represents an injective
245       *          linear map, i.e. has trivial kernel; false otherwise.
246       *
247       * \note This method has to determine which pivots should be considered nonzero.
248       *       For that, it uses the threshold value that you can control by calling
249       *       setThreshold(const RealScalar&).
250       */
isInjective()251     inline bool isInjective() const
252     {
253       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
254       return rank() == cols();
255     }
256 
257     /** \returns true if the matrix of which *this is the QR decomposition represents a surjective
258       *          linear map; false otherwise.
259       *
260       * \note This method has to determine which pivots should be considered nonzero.
261       *       For that, it uses the threshold value that you can control by calling
262       *       setThreshold(const RealScalar&).
263       */
isSurjective()264     inline bool isSurjective() const
265     {
266       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
267       return rank() == rows();
268     }
269 
270     /** \returns true if the matrix of which *this is the QR decomposition is invertible.
271       *
272       * \note This method has to determine which pivots should be considered nonzero.
273       *       For that, it uses the threshold value that you can control by calling
274       *       setThreshold(const RealScalar&).
275       */
isInvertible()276     inline bool isInvertible() const
277     {
278       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
279       return isInjective() && isSurjective();
280     }
281 
282     /** \returns the inverse of the matrix of which *this is the QR decomposition.
283       *
284       * \note If this matrix is not invertible, the returned matrix has undefined coefficients.
285       *       Use isInvertible() to first determine whether this matrix is invertible.
286       */
287     inline const
288     internal::solve_retval<ColPivHouseholderQR, typename MatrixType::IdentityReturnType>
inverse()289     inverse() const
290     {
291       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
292       return internal::solve_retval<ColPivHouseholderQR,typename MatrixType::IdentityReturnType>
293                (*this, MatrixType::Identity(m_qr.rows(), m_qr.cols()));
294     }
295 
rows()296     inline Index rows() const { return m_qr.rows(); }
cols()297     inline Index cols() const { return m_qr.cols(); }
298 
299     /** \returns a const reference to the vector of Householder coefficients used to represent the factor \c Q.
300       *
301       * For advanced uses only.
302       */
hCoeffs()303     const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
304 
305     /** Allows to prescribe a threshold to be used by certain methods, such as rank(),
306       * who need to determine when pivots are to be considered nonzero. This is not used for the
307       * QR decomposition itself.
308       *
309       * When it needs to get the threshold value, Eigen calls threshold(). By default, this
310       * uses a formula to automatically determine a reasonable threshold.
311       * Once you have called the present method setThreshold(const RealScalar&),
312       * your value is used instead.
313       *
314       * \param threshold The new value to use as the threshold.
315       *
316       * A pivot will be considered nonzero if its absolute value is strictly greater than
317       *  \f$ \vert pivot \vert \leqslant threshold \times \vert maxpivot \vert \f$
318       * where maxpivot is the biggest pivot.
319       *
320       * If you want to come back to the default behavior, call setThreshold(Default_t)
321       */
setThreshold(const RealScalar & threshold)322     ColPivHouseholderQR& setThreshold(const RealScalar& threshold)
323     {
324       m_usePrescribedThreshold = true;
325       m_prescribedThreshold = threshold;
326       return *this;
327     }
328 
329     /** Allows to come back to the default behavior, letting Eigen use its default formula for
330       * determining the threshold.
331       *
332       * You should pass the special object Eigen::Default as parameter here.
333       * \code qr.setThreshold(Eigen::Default); \endcode
334       *
335       * See the documentation of setThreshold(const RealScalar&).
336       */
setThreshold(Default_t)337     ColPivHouseholderQR& setThreshold(Default_t)
338     {
339       m_usePrescribedThreshold = false;
340       return *this;
341     }
342 
343     /** Returns the threshold that will be used by certain methods such as rank().
344       *
345       * See the documentation of setThreshold(const RealScalar&).
346       */
threshold()347     RealScalar threshold() const
348     {
349       eigen_assert(m_isInitialized || m_usePrescribedThreshold);
350       return m_usePrescribedThreshold ? m_prescribedThreshold
351       // this formula comes from experimenting (see "LU precision tuning" thread on the list)
352       // and turns out to be identical to Higham's formula used already in LDLt.
353                                       : NumTraits<Scalar>::epsilon() * RealScalar(m_qr.diagonalSize());
354     }
355 
356     /** \returns the number of nonzero pivots in the QR decomposition.
357       * Here nonzero is meant in the exact sense, not in a fuzzy sense.
358       * So that notion isn't really intrinsically interesting, but it is
359       * still useful when implementing algorithms.
360       *
361       * \sa rank()
362       */
nonzeroPivots()363     inline Index nonzeroPivots() const
364     {
365       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
366       return m_nonzero_pivots;
367     }
368 
369     /** \returns the absolute value of the biggest pivot, i.e. the biggest
370       *          diagonal coefficient of R.
371       */
maxPivot()372     RealScalar maxPivot() const { return m_maxpivot; }
373 
374     /** \brief Reports whether the QR factorization was succesful.
375       *
376       * \note This function always returns \c Success. It is provided for compatibility
377       * with other factorization routines.
378       * \returns \c Success
379       */
info()380     ComputationInfo info() const
381     {
382       eigen_assert(m_isInitialized && "Decomposition is not initialized.");
383       return Success;
384     }
385 
386   protected:
387 
check_template_parameters()388     static void check_template_parameters()
389     {
390       EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
391     }
392 
393     MatrixType m_qr;
394     HCoeffsType m_hCoeffs;
395     PermutationType m_colsPermutation;
396     IntRowVectorType m_colsTranspositions;
397     RowVectorType m_temp;
398     RealRowVectorType m_colSqNorms;
399     bool m_isInitialized, m_usePrescribedThreshold;
400     RealScalar m_prescribedThreshold, m_maxpivot;
401     Index m_nonzero_pivots;
402     Index m_det_pq;
403 };
404 
405 template<typename MatrixType>
absDeterminant()406 typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::absDeterminant() const
407 {
408   using std::abs;
409   eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
410   eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
411   return abs(m_qr.diagonal().prod());
412 }
413 
414 template<typename MatrixType>
logAbsDeterminant()415 typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::logAbsDeterminant() const
416 {
417   eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
418   eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
419   return m_qr.diagonal().cwiseAbs().array().log().sum();
420 }
421 
422 /** Performs the QR factorization of the given matrix \a matrix. The result of
423   * the factorization is stored into \c *this, and a reference to \c *this
424   * is returned.
425   *
426   * \sa class ColPivHouseholderQR, ColPivHouseholderQR(const MatrixType&)
427   */
428 template<typename MatrixType>
compute(const MatrixType & matrix)429 ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const MatrixType& matrix)
430 {
431   check_template_parameters();
432 
433   using std::abs;
434   Index rows = matrix.rows();
435   Index cols = matrix.cols();
436   Index size = matrix.diagonalSize();
437 
438   // the column permutation is stored as int indices, so just to be sure:
439   eigen_assert(cols<=NumTraits<int>::highest());
440 
441   m_qr = matrix;
442   m_hCoeffs.resize(size);
443 
444   m_temp.resize(cols);
445 
446   m_colsTranspositions.resize(matrix.cols());
447   Index number_of_transpositions = 0;
448 
449   m_colSqNorms.resize(cols);
450   for(Index k = 0; k < cols; ++k)
451     m_colSqNorms.coeffRef(k) = m_qr.col(k).squaredNorm();
452 
453   RealScalar threshold_helper = m_colSqNorms.maxCoeff() * numext::abs2(NumTraits<Scalar>::epsilon()) / RealScalar(rows);
454 
455   m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
456   m_maxpivot = RealScalar(0);
457 
458   for(Index k = 0; k < size; ++k)
459   {
460     // first, we look up in our table m_colSqNorms which column has the biggest squared norm
461     Index biggest_col_index;
462     RealScalar biggest_col_sq_norm = m_colSqNorms.tail(cols-k).maxCoeff(&biggest_col_index);
463     biggest_col_index += k;
464 
465     // since our table m_colSqNorms accumulates imprecision at every step, we must now recompute
466     // the actual squared norm of the selected column.
467     // Note that not doing so does result in solve() sometimes returning inf/nan values
468     // when running the unit test with 1000 repetitions.
469     biggest_col_sq_norm = m_qr.col(biggest_col_index).tail(rows-k).squaredNorm();
470 
471     // we store that back into our table: it can't hurt to correct our table.
472     m_colSqNorms.coeffRef(biggest_col_index) = biggest_col_sq_norm;
473 
474     // Track the number of meaningful pivots but do not stop the decomposition to make
475     // sure that the initial matrix is properly reproduced. See bug 941.
476     if(m_nonzero_pivots==size && biggest_col_sq_norm < threshold_helper * RealScalar(rows-k))
477       m_nonzero_pivots = k;
478 
479     // apply the transposition to the columns
480     m_colsTranspositions.coeffRef(k) = biggest_col_index;
481     if(k != biggest_col_index) {
482       m_qr.col(k).swap(m_qr.col(biggest_col_index));
483       std::swap(m_colSqNorms.coeffRef(k), m_colSqNorms.coeffRef(biggest_col_index));
484       ++number_of_transpositions;
485     }
486 
487     // generate the householder vector, store it below the diagonal
488     RealScalar beta;
489     m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta);
490 
491     // apply the householder transformation to the diagonal coefficient
492     m_qr.coeffRef(k,k) = beta;
493 
494     // remember the maximum absolute value of diagonal coefficients
495     if(abs(beta) > m_maxpivot) m_maxpivot = abs(beta);
496 
497     // apply the householder transformation
498     m_qr.bottomRightCorner(rows-k, cols-k-1)
499         .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1));
500 
501     // update our table of squared norms of the columns
502     m_colSqNorms.tail(cols-k-1) -= m_qr.row(k).tail(cols-k-1).cwiseAbs2();
503   }
504 
505   m_colsPermutation.setIdentity(PermIndexType(cols));
506   for(PermIndexType k = 0; k < size/*m_nonzero_pivots*/; ++k)
507     m_colsPermutation.applyTranspositionOnTheRight(k, PermIndexType(m_colsTranspositions.coeff(k)));
508 
509   m_det_pq = (number_of_transpositions%2) ? -1 : 1;
510   m_isInitialized = true;
511 
512   return *this;
513 }
514 
515 namespace internal {
516 
517 template<typename _MatrixType, typename Rhs>
518 struct solve_retval<ColPivHouseholderQR<_MatrixType>, Rhs>
519   : solve_retval_base<ColPivHouseholderQR<_MatrixType>, Rhs>
520 {
521   EIGEN_MAKE_SOLVE_HELPERS(ColPivHouseholderQR<_MatrixType>,Rhs)
522 
523   template<typename Dest> void evalTo(Dest& dst) const
524   {
525     eigen_assert(rhs().rows() == dec().rows());
526 
527     const Index cols = dec().cols(),
528 				nonzero_pivots = dec().nonzeroPivots();
529 
530     if(nonzero_pivots == 0)
531     {
532       dst.setZero();
533       return;
534     }
535 
536     typename Rhs::PlainObject c(rhs());
537 
538     // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
539     c.applyOnTheLeft(householderSequence(dec().matrixQR(), dec().hCoeffs())
540                      .setLength(dec().nonzeroPivots())
541 		     .transpose()
542       );
543 
544     dec().matrixR()
545        .topLeftCorner(nonzero_pivots, nonzero_pivots)
546        .template triangularView<Upper>()
547        .solveInPlace(c.topRows(nonzero_pivots));
548 
549     for(Index i = 0; i < nonzero_pivots; ++i) dst.row(dec().colsPermutation().indices().coeff(i)) = c.row(i);
550     for(Index i = nonzero_pivots; i < cols; ++i) dst.row(dec().colsPermutation().indices().coeff(i)).setZero();
551   }
552 };
553 
554 } // end namespace internal
555 
556 /** \returns the matrix Q as a sequence of householder transformations.
557   * You can extract the meaningful part only by using:
558   * \code qr.householderQ().setLength(qr.nonzeroPivots()) \endcode*/
559 template<typename MatrixType>
560 typename ColPivHouseholderQR<MatrixType>::HouseholderSequenceType ColPivHouseholderQR<MatrixType>
561   ::householderQ() const
562 {
563   eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
564   return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate());
565 }
566 
567 /** \return the column-pivoting Householder QR decomposition of \c *this.
568   *
569   * \sa class ColPivHouseholderQR
570   */
571 template<typename Derived>
572 const ColPivHouseholderQR<typename MatrixBase<Derived>::PlainObject>
573 MatrixBase<Derived>::colPivHouseholderQr() const
574 {
575   return ColPivHouseholderQR<PlainObject>(eval());
576 }
577 
578 } // end namespace Eigen
579 
580 #endif // EIGEN_COLPIVOTINGHOUSEHOLDERQR_H
581