1// This file is part of Eigen, a lightweight C++ template library
2// for linear algebra.
3//
4// Copyright (C) 2016 Rasmus Munk Larsen <rmlarsen@google.com>
5//
6// This Source Code Form is subject to the terms of the Mozilla
7// Public License v. 2.0. If a copy of the MPL was not distributed
8// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9
10#ifndef EIGEN_COMPLETEORTHOGONALDECOMPOSITION_H
11#define EIGEN_COMPLETEORTHOGONALDECOMPOSITION_H
12
13namespace Eigen {
14
15namespace internal {
16template <typename _MatrixType>
17struct traits<CompleteOrthogonalDecomposition<_MatrixType> >
18 : traits<_MatrixType> {
19 enum { Flags = 0 };
20};
21
22} // end namespace internal
23
24/** \ingroup QR_Module
25 *
26 * \class CompleteOrthogonalDecomposition
27 *
28 * \brief Complete orthogonal decomposition (COD) of a matrix.
29 *
30 * \param MatrixType the type of the matrix of which we are computing the COD.
31 *
32 * This class performs a rank-revealing complete orthogonal decomposition of a
33 * matrix \b A into matrices \b P, \b Q, \b T, and \b Z such that
34 * \f[
35 * \mathbf{A} \, \mathbf{P} = \mathbf{Q} \,
36 * \begin{bmatrix} \mathbf{T} & \mathbf{0} \\
37 * \mathbf{0} & \mathbf{0} \end{bmatrix} \, \mathbf{Z}
38 * \f]
39 * by using Householder transformations. Here, \b P is a permutation matrix,
40 * \b Q and \b Z are unitary matrices and \b T an upper triangular matrix of
41 * size rank-by-rank. \b A may be rank deficient.
42 *
43 * This class supports the \link InplaceDecomposition inplace decomposition \endlink mechanism.
44 *
45 * \sa MatrixBase::completeOrthogonalDecomposition()
46 */
47template <typename _MatrixType>
48class CompleteOrthogonalDecomposition {
49 public:
50 typedef _MatrixType MatrixType;
51 enum {
52 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
53 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
54 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
55 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
56 };
57 typedef typename MatrixType::Scalar Scalar;
58 typedef typename MatrixType::RealScalar RealScalar;
59 typedef typename MatrixType::StorageIndex StorageIndex;
60 typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
61 typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime>
62 PermutationType;
63 typedef typename internal::plain_row_type<MatrixType, Index>::type
64 IntRowVectorType;
65 typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
66 typedef typename internal::plain_row_type<MatrixType, RealScalar>::type
67 RealRowVectorType;
68 typedef HouseholderSequence<
69 MatrixType, typename internal::remove_all<
70 typename HCoeffsType::ConjugateReturnType>::type>
71 HouseholderSequenceType;
72 typedef typename MatrixType::PlainObject PlainObject;
73
74 private:
75 typedef typename PermutationType::Index PermIndexType;
76
77 public:
78 /**
79 * \brief Default Constructor.
80 *
81 * The default constructor is useful in cases in which the user intends to
82 * perform decompositions via
83 * \c CompleteOrthogonalDecomposition::compute(const* MatrixType&).
84 */
85 CompleteOrthogonalDecomposition() : m_cpqr(), m_zCoeffs(), m_temp() {}
86
87 /** \brief Default Constructor with memory preallocation
88 *
89 * Like the default constructor but with preallocation of the internal data
90 * according to the specified problem \a size.
91 * \sa CompleteOrthogonalDecomposition()
92 */
93 CompleteOrthogonalDecomposition(Index rows, Index cols)
94 : m_cpqr(rows, cols), m_zCoeffs((std::min)(rows, cols)), m_temp(cols) {}
95
96 /** \brief Constructs a complete orthogonal decomposition from a given
97 * matrix.
98 *
99 * This constructor computes the complete orthogonal decomposition of the
100 * matrix \a matrix by calling the method compute(). The default
101 * threshold for rank determination will be used. It is a short cut for:
102 *
103 * \code
104 * CompleteOrthogonalDecomposition<MatrixType> cod(matrix.rows(),
105 * matrix.cols());
106 * cod.setThreshold(Default);
107 * cod.compute(matrix);
108 * \endcode
109 *
110 * \sa compute()
111 */
112 template <typename InputType>
113 explicit CompleteOrthogonalDecomposition(const EigenBase<InputType>& matrix)
114 : m_cpqr(matrix.rows(), matrix.cols()),
115 m_zCoeffs((std::min)(matrix.rows(), matrix.cols())),
116 m_temp(matrix.cols())
117 {
118 compute(matrix.derived());
119 }
120
121 /** \brief Constructs a complete orthogonal decomposition from a given matrix
122 *
123 * This overloaded constructor is provided for \link InplaceDecomposition inplace decomposition \endlink when \c MatrixType is a Eigen::Ref.
124 *
125 * \sa CompleteOrthogonalDecomposition(const EigenBase&)
126 */
127 template<typename InputType>
128 explicit CompleteOrthogonalDecomposition(EigenBase<InputType>& matrix)
129 : m_cpqr(matrix.derived()),
130 m_zCoeffs((std::min)(matrix.rows(), matrix.cols())),
131 m_temp(matrix.cols())
132 {
133 computeInPlace();
134 }
135
136
137 /** This method computes the minimum-norm solution X to a least squares
138 * problem \f[\mathrm{minimize} \|A X - B\|, \f] where \b A is the matrix of
139 * which \c *this is the complete orthogonal decomposition.
140 *
141 * \param b the right-hand sides of the problem to solve.
142 *
143 * \returns a solution.
144 *
145 */
146 template <typename Rhs>
147 inline const Solve<CompleteOrthogonalDecomposition, Rhs> solve(
148 const MatrixBase<Rhs>& b) const {
149 eigen_assert(m_cpqr.m_isInitialized &&
150 "CompleteOrthogonalDecomposition is not initialized.");
151 return Solve<CompleteOrthogonalDecomposition, Rhs>(*this, b.derived());
152 }
153
154 HouseholderSequenceType householderQ(void) const;
155 HouseholderSequenceType matrixQ(void) const { return m_cpqr.householderQ(); }
156
157 /** \returns the matrix \b Z.
158 */
159 MatrixType matrixZ() const {
160 MatrixType Z = MatrixType::Identity(m_cpqr.cols(), m_cpqr.cols());
161 applyZAdjointOnTheLeftInPlace(Z);
162 return Z.adjoint();
163 }
164
165 /** \returns a reference to the matrix where the complete orthogonal
166 * decomposition is stored
167 */
168 const MatrixType& matrixQTZ() const { return m_cpqr.matrixQR(); }
169
170 /** \returns a reference to the matrix where the complete orthogonal
171 * decomposition is stored.
172 * \warning The strict lower part and \code cols() - rank() \endcode right
173 * columns of this matrix contains internal values.
174 * Only the upper triangular part should be referenced. To get it, use
175 * \code matrixT().template triangularView<Upper>() \endcode
176 * For rank-deficient matrices, use
177 * \code
178 * matrixR().topLeftCorner(rank(), rank()).template triangularView<Upper>()
179 * \endcode
180 */
181 const MatrixType& matrixT() const { return m_cpqr.matrixQR(); }
182
183 template <typename InputType>
184 CompleteOrthogonalDecomposition& compute(const EigenBase<InputType>& matrix) {
185 // Compute the column pivoted QR factorization A P = Q R.
186 m_cpqr.compute(matrix);
187 computeInPlace();
188 return *this;
189 }
190
191 /** \returns a const reference to the column permutation matrix */
192 const PermutationType& colsPermutation() const {
193 return m_cpqr.colsPermutation();
194 }
195
196 /** \returns the absolute value of the determinant of the matrix of which
197 * *this is the complete orthogonal decomposition. It has only linear
198 * complexity (that is, O(n) where n is the dimension of the square matrix)
199 * as the complete orthogonal decomposition has already been computed.
200 *
201 * \note This is only for square matrices.
202 *
203 * \warning a determinant can be very big or small, so for matrices
204 * of large enough dimension, there is a risk of overflow/underflow.
205 * One way to work around that is to use logAbsDeterminant() instead.
206 *
207 * \sa logAbsDeterminant(), MatrixBase::determinant()
208 */
209 typename MatrixType::RealScalar absDeterminant() const;
210
211 /** \returns the natural log of the absolute value of the determinant of the
212 * matrix of which *this is the complete orthogonal decomposition. It has
213 * only linear complexity (that is, O(n) where n is the dimension of the
214 * square matrix) as the complete orthogonal decomposition has already been
215 * computed.
216 *
217 * \note This is only for square matrices.
218 *
219 * \note This method is useful to work around the risk of overflow/underflow
220 * that's inherent to determinant computation.
221 *
222 * \sa absDeterminant(), MatrixBase::determinant()
223 */
224 typename MatrixType::RealScalar logAbsDeterminant() const;
225
226 /** \returns the rank of the matrix of which *this is the complete orthogonal
227 * decomposition.
228 *
229 * \note This method has to determine which pivots should be considered
230 * nonzero. For that, it uses the threshold value that you can control by
231 * calling setThreshold(const RealScalar&).
232 */
233 inline Index rank() const { return m_cpqr.rank(); }
234
235 /** \returns the dimension of the kernel of the matrix of which *this is the
236 * complete orthogonal decomposition.
237 *
238 * \note This method has to determine which pivots should be considered
239 * nonzero. For that, it uses the threshold value that you can control by
240 * calling setThreshold(const RealScalar&).
241 */
242 inline Index dimensionOfKernel() const { return m_cpqr.dimensionOfKernel(); }
243
244 /** \returns true if the matrix of which *this is the decomposition represents
245 * an injective linear map, i.e. has trivial kernel; false otherwise.
246 *
247 * \note This method has to determine which pivots should be considered
248 * nonzero. For that, it uses the threshold value that you can control by
249 * calling setThreshold(const RealScalar&).
250 */
251 inline bool isInjective() const { return m_cpqr.isInjective(); }
252
253 /** \returns true if the matrix of which *this is the decomposition represents
254 * a surjective linear map; false otherwise.
255 *
256 * \note This method has to determine which pivots should be considered
257 * nonzero. For that, it uses the threshold value that you can control by
258 * calling setThreshold(const RealScalar&).
259 */
260 inline bool isSurjective() const { return m_cpqr.isSurjective(); }
261
262 /** \returns true if the matrix of which *this is the complete orthogonal
263 * decomposition is invertible.
264 *
265 * \note This method has to determine which pivots should be considered
266 * nonzero. For that, it uses the threshold value that you can control by
267 * calling setThreshold(const RealScalar&).
268 */
269 inline bool isInvertible() const { return m_cpqr.isInvertible(); }
270
271 /** \returns the pseudo-inverse of the matrix of which *this is the complete
272 * orthogonal decomposition.
273 * \warning: Do not compute \c this->pseudoInverse()*rhs to solve a linear systems.
274 * It is more efficient and numerically stable to call \c this->solve(rhs).
275 */
276 inline const Inverse<CompleteOrthogonalDecomposition> pseudoInverse() const
277 {
278 return Inverse<CompleteOrthogonalDecomposition>(*this);
279 }
280
281 inline Index rows() const { return m_cpqr.rows(); }
282 inline Index cols() const { return m_cpqr.cols(); }
283
284 /** \returns a const reference to the vector of Householder coefficients used
285 * to represent the factor \c Q.
286 *
287 * For advanced uses only.
288 */
289 inline const HCoeffsType& hCoeffs() const { return m_cpqr.hCoeffs(); }
290
291 /** \returns a const reference to the vector of Householder coefficients
292 * used to represent the factor \c Z.
293 *
294 * For advanced uses only.
295 */
296 const HCoeffsType& zCoeffs() const { return m_zCoeffs; }
297
298 /** Allows to prescribe a threshold to be used by certain methods, such as
299 * rank(), who need to determine when pivots are to be considered nonzero.
300 * Most be called before calling compute().
301 *
302 * When it needs to get the threshold value, Eigen calls threshold(). By
303 * default, this uses a formula to automatically determine a reasonable
304 * threshold. Once you have called the present method
305 * setThreshold(const RealScalar&), your value is used instead.
306 *
307 * \param threshold The new value to use as the threshold.
308 *
309 * A pivot will be considered nonzero if its absolute value is strictly
310 * greater than
311 * \f$ \vert pivot \vert \leqslant threshold \times \vert maxpivot \vert \f$
312 * where maxpivot is the biggest pivot.
313 *
314 * If you want to come back to the default behavior, call
315 * setThreshold(Default_t)
316 */
317 CompleteOrthogonalDecomposition& setThreshold(const RealScalar& threshold) {
318 m_cpqr.setThreshold(threshold);
319 return *this;
320 }
321
322 /** Allows to come back to the default behavior, letting Eigen use its default
323 * formula for determining the threshold.
324 *
325 * You should pass the special object Eigen::Default as parameter here.
326 * \code qr.setThreshold(Eigen::Default); \endcode
327 *
328 * See the documentation of setThreshold(const RealScalar&).
329 */
330 CompleteOrthogonalDecomposition& setThreshold(Default_t) {
331 m_cpqr.setThreshold(Default);
332 return *this;
333 }
334
335 /** Returns the threshold that will be used by certain methods such as rank().
336 *
337 * See the documentation of setThreshold(const RealScalar&).
338 */
339 RealScalar threshold() const { return m_cpqr.threshold(); }
340
341 /** \returns the number of nonzero pivots in the complete orthogonal
342 * decomposition. Here nonzero is meant in the exact sense, not in a
343 * fuzzy sense. So that notion isn't really intrinsically interesting,
344 * but it is still useful when implementing algorithms.
345 *
346 * \sa rank()
347 */
348 inline Index nonzeroPivots() const { return m_cpqr.nonzeroPivots(); }
349
350 /** \returns the absolute value of the biggest pivot, i.e. the biggest
351 * diagonal coefficient of R.
352 */
353 inline RealScalar maxPivot() const { return m_cpqr.maxPivot(); }
354
355 /** \brief Reports whether the complete orthogonal decomposition was
356 * succesful.
357 *
358 * \note This function always returns \c Success. It is provided for
359 * compatibility
360 * with other factorization routines.
361 * \returns \c Success
362 */
363 ComputationInfo info() const {
364 eigen_assert(m_cpqr.m_isInitialized && "Decomposition is not initialized.");
365 return Success;
366 }
367
368#ifndef EIGEN_PARSED_BY_DOXYGEN
369 template <typename RhsType, typename DstType>
370 EIGEN_DEVICE_FUNC void _solve_impl(const RhsType& rhs, DstType& dst) const;
371#endif
372
373 protected:
374 static void check_template_parameters() {
375 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
376 }
377
378 void computeInPlace();
379
380 /** Overwrites \b rhs with \f$ \mathbf{Z}^* * \mathbf{rhs} \f$.
381 */
382 template <typename Rhs>
383 void applyZAdjointOnTheLeftInPlace(Rhs& rhs) const;
384
385 ColPivHouseholderQR<MatrixType> m_cpqr;
386 HCoeffsType m_zCoeffs;
387 RowVectorType m_temp;
388};
389
390template <typename MatrixType>
391typename MatrixType::RealScalar
392CompleteOrthogonalDecomposition<MatrixType>::absDeterminant() const {
393 return m_cpqr.absDeterminant();
394}
395
396template <typename MatrixType>
397typename MatrixType::RealScalar
398CompleteOrthogonalDecomposition<MatrixType>::logAbsDeterminant() const {
399 return m_cpqr.logAbsDeterminant();
400}
401
402/** Performs the complete orthogonal decomposition of the given matrix \a
403 * matrix. The result of the factorization is stored into \c *this, and a
404 * reference to \c *this is returned.
405 *
406 * \sa class CompleteOrthogonalDecomposition,
407 * CompleteOrthogonalDecomposition(const MatrixType&)
408 */
409template <typename MatrixType>
410void CompleteOrthogonalDecomposition<MatrixType>::computeInPlace()
411{
412 check_template_parameters();
413
414 // the column permutation is stored as int indices, so just to be sure:
415 eigen_assert(m_cpqr.cols() <= NumTraits<int>::highest());
416
417 const Index rank = m_cpqr.rank();
418 const Index cols = m_cpqr.cols();
419 const Index rows = m_cpqr.rows();
420 m_zCoeffs.resize((std::min)(rows, cols));
421 m_temp.resize(cols);
422
423 if (rank < cols) {
424 // We have reduced the (permuted) matrix to the form
425 // [R11 R12]
426 // [ 0 R22]
427 // where R11 is r-by-r (r = rank) upper triangular, R12 is
428 // r-by-(n-r), and R22 is empty or the norm of R22 is negligible.
429 // We now compute the complete orthogonal decomposition by applying
430 // Householder transformations from the right to the upper trapezoidal
431 // matrix X = [R11 R12] to zero out R12 and obtain the factorization
432 // [R11 R12] = [T11 0] * Z, where T11 is r-by-r upper triangular and
433 // Z = Z(0) * Z(1) ... Z(r-1) is an n-by-n orthogonal matrix.
434 // We store the data representing Z in R12 and m_zCoeffs.
435 for (Index k = rank - 1; k >= 0; --k) {
436 if (k != rank - 1) {
437 // Given the API for Householder reflectors, it is more convenient if
438 // we swap the leading parts of columns k and r-1 (zero-based) to form
439 // the matrix X_k = [X(0:k, k), X(0:k, r:n)]
440 m_cpqr.m_qr.col(k).head(k + 1).swap(
441 m_cpqr.m_qr.col(rank - 1).head(k + 1));
442 }
443 // Construct Householder reflector Z(k) to zero out the last row of X_k,
444 // i.e. choose Z(k) such that
445 // [X(k, k), X(k, r:n)] * Z(k) = [beta, 0, .., 0].
446 RealScalar beta;
447 m_cpqr.m_qr.row(k)
448 .tail(cols - rank + 1)
449 .makeHouseholderInPlace(m_zCoeffs(k), beta);
450 m_cpqr.m_qr(k, rank - 1) = beta;
451 if (k > 0) {
452 // Apply Z(k) to the first k rows of X_k
453 m_cpqr.m_qr.topRightCorner(k, cols - rank + 1)
454 .applyHouseholderOnTheRight(
455 m_cpqr.m_qr.row(k).tail(cols - rank).transpose(), m_zCoeffs(k),
456 &m_temp(0));
457 }
458 if (k != rank - 1) {
459 // Swap X(0:k,k) back to its proper location.
460 m_cpqr.m_qr.col(k).head(k + 1).swap(
461 m_cpqr.m_qr.col(rank - 1).head(k + 1));
462 }
463 }
464 }
465}
466
467template <typename MatrixType>
468template <typename Rhs>
469void CompleteOrthogonalDecomposition<MatrixType>::applyZAdjointOnTheLeftInPlace(
470 Rhs& rhs) const {
471 const Index cols = this->cols();
472 const Index nrhs = rhs.cols();
473 const Index rank = this->rank();
474 Matrix<typename MatrixType::Scalar, Dynamic, 1> temp((std::max)(cols, nrhs));
475 for (Index k = 0; k < rank; ++k) {
476 if (k != rank - 1) {
477 rhs.row(k).swap(rhs.row(rank - 1));
478 }
479 rhs.middleRows(rank - 1, cols - rank + 1)
480 .applyHouseholderOnTheLeft(
481 matrixQTZ().row(k).tail(cols - rank).adjoint(), zCoeffs()(k),
482 &temp(0));
483 if (k != rank - 1) {
484 rhs.row(k).swap(rhs.row(rank - 1));
485 }
486 }
487}
488
489#ifndef EIGEN_PARSED_BY_DOXYGEN
490template <typename _MatrixType>
491template <typename RhsType, typename DstType>
492void CompleteOrthogonalDecomposition<_MatrixType>::_solve_impl(
493 const RhsType& rhs, DstType& dst) const {
494 eigen_assert(rhs.rows() == this->rows());
495
496 const Index rank = this->rank();
497 if (rank == 0) {
498 dst.setZero();
499 return;
500 }
501
502 // Compute c = Q^* * rhs
503 // Note that the matrix Q = H_0^* H_1^*... so its inverse is
504 // Q^* = (H_0 H_1 ...)^T
505 typename RhsType::PlainObject c(rhs);
506 c.applyOnTheLeft(
507 householderSequence(matrixQTZ(), hCoeffs()).setLength(rank).transpose());
508
509 // Solve T z = c(1:rank, :)
510 dst.topRows(rank) = matrixT()
511 .topLeftCorner(rank, rank)
512 .template triangularView<Upper>()
513 .solve(c.topRows(rank));
514
515 const Index cols = this->cols();
516 if (rank < cols) {
517 // Compute y = Z^* * [ z ]
518 // [ 0 ]
519 dst.bottomRows(cols - rank).setZero();
520 applyZAdjointOnTheLeftInPlace(dst);
521 }
522
523 // Undo permutation to get x = P^{-1} * y.
524 dst = colsPermutation() * dst;
525}
526#endif
527
528namespace internal {
529
530template<typename DstXprType, typename MatrixType>
531struct Assignment<DstXprType, Inverse<CompleteOrthogonalDecomposition<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename CompleteOrthogonalDecomposition<MatrixType>::Scalar>, Dense2Dense>
532{
533 typedef CompleteOrthogonalDecomposition<MatrixType> CodType;
534 typedef Inverse<CodType> SrcXprType;
535 static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename CodType::Scalar> &)
536 {
537 dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.rows()));
538 }
539};
540
541} // end namespace internal
542
543/** \returns the matrix Q as a sequence of householder transformations */
544template <typename MatrixType>
545typename CompleteOrthogonalDecomposition<MatrixType>::HouseholderSequenceType
546CompleteOrthogonalDecomposition<MatrixType>::householderQ() const {
547 return m_cpqr.householderQ();
548}
549
550/** \return the complete orthogonal decomposition of \c *this.
551 *
552 * \sa class CompleteOrthogonalDecomposition
553 */
554template <typename Derived>
555const CompleteOrthogonalDecomposition<typename MatrixBase<Derived>::PlainObject>
556MatrixBase<Derived>::completeOrthogonalDecomposition() const {
557 return CompleteOrthogonalDecomposition<PlainObject>(eval());
558}
559
560} // end namespace Eigen
561
562#endif // EIGEN_COMPLETEORTHOGONALDECOMPOSITION_H
563