source: pacpussensors/trunk/Vislab/lib3dv/eigen/Eigen/src/QR/ColPivHouseholderQR.h@ 136

Last change on this file since 136 was 136, checked in by ldecherf, 7 years ago

Doc
File size: 21.7 KB
Line 
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 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 */
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 */
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 */
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>
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;
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 */
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 */
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 */
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 */
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 */
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 */
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 */
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 */
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>
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
296 inline Index rows() const { return m_qr.rows(); }
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 */
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 */
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 */
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 */
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 */
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 */
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 */
380 ComputationInfo info() const
381 {
382 eigen_assert(m_isInitialized && "Decomposition is not initialized.");
383 return Success;
384 }
385
386 protected:
387
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
405template<typename MatrixType>
406typename 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
414template<typename MatrixType>
415typename 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 */
428template<typename MatrixType>
429ColPivHouseholderQR<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
515namespace internal {
516
517template<typename _MatrixType, typename Rhs>
518struct 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*/
559template<typename MatrixType>
560typename 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 */
571template<typename Derived>
572const ColPivHouseholderQR<typename MatrixBase<Derived>::PlainObject>
573MatrixBase<Derived>::colPivHouseholderQr() const
574{
575 return ColPivHouseholderQR<PlainObject>(eval());
576}
577
578} // end namespace Eigen
579
580#endif // EIGEN_COLPIVOTINGHOUSEHOLDERQR_H
Note: See TracBrowser for help on using the repository browser.