boost::multiprecision::mpfr_float 的 Eigen3 动态矩阵

Question

我想制作矩阵并使用 Eigen3 库使用它们，我的数字类型是 Boost.Multiprecision 的 mpfr_float 包装器。我可以很好地制作矩阵，但是除了矩阵加法之外，我尝试过的所有操作都失败了。仅将两个单位矩阵相乘会产生垃圾结果！

这是一个 MWE：

#include <eigen3/Eigen/Dense>
#include <eigen3/Eigen/LU>
#include <boost/multiprecision/mpfr.hpp>
#include <iostream>

namespace Eigen{

using boost::multiprecision::mpfr_float;
    template<> struct NumTraits<boost::multiprecision::mpfr_float>
    {

    typedef boost::multiprecision::mpfr_float Real;
    typedef boost::multiprecision::mpfr_float NonInteger;
    typedef boost::multiprecision::mpfr_float Nested;
    enum {
        IsComplex = 0,
        IsInteger = 0,
        IsSigned = 1,
        RequireInitialization = 1,
        ReadCost = 20, //these have no impact
        AddCost = 30,
        MulCost = 40
    };


    inline static Real highest() { // these seem to have no impact
        return (mpfr_float(1) - epsilon()) * pow(mpfr_float(2),mpfr_get_emax()-1);
    }

    inline static Real lowest() {
        return -highest();
    }

    inline static Real dummy_precision(){
        return pow(mpfr_float(10),-int(mpfr_float::default_precision()-3));
    }

    inline static Real epsilon(){
        return pow(mpfr_float(10),-int( mpfr_float::default_precision()));
    }
    //http://www.manpagez.com/info/mpfr/mpfr-2.3.2/mpfr_31.php
};
} // namespace eigen


int main()
{
    int size = 10;

    typedef Eigen::Matrix<boost::multiprecision::mpfr_float, Eigen::Dynamic, Eigen::Dynamic> mp_matrix;
    mp_matrix A = mp_matrix::Identity(size, size);
    std::cout << A * A << std::endl; // produces nan's every other row!!!

    return 0;
}

它生成的单位矩阵很好，但是在我的机器上，使用此代码依赖项的最新自制软件分布式版本（和其他版本）（Boost 1.57，Eigen 3.2.4），我的程序生成 NaN 的每个矩阵中的另一行：

  1   0   0   0   0   0   0   0   0   0
nan nan nan nan nan nan nan nan nan nan
  0   0   1   0   0   0   0   0   0   0
nan nan nan nan nan nan nan nan nan nan
  0   0   0   0   1   0   0   0   0   0
nan nan nan nan nan nan nan nan nan nan
  0   0   0   0   0   0   1   0   0   0
nan nan nan nan nan nan nan nan nan nan
  0   0   0   0   0   0   0   0   1   0
nan nan nan nan nan nan nan nan nan nan

奇数矩阵大小在底部产生两行 nan...

这似乎不依赖于默认精度，也不依赖于我定义的 NumTraits 结构的细节，甚至不依赖于我是否定义了一个。我可以从 GenericTraits<mpfr_float> 继承，也可以不继承；我可以说 RequireInitialization = 1，或 0。我得到了 NaN。如果我尝试通过 LU 求逆来求解系统，则返回的矩阵完全是 NaN。如果矩阵的大小是 1x1，我什至可以从矩阵乘法中得到一个 NaN。更改各种静态函数也没有影响。

我觉得最奇怪的部分是，如果我定义一个自定义复合体 class（不是 std::complex，因为数据丢失原因），将 mpfr_float 作为基础类型对于实部和虚部，我确实得到了函数矩阵。

编辑：这里是复杂类型的 NumTraits：

/**
 \brief this templated struct permits us to use the Float type in Eigen matrices.
 */
template<> struct NumTraits<mynamespace::complex> : NumTraits<boost::multiprecision::mpfr_float> // permits to get the epsilon, dummy_precision, lowest, highest functions
{
    typedef boost::multiprecision::mpfr_float Real;
    typedef boost::multiprecision::mpfr_float NonInteger;
    typedef mynamespace::complex Nested;
    enum {
        IsComplex = 1,
        IsInteger = 0,
        IsSigned = 1,
        RequireInitialization = 1, // yes, require initialization, otherwise get crashes
        ReadCost = 2 * NumTraits<Real>::ReadCost,
        AddCost = 2 * NumTraits<Real>::AddCost,
        MulCost = 4 * NumTraits<Real>::MulCost + 2 * NumTraits<Real>::AddCost
    };
};

这是复杂的class我写的：

#include <boost/multiprecision/mpfr.hpp>
#include <boost/multiprecision/random.hpp>


#include <boost/archive/text_oarchive.hpp>
#include <boost/archive/text_iarchive.hpp>
#include <boost/serialization/split_member.hpp>

#include <eigen3/Eigen/Core>

#include <assert.h>
namespace mynamespace {
    using boost::multiprecision::mpfr_float;

    class complex {

    private:

