HouseholderQR.h
Go to the documentation of this file.
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 // Eigen is free software; you can redistribute it and/or
9 // modify it under the terms of the GNU Lesser General Public
10 // License as published by the Free Software Foundation; either
11 // version 3 of the License, or (at your option) any later version.
12 //
13 // Alternatively, you can redistribute it and/or
14 // modify it under the terms of the GNU General Public License as
15 // published by the Free Software Foundation; either version 2 of
16 // the License, or (at your option) any later version.
17 //
18 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
19 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
20 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
21 // GNU General Public License for more details.
22 //
23 // You should have received a copy of the GNU Lesser General Public
24 // License and a copy of the GNU General Public License along with
25 // Eigen. If not, see <http://www.gnu.org/licenses/>.
26 
27 #ifndef EIGEN_QR_H
28 #define EIGEN_QR_H
29 
30 namespace Eigen {
31 
57 template<typename _MatrixType> class HouseholderQR
58 {
59  public:
60 
61  typedef _MatrixType MatrixType;
62  enum {
67  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
68  };
69  typedef typename MatrixType::Scalar Scalar;
70  typedef typename MatrixType::RealScalar RealScalar;
71  typedef typename MatrixType::Index Index;
73  typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
74  typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
76 
84 
92  : m_qr(rows, cols),
93  m_hCoeffs((std::min)(rows,cols)),
94  m_temp(cols),
95  m_isInitialized(false) {}
96 
97  HouseholderQR(const MatrixType& 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);
104  }
105 
123  template<typename Rhs>
124  inline const internal::solve_retval<HouseholderQR, Rhs>
125  solve(const MatrixBase<Rhs>& b) const
126  {
127  eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
128  return internal::solve_retval<HouseholderQR, Rhs>(*this, b.derived());
129  }
130 
132  {
133  eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
134  return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate());
135  }
136 
140  const MatrixType& matrixQR() const
141  {
142  eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
143  return m_qr;
144  }
145 
146  HouseholderQR& compute(const MatrixType& matrix);
147 
161  typename MatrixType::RealScalar absDeterminant() const;
162 
175  typename MatrixType::RealScalar logAbsDeterminant() const;
176 
177  inline Index rows() const { return m_qr.rows(); }
178  inline Index cols() const { return m_qr.cols(); }
179  const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
180 
181  protected:
186 };
187 
188 template<typename MatrixType>
189 typename MatrixType::RealScalar HouseholderQR<MatrixType>::absDeterminant() const
190 {
191  eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
192  eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
193  return internal::abs(m_qr.diagonal().prod());
194 }
195 
196 template<typename MatrixType>
197 typename MatrixType::RealScalar HouseholderQR<MatrixType>::logAbsDeterminant() const
198 {
199  eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
200  eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
201  return m_qr.diagonal().cwiseAbs().array().log().sum();
202 }
203 
204 namespace internal {
205 
207 template<typename MatrixQR, typename HCoeffs>
208 void householder_qr_inplace_unblocked(MatrixQR& mat, HCoeffs& hCoeffs, typename MatrixQR::Scalar* tempData = 0)
209 {
210  typedef typename MatrixQR::Index Index;
211  typedef typename MatrixQR::Scalar Scalar;
212  typedef typename MatrixQR::RealScalar RealScalar;
213  Index rows = mat.rows();
214  Index cols = mat.cols();
215  Index size = (std::min)(rows,cols);
216 
217  eigen_assert(hCoeffs.size() == size);
218 
220  TempType tempVector;
221  if(tempData==0)
222  {
223  tempVector.resize(cols);
224  tempData = tempVector.data();
225  }
226 
227  for(Index k = 0; k < size; ++k)
228  {
229  Index remainingRows = rows - k;
230  Index remainingCols = cols - k - 1;
231 
232  RealScalar beta;
233  mat.col(k).tail(remainingRows).makeHouseholderInPlace(hCoeffs.coeffRef(k), beta);
234  mat.coeffRef(k,k) = beta;
235 
236  // apply H to remaining part of m_qr from the left
237  mat.bottomRightCorner(remainingRows, remainingCols)
238  .applyHouseholderOnTheLeft(mat.col(k).tail(remainingRows-1), hCoeffs.coeffRef(k), tempData+k+1);
239  }
240 }
241 
243 template<typename MatrixQR, typename HCoeffs>
244 void householder_qr_inplace_blocked(MatrixQR& mat, HCoeffs& hCoeffs,
245  typename MatrixQR::Index maxBlockSize=32,
246  typename MatrixQR::Scalar* tempData = 0)
247 {
248  typedef typename MatrixQR::Index Index;
249  typedef typename MatrixQR::Scalar Scalar;
250  typedef typename MatrixQR::RealScalar RealScalar;
251  typedef Block<MatrixQR,Dynamic,Dynamic> BlockType;
252 
253  Index rows = mat.rows();
254  Index cols = mat.cols();
255  Index size = (std::min)(rows, cols);
256 
258  TempType tempVector;
259  if(tempData==0)
260  {
261  tempVector.resize(cols);
262  tempData = tempVector.data();
263  }
264 
265  Index blockSize = (std::min)(maxBlockSize,size);
266 
267  Index k = 0;
268  for (k = 0; k < size; k += blockSize)
269  {
270  Index bs = (std::min)(size-k,blockSize); // actual size of the block
271  Index tcols = cols - k - bs; // trailing columns
272  Index brows = rows-k; // rows of the block
273 
274  // partition the matrix:
275  // A00 | A01 | A02
276  // mat = A10 | A11 | A12
277  // A20 | A21 | A22
278  // and performs the qr dec of [A11^T A12^T]^T
279  // and update [A21^T A22^T]^T using level 3 operations.
280  // Finally, the algorithm continue on A22
281 
282  BlockType A11_21 = mat.block(k,k,brows,bs);
283  Block<HCoeffs,Dynamic,1> hCoeffsSegment = hCoeffs.segment(k,bs);
284 
285  householder_qr_inplace_unblocked(A11_21, hCoeffsSegment, tempData);
286 
287  if(tcols)
288  {
289  BlockType A21_22 = mat.block(k,k+bs,brows,tcols);
290  apply_block_householder_on_the_left(A21_22,A11_21,hCoeffsSegment.adjoint());
291  }
292  }
293 }
294 
295 template<typename _MatrixType, typename Rhs>
296 struct solve_retval<HouseholderQR<_MatrixType>, Rhs>
297  : solve_retval_base<HouseholderQR<_MatrixType>, Rhs>
298 {
300 
301  template<typename Dest> void evalTo(Dest& dst) const
302  {
303  const Index rows = dec().rows(), cols = dec().cols();
304  const Index rank = (std::min)(rows, cols);
305  eigen_assert(rhs().rows() == rows);
306 
307  typename Rhs::PlainObject c(rhs());
308 
309  // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
310  c.applyOnTheLeft(householderSequence(
311  dec().matrixQR().leftCols(rank),
312  dec().hCoeffs().head(rank)).transpose()
313  );
314 
315  dec().matrixQR()
316  .topLeftCorner(rank, rank)
317  .template triangularView<Upper>()
318  .solveInPlace(c.topRows(rank));
319 
320  dst.topRows(rank) = c.topRows(rank);
321  dst.bottomRows(cols-rank).setZero();
322  }
323 };
324 
325 } // end namespace internal
326 
327 template<typename MatrixType>
329 {
330  Index rows = matrix.rows();
331  Index cols = matrix.cols();
332  Index size = (std::min)(rows,cols);
333 
334  m_qr = matrix;
335  m_hCoeffs.resize(size);
336 
337  m_temp.resize(cols);
338 
339  internal::householder_qr_inplace_blocked(m_qr, m_hCoeffs, 48, m_temp.data());
340 
341  m_isInitialized = true;
342  return *this;
343 }
344 
349 template<typename Derived>
352 {
353  return HouseholderQR<PlainObject>(eval());
354 }
355 
356 } // end namespace Eigen
357 
358 #endif // EIGEN_QR_H