代码之家 › 专栏 › 技术社区 › noirritchandra

superlu导致Rcpp中mlpack出现问题

mlpack rcpparmadillo rcpp sparse-matrix r

1

noirritchandra · 技术社区 · 9 月前

我有以下Rcpp代码

#include <RcppArmadillo.h>
#define _USE_MATH_DEFINES
#include<cmath>
#include <boost/math/special_functions/bessel.hpp> 

#include <mlpack.h>
#include <mlpack/methods/kmeans.hpp>
// [[Rcpp::depends(RcppEnsmallen)]]
// [[Rcpp::depends(mlpack)]]

// [[Rcpp::plugins(cpp17)]]
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::plugins(openmp)]]
// [[Rcpp::depends(BH)]]

using namespace arma;
using namespace Rcpp;

#define ARMA_USE_ARPACK 1
#define ARMA_USE_SUPERLU 1
#define tolr 1e-5
#define crossprod(x) symmatu((x).t() * (x))
#define tcrossprod(x) symmatu((x) * (x).t())
#define pi_sq 9.869604401089358


// [[Rcpp::export]]
Rcpp::List NNGP_fit(const int n){
  sp_mat A = sprandu<sp_mat>(n, n, 0.1);
  sp_mat B = A.t()*A;
  
  vec eigval;
  mat eigvec;
  
  eigs_sym(eigval, eigvec, B, n-1);  // find n-1 eigenvalues/eigenvectors
  Rcpp::Rcout<<"flag eig"<<endl;
  
  superlu_opts opts;
  opts.symmetric  = true;
  opts.equilibrate =true;
  vec mu_f,
   y(n, fill::randn);
  
  bool status = spsolve(mu_f,B, y, "superlu",opts);
  Rcpp::Rcout<<"status= "<<status<<endl;
  
  return Rcpp::List::create(Rcpp::Named("covmat") =wrap(B),
                            Rcpp::Named("eigval") =wrap(eigval),
                            Rcpp::Named("eigvec") =wrap(eigvec),
                            Rcpp::Named("y") =wrap(y),
                            Rcpp::Named("mu_f") =wrap(mu_f)
  );
}

我得到以下输出

flag eig
Error: spsolve(): use of SuperLU must be enabled

也就是说,代码在 spsolve 命令但执行 eigs_sym .

但是,删除 mlpack 通过删除以下内容可以解决问题。

#include <mlpack.h>
#include <mlpack/methods/kmeans.hpp>
// [[Rcpp::depends(RcppEnsmallen)]]
// [[Rcpp::depends(mlpack)]]

所以我相信我已经联系上了 superlu 正确地。

这只是代码的可复制部分,我想使用 mlpack 在我的项目中。我该如何解决这个问题?我还需要来自 boost 项目的库,并保留了标题。

1 回复 | 直到 9 月前

1

0

Ada Lovelace 9 月前

当你使用SuperLU时,你将离开(非常方便的)只使用头文件库的世界,比如你在这里使用的库:

Armadillo通过RcpArmadillo(因为从R使用时不需要链接到LAPACK/BLAS)
ensmallen通过RcppEnsmall,因为它只是标题
通过BH进行增强,因为它专注于仅包含头部的子集。
mlpack,从4.0版本开始,它只与R一起工作。

但是稀疏线性代数需要链接。这可以通过以下方式解释和演示 an Rcpp Gallery article 来自RcppArmadillo的两位作者。

一般来说,通常建议先解决较小的子集问题。你不是真的 mlpack 这里的问题:您在使用时遇到问题 RcppArmadillo 通过稀疏分解,您可以链接到 -lsuperlu 正如Rcpp画廊所描述和展示的那样。我会试着让它发挥作用,然后试着继续 mlpack 问题。

PS 1 你的代码还有另一个问题阻止它运行:行为改变 #define 必须来之前 #include 所以请移动

#define ARMA_USE_ARPACK 1
#define ARMA_USE_SUPERLU 1

到最顶端。然后,当运行时出现正常的编码问题时,它就会工作:

> Rcpp::sourceCpp("soquestion20241227.cpp")
> set.seed(123)
> NNGP_fit(3)
flag eig
status= 0
Error in getClass("dgCMatrix") : âdgCMatrixâ is not a defined class

请注意,为了走到这一步,我必须启用链接通过envoirment变量使用SuperLU,我根据Rcpp Gallery示例代码进行了验证:

> Sys.setenv("PKG_LIBS"="-lsuperlu")
> Rcpp::sourceCpp("gallery20170217.cpp")
> superLuDemo()
Done.
> 
>

设置此项很重要之前你第一次打电话 Rcpp::sourceCpp() 因为其中一些内容似乎是按会话缓存的。

PS 2 这是一个 可重复性 成功。请注意,我加载了 Matrix 打包(用于稀疏矩阵支持)并将种子设置为可重复的。

> library(Matrix)
> Sys.setenv("PKG_LIBS"="-lsuperlu")
> Rcpp::sourceCpp("soquestion20241227.cpp")
> set.seed(123)
> NNGP_fit(5)
flag eig
status= 0
$covmat
5 x 5 sparse Matrix of class "dgCMatrix"
                             
[1,] 0.0827008 . . . 0.226699
[2,] .         . . . .       
[3,] .         . . . .       
[4,] .         . . . .       
[5,] 0.2266988 . . . 0.788687

$eigval
            [,1]
[1,] 7.05517e-33
[2,] 6.13069e-17
[3,] 1.61746e-02
[4,] 8.55213e-01

$eigvec
             [,1]         [,2]         [,3]         [,4]
[1,]  1.97882e-16 -7.23362e-15  9.59537e-01 -2.81582e-01
[2,]  1.64671e-01  7.22721e-01  5.40204e-15  2.43186e-17
[3,] -3.27270e-02 -6.76148e-01 -5.07065e-15 -1.23420e-17
[4,] -9.85805e-01  1.43172e-01  1.05140e-15  1.94107e-18
[5,]  3.29803e-17  2.09006e-15 -2.81582e-01 -9.59537e-01

$y
           [,1]
[1,]  0.0699539
[2,]  0.9767215
[3,]  2.1650415
[4,] -1.8262054
[5,]  0.5804146

$mu_f
     [,1]

>