Optimizing a function in Eigen

Question

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.

Thanks in advance for your help!

#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.

Answer 1

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());
  ...