        mpfr_float real_, imag_; ///< the real and imaginary parts of the complex number


        // let the boost serialization library have access to the private members of this class.
        friend class boost::serialization::access;

        template<class Archive>
        void save(Archive & ar, const unsigned int version) const {
            // note, version is always the latest when saving
            ar & real_;
            ar & imag_;
        }


        /**
         \brief load method for archiving a bertini::complex
         */
        template<class Archive>
        void load(Archive & ar, const unsigned int version) {
            ar & real_;
            ar & imag_;
        }

        BOOST_SERIALIZATION_SPLIT_MEMBER()


    public:

        complex():real_(), imag_(){}

        complex(double re) : real_(re), imag_("0.0"){}

        complex(const mpfr_float & re) : real_(re), imag_("0.0"){}

        complex(const std::string & re) : real_(re), imag_("0.0"){}

        complex(const mpfr_float & re, const mpfr_float & im) : real_(re), imag_(im) {}

        complex(double re, double im) : real_(re), imag_(im) {}

        complex(const std::string & re, const std::string & im) : real_(re), imag_(im) {}

        complex(const mpfr_float & re, const std::string & im) : real_(re), imag_(im) {}

        complex(const std::string & re, const mpfr_float & im) : real_(re), imag_(im) {}

        complex(complex&& other) : complex() {
            swap(*this, other);
        }

        complex(const complex & other) : real_(other.real_), imag_(other.imag_) {}

        friend void swap(complex& first, complex& second)  {
            using std::swap;
            swap(first.real_,second.real_);
            swap(first.imag_,second.imag_);
        }

        complex& operator=(complex other) {
            swap(*this, other);
            return *this;
        }

        mpfr_float real() const {return real_;}

        mpfr_float imag() const {return imag_;}

        void real(const mpfr_float & new_real){real_ = new_real;}

        void imag(const mpfr_float & new_imag){imag_ = new_imag;}

        void real(const std::string & new_real){real_ = mpfr_float(new_real);}

        void imag(const std::string & new_imag){imag_ = mpfr_float(new_imag);}


        complex& operator+=(const complex & rhs) {
            real_+=rhs.real_;
            imag_+=rhs.imag_;
            return *this;
        }

        complex& operator-=(const complex & rhs) {
            real_-=rhs.real_;
            imag_-=rhs.imag_;
            return *this;
        }

        complex& operator*=(const complex & rhs) {
            mpfr_float a = real_*rhs.real_ - imag_*rhs.imag_; // cache the real part of the result
            imag_ = real_*rhs.imag_ + imag_*rhs.real_;
            real_ = a;
            return *this;
        }

        complex& operator/=(const complex & rhs) {
            mpfr_float d = rhs.abs2();
            mpfr_float a = real_*rhs.real_ + imag_*rhs.imag_; // cache the numerator of the real part of the result
            imag_ = imag_*rhs.real_ - real_*rhs.imag_/d;
            real_ = a/d;

            return *this;
        }

        complex operator-() const 
            return complex(-real(), -imag());
        }

        mpfr_float abs2() const {
            return pow(real(),2)+pow(imag(),2);
        }

        mpfr_float abs() const {
            return sqrt(abs2());
        }

        mpfr_float arg() const {
            return boost::multiprecision::atan2(imag(),real());
        }

        mpfr_float norm() const {
            return abs2();
        }

        complex conj() const {
            return complex(real(), -imag());
        }

        void precision(unsigned int prec) {
            real_.precision(prec);
            imag_.precision(prec);
        }

        unsigned int precision() const {
            assert(real_.precision()==imag_.precision());
            return real_.precision();
        }

        friend std::ostream& operator<<(std::ostream& out, const complex & z) {
            out << "(" << z.real() << "," << z.imag() << ")";
            return out;
        }

