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