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