将 boost 多精度与 RcppEigen 结合使用

Question

我使用 RcppArmadillo 编写了一个库。我遇到的问题是，对于某些参数，arma::solve 函数没有给我确切的解决方案，因为 "A" 矩阵的 rcond 接近 0。如果 arma::solve 可以解决线性问题完全等式，那不会有问题。但它给了我一个近似解，这对我来说还不够好。

然后，我想到了使用RcppEigen，使用boost多精度变量。如果我正确理解 Eigen，Eigen 求解器会给我一个多精度的解决方案，这个解决方案很可能对我来说足够好（即使是 float128）。

但是当我尝试执行这个计划时，我遇到了错误，我不知道该怎么办。例如下面的作品：

#include <RcppArmadillo.h>
#include <RcppEigen.h>
#include <cmath>
#include <cstdint>
#include <boost/multiprecision/float128.hpp>
#include <boost/multiprecision/eigen.hpp>

// Correctly setup the build environment 
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::depends(RcppEigen)]]
// [[Rcpp::depends(BH)]]
using namespace Rcpp;
using namespace arma;
using namespace Eigen;
using namespace boost;
using namespace boost::multiprecision;
using Eigen::Map;
using Eigen::MatrixXd;                  // variable size matrix, double precision
using Eigen::VectorXd;
using Eigen::Matrix;
using Eigen::Dynamic;
namespace mp = boost::multiprecision;


    // [[Rcpp::export]]
    Eigen::MatrixXd onesfgh_LPPLS_RcppEigen(int t2, int t1, double tc, double m, double w) {
      //Definitions & pre-computations:
      unsigned int tdim;
      double ttc, powttc, wlog, sinwlog, coswlog;
      tdim = t2 - t1 + 1;
      ttc = tc - t1 + 1;
      Eigen::MatrixXd output(tdim, 4); //first dimension is time, second is onesfgh
      //Main calculations:
      for (unsigned int i = 0; i < tdim; i++) {
        output(i, 0) = 1.0;
        ttc--;
        powttc = pow(ttc, m);
        wlog = w * log(ttc);
        sinwlog = sin(wlog);
        coswlog = cos(wlog);
        //sincos(wlog, &sinwlog, &coswlog);
        output(i, 1) = powttc;
        output(i, 2) = coswlog * powttc;
        output(i, 3) = sinwlog * powttc;
      }
      return output;
    }

但是：

// [[Rcpp::export]]
Eigen::Matrix<boost::multiprecision::number<boost::multiprecision::float128, boost::multiprecision::et_off>, Eigen::Dynamic, 4> onesfgh_LPPLS_RcppEigen(int t2, int t1, double tc, double m, double w) {
  //Definitions & pre-computations:
  unsigned int tdim;
  boost::multiprecision::number<boost::multiprecision::float128, boost::multiprecision::et_off> tc128(tc), m128(m), w128(w);
  boost::multiprecision::number<boost::multiprecision::float128, boost::multiprecision::et_off> ttc, powttc, wlog, sinwlog, coswlog;
  boost::multiprecision::number<boost::multiprecision::float128, boost::multiprecision::et_off> t1128(double(t1));
  tdim = t2 - t1 + 1;
  ttc = tc128 - t1128 + 1.0;
  Eigen::Matrix<boost::multiprecision::number<boost::multiprecision::float128, boost::multiprecision::et_off>, Eigen::Dynamic, 4> output(tdim, 4); //first dimension is time, second is onesfgh
  //Main calculations:
  for (unsigned int i = 0; i < tdim; i++) {
    output(i, 0) = 1.0;
    ttc--;
    powttc = boost::multiprecision::pow(ttc, m128);
    wlog = w128 * boost::multiprecision::log(ttc);
    sinwlog = boost::multiprecision::sin(wlog);
    coswlog = boost::multiprecision::cos(wlog);
    //sincos(wlog, &sinwlog, &coswlog);
    output(i, 1) = powttc;
    output(i, 2) = coswlog * powttc;
    output(i, 3) = sinwlog * powttc;
  }
  return output;
}

