6
$\begingroup$

I want to use the C++ Eigen library in Mathematica's LibraryLink function.

In the past, matheorem called C++ Eigen in Mathematica successfully. However, I am not familiar with the call between different languages, especially the call to a language function package. I can not follow matheorem's method, because I don't even know where to put the Eigen package so that Mathematica can find it.

My goal is inputting two sparse matrices in Mathematica, and call C++ Eigen by LibraryLink to calculate their product.

The reason I want to call C++ Eigen is that the sparse matrices multiplication in Mathematica is very slow. I used LTemplate before. It provides a convenient way to call C++ Armadillo. After my verification, the time Armadillo sparse matrices multiplication cost is five sixth of Mathematica. Some people say that Eigen is faster than Armadillo in matrix multiplication, so I want to try.

This problem has been bothering me for a long time, and I would like to thank everyone for their suggestions.

$\endgroup$
9
  • $\begingroup$ You surely know that MA can deal with sparse matrices? What do you expect from calling an external library? $\endgroup$
    – yarchik
    Jul 20 '21 at 7:57
  • $\begingroup$ I updated the question. $\endgroup$
    – sidy
    Jul 20 '21 at 8:19
  • 2
    $\begingroup$ This question is entirely too broad, and should be closed if it can't be narrowed some more. As it reads you are asking someone to write an entire LibraryLink interface for you. You say you have used LTemplate before, so you know the basics of writing/compiling/linking C++ code. So how far have you gotten in connecting to eigen? $\endgroup$
    – Jason B.
    Jul 20 '21 at 11:26
  • 1
    $\begingroup$ There is just no point in using Eigen from Mathematica. Either Eigen is linked to MKL (what Mathematica already is) or it is very, very slow. In total, I really do not get why people are so fanatic about Eigen. $\endgroup$ Jul 20 '21 at 11:28
  • 1
    $\begingroup$ Btw., what makes Amardillo fast is again a library like MKL or OpenBLAS. Very likely, the speed differences is caused by a superfluous copy made by LibraryLink when returning the result. This is annoying but not so easy to prevent without ruining the user experience. $\endgroup$ Jul 20 '21 at 11:34
13
$\begingroup$

Even if Eigen is probably not the solution for you problem at hand, it contains a couple of other interesting features that one might want to use from within Mathematica.

So here a minimal example with LibraryLink that compiles and load a function that computes the positive-definite part of a symmetric matrix A with respect to given positive-definite matrix B.

First you have to download Eigen and place it somewhere on your harddrive. You can do it with a package manager (which might make it automatically visible to you C++ compiler). But you can also simply store anywhere you want it to be and write a string with the path to the eigen3 directory into the variable $EigenIncludeDirectory; that's what I assume in the example below:

srcpath = "~";
outpath = "~";

Needs["CCompilerDriver`"];
Module[{opts, path, file, lib},
 If[! FileExistsQ[srcpath], CreateDirectory[srcpath]];
 If[! FileExistsQ[outpath], CreateDirectory[outpath]];
 file = Export[FileNameJoin[{srcpath, "cClipGeneralizedEigenvalues.cpp"}],
   "
#include\"WolframLibrary.h\"
#include<Eigen/Eigenvalues>

EXTERN_C DLLEXPORT int cClipGeneralizedEigenvalues(WolframLibraryData libData, mint Argc, MArgument *Args, MArgument Res)
{
    MTensor Ain = MArgument_getMTensor(Args[0]);
    MTensor Bin = MArgument_getMTensor(Args[1]);
    MTensor Aout;
    mint n = libData->MTensor_getDimensions(Ain)[0];

    int err = libData->MTensor_new(MType_Real, 2, libData->MTensor_getDimensions(Ain), &Aout);
    
    Eigen::GeneralizedSelfAdjointEigenSolver<Eigen::MatrixXd>eigs(n);
    Eigen::Map<Eigen::MatrixXd>A(libData->MTensor_getRealData(Ain),n,n);
    Eigen::Map<Eigen::MatrixXd>B(libData->MTensor_getRealData(Bin),n,n);
    Eigen::Map<Eigen::MatrixXd>C(libData->MTensor_getRealData(Aout),n,n);

    eigs.compute(A,B);

    C=B*eigs.eigenvectors()*eigs.eigenvalues().cwiseMax(0).asDiagonal()*eigs.eigenvectors().inverse();

    MArgument_setMTensor(Res,Aout);
    return LIBRARY_NO_ERROR;
}"
   ,
   "Text"
   ];
 lib = CreateLibrary[{file}, "cClipGeneralizedEigenvalues",
   "TargetDirectory" -> outpath,
   (*"ShellCommandFunction"\[Rule]Print,*)
   "ShellOutputFunction" -> Print,
   "IncludeDirectories" -> {$EigenIncludeDirectory}
   ];
 
 With[{libfile = lib},
  
  cClipGeneralizedEigenvalues::usage = "";
  
  cClipGeneralizedEigenvalues := cClipGeneralizedEigenvalues = 
    LibraryFunctionLoad[libfile, "cClipGeneralizedEigenvalues", {{Real, 2}, {Real, 2}}, {Real, 2}];
  ]
 ]

Here is a simple usage example with timing comparison:

n = 1200;
A = Transpose[#].# &@RandomReal[{-1, 1}, {n, n}];
B = Transpose[#].# &@RandomReal[{-1, 1}, {n, n}];
resultEigen = cClipGeneralizedEigenvalues[A, B]; // AbsoluteTiming // First
resultEigen = cClipGeneralizedEigenvalues[A, B]; // AbsoluteTiming // First

First@AbsoluteTiming[
  {\[CapitalLambda], U} = Eigensystem[{A, B}];
  result = LinearSolve[U, Ramp[\[CapitalLambda]] (U . B)];
  ]

Max[Abs[resultEigen - result]]/Max[Abs[resultEigen]]

2.8983

2.75963

0.242256

1.60978*10^-9

When you run this the first time after compilation, you see that the first timing (2.8983) is a bit greater than the second (2.75963). This is due to the call to LibraryFunctionLoad. You also see that Eigen is very sluggish when not linked to any good linear algebra libraries. Eigen and Armadillo merely provide interfaces to dedicated linear algebra libraries -- something that Mathematica does, too. So when working with the latter, using the former should have little value in general.

$\endgroup$
2
  • $\begingroup$ Through above method, I successfully called Eigen. Thank you very much. $\endgroup$
    – sidy
    Jul 20 '21 at 15:34
  • $\begingroup$ You're welcome! $\endgroup$ Jul 20 '21 at 15:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.