src/SPQRSupport/SuiteSparseQRSupport.h

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2012 Desire Nuentsa <desire.nuentsa_wakam@inria.fr>
0005 // Copyright (C) 2014 Gael Guennebaud <gael.guennebaud@inria.fr>
0006 //
0007 // This Source Code Form is subject to the terms of the Mozilla
0008 // Public License v. 2.0. If a copy of the MPL was not distributed
0009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0010
0011 #ifndef EIGEN_SUITESPARSEQRSUPPORT_H
0012 #define EIGEN_SUITESPARSEQRSUPPORT_H
0013
0014 namespace Eigen {
0015
0016   template<typename MatrixType> class SPQR;
0017   template<typename SPQRType> struct SPQRMatrixQReturnType;
0018   template<typename SPQRType> struct SPQRMatrixQTransposeReturnType;
0019   template <typename SPQRType, typename Derived> struct SPQR_QProduct;
0020   namespace internal {
0021     template <typename SPQRType> struct traits<SPQRMatrixQReturnType<SPQRType> >
0022     {
0023       typedef typename SPQRType::MatrixType ReturnType;
0024     };
0025     template <typename SPQRType> struct traits<SPQRMatrixQTransposeReturnType<SPQRType> >
0026     {
0027       typedef typename SPQRType::MatrixType ReturnType;
0028     };
0029     template <typename SPQRType, typename Derived> struct traits<SPQR_QProduct<SPQRType, Derived> >
0030     {
0031       typedef typename Derived::PlainObject ReturnType;
0032     };
0033   } // End namespace internal
0034
0035 /**
0036   * \ingroup SPQRSupport_Module
0037   * \class SPQR
0038   * \brief Sparse QR factorization based on SuiteSparseQR library
0039   *
0040   * This class is used to perform a multithreaded and multifrontal rank-revealing QR decomposition
0041   * of sparse matrices. The result is then used to solve linear leasts_square systems.
0042   * Clearly, a QR factorization is returned such that A*P = Q*R where :
0043   *
0044   * P is the column permutation. Use colsPermutation() to get it.
0045   *
0046   * Q is the orthogonal matrix represented as Householder reflectors.
0047   * Use matrixQ() to get an expression and matrixQ().transpose() to get the transpose.
0048   * You can then apply it to a vector.
0049   *
0050   * R is the sparse triangular factor. Use matrixQR() to get it as SparseMatrix.
0051   * NOTE : The Index type of R is always SuiteSparse_long. You can get it with SPQR::Index
0052   *
0053   * \tparam _MatrixType The type of the sparse matrix A, must be a column-major SparseMatrix<>
0054   *
0055   * \implsparsesolverconcept
0056   *
0057   *
0058   */
0059 template<typename _MatrixType>
0060 class SPQR : public SparseSolverBase<SPQR<_MatrixType> >
0061 {
0062   protected:
0063     typedef SparseSolverBase<SPQR<_MatrixType> > Base;
0064     using Base::m_isInitialized;
0065   public:
0066     typedef typename _MatrixType::Scalar Scalar;
0067     typedef typename _MatrixType::RealScalar RealScalar;
0068     typedef SuiteSparse_long StorageIndex ;
0069     typedef SparseMatrix<Scalar, ColMajor, StorageIndex> MatrixType;
0070     typedef Map<PermutationMatrix<Dynamic, Dynamic, StorageIndex> > PermutationType;
0071     enum {
0072       ColsAtCompileTime = Dynamic,
0073       MaxColsAtCompileTime = Dynamic
0074     };
0075   public:
0076     SPQR()
0077       : m_analysisIsOk(false),
0078         m_factorizationIsOk(false),
0079         m_isRUpToDate(false),
0080         m_ordering(SPQR_ORDERING_DEFAULT),
0081         m_allow_tol(SPQR_DEFAULT_TOL),
0082         m_tolerance (NumTraits<Scalar>::epsilon()),
0083         m_cR(0),
0084         m_E(0),
0085         m_H(0),
0086         m_HPinv(0),
0087         m_HTau(0),
0088         m_useDefaultThreshold(true)
0089     {
0090       cholmod_l_start(&m_cc);
0091     }
0092
0093     explicit SPQR(const _MatrixType& matrix)
0094       : m_analysisIsOk(false),
0095         m_factorizationIsOk(false),
0096         m_isRUpToDate(false),
0097         m_ordering(SPQR_ORDERING_DEFAULT),
0098         m_allow_tol(SPQR_DEFAULT_TOL),
0099         m_tolerance (NumTraits<Scalar>::epsilon()),
0100         m_cR(0),
0101         m_E(0),
0102         m_H(0),
0103         m_HPinv(0),
0104         m_HTau(0),
0105         m_useDefaultThreshold(true)
0106     {
0107       cholmod_l_start(&m_cc);
0108       compute(matrix);
0109     }
0110
0111     ~SPQR()
0112     {
0113       SPQR_free();
0114       cholmod_l_finish(&m_cc);
0115     }
0116     void SPQR_free()
0117     {
0118       cholmod_l_free_sparse(&m_H, &m_cc);
0119       cholmod_l_free_sparse(&m_cR, &m_cc);
0120       cholmod_l_free_dense(&m_HTau, &m_cc);
0121       std::free(m_E);
0122       std::free(m_HPinv);
0123     }
0124
0125     void compute(const _MatrixType& matrix)
0126     {
0127       if(m_isInitialized) SPQR_free();
0128
0129       MatrixType mat(matrix);
0130
0131       /* Compute the default threshold as in MatLab, see:
0132        * Tim Davis, "Algorithm 915, SuiteSparseQR: Multifrontal Multithreaded Rank-Revealing
0133        * Sparse QR Factorization, ACM Trans. on Math. Soft. 38(1), 2011, Page 8:3
0134        */
0135       RealScalar pivotThreshold = m_tolerance;
0136       if(m_useDefaultThreshold)
0137       {
0138         RealScalar max2Norm = 0.0;
0139         for (int j = 0; j < mat.cols(); j++) max2Norm = numext::maxi(max2Norm, mat.col(j).norm());
0140         if(max2Norm==RealScalar(0))
0141           max2Norm = RealScalar(1);
0142         pivotThreshold = 20 * (mat.rows() + mat.cols()) * max2Norm * NumTraits<RealScalar>::epsilon();
0143       }
0144       cholmod_sparse A;
0145       A = viewAsCholmod(mat);
0146       m_rows = matrix.rows();
0147       Index col = matrix.cols();
0148       m_rank = SuiteSparseQR<Scalar>(m_ordering, pivotThreshold, col, &A,
0149                              &m_cR, &m_E, &m_H, &m_HPinv, &m_HTau, &m_cc);
0150
0151       if (!m_cR)
0152       {
0153         m_info = NumericalIssue;
0154         m_isInitialized = false;
0155         return;
0156       }
0157       m_info = Success;
0158       m_isInitialized = true;
0159       m_isRUpToDate = false;
0160     }
0161     /**
0162      * Get the number of rows of the input matrix and the Q matrix
0163      */
0164     inline Index rows() const {return m_rows; }
0165
0166     /**
0167      * Get the number of columns of the input matrix.
0168      */
0169     inline Index cols() const { return m_cR->ncol; }
0170
0171     template<typename Rhs, typename Dest>
0172     void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
0173     {
0174       eigen_assert(m_isInitialized && " The QR factorization should be computed first, call compute()");
0175       eigen_assert(b.cols()==1 && "This method is for vectors only");
0176
0177       //Compute Q^T * b
0178       typename Dest::PlainObject y, y2;
0179       y = matrixQ().transpose() * b;
0180
0181       // Solves with the triangular matrix R
0182       Index rk = this->rank();
0183       y2 = y;
0184       y.resize((std::max)(cols(),Index(y.rows())),y.cols());
0185       y.topRows(rk) = this->matrixR().topLeftCorner(rk, rk).template triangularView<Upper>().solve(y2.topRows(rk));
0186
0187       // Apply the column permutation
0188       // colsPermutation() performs a copy of the permutation,
0189       // so let's apply it manually:
0190       for(Index i = 0; i < rk; ++i) dest.row(m_E[i]) = y.row(i);
0191       for(Index i = rk; i < cols(); ++i) dest.row(m_E[i]).setZero();
0192
0193 //       y.bottomRows(y.rows()-rk).setZero();
0194 //       dest = colsPermutation() * y.topRows(cols());
0195
0196       m_info = Success;
0197     }
0198
0199     /** \returns the sparse triangular factor R. It is a sparse matrix
0200      */
0201     const MatrixType matrixR() const
0202     {
0203       eigen_assert(m_isInitialized && " The QR factorization should be computed first, call compute()");
0204       if(!m_isRUpToDate) {
0205         m_R = viewAsEigen<Scalar,ColMajor, typename MatrixType::StorageIndex>(*m_cR);
0206         m_isRUpToDate = true;
0207       }
0208       return m_R;
0209     }
0210     /// Get an expression of the matrix Q
0211     SPQRMatrixQReturnType<SPQR> matrixQ() const
0212     {
0213       return SPQRMatrixQReturnType<SPQR>(*this);
0214     }
0215     /// Get the permutation that was applied to columns of A
0216     PermutationType colsPermutation() const
0217     {
0218       eigen_assert(m_isInitialized && "Decomposition is not initialized.");
0219       return PermutationType(m_E, m_cR->ncol);
0220     }
0221     /**
0222      * Gets the rank of the matrix.
0223      * It should be equal to matrixQR().cols if the matrix is full-rank
0224      */
0225     Index rank() const
0226     {
0227       eigen_assert(m_isInitialized && "Decomposition is not initialized.");
0228       return m_cc.SPQR_istat[4];
0229     }
0230     /// Set the fill-reducing ordering method to be used
0231     void setSPQROrdering(int ord) { m_ordering = ord;}
0232     /// Set the tolerance tol to treat columns with 2-norm < =tol as zero
0233     void setPivotThreshold(const RealScalar& tol)
0234     {
0235       m_useDefaultThreshold = false;
0236       m_tolerance = tol;
0237     }
0238
0239     /** \returns a pointer to the SPQR workspace */
0240     cholmod_common *cholmodCommon() const { return &m_cc; }
0241
0242
0243     /** \brief Reports whether previous computation was successful.
0244       *
0245       * \returns \c Success if computation was successful,
0246       *          \c NumericalIssue if the sparse QR can not be computed
0247       */
0248     ComputationInfo info() const
0249     {
0250       eigen_assert(m_isInitialized && "Decomposition is not initialized.");
0251       return m_info;
0252     }
0253   protected:
0254     bool m_analysisIsOk;
0255     bool m_factorizationIsOk;
0256     mutable bool m_isRUpToDate;
0257     mutable ComputationInfo m_info;
0258     int m_ordering; // Ordering method to use, see SPQR's manual
0259     int m_allow_tol; // Allow to use some tolerance during numerical factorization.
0260     RealScalar m_tolerance; // treat columns with 2-norm below this tolerance as zero
0261     mutable cholmod_sparse *m_cR; // The sparse R factor in cholmod format
0262     mutable MatrixType m_R; // The sparse matrix R in Eigen format
0263     mutable StorageIndex *m_E; // The permutation applied to columns
0264     mutable cholmod_sparse *m_H;  //The householder vectors
0265     mutable StorageIndex *m_HPinv; // The row permutation of H
0266     mutable cholmod_dense *m_HTau; // The Householder coefficients
0267     mutable Index m_rank; // The rank of the matrix
0268     mutable cholmod_common m_cc; // Workspace and parameters
0269     bool m_useDefaultThreshold;     // Use default threshold
0270     Index m_rows;
0271     template<typename ,typename > friend struct SPQR_QProduct;
0272 };
0273
0274 template <typename SPQRType, typename Derived>
0275 struct SPQR_QProduct : ReturnByValue<SPQR_QProduct<SPQRType,Derived> >
0276 {
0277   typedef typename SPQRType::Scalar Scalar;
0278   typedef typename SPQRType::StorageIndex StorageIndex;
0279   //Define the constructor to get reference to argument types
0280   SPQR_QProduct(const SPQRType& spqr, const Derived& other, bool transpose) : m_spqr(spqr),m_other(other),m_transpose(transpose) {}
0281
0282   inline Index rows() const { return m_transpose ? m_spqr.rows() : m_spqr.cols(); }
0283   inline Index cols() const { return m_other.cols(); }
0284   // Assign to a vector
0285   template<typename ResType>
0286   void evalTo(ResType& res) const
0287   {
0288     cholmod_dense y_cd;
0289     cholmod_dense *x_cd;
0290     int method = m_transpose ? SPQR_QTX : SPQR_QX;
0291     cholmod_common *cc = m_spqr.cholmodCommon();
0292     y_cd = viewAsCholmod(m_other.const_cast_derived());
0293     x_cd = SuiteSparseQR_qmult<Scalar>(method, m_spqr.m_H, m_spqr.m_HTau, m_spqr.m_HPinv, &y_cd, cc);
0294     res = Matrix<Scalar,ResType::RowsAtCompileTime,ResType::ColsAtCompileTime>::Map(reinterpret_cast<Scalar*>(x_cd->x), x_cd->nrow, x_cd->ncol);
0295     cholmod_l_free_dense(&x_cd, cc);
0296   }
0297   const SPQRType& m_spqr;
0298   const Derived& m_other;
0299   bool m_transpose;
0300
0301 };
0302 template<typename SPQRType>
0303 struct SPQRMatrixQReturnType{
0304
0305   SPQRMatrixQReturnType(const SPQRType& spqr) : m_spqr(spqr) {}
0306   template<typename Derived>
0307   SPQR_QProduct<SPQRType, Derived> operator*(const MatrixBase<Derived>& other)
0308   {
0309     return SPQR_QProduct<SPQRType,Derived>(m_spqr,other.derived(),false);
0310   }
0311   SPQRMatrixQTransposeReturnType<SPQRType> adjoint() const
0312   {
0313     return SPQRMatrixQTransposeReturnType<SPQRType>(m_spqr);
0314   }
0315   // To use for operations with the transpose of Q
0316   SPQRMatrixQTransposeReturnType<SPQRType> transpose() const
0317   {
0318     return SPQRMatrixQTransposeReturnType<SPQRType>(m_spqr);
0319   }
0320   const SPQRType& m_spqr;
0321 };
0322
0323 template<typename SPQRType>
0324 struct SPQRMatrixQTransposeReturnType{
0325   SPQRMatrixQTransposeReturnType(const SPQRType& spqr) : m_spqr(spqr) {}
0326   template<typename Derived>
0327   SPQR_QProduct<SPQRType,Derived> operator*(const MatrixBase<Derived>& other)
0328   {
0329     return SPQR_QProduct<SPQRType,Derived>(m_spqr,other.derived(), true);
0330   }
0331   const SPQRType& m_spqr;
0332 };
0333
0334 }// End namespace Eigen
0335 #endif