1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2008-2010 Gael Guennebaud <gael.guennebaud (at) inria.fr> 5 // Copyright (C) 2009 Benoit Jacob <jacob.benoit.1 (at) gmail.com> 6 // Copyright (C) 2010 Vincent Lejeune 7 // 8 // This Source Code Form is subject to the terms of the Mozilla 9 // Public License v. 2.0. If a copy of the MPL was not distributed 10 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 11 12 #ifndef EIGEN_QR_H 13 #define EIGEN_QR_H 14 15 namespace Eigen { 16 17 /** \ingroup QR_Module 18 * 19 * 20 * \class HouseholderQR 21 * 22 * \brief Householder QR decomposition of a matrix 23 * 24 * \param MatrixType the type of the matrix of which we are computing the QR decomposition 25 * 26 * This class performs a QR decomposition of a matrix \b A into matrices \b Q and \b R 27 * such that 28 * \f[ 29 * \mathbf{A} = \mathbf{Q} \, \mathbf{R} 30 * \f] 31 * by using Householder transformations. Here, \b Q a unitary matrix and \b R an upper triangular matrix. 32 * The result is stored in a compact way compatible with LAPACK. 33 * 34 * Note that no pivoting is performed. This is \b not a rank-revealing decomposition. 35 * If you want that feature, use FullPivHouseholderQR or ColPivHouseholderQR instead. 36 * 37 * This Householder QR decomposition is faster, but less numerically stable and less feature-full than 38 * FullPivHouseholderQR or ColPivHouseholderQR. 39 * 40 * \sa MatrixBase::householderQr() 41 */ 42 template<typename _MatrixType> class HouseholderQR 43 { 44 public: 45 46 typedef _MatrixType MatrixType; 47 enum { 48 RowsAtCompileTime = MatrixType::RowsAtCompileTime, 49 ColsAtCompileTime = MatrixType::ColsAtCompileTime, 50 Options = MatrixType::Options, 51 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 52 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime 53 }; 54 typedef typename MatrixType::Scalar Scalar; 55 typedef typename MatrixType::RealScalar RealScalar; 56 typedef typename MatrixType::Index Index; 57 typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, (MatrixType::Flags&RowMajorBit) ? RowMajor : ColMajor, MaxRowsAtCompileTime, MaxRowsAtCompileTime> MatrixQType; 58 typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType; 59 typedef typename internal::plain_row_type<MatrixType>::type RowVectorType; 60 typedef HouseholderSequence<MatrixType,typename internal::remove_all<typename HCoeffsType::ConjugateReturnType>::type> HouseholderSequenceType; 61 62 /** 63 * \brief Default Constructor. 64 * 65 * The default constructor is useful in cases in which the user intends to 66 * perform decompositions via HouseholderQR::compute(const MatrixType&). 67 */ 68 HouseholderQR() : m_qr(), m_hCoeffs(), m_temp(), m_isInitialized(false) {} 69 70 /** \brief Default Constructor with memory preallocation 71 * 72 * Like the default constructor but with preallocation of the internal data 73 * according to the specified problem \a size. 74 * \sa HouseholderQR() 75 */ 76 HouseholderQR(Index rows, Index cols) 77 : m_qr(rows, cols), 78 m_hCoeffs((std::min)(rows,cols)), 79 m_temp(cols), 80 m_isInitialized(false) {} 81 82 /** \brief Constructs a QR factorization from a given matrix 83 * 84 * This constructor computes the QR factorization of the matrix \a matrix by calling 85 * the method compute(). It is a short cut for: 86 * 87 * \code 88 * HouseholderQR<MatrixType> qr(matrix.rows(), matrix.cols()); 89 * qr.compute(matrix); 90 * \endcode 91 * 92 * \sa compute() 93 */ 94 HouseholderQR(const MatrixType& matrix) 95 : m_qr(matrix.rows(), matrix.cols()), 96 m_hCoeffs((std::min)(matrix.rows(),matrix.cols())), 97 m_temp(matrix.cols()), 98 m_isInitialized(false) 99 { 100 compute(matrix); 101 } 102 103 /** This method finds a solution x to the equation Ax=b, where A is the matrix of which 104 * *this is the QR decomposition, if any exists. 105 * 106 * \param b the right-hand-side of the equation to solve. 107 * 108 * \returns a solution. 109 * 110 * \note The case where b is a matrix is not yet implemented. Also, this 111 * code is space inefficient. 112 * 113 * \note_about_checking_solutions 114 * 115 * \note_about_arbitrary_choice_of_solution 116 * 117 * Example: \include HouseholderQR_solve.cpp 118 * Output: \verbinclude HouseholderQR_solve.out 119 */ 120 template<typename Rhs> 121 inline const internal::solve_retval<HouseholderQR, Rhs> 122 solve(const MatrixBase<Rhs>& b) const 123 { 124 eigen_assert(m_isInitialized && "HouseholderQR is not initialized."); 125 return internal::solve_retval<HouseholderQR, Rhs>(*this, b.derived()); 126 } 127 128 /** This method returns an expression of the unitary matrix Q as a sequence of Householder transformations. 129 * 130 * The returned expression can directly be used to perform matrix products. It can also be assigned to a dense Matrix object. 131 * Here is an example showing how to recover the full or thin matrix Q, as well as how to perform matrix products using operator*: 132 * 133 * Example: \include HouseholderQR_householderQ.cpp 134 * Output: \verbinclude HouseholderQR_householderQ.out 135 */ 136 HouseholderSequenceType householderQ() const 137 { 138 eigen_assert(m_isInitialized && "HouseholderQR is not initialized."); 139 return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate()); 140 } 141 142 /** \returns a reference to the matrix where the Householder QR decomposition is stored 143 * in a LAPACK-compatible way. 144 */ 145 const MatrixType& matrixQR() const 146 { 147 eigen_assert(m_isInitialized && "HouseholderQR is not initialized."); 148 return m_qr; 149 } 150 151 HouseholderQR& compute(const MatrixType& matrix); 152 153 /** \returns the absolute value of the determinant of the matrix of which 154 * *this is the QR decomposition. It has only linear complexity 155 * (that is, O(n) where n is the dimension of the square matrix) 156 * as the QR decomposition has already been computed. 157 * 158 * \note This is only for square matrices. 159 * 160 * \warning a determinant can be very big or small, so for matrices 161 * of large enough dimension, there is a risk of overflow/underflow. 162 * One way to work around that is to use logAbsDeterminant() instead. 163 * 164 * \sa logAbsDeterminant(), MatrixBase::determinant() 165 */ 166 typename MatrixType::RealScalar absDeterminant() const; 167 168 /** \returns the natural log of the absolute value of the determinant of the matrix of which 169 * *this is the QR decomposition. It has only linear complexity 170 * (that is, O(n) where n is the dimension of the square matrix) 171 * as the QR decomposition has already been computed. 172 * 173 * \note This is only for square matrices. 174 * 175 * \note This method is useful to work around the risk of overflow/underflow that's inherent 176 * to determinant computation. 177 * 178 * \sa absDeterminant(), MatrixBase::determinant() 179 */ 180 typename MatrixType::RealScalar logAbsDeterminant() const; 181 182 inline Index rows() const { return m_qr.rows(); } 183 inline Index cols() const { return m_qr.cols(); } 184 185 /** \returns a const reference to the vector of Householder coefficients used to represent the factor \c Q. 186 * 187 * For advanced uses only. 188 */ 189 const HCoeffsType& hCoeffs() const { return m_hCoeffs; } 190 191 protected: 192 MatrixType m_qr; 193 HCoeffsType m_hCoeffs; 194 RowVectorType m_temp; 195 bool m_isInitialized; 196 }; 197 198 template<typename MatrixType> 199 typename MatrixType::RealScalar HouseholderQR<MatrixType>::absDeterminant() const 200 { 201 using std::abs; 202 eigen_assert(m_isInitialized && "HouseholderQR is not initialized."); 203 eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!"); 204 return abs(m_qr.diagonal().prod()); 205 } 206 207 template<typename MatrixType> 208 typename MatrixType::RealScalar HouseholderQR<MatrixType>::logAbsDeterminant() const 209 { 210 eigen_assert(m_isInitialized && "HouseholderQR is not initialized."); 211 eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!"); 212 return m_qr.diagonal().cwiseAbs().array().log().sum(); 213 } 214 215 namespace internal { 216 217 /** \internal */ 218 template<typename MatrixQR, typename HCoeffs> 219 void householder_qr_inplace_unblocked(MatrixQR& mat, HCoeffs& hCoeffs, typename MatrixQR::Scalar* tempData = 0) 220 { 221 typedef typename MatrixQR::Index Index; 222 typedef typename MatrixQR::Scalar Scalar; 223 typedef typename MatrixQR::RealScalar RealScalar; 224 Index rows = mat.rows(); 225 Index cols = mat.cols(); 226 Index size = (std::min)(rows,cols); 227 228 eigen_assert(hCoeffs.size() == size); 229 230 typedef Matrix<Scalar,MatrixQR::ColsAtCompileTime,1> TempType; 231 TempType tempVector; 232 if(tempData==0) 233 { 234 tempVector.resize(cols); 235 tempData = tempVector.data(); 236 } 237 238 for(Index k = 0; k < size; ++k) 239 { 240 Index remainingRows = rows - k; 241 Index remainingCols = cols - k - 1; 242 243 RealScalar beta; 244 mat.col(k).tail(remainingRows).makeHouseholderInPlace(hCoeffs.coeffRef(k), beta); 245 mat.coeffRef(k,k) = beta; 246 247 // apply H to remaining part of m_qr from the left 248 mat.bottomRightCorner(remainingRows, remainingCols) 249 .applyHouseholderOnTheLeft(mat.col(k).tail(remainingRows-1), hCoeffs.coeffRef(k), tempData+k+1); 250 } 251 } 252 253 /** \internal */ 254 template<typename MatrixQR, typename HCoeffs> 255 void householder_qr_inplace_blocked(MatrixQR& mat, HCoeffs& hCoeffs, 256 typename MatrixQR::Index maxBlockSize=32, 257 typename MatrixQR::Scalar* tempData = 0) 258 { 259 typedef typename MatrixQR::Index Index; 260 typedef typename MatrixQR::Scalar Scalar; 261 typedef Block<MatrixQR,Dynamic,Dynamic> BlockType; 262 263 Index rows = mat.rows(); 264 Index cols = mat.cols(); 265 Index size = (std::min)(rows, cols); 266 267 typedef Matrix<Scalar,Dynamic,1,ColMajor,MatrixQR::MaxColsAtCompileTime,1> TempType; 268 TempType tempVector; 269 if(tempData==0) 270 { 271 tempVector.resize(cols); 272 tempData = tempVector.data(); 273 } 274 275 Index blockSize = (std::min)(maxBlockSize,size); 276 277 Index k = 0; 278 for (k = 0; k < size; k += blockSize) 279 { 280 Index bs = (std::min)(size-k,blockSize); // actual size of the block 281 Index tcols = cols - k - bs; // trailing columns 282 Index brows = rows-k; // rows of the block 283 284 // partition the matrix: 285 // A00 | A01 | A02 286 // mat = A10 | A11 | A12 287 // A20 | A21 | A22 288 // and performs the qr dec of [A11^T A12^T]^T 289 // and update [A21^T A22^T]^T using level 3 operations. 290 // Finally, the algorithm continue on A22 291 292 BlockType A11_21 = mat.block(k,k,brows,bs); 293 Block<HCoeffs,Dynamic,1> hCoeffsSegment = hCoeffs.segment(k,bs); 294 295 householder_qr_inplace_unblocked(A11_21, hCoeffsSegment, tempData); 296 297 if(tcols) 298 { 299 BlockType A21_22 = mat.block(k,k+bs,brows,tcols); 300 apply_block_householder_on_the_left(A21_22,A11_21,hCoeffsSegment.adjoint()); 301 } 302 } 303 } 304 305 template<typename _MatrixType, typename Rhs> 306 struct solve_retval<HouseholderQR<_MatrixType>, Rhs> 307 : solve_retval_base<HouseholderQR<_MatrixType>, Rhs> 308 { 309 EIGEN_MAKE_SOLVE_HELPERS(HouseholderQR<_MatrixType>,Rhs) 310 311 template<typename Dest> void evalTo(Dest& dst) const 312 { 313 const Index rows = dec().rows(), cols = dec().cols(); 314 const Index rank = (std::min)(rows, cols); 315 eigen_assert(rhs().rows() == rows); 316 317 typename Rhs::PlainObject c(rhs()); 318 319 // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T 320 c.applyOnTheLeft(householderSequence( 321 dec().matrixQR().leftCols(rank), 322 dec().hCoeffs().head(rank)).transpose() 323 ); 324 325 dec().matrixQR() 326 .topLeftCorner(rank, rank) 327 .template triangularView<Upper>() 328 .solveInPlace(c.topRows(rank)); 329 330 dst.topRows(rank) = c.topRows(rank); 331 dst.bottomRows(cols-rank).setZero(); 332 } 333 }; 334 335 } // end namespace internal 336 337 /** Performs the QR factorization of the given matrix \a matrix. The result of 338 * the factorization is stored into \c *this, and a reference to \c *this 339 * is returned. 340 * 341 * \sa class HouseholderQR, HouseholderQR(const MatrixType&) 342 */ 343 template<typename MatrixType> 344 HouseholderQR<MatrixType>& HouseholderQR<MatrixType>::compute(const MatrixType& matrix) 345 { 346 Index rows = matrix.rows(); 347 Index cols = matrix.cols(); 348 Index size = (std::min)(rows,cols); 349 350 m_qr = matrix; 351 m_hCoeffs.resize(size); 352 353 m_temp.resize(cols); 354 355 internal::householder_qr_inplace_blocked(m_qr, m_hCoeffs, 48, m_temp.data()); 356 357 m_isInitialized = true; 358 return *this; 359 } 360 361 /** \return the Householder QR decomposition of \c *this. 362 * 363 * \sa class HouseholderQR 364 */ 365 template<typename Derived> 366 const HouseholderQR<typename MatrixBase<Derived>::PlainObject> 367 MatrixBase<Derived>::householderQr() const 368 { 369 return HouseholderQR<PlainObject>(eval()); 370 } 371 372 } // end namespace Eigen 373 374 #endif // EIGEN_QR_H 375