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 typename HouseholderSequence<MatrixType,HCoeffsType>::ConjugateReturnType 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 
82  HouseholderQR(const MatrixType& matrix)
83  : m_qr(matrix.rows(), matrix.cols()),
84  m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
85  m_temp(matrix.cols()),
86  m_isInitialized(false)
87  {
88  compute(matrix);
89  }
90 
108  template<typename Rhs>
109  inline const internal::solve_retval<HouseholderQR, Rhs>
110  solve(const MatrixBase<Rhs>& b) const
111  {
112  eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
113  return internal::solve_retval<HouseholderQR, Rhs>(*this, b.derived());
114  }
115 
116  HouseholderSequenceType householderQ() const
117  {
118  eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
119  return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate());
120  }
121 
125  const MatrixType& matrixQR() const
126  {
127  eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
128  return m_qr;
129  }
130 
131  HouseholderQR& compute(const MatrixType& matrix);
132 
146  typename MatrixType::RealScalar absDeterminant() const;
147 
160  typename MatrixType::RealScalar logAbsDeterminant() const;
161 
162  inline Index rows() const { return m_qr.rows(); }
163  inline Index cols() const { return m_qr.cols(); }
164  const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
165 
166  protected:
167  MatrixType m_qr;
168  HCoeffsType m_hCoeffs;
169  RowVectorType m_temp;
170  bool m_isInitialized;
171 };
172 
173 template<typename MatrixType>
174 typename MatrixType::RealScalar HouseholderQR<MatrixType>::absDeterminant() const
175 {
176  eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
177  eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
178  return internal::abs(m_qr.diagonal().prod());
179 }
180 
181 template<typename MatrixType>
182 typename MatrixType::RealScalar HouseholderQR<MatrixType>::logAbsDeterminant() const
183 {
184  eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
185  eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
186  return m_qr.diagonal().cwiseAbs().array().log().sum();
187 }
188 
189 namespace internal {
190 
192 template<typename MatrixQR, typename HCoeffs>
193 void householder_qr_inplace_unblocked(MatrixQR& mat, HCoeffs& hCoeffs, typename MatrixQR::Scalar* tempData = 0)
194 {
195  typedef typename MatrixQR::Index Index;
196  typedef typename MatrixQR::Scalar Scalar;
197  typedef typename MatrixQR::RealScalar RealScalar;
198  Index rows = mat.rows();
199  Index cols = mat.cols();
200  Index size = (std::min)(rows,cols);
201 
202  eigen_assert(hCoeffs.size() == size);
203 
205  TempType tempVector;
206  if(tempData==0)
207  {
208  tempVector.resize(cols);
209  tempData = tempVector.data();
210  }
211 
212  for(Index k = 0; k < size; ++k)
213  {
214  Index remainingRows = rows - k;
215  Index remainingCols = cols - k - 1;
216 
217  RealScalar beta;
218  mat.col(k).tail(remainingRows).makeHouseholderInPlace(hCoeffs.coeffRef(k), beta);
219  mat.coeffRef(k,k) = beta;
220 
221  // apply H to remaining part of m_qr from the left
222  mat.bottomRightCorner(remainingRows, remainingCols)
223  .applyHouseholderOnTheLeft(mat.col(k).tail(remainingRows-1), hCoeffs.coeffRef(k), tempData+k+1);
224  }
225 }
226 
228 template<typename MatrixQR, typename HCoeffs>
229 void householder_qr_inplace_blocked(MatrixQR& mat, HCoeffs& hCoeffs,
230  typename MatrixQR::Index maxBlockSize=32,
231  typename MatrixQR::Scalar* tempData = 0)
232 {
233  typedef typename MatrixQR::Index Index;
234  typedef typename MatrixQR::Scalar Scalar;
235  // typedef typename MatrixQR::RealScalar RealScalar;
236  typedef Block<MatrixQR,Dynamic,Dynamic> BlockType;
237 
238  Index rows = mat.rows();
239  Index cols = mat.cols();
240  Index size = (std::min)(rows, cols);
241 
242  typedef Matrix<Scalar,Dynamic,1,ColMajor,MatrixQR::MaxColsAtCompileTime,1> TempType;
243  TempType tempVector;
244  if(tempData==0)
245  {
246  tempVector.resize(cols);
247  tempData = tempVector.data();
248  }
249 
250  Index blockSize = (std::min)(maxBlockSize,size);
251 
252  Index k = 0;
253  for (k = 0; k < size; k += blockSize)
254  {
255  Index bs = (std::min)(size-k,blockSize); // actual size of the block
256  Index tcols = cols - k - bs; // trailing columns
257  Index brows = rows-k; // rows of the block
258 
259  // partition the matrix:
260  // A00 | A01 | A02
261  // mat = A10 | A11 | A12
262  // A20 | A21 | A22
263  // and performs the qr dec of [A11^T A12^T]^T
264  // and update [A21^T A22^T]^T using level 3 operations.
265  // Finally, the algorithm continue on A22
266 
267  BlockType A11_21 = mat.block(k,k,brows,bs);
268  Block<HCoeffs,Dynamic,1> hCoeffsSegment = hCoeffs.segment(k,bs);
269 
270  householder_qr_inplace_unblocked(A11_21, hCoeffsSegment, tempData);
271 
272  if(tcols)
273  {
274  BlockType A21_22 = mat.block(k,k+bs,brows,tcols);
275  apply_block_householder_on_the_left(A21_22,A11_21,hCoeffsSegment.adjoint());
276  }
277  }
278 }
279 
280 template<typename _MatrixType, typename Rhs>
281 struct solve_retval<HouseholderQR<_MatrixType>, Rhs>
282  : solve_retval_base<HouseholderQR<_MatrixType>, Rhs>
283 {
284  EIGEN_MAKE_SOLVE_HELPERS(HouseholderQR<_MatrixType>,Rhs)
285 
286  template<typename Dest> void evalTo(Dest& dst) const
287  {
288  const Index rows = dec().rows(), cols = dec().cols();
289  const Index rank = (std::min)(rows, cols);
290  eigen_assert(rhs().rows() == rows);
291 
292  typename Rhs::PlainObject c(rhs());
293 
294  // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
295  c.applyOnTheLeft(householderSequence(
296  dec().matrixQR().leftCols(rank),
297  dec().hCoeffs().head(rank)).transpose()
298  );
299 
300  dec().matrixQR()
301  .topLeftCorner(rank, rank)
302  .template triangularView<Upper>()
303  .solveInPlace(c.topRows(rank));
304 
305  dst.topRows(rank) = c.topRows(rank);
306  dst.bottomRows(cols-rank).setZero();
307  }
308 };
309 
310 } // end namespace internal
311 
312 template<typename MatrixType>
313 HouseholderQR<MatrixType>& HouseholderQR<MatrixType>::compute(const MatrixType& matrix)
314 {
315  Index rows = matrix.rows();
316  Index cols = matrix.cols();
317  Index size = (std::min)(rows,cols);
318 
319  m_qr = matrix;
320  m_hCoeffs.resize(size);
321 
322  m_temp.resize(cols);
323 
324  internal::householder_qr_inplace_blocked(m_qr, m_hCoeffs, 48, m_temp.data());
325 
326  m_isInitialized = true;
327  return *this;
328 }
329 
334 template<typename Derived>
335 const HouseholderQR<typename MatrixBase<Derived>::PlainObject>
337 {
338  return HouseholderQR<PlainObject>(eval());
339 }
340 
341 } // end namespace Eigen
342 
343 #endif // EIGEN_QR_H