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