Search code examples
eigeneigen3

Solving linear systems in Eigen3 with GMP types

[I know there are two questions related to Eigen3 and GMP, but they are not addressing my problem.]

I'm trying to do arbitrary-precision linear algebra in Eigen3. Therefore, I'd like to use GMP's mpq_class as scalar type in Eigen's matrix class. Luckily, the documentation provides an example on how to do this. Less luckily, it fails as soon as I want to compute a Householder decomposition. In particular, consider the following example (almost identical to the example linked to above):

#include <eigen3/Eigen/Dense>
#include <eigen3/Eigen/Core>
#include <gmpxx.h>

typedef mpq_class entry_t;

namespace Eigen {
template<> struct NumTraits<entry_t> : GenericNumTraits<entry_t>
{
        typedef entry_t Real;
        typedef entry_t NonInteger;
        typedef entry_t Nested;
        static inline Real epsilon() { return 0; }
        static inline Real dummy_precision() { return 0; }
    static inline Real digits10() { return 0; }

        enum {
        IsInteger = 0,
        IsSigned = 1,
        IsComplex = 0,
        RequireInitialization = 1,
        ReadCost = 6,
        AddCost = 150,
        MulCost = 100
    };
};
}

int main()
{
    Eigen::Matrix<entry_t,1,1> mat;
    mat.householderQr();
}

What happens now is that I am bombarded with one of C++'s infamous template error messages, to wit (with eigentest.cpp being the file I compiled this from):

