// This file is part of Eigen, a lightweight C++ template library // for linear algebra. // // Copyright (C) 2008-2009 Gael Guennebaud // Copyright (C) 2009 Benoit Jacob // // This Source Code Form is subject to the terms of the Mozilla // Public License v. 2.0. If a copy of the MPL was not distributed // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. #ifndef EIGEN_COLPIVOTINGHOUSEHOLDERQR_H #define EIGEN_COLPIVOTINGHOUSEHOLDERQR_H namespace Eigen { /** \ingroup QR_Module * * \class ColPivHouseholderQR * * \brief Householder rank-revealing QR decomposition of a matrix with column-pivoting * * \param MatrixType the type of the matrix of which we are computing the QR decomposition * * This class performs a rank-revealing QR decomposition of a matrix \b A into matrices \b P, \b Q and \b R * such that * \f[ * \mathbf{A} \, \mathbf{P} = \mathbf{Q} \, \mathbf{R} * \f] * by using Householder transformations. Here, \b P is a permutation matrix, \b Q a unitary matrix and \b R an * upper triangular matrix. * * This decomposition performs column pivoting in order to be rank-revealing and improve * numerical stability. It is slower than HouseholderQR, and faster than FullPivHouseholderQR. * * \sa MatrixBase::colPivHouseholderQr() */ template class ColPivHouseholderQR { public: typedef _MatrixType MatrixType; enum { RowsAtCompileTime = MatrixType::RowsAtCompileTime, ColsAtCompileTime = MatrixType::ColsAtCompileTime, Options = MatrixType::Options, MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime }; typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::RealScalar RealScalar; typedef typename MatrixType::Index Index; typedef Matrix MatrixQType; typedef typename internal::plain_diag_type::type HCoeffsType; typedef PermutationMatrix PermutationType; typedef typename internal::plain_row_type::type IntRowVectorType; typedef typename internal::plain_row_type::type RowVectorType; typedef typename internal::plain_row_type::type RealRowVectorType; typedef HouseholderSequence::type> HouseholderSequenceType; private: typedef typename PermutationType::Index PermIndexType; public: /** * \brief Default Constructor. * * The default constructor is useful in cases in which the user intends to * perform decompositions via ColPivHouseholderQR::compute(const MatrixType&). */ ColPivHouseholderQR() : m_qr(), m_hCoeffs(), m_colsPermutation(), m_colsTranspositions(), m_temp(), m_colSqNorms(), m_isInitialized(false), m_usePrescribedThreshold(false) {} /** \brief Default Constructor with memory preallocation * * Like the default constructor but with preallocation of the internal data * according to the specified problem \a size. * \sa ColPivHouseholderQR() */ ColPivHouseholderQR(Index rows, Index cols) : m_qr(rows, cols), m_hCoeffs((std::min)(rows,cols)), m_colsPermutation(PermIndexType(cols)), m_colsTranspositions(cols), m_temp(cols), m_colSqNorms(cols), m_isInitialized(false), m_usePrescribedThreshold(false) {} /** \brief Constructs a QR factorization from a given matrix * * This constructor computes the QR factorization of the matrix \a matrix by calling * the method compute(). It is a short cut for: * * \code * ColPivHouseholderQR qr(matrix.rows(), matrix.cols()); * qr.compute(matrix); * \endcode * * \sa compute() */ ColPivHouseholderQR(const MatrixType& matrix) : m_qr(matrix.rows(), matrix.cols()), m_hCoeffs((std::min)(matrix.rows(),matrix.cols())), m_colsPermutation(PermIndexType(matrix.cols())), m_colsTranspositions(matrix.cols()), m_temp(matrix.cols()), m_colSqNorms(matrix.cols()), m_isInitialized(false), m_usePrescribedThreshold(false) { compute(matrix); } /** This method finds a solution x to the equation Ax=b, where A is the matrix of which * *this is the QR decomposition, if any exists. * * \param b the right-hand-side of the equation to solve. * * \returns a solution. * * \note The case where b is a matrix is not yet implemented. Also, this * code is space inefficient. * * \note_about_checking_solutions * * \note_about_arbitrary_choice_of_solution * * Example: \include ColPivHouseholderQR_solve.cpp * Output: \verbinclude ColPivHouseholderQR_solve.out */ template inline const internal::solve_retval solve(const MatrixBase& b) const { eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); return internal::solve_retval(*this, b.derived()); } HouseholderSequenceType householderQ(void) const; HouseholderSequenceType matrixQ(void) const { return householderQ(); } /** \returns a reference to the matrix where the Householder QR decomposition is stored */ const MatrixType& matrixQR() const { eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); return m_qr; } /** \returns a reference to the matrix where the result Householder QR is stored * \warning The strict lower part of this matrix contains internal values. * Only the upper triangular part should be referenced. To get it, use * \code matrixR().template triangularView() \endcode * For rank-deficient matrices, use * \code * matrixR().topLeftCorner(rank(), rank()).template triangularView() * \endcode */ const MatrixType& matrixR() const { eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); return m_qr; } ColPivHouseholderQR& compute(const MatrixType& matrix); /** \returns a const reference to the column permutation matrix */ const PermutationType& colsPermutation() const { eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); return m_colsPermutation; } /** \returns the absolute value of the determinant of the matrix of which * *this is the QR decomposition. It has only linear complexity * (that is, O(n) where n is the dimension of the square matrix) * as the QR decomposition has already been computed. * * \note This is only for square matrices. * * \warning a determinant can be very big or small, so for matrices * of large enough dimension, there is a risk of overflow/underflow. * One way to work around that is to use logAbsDeterminant() instead. * * \sa logAbsDeterminant(), MatrixBase::determinant() */ typename MatrixType::RealScalar absDeterminant() const; /** \returns the natural log of the absolute value of the determinant of the matrix of which * *this is the QR decomposition. It has only linear complexity * (that is, O(n) where n is the dimension of the square matrix) * as the QR decomposition has already been computed. * * \note This is only for square matrices. * * \note This method is useful to work around the risk of overflow/underflow that's inherent * to determinant computation. * * \sa absDeterminant(), MatrixBase::determinant() */ typename MatrixType::RealScalar logAbsDeterminant() const; /** \returns the rank of the matrix of which *this is the QR decomposition. * * \note This method has to determine which pivots should be considered nonzero. * For that, it uses the threshold value that you can control by calling * setThreshold(const RealScalar&). */ inline Index rank() const { using std::abs; eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); RealScalar premultiplied_threshold = abs(m_maxpivot) * threshold(); Index result = 0; for(Index i = 0; i < m_nonzero_pivots; ++i) result += (abs(m_qr.coeff(i,i)) > premultiplied_threshold); return result; } /** \returns the dimension of the kernel of the matrix of which *this is the QR decomposition. * * \note This method has to determine which pivots should be considered nonzero. * For that, it uses the threshold value that you can control by calling * setThreshold(const RealScalar&). */ inline Index dimensionOfKernel() const { eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); return cols() - rank(); } /** \returns true if the matrix of which *this is the QR decomposition represents an injective * linear map, i.e. has trivial kernel; false otherwise. * * \note This method has to determine which pivots should be considered nonzero. * For that, it uses the threshold value that you can control by calling * setThreshold(const RealScalar&). */ inline bool isInjective() const { eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); return rank() == cols(); } /** \returns true if the matrix of which *this is the QR decomposition represents a surjective * linear map; false otherwise. * * \note This method has to determine which pivots should be considered nonzero. * For that, it uses the threshold value that you can control by calling * setThreshold(const RealScalar&). */ inline bool isSurjective() const { eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); return rank() == rows(); } /** \returns true if the matrix of which *this is the QR decomposition is invertible. * * \note This method has to determine which pivots should be considered nonzero. * For that, it uses the threshold value that you can control by calling * setThreshold(const RealScalar&). */ inline bool isInvertible() const { eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); return isInjective() && isSurjective(); } /** \returns the inverse of the matrix of which *this is the QR decomposition. * * \note If this matrix is not invertible, the returned matrix has undefined coefficients. * Use isInvertible() to first determine whether this matrix is invertible. */ inline const internal::solve_retval inverse() const { eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); return internal::solve_retval (*this, MatrixType::Identity(m_qr.rows(), m_qr.cols())); } inline Index rows() const { return m_qr.rows(); } inline Index cols() const { return m_qr.cols(); } /** \returns a const reference to the vector of Householder coefficients used to represent the factor \c Q. * * For advanced uses only. */ const HCoeffsType& hCoeffs() const { return m_hCoeffs; } /** Allows to prescribe a threshold to be used by certain methods, such as rank(), * who need to determine when pivots are to be considered nonzero. This is not used for the * QR decomposition itself. * * When it needs to get the threshold value, Eigen calls threshold(). By default, this * uses a formula to automatically determine a reasonable threshold. * Once you have called the present method setThreshold(const RealScalar&), * your value is used instead. * * \param threshold The new value to use as the threshold. * * A pivot will be considered nonzero if its absolute value is strictly greater than * \f$\vert pivot \vert \leqslant threshold \times \vert maxpivot \vert \f$ * where maxpivot is the biggest pivot. * * If you want to come back to the default behavior, call setThreshold(Default_t) */ ColPivHouseholderQR& setThreshold(const RealScalar& threshold) { m_usePrescribedThreshold = true; m_prescribedThreshold = threshold; return *this; } /** Allows to come back to the default behavior, letting Eigen use its default formula for * determining the threshold. * * You should pass the special object Eigen::Default as parameter here. * \code qr.setThreshold(Eigen::Default); \endcode * * See the documentation of setThreshold(const RealScalar&). */ ColPivHouseholderQR& setThreshold(Default_t) { m_usePrescribedThreshold = false; return *this; } /** Returns the threshold that will be used by certain methods such as rank(). * * See the documentation of setThreshold(const RealScalar&). */ RealScalar threshold() const { eigen_assert(m_isInitialized || m_usePrescribedThreshold); return m_usePrescribedThreshold ? m_prescribedThreshold // this formula comes from experimenting (see "LU precision tuning" thread on the list) // and turns out to be identical to Higham's formula used already in LDLt. : NumTraits::epsilon() * RealScalar(m_qr.diagonalSize()); } /** \returns the number of nonzero pivots in the QR decomposition. * Here nonzero is meant in the exact sense, not in a fuzzy sense. * So that notion isn't really intrinsically interesting, but it is * still useful when implementing algorithms. * * \sa rank() */ inline Index nonzeroPivots() const { eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); return m_nonzero_pivots; } /** \returns the absolute value of the biggest pivot, i.e. the biggest * diagonal coefficient of R. */ RealScalar maxPivot() const { return m_maxpivot; } /** \brief Reports whether the QR factorization was succesful. * * \note This function always returns \c Success. It is provided for compatibility * with other factorization routines. * \returns \c Success */ ComputationInfo info() const { eigen_assert(m_isInitialized && "Decomposition is not initialized."); return Success; } protected: static void check_template_parameters() { EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar); } MatrixType m_qr; HCoeffsType m_hCoeffs; PermutationType m_colsPermutation; IntRowVectorType m_colsTranspositions; RowVectorType m_temp; RealRowVectorType m_colSqNorms; bool m_isInitialized, m_usePrescribedThreshold; RealScalar m_prescribedThreshold, m_maxpivot; Index m_nonzero_pivots; Index m_det_pq; }; template typename MatrixType::RealScalar ColPivHouseholderQR::absDeterminant() const { using std::abs; eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!"); return abs(m_qr.diagonal().prod()); } template typename MatrixType::RealScalar ColPivHouseholderQR::logAbsDeterminant() const { eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!"); return m_qr.diagonal().cwiseAbs().array().log().sum(); } /** Performs the QR factorization of the given matrix \a matrix. The result of * the factorization is stored into \c *this, and a reference to \c *this * is returned. * * \sa class ColPivHouseholderQR, ColPivHouseholderQR(const MatrixType&) */ template ColPivHouseholderQR& ColPivHouseholderQR::compute(const MatrixType& matrix) { check_template_parameters(); using std::abs; Index rows = matrix.rows(); Index cols = matrix.cols(); Index size = matrix.diagonalSize(); // the column permutation is stored as int indices, so just to be sure: eigen_assert(cols<=NumTraits::highest()); m_qr = matrix; m_hCoeffs.resize(size); m_temp.resize(cols); m_colsTranspositions.resize(matrix.cols()); Index number_of_transpositions = 0; m_colSqNorms.resize(cols); for(Index k = 0; k < cols; ++k) m_colSqNorms.coeffRef(k) = m_qr.col(k).squaredNorm(); RealScalar threshold_helper = m_colSqNorms.maxCoeff() * numext::abs2(NumTraits::epsilon()) / RealScalar(rows); m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case) m_maxpivot = RealScalar(0); for(Index k = 0; k < size; ++k) { // first, we look up in our table m_colSqNorms which column has the biggest squared norm Index biggest_col_index; RealScalar biggest_col_sq_norm = m_colSqNorms.tail(cols-k).maxCoeff(&biggest_col_index); biggest_col_index += k; // since our table m_colSqNorms accumulates imprecision at every step, we must now recompute // the actual squared norm of the selected column. // Note that not doing so does result in solve() sometimes returning inf/nan values // when running the unit test with 1000 repetitions. biggest_col_sq_norm = m_qr.col(biggest_col_index).tail(rows-k).squaredNorm(); // we store that back into our table: it can't hurt to correct our table. m_colSqNorms.coeffRef(biggest_col_index) = biggest_col_sq_norm; // Track the number of meaningful pivots but do not stop the decomposition to make // sure that the initial matrix is properly reproduced. See bug 941. if(m_nonzero_pivots==size && biggest_col_sq_norm < threshold_helper * RealScalar(rows-k)) m_nonzero_pivots = k; // apply the transposition to the columns m_colsTranspositions.coeffRef(k) = biggest_col_index; if(k != biggest_col_index) { m_qr.col(k).swap(m_qr.col(biggest_col_index)); std::swap(m_colSqNorms.coeffRef(k), m_colSqNorms.coeffRef(biggest_col_index)); ++number_of_transpositions; } // generate the householder vector, store it below the diagonal RealScalar beta; m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta); // apply the householder transformation to the diagonal coefficient m_qr.coeffRef(k,k) = beta; // remember the maximum absolute value of diagonal coefficients if(abs(beta) > m_maxpivot) m_maxpivot = abs(beta); // apply the householder transformation m_qr.bottomRightCorner(rows-k, cols-k-1) .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1)); // update our table of squared norms of the columns m_colSqNorms.tail(cols-k-1) -= m_qr.row(k).tail(cols-k-1).cwiseAbs2(); } m_colsPermutation.setIdentity(PermIndexType(cols)); for(PermIndexType k = 0; k < size/*m_nonzero_pivots*/; ++k) m_colsPermutation.applyTranspositionOnTheRight(k, PermIndexType(m_colsTranspositions.coeff(k))); m_det_pq = (number_of_transpositions%2) ? -1 : 1; m_isInitialized = true; return *this; } namespace internal { template struct solve_retval, Rhs> : solve_retval_base, Rhs> { EIGEN_MAKE_SOLVE_HELPERS(ColPivHouseholderQR<_MatrixType>,Rhs) template void evalTo(Dest& dst) const { eigen_assert(rhs().rows() == dec().rows()); const Index cols = dec().cols(), nonzero_pivots = dec().nonzeroPivots(); if(nonzero_pivots == 0) { dst.setZero(); return; } typename Rhs::PlainObject c(rhs()); // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T c.applyOnTheLeft(householderSequence(dec().matrixQR(), dec().hCoeffs()) .setLength(dec().nonzeroPivots()) .transpose() ); dec().matrixR() .topLeftCorner(nonzero_pivots, nonzero_pivots) .template triangularView() .solveInPlace(c.topRows(nonzero_pivots)); for(Index i = 0; i < nonzero_pivots; ++i) dst.row(dec().colsPermutation().indices().coeff(i)) = c.row(i); for(Index i = nonzero_pivots; i < cols; ++i) dst.row(dec().colsPermutation().indices().coeff(i)).setZero(); } }; } // end namespace internal /** \returns the matrix Q as a sequence of householder transformations. * You can extract the meaningful part only by using: * \code qr.householderQ().setLength(qr.nonzeroPivots()) \endcode*/ template typename ColPivHouseholderQR::HouseholderSequenceType ColPivHouseholderQR ::householderQ() const { eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate()); } /** \return the column-pivoting Householder QR decomposition of \c *this. * * \sa class ColPivHouseholderQR */ template const ColPivHouseholderQR::PlainObject> MatrixBase::colPivHouseholderQr() const { return ColPivHouseholderQR(eval()); } } // end namespace Eigen #endif // EIGEN_COLPIVOTINGHOUSEHOLDERQR_H