使用 GMP 类型求解 Eigen3 中的线性系统

Question

[我知道有两个问题与 Eigen3 和 GMP 相关，但它们没有解决我的问题。]

我想在 Eigen3 中做任意精度的线性代数。因此，我想使用 GMP 的 mpq_class 作为 Eigen 矩阵 class 中的标量类型。幸运的是，文档提供了有关如何执行此操作的 an example。不幸的是，一旦我想计算 Householder 分解，它就会失败。特别是考虑以下示例（几乎与上面链接的示例相同）：

#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();
}

现在发生的事情是，我被 C++ 臭名昭著的模板错误消息之一轰炸了，也就是说（eigentest.cpp 是我从中编译的文件）：

/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); }

我不知道这里发生了什么，更不知道如何解决它。非常感谢任何帮助（注意：帮助还可能包括以任意精度将我指向线性代数库的 C++ 接口）。

Answer 1

解密编译器的错误信息你会看到它试图对有理数取平方根（并将其存储为有理数），这是不可能的。

如果您想使用有理数，则需要进行不需要平方根的分解，例如任何类型的 LU 分解或 LDLt 分解（如果您的矩阵是对称的和半定的） .

使用 GMP 类型求解 Eigen3 中的线性系统

Solving linear systems in Eigen3 with GMP types

eigen

eigen3