Eigen  3.2.10
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 namespace Eigen {
16 
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 
68  HouseholderQR() : m_qr(), m_hCoeffs(), m_temp(), m_isInitialized(false) {}
69 
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 
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 
117  template<typename Rhs>
118  inline const internal::solve_retval<HouseholderQR, Rhs>
119  solve(const MatrixBase<Rhs>& b) const
120  {
121  eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
122  return internal::solve_retval<HouseholderQR, Rhs>(*this, b.derived());
123  }
124 
133  HouseholderSequenceType householderQ() const
134  {
135  eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
136  return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate());
137  }
138 
142  const MatrixType& matrixQR() const
143  {
144  eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
145  return m_qr;
146  }
147 
148  HouseholderQR& compute(const MatrixType& matrix);
149 
163  typename MatrixType::RealScalar absDeterminant() const;
164 
177  typename MatrixType::RealScalar logAbsDeterminant() const;
178 
179  inline Index rows() const { return m_qr.rows(); }
180  inline Index cols() const { return m_qr.cols(); }
181 
186  const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
187 
188  protected:
189 
190  static void check_template_parameters()
191  {
192  EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
193  }
194 
195  MatrixType m_qr;
196  HCoeffsType m_hCoeffs;
197  RowVectorType m_temp;
198  bool m_isInitialized;
199 };
200 
201 template<typename MatrixType>
202 typename MatrixType::RealScalar HouseholderQR<MatrixType>::absDeterminant() const
203 {
204  using std::abs;
205  eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
206  eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
207  return abs(m_qr.diagonal().prod());
208 }
209 
210 template<typename MatrixType>
211 typename MatrixType::RealScalar HouseholderQR<MatrixType>::logAbsDeterminant() const
212 {
213  eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
214  eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
215  return m_qr.diagonal().cwiseAbs().array().log().sum();
216 }
217 
218 namespace internal {
219 
221 template<typename MatrixQR, typename HCoeffs>
222 void householder_qr_inplace_unblocked(MatrixQR& mat, HCoeffs& hCoeffs, typename MatrixQR::Scalar* tempData = 0)
223 {
224  typedef typename MatrixQR::Index Index;
225  typedef typename MatrixQR::Scalar Scalar;
226  typedef typename MatrixQR::RealScalar RealScalar;
227  Index rows = mat.rows();
228  Index cols = mat.cols();
229  Index size = (std::min)(rows,cols);
230 
231  eigen_assert(hCoeffs.size() == size);
232 
234  TempType tempVector;
235  if(tempData==0)
236  {
237  tempVector.resize(cols);
238  tempData = tempVector.data();
239  }
240 
241  for(Index k = 0; k < size; ++k)
242  {
243  Index remainingRows = rows - k;
244  Index remainingCols = cols - k - 1;
245 
246  RealScalar beta;
247  mat.col(k).tail(remainingRows).makeHouseholderInPlace(hCoeffs.coeffRef(k), beta);
248  mat.coeffRef(k,k) = beta;
249 
250  // apply H to remaining part of m_qr from the left
251  mat.bottomRightCorner(remainingRows, remainingCols)
252  .applyHouseholderOnTheLeft(mat.col(k).tail(remainingRows-1), hCoeffs.coeffRef(k), tempData+k+1);
253  }
254 }
255 
257 template<typename MatrixQR, typename HCoeffs,
258  typename MatrixQRScalar = typename MatrixQR::Scalar,
259  bool InnerStrideIsOne = (MatrixQR::InnerStrideAtCompileTime == 1 && HCoeffs::InnerStrideAtCompileTime == 1)>
260 struct householder_qr_inplace_blocked
261 {
262  // This is specialized for MKL-supported Scalar types in HouseholderQR_MKL.h
263  static void run(MatrixQR& mat, HCoeffs& hCoeffs,
264  typename MatrixQR::Index maxBlockSize=32,
265  typename MatrixQR::Scalar* tempData = 0)
266  {
267  typedef typename MatrixQR::Index Index;
268  typedef typename MatrixQR::Scalar Scalar;
269  typedef Block<MatrixQR,Dynamic,Dynamic> BlockType;
270 
271  Index rows = mat.rows();
272  Index cols = mat.cols();
273  Index size = (std::min)(rows, cols);
274 
276  TempType tempVector;
277  if(tempData==0)
278  {
279  tempVector.resize(cols);
280  tempData = tempVector.data();
281  }
282 
283  Index blockSize = (std::min)(maxBlockSize,size);
284 
285  Index k = 0;
286  for (k = 0; k < size; k += blockSize)
287  {
288  Index bs = (std::min)(size-k,blockSize); // actual size of the block
289  Index tcols = cols - k - bs; // trailing columns
290  Index brows = rows-k; // rows of the block
291 
292  // partition the matrix:
293  // A00 | A01 | A02
294  // mat = A10 | A11 | A12
295  // A20 | A21 | A22
296  // and performs the qr dec of [A11^T A12^T]^T
297  // and update [A21^T A22^T]^T using level 3 operations.
298  // Finally, the algorithm continue on A22
299 
300  BlockType A11_21 = mat.block(k,k,brows,bs);
301  Block<HCoeffs,Dynamic,1> hCoeffsSegment = hCoeffs.segment(k,bs);
302 
303  householder_qr_inplace_unblocked(A11_21, hCoeffsSegment, tempData);
304 
305  if(tcols)
306  {
307  BlockType A21_22 = mat.block(k,k+bs,brows,tcols);
308  apply_block_householder_on_the_left(A21_22,A11_21,hCoeffsSegment.adjoint());
309  }
310  }
311  }
312 };
313 
314 template<typename _MatrixType, typename Rhs>
315 struct solve_retval<HouseholderQR<_MatrixType>, Rhs>
316  : solve_retval_base<HouseholderQR<_MatrixType>, Rhs>
317 {
318  EIGEN_MAKE_SOLVE_HELPERS(HouseholderQR<_MatrixType>,Rhs)
319 
320  template<typename Dest> void evalTo(Dest& dst) const
321  {
322  const Index rows = dec().rows(), cols = dec().cols();
323  const Index rank = (std::min)(rows, cols);
324  eigen_assert(rhs().rows() == rows);
325 
326  typename Rhs::PlainObject c(rhs());
327 
328  // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
329  c.applyOnTheLeft(householderSequence(
330  dec().matrixQR().leftCols(rank),
331  dec().hCoeffs().head(rank)).transpose()
332  );
333 
334  dec().matrixQR()
335  .topLeftCorner(rank, rank)
336  .template triangularView<Upper>()
337  .solveInPlace(c.topRows(rank));
338 
339  dst.topRows(rank) = c.topRows(rank);
340  dst.bottomRows(cols-rank).setZero();
341  }
342 };
343 
344 } // end namespace internal
345 
352 template<typename MatrixType>
354 {
355  check_template_parameters();
356 
357  Index rows = matrix.rows();
358  Index cols = matrix.cols();
359  Index size = (std::min)(rows,cols);
360 
361  m_qr = matrix;
362  m_hCoeffs.resize(size);
363 
364  m_temp.resize(cols);
365 
366  internal::householder_qr_inplace_blocked<MatrixType, HCoeffsType>::run(m_qr, m_hCoeffs, 48, m_temp.data());
367 
368  m_isInitialized = true;
369  return *this;
370 }
371 
376 template<typename Derived>
379 {
380  return HouseholderQR<PlainObject>(eval());
381 }
382 
383 } // end namespace Eigen
384 
385 #endif // EIGEN_QR_H
HouseholderQR(const MatrixType &matrix)
Constructs a QR factorization from a given matrix.
Definition: HouseholderQR.h:94
MatrixType::RealScalar logAbsDeterminant() const
Definition: HouseholderQR.h:211
Definition: LDLT.h:16
Definition: StdDeque.h:50
HouseholderSequence< VectorsType, CoeffsType > householderSequence(const VectorsType &v, const CoeffsType &h)
Convenience function for constructing a Householder sequence.
Definition: HouseholderSequence.h:423
HouseholderSequenceType householderQ() const
Definition: HouseholderQR.h:133
const HouseholderQR< PlainObject > householderQr() const
Definition: HouseholderQR.h:378
const internal::solve_retval< HouseholderQR, Rhs > solve(const MatrixBase< Rhs > &b) const
Definition: HouseholderQR.h:119
HouseholderQR()
Default Constructor.
Definition: HouseholderQR.h:68
MatrixType::RealScalar absDeterminant() const
Definition: HouseholderQR.h:202
const HCoeffsType & hCoeffs() const
Definition: HouseholderQR.h:186
HouseholderQR & compute(const MatrixType &matrix)
Definition: HouseholderQR.h:353
Definition: Eigen_Colamd.h:50
Expression of a fixed-size or dynamic-size block.
Definition: Block.h:102
Householder QR decomposition of a matrix.
Definition: ForwardDeclarations.h:221
HouseholderQR(Index rows, Index cols)
Default Constructor with memory preallocation.
Definition: HouseholderQR.h:76
void resize(Index nbRows, Index nbCols)
Definition: PlainObjectBase.h:235
const MatrixType & matrixQR() const
Definition: HouseholderQR.h:142
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:127
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48