它不起作用。我收到错误消息：

RcppExports.cpp:11:15: error: 'boost' was not declared in this scope

特别是，我尝试了 https://www.boost.org/doc/libs/1_60_0/libs/math/doc/html/math_toolkit/high_precision/using_test.html 中给出的建议，但没有使用 typedef，因为在另一个 post 中，Dirk Eddelbuettel 表示 typedef 不容易处理（它们应该是放在另一个 .h 文件中）。

有什么建议可以继续吗？

编辑：

我修改了函数，不导出。现在代码是：

  #include <RcppArmadillo.h>
    #include <RcppEigen.h>
    #include <cmath>
    #include <cstdint>
    #include <boost/multiprecision/float128.hpp>
    #include <boost/multiprecision/eigen.hpp>

    // Correctly setup the build environment 
    // [[Rcpp::depends(RcppArmadillo)]]
    // [[Rcpp::depends(RcppEigen)]]
    // [[Rcpp::depends(BH)]]
    using namespace Rcpp;
    using namespace arma;
    using namespace Eigen;
    using namespace boost;
    using namespace boost::multiprecision;
    using Eigen::Map;
    using Eigen::MatrixXd;                  
    using Eigen::VectorXd;
    using Eigen::Matrix;
    using Eigen::Dynamic;
    namespace mp = boost::multiprecision;

    Eigen::Matrix<mp::float128, Eigen::Dynamic, 4> onesfgh_LPPLS_RcppEigen(int `t2, int t1, double tc, double m, double w) {`
      //Definitions & pre-computations:
      unsigned int tdim;
      mp::float128 tc128(tc), m128(m), w128(w);
      mp::float128 ttc, powttc, wlog, sinwlog, coswlog;
      mp::float128 t1128(double(t1));
      tdim = t2 - t1 + 1;
      ttc = tc128 - t1128 + 1.0;
      Eigen::Matrix<mp::float128, Eigen::Dynamic, 4> output(tdim, 4); 
      //Main calculations:
      for (unsigned int i = 0; i < tdim; i++) {
        output(i, 0) = 1.0;
        ttc--;
        powttc = mp::pow(ttc, m128);
        wlog = w128 * mp::log(ttc);
        sinwlog = mp::sin(wlog);
        coswlog = mp::cos(wlog);
        //sincos(wlog, &sinwlog, &coswlog);
        output(i, 1) = powttc;
        output(i, 2) = coswlog * powttc;
        output(i, 3) = sinwlog * powttc;
      }
      return output;
    }

现在我收到以下错误消息：

babel_RcppArmadillo.cpp:6:42: fatal error: boost/multiprecision/eigen.hpp: 
No such file or directory
#include <boost/multiprecision/eigen.hpp>                                          
compilation terminated.

这很奇怪，因为我假设 BH 包含 boost 的子文件夹。

Answer 1

更新： 现在应该不需要了，因为自 2019 年 1 月以来 BH package 包含所需的 header。

boost/multiprecision/eigen.hpp header 是在版本 1.68 中添加的，而 BH package 目前提供 boost 1.66。您必须单独安装更新的 boost。以下应该有效但未经测试：

下载 boost 1.68 并将其解压到合适的目录中。在 Linux 和其他 Unix-like 系统上，我可能会使用 /usr/local/include。否则我会使用名称中没有空格的任何路径。
不要依赖 BH 包，即为 sourceCpp 等删除 // [[Rcpp::depends(BH)]]。或 Imports: BH 包使用。
如果在第 1 步中将 boost 安装在 non-standard 位置，则必须告诉编译器到那里查找。这可以用
来完成
```
PKG_CPPFLAGS=-I<path-to-boost>
```
通过 Sys.setenv 为 sourceCpp 等人。或在 src/Makevars(.win) 内用于包使用。

将 boost 多精度与 RcppEigen 结合使用

Using boost multiprecision with RcppEigen

boost

r

armadillo

rcpp

eigen