c++ - boost 间隔矩阵

Question

我正在尝试使用区间算法为 R^n->R^m 函数编写一个“分而治之”的根估计算法。我以前在 python 中做过这个，但是性能太慢所以我决定试一试 C++，我是一个初学者。

我搜索了一段时间，找到了 Boosts interval library，它似乎可以方便地处理我的问题。但是，我需要用这种间隔的矩阵进行计算。因为稍后我需要使用区间中点矩阵并(伪)反转它们，所以我认为 Eigen 可能是表示我的矩阵的好方法。

我的问题:这种方法真的是个好主意，尤其是在性能方面？如果是这样，我将如何继续使这样的矩阵能够包含区间作为入口并处理基本操作(矩阵乘法等)

我希望能够做的事情的一个例子:

#include <eigen3/Eigen/Dense>
#include <boost/numeric/interval.hpp>
using Eigen::MatrixXd;
using Eigen::VectorXd;
int main()
{
  MatrixXd m(2,2);
  VectorXd v(2);
  Interval i1(0.0, 1.0);
  Interval i2(1.0, 2.0);
  Interval i3(0.0, 0.0);
  m(0,0) = i1;
  m(1,0) = i3;
  m(0,1) = i3;
  m(1,1) = i1;

  v(0) = i1;
  v(1) = i2;

  std::cout << m*v << std::endl;
}

我对区间类型做的规范是

 typedef boost::numeric::interval<double, boost::numeric::interval_lib::policies
                            <boost::numeric::interval_lib::save_state
                            <boost::numeric::interval_lib::rounded_transc_std<double> >,
                            boost::numeric::interval_lib::checking_base<double> > > Interval;

这最终应该输出类似 (0.0,1.0)、(0.0,2.0) 的内容。现在这会引发错误“将标量分配给不可比较的类型间隔”，我不知道如何解决它，因为 Eigens 矩阵类型自然不支持间隔作为条目。

assigning to 'Scalar' (aka 'double') from incompatible type 'Interval'
(aka 'interval < double, boost::numeric::interval_lib::policies <
boost::numeric::interval_lib::save_state <
boost::numeric::interval_lib::rounded_transc_std < double > >,
boost::numeric::interval_lib::checking_base<double> > >')

欢迎任何(新手级!)引用和建议。

Answer 1

Docs say那个

typedef Matrix<double, Dynamic, Dynamic> MatrixXd;

这里有一些修复:

Live On Coliru

#include <boost/numeric/interval.hpp>
#include <boost/numeric/interval/utility.hpp>
#include <boost/numeric/interval/io.hpp>

namespace bn = boost::numeric;
namespace bi = bn::interval_lib;

using Interval = bn::interval<
        double, 
        bi::policies<
            bi::save_state<bi::rounded_transc_std<double> >,
            bi::checking_base<double>
        > 
    >;

#include <eigen3/Eigen/Dense>
using Matrix = Eigen::Matrix<Interval, 2, 2>;
using Vector = Eigen::Matrix<Interval, 2, 1>;

#include <iostream>

int main() {
    Matrix m;
    m << 
       Interval {0.0, 1.0},
       Interval {0.0, 0.0},
       Interval {0.0, 0.0},
       Interval {0.0, 1.0};

    Vector v;
    v << 
       Interval {0.0, 1.0},
       Interval {1.0, 2.0};

    Vector prod = (m*v).eval();

    std::cout << prod(0,0) << std::endl;
    std::cout << prod(1,0) << std::endl;
}

打印:

[0,1]
[0,2]

当然你可以使用Dynamic extents。