        friend std::istream& operator>>(std::istream& in, complex & z) {
            std::string gotten;
            in >> gotten;

            if (gotten[0]=='(') {
                if (*(gotten.end()-1)!=')') {
                    in.setstate(std::ios::failbit);
                    z.real("NaN");
                    z.imag("NaN");
                    return in;
                }
                else{
                    // try to find a comma in the string.
                    size_t comma_pos = gotten.find(",");

                    // if the second character, have no numbers in the real part.
                    // if the second to last character, have no numbers in the imag part.

                    if (comma_pos!=std::string::npos){
                        if (comma_pos==1 || comma_pos==gotten.size()-2) {
                            in.setstate(std::ios::failbit);
                            z.real("NaN");
                            z.imag("NaN");
                            return in;
                        }
                        else{
                            z.real(gotten.substr(1, comma_pos-1));
                            z.imag(gotten.substr(comma_pos+1, gotten.size()-2 - (comma_pos)));
                            return in;
                        }
                    }
                    // did not find a comma
                    else{
                        z.real(gotten.substr(1,gotten.size()-2));
                        z.imag("0.0");
                        return in;
                    }

                }
            }
            else{
                z.real(gotten);
                z.imag("0.0");
                return in;
            }
        }
    }; // end declaration of the mynamespace::complex number class

    inline complex operator+(complex lhs, const complex & rhs){
        lhs += rhs;
        return lhs;
    }

    inline complex operator+(complex lhs, const mpfr_float & rhs) {
        lhs.real(lhs.real()+rhs);
        return lhs;
    }

    inline complex operator+(const mpfr_float & lhs, complex rhs) {
        return rhs+lhs;
    }

    inline complex operator-(complex lhs, const complex & rhs){
        lhs -= rhs;
        return lhs;
    }

    inline complex operator-(complex lhs, const mpfr_float & rhs) {
        lhs.real(lhs.real()-rhs);
        return lhs;
    }

    inline complex operator-(const mpfr_float & lhs, complex rhs) {
        rhs.real(lhs - rhs.real());
        return rhs;
    }

    inline complex operator*(complex lhs, const complex & rhs){
        lhs *= rhs;
        return lhs;
    }

    inline complex operator*(complex lhs, const mpfr_float & rhs) {
        lhs.real(lhs.real()*rhs);
        lhs.imag(lhs.imag()*rhs);
        return lhs;
    }

    inline complex operator*(const mpfr_float & lhs, complex rhs) {
        return rhs*lhs; // it commutes!
    }

    inline complex operator/(complex lhs, const complex & rhs){
        lhs /= rhs;
        return lhs;
    }

    inline complex operator/(complex lhs, const mpfr_float & rhs) {
        lhs.real(lhs.real()/rhs);
        lhs.imag(lhs.imag()/rhs);
        return lhs;
    }

    inline complex operator/(const mpfr_float & lhs, const complex & rhs) {
        mpfr_float d = rhs.abs2();
        return complex(lhs*rhs.real()/d, -lhs*rhs.imag()/d);
    }

    inline mpfr_float real(const complex & z) {
        return z.real();
    }

    inline mpfr_float imag(const complex & z) {
        return z.imag();
    }

    inline complex conj(const complex & z) {
        return z.conj();
    }

    inline mpfr_float abs2(const complex & z) {
        return z.abs2();
    }

    inline mpfr_float abs(const complex & z) {
        return boost::multiprecision::sqrt(abs2(z));
    }

    inline mpfr_float arg(const complex & z) {
        return boost::multiprecision::atan2(z.imag(),z.real());
    }

    inline complex inverse(const complex & z) {
        mpfr_float d = z.abs2();

        return complex(z.real()/d, -z.imag()/d);
    }

    inline complex square(const complex & z) {
        return complex(z.real()*z.real() - z.imag()*z.imag(), mpfr_float("2.0")*z.real()*z.imag());
    }

    inline complex pow(const complex & z, int power) {
        if (power < 0) {
            return pow(inverse(z), -power);
        }
        else if (power==0)
            return complex("1.0","0.0");
        else if(power==1)
            return z;
        else if(power==2)
            return z*z;
        else if(power==3)
            return z*z*z;
        else {
            unsigned int p(power);
            complex result("1.0","0.0"), z_to_the_current_power_of_two = z;
            // have copy of p in memory, can freely modify it.
            do {
                if ( (p & 1) == 1 ) { // get the lowest bit of the number
                    result *= z_to_the_current_power_of_two;
                }
                z_to_the_current_power_of_two *= z_to_the_current_power_of_two; // square z_to_the_current_power_of_two
            } while (p  >>= 1);

            return result;
        }
    }

    inline complex polar(const mpfr_float & rho, const mpfr_float & theta) {
        return complex(rho*cos(theta), rho*sin(theta));
    }

    inline complex sqrt(const complex & z) {
        return polar(sqrt(abs(z)), arg(z)/2);
    }
} // re: namespace

我做错了什么？我可以对 Eigen / NumTraits 等做些什么才能使矩阵运算正常工作？

Answer 1

这种奇怪的 nan 行为是由 boost/multiprecision/mpfr.hpp 文件中的错误引起的，该错误出现在 Boost 版本 1.56、1.57 和可能更早的版本中。赋值运算符未测试自赋值——因此数字设置为 nan。

为什么这每隔一行发生一次，如果矩阵大小为奇数，则在最后一行发生，我仍然不清楚。然而，添加一个自我分配测试，就像 GitHub here 上的提交一样，解决了这个问题。维护者即将发布官方补丁，但 will/did 不会进入 1.58。感谢 Boost 用户邮件列表中的 John 发现了这个错误。

如果您想在 boost/multiprecision/mpfr.hpp 中修复 Boost 安装中的这个错误，只需将引用赋值运算符的内容（不包括 return）包装在 if(this != &o){...} 中即可。

我有一个关于 Eigen 实现的挥之不去的问题 -- 是什么导致了自赋值？是一般问题吗？

boost::multiprecision::mpfr_float 的 Eigen3 动态矩阵

Eigen3 dynamic matrix with boost::multiprecision::mpfr_float

c++

boost

matrix

eigen

multiprecision