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 * \tparam _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 * This class supports the \link InplaceDecomposition inplace decomposition \endlink mechanism. 41 * 42 * \sa MatrixBase::householderQr() 43 */ 44 template<typename _MatrixType> class HouseholderQR 45 { 46 public: 47 48 typedef _MatrixType MatrixType; 49 enum { 50 RowsAtCompileTime = MatrixType::RowsAtCompileTime, 51 ColsAtCompileTime = MatrixType::ColsAtCompileTime, 52 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 53 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime 54 }; 55 typedef typename MatrixType::Scalar Scalar; 56 typedef typename MatrixType::RealScalar RealScalar; 57 // FIXME should be int 58 typedef typename MatrixType::StorageIndex StorageIndex; 59 typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, (MatrixType::Flags&RowMajorBit) ? RowMajor : ColMajor, MaxRowsAtCompileTime, MaxRowsAtCompileTime> MatrixQType; 60 typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType; 61 typedef typename internal::plain_row_type<MatrixType>::type RowVectorType; 62 typedef HouseholderSequence<MatrixType,typename internal::remove_all<typename HCoeffsType::ConjugateReturnType>::type> HouseholderSequenceType; 63 64 /** 65 * \brief Default Constructor. 66 * 67 * The default constructor is useful in cases in which the user intends to 68 * perform decompositions via HouseholderQR::compute(const MatrixType&). 69 */ 70 HouseholderQR() : m_qr(), m_hCoeffs(), m_temp(), m_isInitialized(false) {} 71 72 /** \brief Default Constructor with memory preallocation 73 * 74 * Like the default constructor but with preallocation of the internal data 75 * according to the specified problem \a size. 76 * \sa HouseholderQR() 77 */ 78 HouseholderQR(Index rows, Index cols) 79 : m_qr(rows, cols), 80 m_hCoeffs((std::min)(rows,cols)), 81 m_temp(cols), 82 m_isInitialized(false) {} 83 84 /** \brief Constructs a QR factorization from a given matrix 85 * 86 * This constructor computes the QR factorization of the matrix \a matrix by calling 87 * the method compute(). It is a short cut for: 88 * 89 * \code 90 * HouseholderQR<MatrixType> qr(matrix.rows(), matrix.cols()); 91 * qr.compute(matrix); 92 * \endcode 93 * 94 * \sa compute() 95 */ 96 template<typename InputType> 97 explicit HouseholderQR(const EigenBase<InputType>& matrix) 98 : m_qr(matrix.rows(), matrix.cols()), 99 m_hCoeffs((std::min)(matrix.rows(),matrix.cols())), 100 m_temp(matrix.cols()), 101 m_isInitialized(false) 102 { 103 compute(matrix.derived()); 104 } 105 106 107 /** \brief Constructs a QR factorization from a given matrix 108 * 109 * This overloaded constructor is provided for \link InplaceDecomposition inplace decomposition \endlink when 110 * \c MatrixType is a Eigen::Ref. 111 * 112 * \sa HouseholderQR(const EigenBase&) 113 */ 114 template<typename InputType> 115 explicit HouseholderQR(EigenBase<InputType>& matrix) 116 : m_qr(matrix.derived()), 117 m_hCoeffs((std::min)(matrix.rows(),matrix.cols())), 118 m_temp(matrix.cols()), 119 m_isInitialized(false) 120 { 121 computeInPlace(); 122 } 123 124 /** This method finds a solution x to the equation Ax=b, where A is the matrix of which 125 * *this is the QR decomposition, if any exists. 126 * 127 * \param b the right-hand-side of the equation to solve. 128 * 129 * \returns a solution. 130 * 131 * \note_about_checking_solutions 132 * 133 * \note_about_arbitrary_choice_of_solution 134 * 135 * Example: \include HouseholderQR_solve.cpp 136 * Output: \verbinclude HouseholderQR_solve.out 137 */ 138 template<typename Rhs> 139 inline const Solve<HouseholderQR, Rhs> 140 solve(const MatrixBase<Rhs>& b) const 141 { 142 eigen_assert(m_isInitialized && "HouseholderQR is not initialized."); 143 return Solve<HouseholderQR, Rhs>(*this, b.derived()); 144 } 145 146 /** This method returns an expression of the unitary matrix Q as a sequence of Householder transformations. 147 * 148 * The returned expression can directly be used to perform matrix products. It can also be assigned to a dense Matrix object. 149 * Here is an example showing how to recover the full or thin matrix Q, as well as how to perform matrix products using operator*: 150 * 151 * Example: \include HouseholderQR_householderQ.cpp 152 * Output: \verbinclude HouseholderQR_householderQ.out 153 */ 154 HouseholderSequenceType householderQ() const 155 { 156 eigen_assert(m_isInitialized && "HouseholderQR is not initialized."); 157 return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate()); 158 } 159 160 /** \returns a reference to the matrix where the Householder QR decomposition is stored 161 * in a LAPACK-compatible way. 162 */ 163 const MatrixType& matrixQR() const 164 { 165 eigen_assert(m_isInitialized && "HouseholderQR is not initialized."); 166 return m_qr; 167 } 168 169 template<typename InputType> 170 HouseholderQR& compute(const EigenBase<InputType>& matrix) { 171 m_qr = matrix.derived(); 172 computeInPlace(); 173 return *this; 174 } 175 176 /** \returns the absolute value of the determinant of the matrix of which 177 * *this is the QR decomposition. It has only linear complexity 178 * (that is, O(n) where n is the dimension of the square matrix) 179 * as the QR decomposition has already been computed. 180 * 181 * \note This is only for square matrices. 182 * 183 * \warning a determinant can be very big or small, so for matrices 184 * of large enough dimension, there is a risk of overflow/underflow. 185 * One way to work around that is to use logAbsDeterminant() instead. 186 * 187 * \sa logAbsDeterminant(), MatrixBase::determinant() 188 */ 189 typename MatrixType::RealScalar absDeterminant() const; 190 191 /** \returns the natural log of the absolute value of the determinant of the matrix of which 192 * *this is the QR decomposition. It has only linear complexity 193 * (that is, O(n) where n is the dimension of the square matrix) 194 * as the QR decomposition has already been computed. 195 * 196 * \note This is only for square matrices. 197 * 198 * \note This method is useful to work around the risk of overflow/underflow that's inherent 199 * to determinant computation. 200 * 201 * \sa absDeterminant(), MatrixBase::determinant() 202 */ 203 typename MatrixType::RealScalar logAbsDeterminant() const; 204 205 inline Index rows() const { return m_qr.rows(); } 206 inline Index cols() const { return m_qr.cols(); } 207 208 /** \returns a const reference to the vector of Householder coefficients used to represent the factor \c Q. 209 * 210 * For advanced uses only. 211 */ 212 const HCoeffsType& hCoeffs() const { return m_hCoeffs; } 213 214 #ifndef EIGEN_PARSED_BY_DOXYGEN 215 template<typename RhsType, typename DstType> 216 EIGEN_DEVICE_FUNC 217 void _solve_impl(const RhsType &rhs, DstType &dst) const; 218 #endif 219 220 protected: 221 222 static void check_template_parameters() 223 { 224 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar); 225 } 226 227 void computeInPlace(); 228 229 MatrixType m_qr; 230 HCoeffsType m_hCoeffs; 231 RowVectorType m_temp; 232 bool m_isInitialized; 233 }; 234 235 template<typename MatrixType> 236 typename MatrixType::RealScalar HouseholderQR<MatrixType>::absDeterminant() const 237 { 238 using std::abs; 239 eigen_assert(m_isInitialized && "HouseholderQR is not initialized."); 240 eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!"); 241 return abs(m_qr.diagonal().prod()); 242 } 243 244 template<typename MatrixType> 245 typename MatrixType::RealScalar HouseholderQR<MatrixType>::logAbsDeterminant() const 246 { 247 eigen_assert(m_isInitialized && "HouseholderQR is not initialized."); 248 eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!"); 249 return m_qr.diagonal().cwiseAbs().array().log().sum(); 250 } 251 252 namespace internal { 253 254 /** \internal */ 255 template<typename MatrixQR, typename HCoeffs> 256 void householder_qr_inplace_unblocked(MatrixQR& mat, HCoeffs& hCoeffs, typename MatrixQR::Scalar* tempData = 0) 257 { 258 typedef typename MatrixQR::Scalar Scalar; 259 typedef typename MatrixQR::RealScalar RealScalar; 260 Index rows = mat.rows(); 261 Index cols = mat.cols(); 262 Index size = (std::min)(rows,cols); 263 264 eigen_assert(hCoeffs.size() == size); 265 266 typedef Matrix<Scalar,MatrixQR::ColsAtCompileTime,1> TempType; 267 TempType tempVector; 268 if(tempData==0) 269 { 270 tempVector.resize(cols); 271 tempData = tempVector.data(); 272 } 273 274 for(Index k = 0; k < size; ++k) 275 { 276 Index remainingRows = rows - k; 277 Index remainingCols = cols - k - 1; 278 279 RealScalar beta; 280 mat.col(k).tail(remainingRows).makeHouseholderInPlace(hCoeffs.coeffRef(k), beta); 281 mat.coeffRef(k,k) = beta; 282 283 // apply H to remaining part of m_qr from the left 284 mat.bottomRightCorner(remainingRows, remainingCols) 285 .applyHouseholderOnTheLeft(mat.col(k).tail(remainingRows-1), hCoeffs.coeffRef(k), tempData+k+1); 286 } 287 } 288 289 /** \internal */ 290 template<typename MatrixQR, typename HCoeffs, 291 typename MatrixQRScalar = typename MatrixQR::Scalar, 292 bool InnerStrideIsOne = (MatrixQR::InnerStrideAtCompileTime == 1 && HCoeffs::InnerStrideAtCompileTime == 1)> 293 struct householder_qr_inplace_blocked 294 { 295 // This is specialized for MKL-supported Scalar types in HouseholderQR_MKL.h 296 static void run(MatrixQR& mat, HCoeffs& hCoeffs, Index maxBlockSize=32, 297 typename MatrixQR::Scalar* tempData = 0) 298 { 299 typedef typename MatrixQR::Scalar Scalar; 300 typedef Block<MatrixQR,Dynamic,Dynamic> BlockType; 301 302 Index rows = mat.rows(); 303 Index cols = mat.cols(); 304 Index size = (std::min)(rows, cols); 305 306 typedef Matrix<Scalar,Dynamic,1,ColMajor,MatrixQR::MaxColsAtCompileTime,1> TempType; 307 TempType tempVector; 308 if(tempData==0) 309 { 310 tempVector.resize(cols); 311 tempData = tempVector.data(); 312 } 313 314 Index blockSize = (std::min)(maxBlockSize,size); 315 316 Index k = 0; 317 for (k = 0; k < size; k += blockSize) 318 { 319 Index bs = (std::min)(size-k,blockSize); // actual size of the block 320 Index tcols = cols - k - bs; // trailing columns 321 Index brows = rows-k; // rows of the block 322 323 // partition the matrix: 324 // A00 | A01 | A02 325 // mat = A10 | A11 | A12 326 // A20 | A21 | A22 327 // and performs the qr dec of [A11^T A12^T]^T 328 // and update [A21^T A22^T]^T using level 3 operations. 329 // Finally, the algorithm continue on A22 330 331 BlockType A11_21 = mat.block(k,k,brows,bs); 332 Block<HCoeffs,Dynamic,1> hCoeffsSegment = hCoeffs.segment(k,bs); 333 334 householder_qr_inplace_unblocked(A11_21, hCoeffsSegment, tempData); 335 336 if(tcols) 337 { 338 BlockType A21_22 = mat.block(k,k+bs,brows,tcols); 339 apply_block_householder_on_the_left(A21_22,A11_21,hCoeffsSegment, false); // false == backward 340 } 341 } 342 } 343 }; 344 345 } // end namespace internal 346 347 #ifndef EIGEN_PARSED_BY_DOXYGEN 348 template<typename _MatrixType> 349 template<typename RhsType, typename DstType> 350 void HouseholderQR<_MatrixType>::_solve_impl(const RhsType &rhs, DstType &dst) const 351 { 352 const Index rank = (std::min)(rows(), cols()); 353 eigen_assert(rhs.rows() == rows()); 354 355 typename RhsType::PlainObject c(rhs); 356 357 // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T 358 c.applyOnTheLeft(householderSequence( 359 m_qr.leftCols(rank), 360 m_hCoeffs.head(rank)).transpose() 361 ); 362 363 m_qr.topLeftCorner(rank, rank) 364 .template triangularView<Upper>() 365 .solveInPlace(c.topRows(rank)); 366 367 dst.topRows(rank) = c.topRows(rank); 368 dst.bottomRows(cols()-rank).setZero(); 369 } 370 #endif 371 372 /** Performs the QR factorization of the given matrix \a matrix. The result of 373 * the factorization is stored into \c *this, and a reference to \c *this 374 * is returned. 375 * 376 * \sa class HouseholderQR, HouseholderQR(const MatrixType&) 377 */ 378 template<typename MatrixType> 379 void HouseholderQR<MatrixType>::computeInPlace() 380 { 381 check_template_parameters(); 382 383 Index rows = m_qr.rows(); 384 Index cols = m_qr.cols(); 385 Index size = (std::min)(rows,cols); 386 387 m_hCoeffs.resize(size); 388 389 m_temp.resize(cols); 390 391 internal::householder_qr_inplace_blocked<MatrixType, HCoeffsType>::run(m_qr, m_hCoeffs, 48, m_temp.data()); 392 393 m_isInitialized = true; 394 } 395 396 /** \return the Householder QR decomposition of \c *this. 397 * 398 * \sa class HouseholderQR 399 */ 400 template<typename Derived> 401 const HouseholderQR<typename MatrixBase<Derived>::PlainObject> 402 MatrixBase<Derived>::householderQr() const 403 { 404 return HouseholderQR<PlainObject>(eval()); 405 } 406 407 } // end namespace Eigen 408 409 #endif // EIGEN_QR_H 410