/usr/include/gmpxx.h: In instantiation of ‘void __gmp_expr<T, __gmp_unary_expr<__gmp_expr<T, U>, Op> >::eval(typename __gmp_resolve_expr<T>::ptr_type) const [with T = __mpq_struct [1]; U = __gmp_binary_expr<__gmp_expr<__mpq_struct [1], __mpq_struct [1]>, __gmp_expr<__mpq_struct [1], __mpq_struct [1]>, __gmp_binary_plus>; Op = __gmp_sqrt_function; typename __gmp_resolve_expr<T>::ptr_type = __mpq_struct*]’:
/usr/include/gmpxx.h:2130:3:   required from ‘void __gmp_set_expr(mpq_ptr, const __gmp_expr<__mpq_struct [1], T>&) [with T = __gmp_unary_expr<__gmp_expr<__mpq_struct [1], __gmp_binary_expr<__gmp_expr<__mpq_struct [1], __mpq_struct [1]>, __gmp_expr<__mpq_struct [1], __mpq_struct [1]>, __gmp_binary_plus> >, __gmp_sqrt_function>; mpq_ptr = __mpq_struct*]’
/usr/include/gmpxx.h:1724:28:   required from ‘__gmp_expr<__mpq_struct [1], __mpq_struct [1]>& __gmp_expr<__mpq_struct [1], __mpq_struct [1]>::operator=(const __gmp_expr<T, U>&) [with T = __mpq_struct [1]; U = __gmp_unary_expr<__gmp_expr<__mpq_struct [1], __gmp_binary_expr<__gmp_expr<__mpq_struct [1], __mpq_struct [1]>, __gmp_expr<__mpq_struct [1], __mpq_struct [1]>, __gmp_binary_plus> >, __gmp_sqrt_function>]’
/usr/include/eigen3/Eigen/src/Householder/Householder.h:87:10:   required from ‘void Eigen::MatrixBase<Derived>::makeHouseholder(EssentialPart&, Eigen::MatrixBase<Derived>::Scalar&, Eigen::MatrixBase<Derived>::RealScalar&) const [with EssentialPart = Eigen::VectorBlock<Eigen::Block<Eigen::Block<Eigen::Block<Eigen::Matrix<__gmp_expr<__mpq_struct [1], __mpq_struct [1]>, 4, 4>, -1, -1, false>, -1, 1, true>, -1, 1, false>, -1>; Derived = Eigen::Block<Eigen::Block<Eigen::Block<Eigen::Matrix<__gmp_expr<__mpq_struct [1], __mpq_struct [1]>, 4, 4>, -1, -1, false>, -1, 1, true>, -1, 1, false>; Eigen::MatrixBase<Derived>::Scalar = __gmp_expr<__mpq_struct [1], __mpq_struct [1]>; Eigen::MatrixBase<Derived>::RealScalar = __gmp_expr<__mpq_struct [1], __mpq_struct [1]>]’
/usr/include/eigen3/Eigen/src/Householder/Householder.h:45:43:   required from ‘void Eigen::MatrixBase<Derived>::makeHouseholderInPlace(Eigen::MatrixBase<Derived>::Scalar&, Eigen::MatrixBase<Derived>::RealScalar&) [with Derived = Eigen::Block<Eigen::Block<Eigen::Block<Eigen::Matrix<__gmp_expr<__mpq_struct [1], __mpq_struct [1]>, 4, 4>, -1, -1, false>, -1, 1, true>, -1, 1, false>; Eigen::MatrixBase<Derived>::Scalar = __gmp_expr<__mpq_struct [1], __mpq_struct [1]>; Eigen::MatrixBase<Derived>::RealScalar = __gmp_expr<__mpq_struct [1], __mpq_struct [1]>]’
/usr/include/eigen3/Eigen/src/QR/HouseholderQR.h:244:5:   required from ‘void Eigen::internal::householder_qr_inplace_unblocked(MatrixQR&, HCoeffs&, typename MatrixQR::Scalar*) [with MatrixQR = Eigen::Block<Eigen::Matrix<__gmp_expr<__mpq_struct [1], __mpq_struct [1]>, 4, 4>, -1, -1, false>; HCoeffs = Eigen::Block<Eigen::Matrix<__gmp_expr<__mpq_struct [1], __mpq_struct [1]>, 4, 1, 0, 4, 1>, -1, 1, false>; typename MatrixQR::Scalar = __gmp_expr<__mpq_struct [1], __mpq_struct [1]>]’
/usr/include/eigen3/Eigen/src/QR/HouseholderQR.h:295:70:   required from ‘void Eigen::internal::householder_qr_inplace_blocked(MatrixQR&, HCoeffs&, typename MatrixQR::Index, typename MatrixQR::Scalar*) [with MatrixQR = Eigen::Matrix<__gmp_expr<__mpq_struct [1], __mpq_struct [1]>, 4, 4>; HCoeffs = Eigen::Matrix<__gmp_expr<__mpq_struct [1], __mpq_struct [1]>, 4, 1, 0, 4, 1>; typename MatrixQR::Index = long int; typename MatrixQR::Scalar = __gmp_expr<__mpq_struct [1], __mpq_struct [1]>]’
/usr/include/eigen3/Eigen/src/QR/HouseholderQR.h:355:78:   required from ‘Eigen::HouseholderQR<_MatrixType>& Eigen::HouseholderQR<MatrixType>::compute(const MatrixType&) [with _MatrixType = Eigen::Matrix<__gmp_expr<__mpq_struct [1], __mpq_struct [1]>, 4, 4>; Eigen::HouseholderQR<MatrixType>::MatrixType = Eigen::Matrix<__gmp_expr<__mpq_struct [1], __mpq_struct [1]>, 4, 4>]’
/usr/include/eigen3/Eigen/src/QR/HouseholderQR.h:100:21:   required from ‘Eigen::HouseholderQR<MatrixType>::HouseholderQR(const MatrixType&) [with _MatrixType = Eigen::Matrix<__gmp_expr<__mpq_struct [1], __mpq_struct [1]>, 4, 4>; Eigen::HouseholderQR<MatrixType>::MatrixType = Eigen::Matrix<__gmp_expr<__mpq_struct [1], __mpq_struct [1]>, 4, 4>]’
/usr/include/eigen3/Eigen/src/QR/HouseholderQR.h:369:43:   required from ‘const Eigen::HouseholderQR<Eigen::Matrix<typename Eigen::internal::traits<T>::Scalar, Eigen::internal::traits<T>::RowsAtCompileTime, Eigen::internal::traits<T>::ColsAtCompileTime, (AutoAlign | ((Eigen::internal::traits<T>::Flags & Eigen::RowMajorBit) ? RowMajor :  ColMajor)), Eigen::internal::traits<T>::MaxRowsAtCompileTime, Eigen::internal::traits<T>::MaxColsAtCompileTime> > Eigen::MatrixBase<Derived>::householderQr() const [with Derived = Eigen::Matrix<__gmp_expr<__mpq_struct [1], __mpq_struct [1]>, 4, 4>; typename Eigen::internal::traits<T>::Scalar = __gmp_expr<__mpq_struct [1], __mpq_struct [1]>]’
eigentest.cpp:33:20:   required from here
/usr/include/gmpxx.h:2236:36: error: no matching function for call to ‘__gmp_sqrt_function::eval(__mpq_struct*&, __mpq_struct*&)’
   { expr.val.eval(p); Op::eval(p, p); }
                                    ^
/usr/include/gmpxx.h:2236:36: note: candidates are:
/usr/include/gmpxx.h:1109:15: note: static void __gmp_sqrt_function::eval(mpz_ptr, mpz_srcptr)
   static void eval(mpz_ptr z, mpz_srcptr w) { mpz_sqrt(z, w); }
               ^
/usr/include/gmpxx.h:1109:15: note:   no known conversion for argument 1 from ‘__gmp_resolve_expr<__mpq_struct [1]>::ptr_type {aka __mpq_struct*}’ to ‘mpz_ptr {aka __mpz_struct*}’
/usr/include/gmpxx.h:1110:15: note: static void __gmp_sqrt_function::eval(mpf_ptr, mpf_srcptr)
   static void eval(mpf_ptr f, mpf_srcptr g) { mpf_sqrt(f, g); }

I have no idea what's going on here, much less how to solve it. Any help is greatly appreciated (note: help may also include pointing me to a C++ interface of a linear algebra library with arbitrary precision).

Solution

  • Decrypting the error message of the compiler you will see that it tries to take a square root of a rational number (and store it into a rational number), which is not possible.

    If you want to work with rational numbers, you need to take a decomposition which does not require square-roots, like any kind of LU decomposition, or the LDLt decomposition (if your matrix is symmetric and semi-definite).