# Optimizing a function in Eigen

I'm a beginner in C++ and I would appreciate advices to optimize the following function I wrote with Eigen (in fact, to be used with RcppEigen).
So far, I observe a 3.5x speed-up compared to the corresponding function written in R, and I was wondering whether I could gain more.

Note that I'm working with very large matrices, so I rely on maps to avoid copies from the corresponding objects in R.

#include <RcppEigen.h>

using namespace Rcpp;
using namespace Eigen;

typedef Map<ArrayXd> MapArr1D;
typedef Map<ArrayXXd> MapArr2D;
typedef Map<MatrixXd> MapMat;
typedef Map<VectorXd> MapVec;

// [[Rcpp::depends(RcppEigen)]]

// [[Rcpp::export]]
void myFct(const MapMat M1, const MapMat M2, MapMat M3, MapMat M4, MapArr2D A1,
MapArr2D A2, const MapArr1D a1, const MapArr1D a2, const MapArr1D a3,
const MapArr1D a4, const MapArr1D a5, const double d1) {

for (int j = 0; j < M1.cols(); ++j) {

M4.noalias() -= M1.col(j) * M3.row(j);

A1.row(j) = a1 * a2 * ((M2 - M4).transpose() * M1.col(j)).array();

A2.row(j) = exp(-Fct(a3(j) - a4(j) - a5 / 2 - d1 / 2 -
pow(A1.row(j).transpose(), 2) / (2 * a1) - log(a1) / 2));

M3.row(j) = A1.row(j) * A2.row(j);

M4.noalias() += M1.col(j) * M3.row(j);

}
}


where Fct is some other function.

• Please provide some information about what myFct is calculating. MathJax is available, if you need to write mathematical text. – 200_success Mar 14 '17 at 8:41
• Hi, myFct(x) computes log(1 + exp(x)) using the "log-sum-exp" trick to avoid overflow. – user79097 Mar 14 '17 at 8:50
• Then why does it take 12 parameters? – 200_success Mar 14 '17 at 8:55
• Cross-posted from Stack Overflow. Please declare your cross-posts, to help users triage your question. – 200_success Mar 14 '17 at 8:58
• OK, I will. Sorry about that. – user79097 Mar 14 '17 at 10:23

It's hard to help without knowing the sizes, but:

1. The first thing to do would be to measure the relative cost of each of the 5 statements to see where the bottleneck is.

2. -(a5/2 + d1/2 + log(a1)) could be precomputed into a temporary outside the loop.

3. Replace pow(A1.row(j).transpose(), 2) by A1.row(j).transpose().square().

4. Compilation flags might also help: make sure you enabled AVX and if supported FMA using e.g., -march=native.

5. Depending on the sizes, you might want to rewrite the expressions to benefits from faster matrix-matrix products:

'

T = (M2-M4).transpose() * M1;
for(j...)
...
A1.row(j) = a1*a2 * (T.col(j) - M1.col(j).squaredNorm() * M3.row(j).transpose());
...

• Thank you very much! I've implemented 2 and 3, and the performance is slightly improved. I'll check the rest. – user79097 Mar 14 '17 at 10:20
• I've tried now using -march=nativebut I don't see any improvement. Regarding your point 5., I'm not sure about this because M3 gets updated later in the loop. Thanks for your help! – user79097 Mar 14 '17 at 22:13
• Regarding the sizes of the matrices, M1has a a very large number, n, of columns (n is of the order of 500K) and so M3, A1and A3 are also large, with n rows. Do you see any further improvements I could implement using this information? Thanks! – user79097 Mar 14 '17 at 22:19