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