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 193 static void check_template_parameters() 194 { 195 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar); 196 } 197 198 MatrixType m_qr; 199 HCoeffsType m_hCoeffs; 200 RowVectorType m_temp; 201 bool m_isInitialized; 202 }; 203 204 template<typename MatrixType> 205 typename MatrixType::RealScalar HouseholderQR<MatrixType>::absDeterminant() const 206 { 207 using std::abs; 208 eigen_assert(m_isInitialized && "HouseholderQR is not initialized."); 209 eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!"); 210 return abs(m_qr.diagonal().prod()); 211 } 212 213 template<typename MatrixType> 214 typename MatrixType::RealScalar HouseholderQR<MatrixType>::logAbsDeterminant() const 215 { 216 eigen_assert(m_isInitialized && "HouseholderQR is not initialized."); 217 eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!"); 218 return m_qr.diagonal().cwiseAbs().array().log().sum(); 219 } 220 221 namespace internal { 222 223 /** \internal */ 224 template<typename MatrixQR, typename HCoeffs> 225 void householder_qr_inplace_unblocked(MatrixQR& mat, HCoeffs& hCoeffs, typename MatrixQR::Scalar* tempData = 0) 226 { 227 typedef typename MatrixQR::Index Index; 228 typedef typename MatrixQR::Scalar Scalar; 229 typedef typename MatrixQR::RealScalar RealScalar; 230 Index rows = mat.rows(); 231 Index cols = mat.cols(); 232 Index size = (std::min)(rows,cols); 233 234 eigen_assert(hCoeffs.size() == size); 235 236 typedef Matrix<Scalar,MatrixQR::ColsAtCompileTime,1> TempType; 237 TempType tempVector; 238 if(tempData==0) 239 { 240 tempVector.resize(cols); 241 tempData = tempVector.data(); 242 } 243 244 for(Index k = 0; k < size; ++k) 245 { 246 Index remainingRows = rows - k; 247 Index remainingCols = cols - k - 1; 248 249 RealScalar beta; 250 mat.col(k).tail(remainingRows).makeHouseholderInPlace(hCoeffs.coeffRef(k), beta); 251 mat.coeffRef(k,k) = beta; 252 253 // apply H to remaining part of m_qr from the left 254 mat.bottomRightCorner(remainingRows, remainingCols) 255 .applyHouseholderOnTheLeft(mat.col(k).tail(remainingRows-1), hCoeffs.coeffRef(k), tempData+k+1); 256 } 257 } 258 259 /** \internal */ 260 template<typename MatrixQR, typename HCoeffs, 261 typename MatrixQRScalar = typename MatrixQR::Scalar, 262 bool InnerStrideIsOne = (MatrixQR::InnerStrideAtCompileTime == 1 && HCoeffs::InnerStrideAtCompileTime == 1)> 263 struct householder_qr_inplace_blocked 264 { 265 // This is specialized for MKL-supported Scalar types in HouseholderQR_MKL.h 266 static void run(MatrixQR& mat, HCoeffs& hCoeffs, 267 typename MatrixQR::Index maxBlockSize=32, 268 typename MatrixQR::Scalar* tempData = 0) 269 { 270 typedef typename MatrixQR::Index Index; 271 typedef typename MatrixQR::Scalar Scalar; 272 typedef Block<MatrixQR,Dynamic,Dynamic> BlockType; 273 274 Index rows = mat.rows(); 275 Index cols = mat.cols(); 276 Index size = (std::min)(rows, cols); 277 278 typedef Matrix<Scalar,Dynamic,1,ColMajor,MatrixQR::MaxColsAtCompileTime,1> TempType; 279 TempType tempVector; 280 if(tempData==0) 281 { 282 tempVector.resize(cols); 283 tempData = tempVector.data(); 284 } 285 286 Index blockSize = (std::min)(maxBlockSize,size); 287 288 Index k = 0; 289 for (k = 0; k < size; k += blockSize) 290 { 291 Index bs = (std::min)(size-k,blockSize); // actual size of the block 292 Index tcols = cols - k - bs; // trailing columns 293 Index brows = rows-k; // rows of the block 294 295 // partition the matrix: 296 // A00 | A01 | A02 297 // mat = A10 | A11 | A12 298 // A20 | A21 | A22 299 // and performs the qr dec of [A11^T A12^T]^T 300 // and update [A21^T A22^T]^T using level 3 operations. 301 // Finally, the algorithm continue on A22 302 303 BlockType A11_21 = mat.block(k,k,brows,bs); 304 Block<HCoeffs,Dynamic,1> hCoeffsSegment = hCoeffs.segment(k,bs); 305 306 householder_qr_inplace_unblocked(A11_21, hCoeffsSegment, tempData); 307 308 if(tcols) 309 { 310 BlockType A21_22 = mat.block(k,k+bs,brows,tcols); 311 apply_block_householder_on_the_left(A21_22,A11_21,hCoeffsSegment.adjoint()); 312 } 313 } 314 } 315 }; 316 317 template<typename _MatrixType, typename Rhs> 318 struct solve_retval<HouseholderQR<_MatrixType>, Rhs> 319 : solve_retval_base<HouseholderQR<_MatrixType>, Rhs> 320 { 321 EIGEN_MAKE_SOLVE_HELPERS(HouseholderQR<_MatrixType>,Rhs) 322 323 template<typename Dest> void evalTo(Dest& dst) const 324 { 325 const Index rows = dec().rows(), cols = dec().cols(); 326 const Index rank = (std::min)(rows, cols); 327 eigen_assert(rhs().rows() == rows); 328 329 typename Rhs::PlainObject c(rhs()); 330 331 // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T 332 c.applyOnTheLeft(householderSequence( 333 dec().matrixQR().leftCols(rank), 334 dec().hCoeffs().head(rank)).transpose() 335 ); 336 337 dec().matrixQR() 338 .topLeftCorner(rank, rank) 339 .template triangularView<Upper>() 340 .solveInPlace(c.topRows(rank)); 341 342 dst.topRows(rank) = c.topRows(rank); 343 dst.bottomRows(cols-rank).setZero(); 344 } 345 }; 346 347 } // end namespace internal 348 349 /** Performs the QR factorization of the given matrix \a matrix. The result of 350 * the factorization is stored into \c *this, and a reference to \c *this 351 * is returned. 352 * 353 * \sa class HouseholderQR, HouseholderQR(const MatrixType&) 354 */ 355 template<typename MatrixType> 356 HouseholderQR<MatrixType>& HouseholderQR<MatrixType>::compute(const MatrixType& matrix) 357 { 358 check_template_parameters(); 359 360 Index rows = matrix.rows(); 361 Index cols = matrix.cols(); 362 Index size = (std::min)(rows,cols); 363 364 m_qr = matrix; 365 m_hCoeffs.resize(size); 366 367 m_temp.resize(cols); 368 369 internal::householder_qr_inplace_blocked<MatrixType, HCoeffsType>::run(m_qr, m_hCoeffs, 48, m_temp.data()); 370 371 m_isInitialized = true; 372 return *this; 373 } 374 375 /** \return the Householder QR decomposition of \c *this. 376 * 377 * \sa class HouseholderQR 378 */ 379 template<typename Derived> 380 const HouseholderQR<typename MatrixBase<Derived>::PlainObject> 381 MatrixBase<Derived>::householderQr() const 382 { 383 return HouseholderQR<PlainObject>(eval()); 384 } 385 386 } // end namespace Eigen 387 388 #endif // EIGEN_QR_H 389
