6

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.

flinty
  • 25,147
  • 2
  • 20
  • 86
sidy
  • 141
  • 4
  • You surely know that MA can deal with sparse matrices? What do you expect from calling an external library? – yarchik Jul 20 '21 at 07:57
  • I updated the question. – sidy Jul 20 '21 at 08:19
  • 2
    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? – Jason B. Jul 20 '21 at 11:26
  • 1
    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. – Henrik Schumacher Jul 20 '21 at 11:28
  • "The reason I want to call C++ Eigen is that the sparse matrices multiplication in Mathematica is very slow." Example please. To my knowledge, Mathematica uses Intel MKL for that. And my experience is that Eigen without linking to MKL can be as much as 10 times slower than MKL. – Henrik Schumacher Jul 20 '21 at 11:29
  • 1
    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. – Henrik Schumacher Jul 20 '21 at 11:34
  • Look up mkl_sparse_d_mm in the MKL docs (https://software.intel.com/content/www/us/en/develop/documentation/onemkl-developer-reference-c/top/blas-and-sparse-blas-routines/inspector-executor-sparse-blas-routines/inspector-executor-sparse-blas-execution-routines/mkl-sparse-mm.html); grab the function's prototype from there and just use it in your "LTemplate`" workflow. The MKL should already be linked by every library created by "LTemplate`". The interface is not as fancy as in Armadillo or Eigen, but it gets the job done. – Henrik Schumacher Jul 20 '21 at 11:37
  • Btw.: In general multiplying two sparse matrices is not a good idea anyways, so you should try to avoid that (it's not always possible, though). – Henrik Schumacher Jul 20 '21 at 11:40
  • Welcome to Mathematica.SE! I hope you will become a regular contributor. To get started, 1) take the introductory [tour] now, 2) when you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge, 3) remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign, and 4) give help too, by answering questions in your areas of expertise. – bbgodfrey Jul 20 '21 at 12:39

1 Answers1

13

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-&gt;MTensor_new(MType_Real, 2, libData-&gt;MTensor_getDimensions(Ain), &amp;Aout);

Eigen::GeneralizedSelfAdjointEigenSolver&lt;Eigen::MatrixXd&gt;eigs(n);
Eigen::Map&lt;Eigen::MatrixXd&gt;A(libData-&gt;MTensor_getRealData(Ain),n,n);
Eigen::Map&lt;Eigen::MatrixXd&gt;B(libData-&gt;MTensor_getRealData(Bin),n,n);
Eigen::Map&lt;Eigen::MatrixXd&gt;C(libData-&gt;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.

Henrik Schumacher
  • 106,770
  • 7
  • 179
  • 309