Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 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; |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 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: |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 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(), |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 79 | m_isInitialized(false), |
| 80 | m_usePrescribedThreshold(false) {} |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 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)), |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 91 | m_colsPermutation(PermIndexType(cols)), |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 92 | m_colsTranspositions(cols), |
| 93 | m_temp(cols), |
| 94 | m_colSqNorms(cols), |
| 95 | m_isInitialized(false), |
| 96 | m_usePrescribedThreshold(false) {} |
| 97 | |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 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 | */ |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 110 | ColPivHouseholderQR(const MatrixType& matrix) |
| 111 | : m_qr(matrix.rows(), matrix.cols()), |
| 112 | m_hCoeffs((std::min)(matrix.rows(),matrix.cols())), |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 113 | m_colsPermutation(PermIndexType(matrix.cols())), |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 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; |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 149 | HouseholderSequenceType matrixQ(void) const |
| 150 | { |
| 151 | return householderQ(); |
| 152 | } |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 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 | } |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 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 | |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 177 | ColPivHouseholderQR& compute(const MatrixType& matrix); |
| 178 | |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 179 | /** \returns a const reference to the column permutation matrix */ |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 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 | { |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 223 | using std::abs; |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 224 | eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 225 | RealScalar premultiplied_threshold = abs(m_maxpivot) * threshold(); |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 226 | Index result = 0; |
| 227 | for(Index i = 0; i < m_nonzero_pivots; ++i) |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 228 | result += (abs(m_qr.coeff(i,i)) > premultiplied_threshold); |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 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(); } |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 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 | */ |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 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. |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 353 | : NumTraits<Scalar>::epsilon() * RealScalar(m_qr.diagonalSize()); |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 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; } |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 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 | } |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 385 | |
| 386 | protected: |
| 387 | MatrixType m_qr; |
| 388 | HCoeffsType m_hCoeffs; |
| 389 | PermutationType m_colsPermutation; |
| 390 | IntRowVectorType m_colsTranspositions; |
| 391 | RowVectorType m_temp; |
| 392 | RealRowVectorType m_colSqNorms; |
| 393 | bool m_isInitialized, m_usePrescribedThreshold; |
| 394 | RealScalar m_prescribedThreshold, m_maxpivot; |
| 395 | Index m_nonzero_pivots; |
| 396 | Index m_det_pq; |
| 397 | }; |
| 398 | |
| 399 | template<typename MatrixType> |
| 400 | typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::absDeterminant() const |
| 401 | { |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 402 | using std::abs; |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 403 | eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); |
| 404 | eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!"); |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 405 | return abs(m_qr.diagonal().prod()); |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 406 | } |
| 407 | |
| 408 | template<typename MatrixType> |
| 409 | typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::logAbsDeterminant() const |
| 410 | { |
| 411 | eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); |
| 412 | eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!"); |
| 413 | return m_qr.diagonal().cwiseAbs().array().log().sum(); |
| 414 | } |
| 415 | |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 416 | /** Performs the QR factorization of the given matrix \a matrix. The result of |
| 417 | * the factorization is stored into \c *this, and a reference to \c *this |
| 418 | * is returned. |
| 419 | * |
| 420 | * \sa class ColPivHouseholderQR, ColPivHouseholderQR(const MatrixType&) |
| 421 | */ |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 422 | template<typename MatrixType> |
| 423 | ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const MatrixType& matrix) |
| 424 | { |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 425 | using std::abs; |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 426 | Index rows = matrix.rows(); |
| 427 | Index cols = matrix.cols(); |
| 428 | Index size = matrix.diagonalSize(); |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 429 | |
| 430 | // the column permutation is stored as int indices, so just to be sure: |
| 431 | eigen_assert(cols<=NumTraits<int>::highest()); |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 432 | |
| 433 | m_qr = matrix; |
| 434 | m_hCoeffs.resize(size); |
| 435 | |
| 436 | m_temp.resize(cols); |
| 437 | |
| 438 | m_colsTranspositions.resize(matrix.cols()); |
| 439 | Index number_of_transpositions = 0; |
| 440 | |
| 441 | m_colSqNorms.resize(cols); |
| 442 | for(Index k = 0; k < cols; ++k) |
| 443 | m_colSqNorms.coeffRef(k) = m_qr.col(k).squaredNorm(); |
| 444 | |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 445 | RealScalar threshold_helper = m_colSqNorms.maxCoeff() * numext::abs2(NumTraits<Scalar>::epsilon()) / RealScalar(rows); |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 446 | |
| 447 | m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case) |
| 448 | m_maxpivot = RealScalar(0); |
| 449 | |
| 450 | for(Index k = 0; k < size; ++k) |
| 451 | { |
| 452 | // first, we look up in our table m_colSqNorms which column has the biggest squared norm |
| 453 | Index biggest_col_index; |
| 454 | RealScalar biggest_col_sq_norm = m_colSqNorms.tail(cols-k).maxCoeff(&biggest_col_index); |
| 455 | biggest_col_index += k; |
| 456 | |
| 457 | // since our table m_colSqNorms accumulates imprecision at every step, we must now recompute |
| 458 | // the actual squared norm of the selected column. |
| 459 | // Note that not doing so does result in solve() sometimes returning inf/nan values |
| 460 | // when running the unit test with 1000 repetitions. |
| 461 | biggest_col_sq_norm = m_qr.col(biggest_col_index).tail(rows-k).squaredNorm(); |
| 462 | |
| 463 | // we store that back into our table: it can't hurt to correct our table. |
| 464 | m_colSqNorms.coeffRef(biggest_col_index) = biggest_col_sq_norm; |
| 465 | |
| 466 | // if the current biggest column is smaller than epsilon times the initial biggest column, |
| 467 | // terminate to avoid generating nan/inf values. |
| 468 | // Note that here, if we test instead for "biggest == 0", we get a failure every 1000 (or so) |
| 469 | // repetitions of the unit test, with the result of solve() filled with large values of the order |
| 470 | // of 1/(size*epsilon). |
| 471 | if(biggest_col_sq_norm < threshold_helper * RealScalar(rows-k)) |
| 472 | { |
| 473 | m_nonzero_pivots = k; |
| 474 | m_hCoeffs.tail(size-k).setZero(); |
| 475 | m_qr.bottomRightCorner(rows-k,cols-k) |
| 476 | .template triangularView<StrictlyLower>() |
| 477 | .setZero(); |
| 478 | break; |
| 479 | } |
| 480 | |
| 481 | // apply the transposition to the columns |
| 482 | m_colsTranspositions.coeffRef(k) = biggest_col_index; |
| 483 | if(k != biggest_col_index) { |
| 484 | m_qr.col(k).swap(m_qr.col(biggest_col_index)); |
| 485 | std::swap(m_colSqNorms.coeffRef(k), m_colSqNorms.coeffRef(biggest_col_index)); |
| 486 | ++number_of_transpositions; |
| 487 | } |
| 488 | |
| 489 | // generate the householder vector, store it below the diagonal |
| 490 | RealScalar beta; |
| 491 | m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta); |
| 492 | |
| 493 | // apply the householder transformation to the diagonal coefficient |
| 494 | m_qr.coeffRef(k,k) = beta; |
| 495 | |
| 496 | // remember the maximum absolute value of diagonal coefficients |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 497 | if(abs(beta) > m_maxpivot) m_maxpivot = abs(beta); |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 498 | |
| 499 | // apply the householder transformation |
| 500 | m_qr.bottomRightCorner(rows-k, cols-k-1) |
| 501 | .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1)); |
| 502 | |
| 503 | // update our table of squared norms of the columns |
| 504 | m_colSqNorms.tail(cols-k-1) -= m_qr.row(k).tail(cols-k-1).cwiseAbs2(); |
| 505 | } |
| 506 | |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 507 | m_colsPermutation.setIdentity(PermIndexType(cols)); |
| 508 | for(PermIndexType k = 0; k < m_nonzero_pivots; ++k) |
| 509 | m_colsPermutation.applyTranspositionOnTheRight(k, PermIndexType(m_colsTranspositions.coeff(k))); |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 510 | |
| 511 | m_det_pq = (number_of_transpositions%2) ? -1 : 1; |
| 512 | m_isInitialized = true; |
| 513 | |
| 514 | return *this; |
| 515 | } |
| 516 | |
| 517 | namespace internal { |
| 518 | |
| 519 | template<typename _MatrixType, typename Rhs> |
| 520 | struct solve_retval<ColPivHouseholderQR<_MatrixType>, Rhs> |
| 521 | : solve_retval_base<ColPivHouseholderQR<_MatrixType>, Rhs> |
| 522 | { |
| 523 | EIGEN_MAKE_SOLVE_HELPERS(ColPivHouseholderQR<_MatrixType>,Rhs) |
| 524 | |
| 525 | template<typename Dest> void evalTo(Dest& dst) const |
| 526 | { |
| 527 | eigen_assert(rhs().rows() == dec().rows()); |
| 528 | |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 529 | const Index cols = dec().cols(), |
| 530 | nonzero_pivots = dec().nonzeroPivots(); |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 531 | |
| 532 | if(nonzero_pivots == 0) |
| 533 | { |
| 534 | dst.setZero(); |
| 535 | return; |
| 536 | } |
| 537 | |
| 538 | typename Rhs::PlainObject c(rhs()); |
| 539 | |
| 540 | // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T |
| 541 | c.applyOnTheLeft(householderSequence(dec().matrixQR(), dec().hCoeffs()) |
| 542 | .setLength(dec().nonzeroPivots()) |
| 543 | .transpose() |
| 544 | ); |
| 545 | |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 546 | dec().matrixR() |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 547 | .topLeftCorner(nonzero_pivots, nonzero_pivots) |
| 548 | .template triangularView<Upper>() |
| 549 | .solveInPlace(c.topRows(nonzero_pivots)); |
| 550 | |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 551 | for(Index i = 0; i < nonzero_pivots; ++i) dst.row(dec().colsPermutation().indices().coeff(i)) = c.row(i); |
| 552 | for(Index i = nonzero_pivots; i < cols; ++i) dst.row(dec().colsPermutation().indices().coeff(i)).setZero(); |
| 553 | } |
| 554 | }; |
| 555 | |
| 556 | } // end namespace internal |
| 557 | |
| 558 | /** \returns the matrix Q as a sequence of householder transformations */ |
| 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()).setLength(m_nonzero_pivots); |
| 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 |