rsparse-matrixrcpprcpparmadillomlpack

superlu causing problems with mlpack in Rcpp


I have the following Rcpp code

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

I am getting the following output

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

That is, the code stops working at the spsolve command but executes eigs_sym.

However, removing mlpack by deleting the following solves the problem.

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

So I believe I have linked superlu correctly.

This is only a reproducible part of the code and I want to use mlpack in my project. How can I solve this issue? I also need functions from the boost library for the project and kept that header as well.


Solution

  • When you work with SuperLU you are leaving the (very convenient) world of header-only libraries such as the ones you use here:

    But sparse linear algebra needs linking. This is explained and demonstrated over by an Rcpp Gallery article from two of the RcppArmadillo authors.

    In general it is often advisable to solve smaller subset problems first. You have not real mlpack issue here: you have an issue using RcppArmadillo with sparse decomposition and for that you to link with -lsuperlu as the Rcpp Gallery describes and shows. I would try to get that to work, and then try to proceed to mlpack issues.

    PS 1 Your code has another issue preventing it from running: the behavior-altering #define must come before the #include so please move

    #define ARMA_USE_ARPACK 1
    #define ARMA_USE_SUPERLU 1
    

    to the very top. It then works as the run-time has a normal coding issue:

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

    Note that in order to get this far I had to enable linking with SuperLU via the envoirnment variable, and I validated against the Rcpp Gallery example code:

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

    It is important to set this before the first time you call Rcpp::sourceCpp() as some of these things appear to be cached per session.

    PS 2 Here is a reproducible success. Note that I load the Matrix package (for sparse matrix support) and set the seed to be reproducible.

    > 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]
    
    >