Eigen  3.4.90 (git rev 67eeba6e720c5745abc77ae6c92ce0a44aa7b7ae)
HouseholderQR.h
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 #include "./InternalHeaderCheck.h"
16 
17 namespace Eigen {
18 
19 namespace internal {
20 template<typename MatrixType_> struct traits<HouseholderQR<MatrixType_> >
21  : traits<MatrixType_>
22 {
23  typedef MatrixXpr XprKind;
24  typedef SolverStorage StorageKind;
25  typedef int StorageIndex;
26  enum { Flags = 0 };
27 };
28 
29 } // end namespace internal
30 
58 template<typename MatrixType_> class HouseholderQR
59  : public SolverBase<HouseholderQR<MatrixType_> >
60 {
61  public:
62 
63  typedef MatrixType_ MatrixType;
65  friend class SolverBase<HouseholderQR>;
66 
67  EIGEN_GENERIC_PUBLIC_INTERFACE(HouseholderQR)
68  enum {
69  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
70  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
71  };
72  typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, (MatrixType::Flags&RowMajorBit) ? RowMajor : ColMajor, MaxRowsAtCompileTime, MaxRowsAtCompileTime> MatrixQType;
73  typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
74  typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
76 
83  HouseholderQR() : m_qr(), m_hCoeffs(), m_temp(), m_isInitialized(false) {}
84 
92  : m_qr(rows, cols),
93  m_hCoeffs((std::min)(rows,cols)),
94  m_temp(cols),
95  m_isInitialized(false) {}
96 
109  template<typename InputType>
110  explicit HouseholderQR(const EigenBase<InputType>& matrix)
111  : m_qr(matrix.rows(), matrix.cols()),
112  m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
113  m_temp(matrix.cols()),
114  m_isInitialized(false)
115  {
116  compute(matrix.derived());
117  }
118 
119 
127  template<typename InputType>
129  : m_qr(matrix.derived()),
130  m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
131  m_temp(matrix.cols()),
132  m_isInitialized(false)
133  {
134  computeInPlace();
135  }
136 
137  #ifdef EIGEN_PARSED_BY_DOXYGEN
152  template<typename Rhs>
153  inline const Solve<HouseholderQR, Rhs>
154  solve(const MatrixBase<Rhs>& b) const;
155  #endif
156 
166  {
167  eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
168  return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate());
169  }
170 
174  const MatrixType& matrixQR() const
175  {
176  eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
177  return m_qr;
178  }
179 
180  template<typename InputType>
181  HouseholderQR& compute(const EigenBase<InputType>& matrix) {
182  m_qr = matrix.derived();
183  computeInPlace();
184  return *this;
185  }
186 
200  typename MatrixType::RealScalar absDeterminant() const;
201 
214  typename MatrixType::RealScalar logAbsDeterminant() const;
215 
216  inline Index rows() const { return m_qr.rows(); }
217  inline Index cols() const { return m_qr.cols(); }
218 
223  const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
224 
225  #ifndef EIGEN_PARSED_BY_DOXYGEN
226  template<typename RhsType, typename DstType>
227  void _solve_impl(const RhsType &rhs, DstType &dst) const;
228 
229  template<bool Conjugate, typename RhsType, typename DstType>
230  void _solve_impl_transposed(const RhsType &rhs, DstType &dst) const;
231  #endif
232 
233  protected:
234 
235  EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar)
236 
237  void computeInPlace();
238 
239  MatrixType m_qr;
240  HCoeffsType m_hCoeffs;
241  RowVectorType m_temp;
242  bool m_isInitialized;
243 };
244 
245 template<typename MatrixType>
246 typename MatrixType::RealScalar HouseholderQR<MatrixType>::absDeterminant() const
247 {
248  using std::abs;
249  eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
250  eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
251  return abs(m_qr.diagonal().prod());
252 }
253 
254 template<typename MatrixType>
255 typename MatrixType::RealScalar HouseholderQR<MatrixType>::logAbsDeterminant() const
256 {
257  eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
258  eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
259  return m_qr.diagonal().cwiseAbs().array().log().sum();
260 }
261 
262 namespace internal {
263 
265 template<typename MatrixQR, typename HCoeffs>
266 void householder_qr_inplace_unblocked(MatrixQR& mat, HCoeffs& hCoeffs, typename MatrixQR::Scalar* tempData = 0)
267 {
268  typedef typename MatrixQR::Scalar Scalar;
269  typedef typename MatrixQR::RealScalar RealScalar;
270  Index rows = mat.rows();
271  Index cols = mat.cols();
272  Index size = (std::min)(rows,cols);
273 
274  eigen_assert(hCoeffs.size() == size);
275 
277  TempType tempVector;
278  if(tempData==0)
279  {
280  tempVector.resize(cols);
281  tempData = tempVector.data();
282  }
283 
284  for(Index k = 0; k < size; ++k)
285  {
286  Index remainingRows = rows - k;
287  Index remainingCols = cols - k - 1;
288 
289  RealScalar beta;
290  mat.col(k).tail(remainingRows).makeHouseholderInPlace(hCoeffs.coeffRef(k), beta);
291  mat.coeffRef(k,k) = beta;
292 
293  // apply H to remaining part of m_qr from the left
294  mat.bottomRightCorner(remainingRows, remainingCols)
295  .applyHouseholderOnTheLeft(mat.col(k).tail(remainingRows-1), hCoeffs.coeffRef(k), tempData+k+1);
296  }
297 }
298 
299 // TODO: add a corresponding public API for updating a QR factorization
308 template <typename MatrixQR, typename HCoeffs, typename VectorQR>
309 void householder_qr_inplace_update(MatrixQR& mat, HCoeffs& hCoeffs, const VectorQR& newColumn,
310  typename MatrixQR::Index k, typename MatrixQR::Scalar* tempData) {
311  typedef typename MatrixQR::Index Index;
312  typedef typename MatrixQR::RealScalar RealScalar;
313  Index rows = mat.rows();
314 
315  eigen_assert(k < mat.cols());
316  eigen_assert(k < rows);
317  eigen_assert(hCoeffs.size() == mat.cols());
318  eigen_assert(newColumn.size() == rows);
319  eigen_assert(tempData);
320 
321  // Store new column in mat at column k
322  mat.col(k) = newColumn;
323  // Apply H = H_1...H_{k-1} on newColumn (skip if k=0)
324  for (Index i = 0; i < k; ++i) {
325  Index remainingRows = rows - i;
326  mat.col(k)
327  .tail(remainingRows)
328  .applyHouseholderOnTheLeft(mat.col(i).tail(remainingRows - 1), hCoeffs.coeffRef(i), tempData + i + 1);
329  }
330  // Construct Householder projector in-place in column k
331  RealScalar beta;
332  mat.col(k).tail(rows - k).makeHouseholderInPlace(hCoeffs.coeffRef(k), beta);
333  mat.coeffRef(k, k) = beta;
334 }
335 
337 template<typename MatrixQR, typename HCoeffs,
338  typename MatrixQRScalar = typename MatrixQR::Scalar,
339  bool InnerStrideIsOne = (MatrixQR::InnerStrideAtCompileTime == 1 && HCoeffs::InnerStrideAtCompileTime == 1)>
340 struct householder_qr_inplace_blocked
341 {
342  // This is specialized for LAPACK-supported Scalar types in HouseholderQR_LAPACKE.h
343  static void run(MatrixQR& mat, HCoeffs& hCoeffs, Index maxBlockSize=32,
344  typename MatrixQR::Scalar* tempData = 0)
345  {
346  typedef typename MatrixQR::Scalar Scalar;
347  typedef Block<MatrixQR,Dynamic,Dynamic> BlockType;
348 
349  Index rows = mat.rows();
350  Index cols = mat.cols();
351  Index size = (std::min)(rows, cols);
352 
353  typedef Matrix<Scalar,Dynamic,1,ColMajor,MatrixQR::MaxColsAtCompileTime,1> TempType;
354  TempType tempVector;
355  if(tempData==0)
356  {
357  tempVector.resize(cols);
358  tempData = tempVector.data();
359  }
360 
361  Index blockSize = (std::min)(maxBlockSize,size);
362 
363  Index k = 0;
364  for (k = 0; k < size; k += blockSize)
365  {
366  Index bs = (std::min)(size-k,blockSize); // actual size of the block
367  Index tcols = cols - k - bs; // trailing columns
368  Index brows = rows-k; // rows of the block
369 
370  // partition the matrix:
371  // A00 | A01 | A02
372  // mat = A10 | A11 | A12
373  // A20 | A21 | A22
374  // and performs the qr dec of [A11^T A12^T]^T
375  // and update [A21^T A22^T]^T using level 3 operations.
376  // Finally, the algorithm continue on A22
377 
378  BlockType A11_21 = mat.block(k,k,brows,bs);
379  Block<HCoeffs,Dynamic,1> hCoeffsSegment = hCoeffs.segment(k,bs);
380 
381  householder_qr_inplace_unblocked(A11_21, hCoeffsSegment, tempData);
382 
383  if(tcols)
384  {
385  BlockType A21_22 = mat.block(k,k+bs,brows,tcols);
386  apply_block_householder_on_the_left(A21_22,A11_21,hCoeffsSegment, false); // false == backward
387  }
388  }
389  }
390 };
391 
392 } // end namespace internal
393 
394 #ifndef EIGEN_PARSED_BY_DOXYGEN
395 template<typename MatrixType_>
396 template<typename RhsType, typename DstType>
397 void HouseholderQR<MatrixType_>::_solve_impl(const RhsType &rhs, DstType &dst) const
398 {
399  const Index rank = (std::min)(rows(), cols());
400 
401  typename RhsType::PlainObject c(rhs);
402 
403  c.applyOnTheLeft(householderQ().setLength(rank).adjoint() );
404 
405  m_qr.topLeftCorner(rank, rank)
406  .template triangularView<Upper>()
407  .solveInPlace(c.topRows(rank));
408 
409  dst.topRows(rank) = c.topRows(rank);
410  dst.bottomRows(cols()-rank).setZero();
411 }
412 
413 template<typename MatrixType_>
414 template<bool Conjugate, typename RhsType, typename DstType>
415 void HouseholderQR<MatrixType_>::_solve_impl_transposed(const RhsType &rhs, DstType &dst) const
416 {
417  const Index rank = (std::min)(rows(), cols());
418 
419  typename RhsType::PlainObject c(rhs);
420 
421  m_qr.topLeftCorner(rank, rank)
422  .template triangularView<Upper>()
423  .transpose().template conjugateIf<Conjugate>()
424  .solveInPlace(c.topRows(rank));
425 
426  dst.topRows(rank) = c.topRows(rank);
427  dst.bottomRows(rows()-rank).setZero();
428 
429  dst.applyOnTheLeft(householderQ().setLength(rank).template conjugateIf<!Conjugate>() );
430 }
431 #endif
432 
439 template<typename MatrixType>
441 {
442  Index rows = m_qr.rows();
443  Index cols = m_qr.cols();
444  Index size = (std::min)(rows,cols);
445 
446  m_hCoeffs.resize(size);
447 
448  m_temp.resize(cols);
449 
450  internal::householder_qr_inplace_blocked<MatrixType, HCoeffsType>::run(m_qr, m_hCoeffs, 48, m_temp.data());
451 
452  m_isInitialized = true;
453 }
454 
459 template<typename Derived>
462 {
463  return HouseholderQR<PlainObject>(eval());
464 }
465 
466 } // end namespace Eigen
467 
468 #endif // EIGEN_QR_H
Expression of a fixed-size or dynamic-size block.
Definition: Block.h:107
Householder QR decomposition of a matrix.
Definition: HouseholderQR.h:60
HouseholderQR(Index rows, Index cols)
Default Constructor with memory preallocation.
Definition: HouseholderQR.h:91
HouseholderQR(EigenBase< InputType > &matrix)
Constructs a QR factorization from a given matrix.
Definition: HouseholderQR.h:128
const Solve< HouseholderQR, Rhs > solve(const MatrixBase< Rhs > &b) const
void computeInPlace()
Definition: HouseholderQR.h:440
HouseholderQR(const EigenBase< InputType > &matrix)
Constructs a QR factorization from a given matrix.
Definition: HouseholderQR.h:110
const HCoeffsType & hCoeffs() const
Definition: HouseholderQR.h:223
HouseholderSequenceType householderQ() const
Definition: HouseholderQR.h:165
HouseholderQR()
Default Constructor.
Definition: HouseholderQR.h:83
MatrixType::RealScalar absDeterminant() const
Definition: HouseholderQR.h:246
const MatrixType & matrixQR() const
Definition: HouseholderQR.h:174
MatrixType::RealScalar logAbsDeterminant() const
Definition: HouseholderQR.h:255
Sequence of Householder reflections acting on subspaces with decreasing size.
Definition: HouseholderSequence.h:123
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:52
const HouseholderQR< PlainObject > householderQr() const
Definition: HouseholderQR.h:461
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:182
void resize(Index rows, Index cols)
Definition: PlainObjectBase.h:283
Pseudo expression representing a solving operation.
Definition: Solve.h:65
A base class for matrix decomposition and solvers.
Definition: SolverBase.h:71
HouseholderQR< MatrixType_ > & derived()
Definition: EigenBase.h:48
@ ColMajor
Definition: Constants.h:321
@ RowMajor
Definition: Constants.h:323
const unsigned int RowMajorBit
Definition: Constants.h:68
Namespace containing all symbols from the Eigen library.
Definition: Core:139
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:59
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_abs_op< typename Derived::Scalar >, const Derived > abs(const Eigen::ArrayBase< Derived > &x)
Definition: EigenBase.h:32
Derived & derived()
Definition: EigenBase.h:48
Eigen::Index Index
The interface type of indices.
Definition: EigenBase.h:41