Eigen C++ Library Over 4x Slower than Fortran

Question

I am experimenting with ways to call compiled programs from Python. My primary interest is iterative optimization methods, so I'm testing an implementation of Newton's method to solve a logistic regression problem.

In my tests, a Fortran implementation with f2py is more than 4x faster than an implementation using C++, Eigen, and Pybind11. I've made sure to include common compiler flages like -DNDEBUG, -mtune=My_Architecture, etc. I haven't included ffast-math because my understanding is that it is unsafe.

This performance gap between implementations is larger than I expected. Is this simply due to Fortran's strength in scientific computing, or am I missing something in my C++ implementation?

Update 1:

I tried some suggestions from the comments. This question is still open so I appreciate additional suggestions.

• Per Spencer's suggestion, I implemented a second timer within the C++ code to make sure that the issue wasn't related to wrapping. It isn't. The time within C++ is on the order of one thousandth of a second less than the time reported from the Python call.

• Per Charlie's suggestion, I made some minor modifications to the code. The code is not discernibly faster.

• I double checked that the two implementations perform the same number of Newton iterations for the same input data. They do.

• I forced eigen to make calls to lapack/blas instead instead of using its own linear algebra implementations. Surprisingly, the code is still much slower than Fortran.

Update 2:

The guilty line is in the Hessian computation. The line

H = X.transpose() * (S*(1.0 - S)/double(y.rows())).matrix().asDiagonal() * X;

takes over half a second each time it is run. Charlie S's suggestion from the comments of changing H to H.no_alias() did not improve the runtime, but his intuition was correct. The comparable line in Fortran is

call dgemm('t', 'n', m, m, n, 1d0, x, n, spread(s*(1d0-s)/float(n), 2, m)*x, n, 0d0, h, m)

where the spread function replaces eigen's use of asDiagonal(). I tried the follow two alternative methods of computing H, but neither improved the speed.

H = (X.transpose().array().rowwise() * (S*(1.0 - S)/double(y.rows())).transpose()).matrix() * X;
H = X.transpose() * (X.array().colwise() * (S*(1.0 - S)/double(y.rows()))).matrix();

It's unclear how I can fix this bottleneck using eigen's syntax.

Original code and run times.

Here are the run times for 10 runs w/ design matrices of size 20000 x 500:

Mean/Std of Python Run times (s): 20.531863737106324 +/- 1.4717737632071584
Mean/Std of Fortran Run times (s): 1.5307265996932984 +/- 0.16236569639705614
Mean/Std of C++ Run times (s): 8.926126623153687 +/- 0.585820111836348

The C++ implementation:

#include <pybind11/pybind11.h>
#include <pybind11/eigen.h>
#include <iostream>
#include <Eigen/Dense>

namespace py = pybind11;
using namespace Eigen;

void logistic(Ref<ArrayXd> Z, Ref<ArrayXd> W){
    W = Z.exp()/(1.0 + Z.exp());
}

void gradf(Ref<MatrixXd> X, Ref<VectorXd> y, Ref<VectorXd> beta, Ref<VectorXd> G, Ref<ArrayXd> S){
    G = X.transpose() * (S.matrix()-y)/double(y.rows());  
}

void hessf(Ref<MatrixXd> X, Ref<VectorXd> y, Ref<VectorXd> beta, Ref<MatrixXd> H, Ref<ArrayXd> S, Ref<ArrayXd> Xbeta){
    Xbeta = X*beta;
    logistic(Xbeta, S);
    H = X.transpose() * (S*(1.0 - S)/double(y.rows())).matrix().asDiagonal() * X;
}

VectorXd Newtons_logreg(Ref<MatrixXd> X, Ref<VectorXd> y, Ref<VectorXd> beta_i, const int max_iter, const double tol){
    int it = 0;
    ArrayXd S(y.rows());
    ArrayXd Xbeta(y.rows());
    VectorXd G(beta_i.rows());
    MatrixXd H(X.cols(), X.cols());
    do{
        hessf(X, y, beta_i, H, S, Xbeta);
        gradf(X, y, beta_i, G, S);
        G = H.ldlt().solve(G);
        beta_i -= G;
        if (G.norm() <= tol){
            return beta_i;
        }
        it += 1;
    }
    while (it <= max_iter);
    std::cout << "Early termination in cpp" << std::endl;
    return beta_i;
}


PYBIND11_MODULE(cpp_funcs, m){
    m.def("Newtons_logreg", &Newtons_logreg);
}

and the Fortran implementation:

! file: fortran_version.f
! compile with f2py -c -m fortran_v fortran_version.f90

subroutine logistic(z, n, w)
    ! compute the logistic function of an array z of length n
    implicit none
    integer n
    !f2py integer intent(hide),depend(z) :: n=shape(z,0)
    real*8, intent(in):: z(n)
    real*8, intent(out):: w(n)
    w = exp(z)/(1d0+exp(z))
end

subroutine gradf(x, y, n, m, g, s)
    ! compute the gradient of the logistic regression likelihood 
    implicit none
    integer n,m
    !f2py intent(hide),depend(x) :: n=shape(x,0), m = shape(x,1)
    real*8, intent(in):: s(n), y(n), x(n,m)
    real*8, intent(out):: g(m)
    call dgemm('n', 'n', 1, m, n, 1d0, (s-y)/float(n), 1, x, n, 0d0, g, 1)
    !g = matmul((s-y)/float(n), x) !i've tried this option and above; unclear which is faster
end

subroutine hessf(x, beta, n, m, h, xbeta, s)
    ! compute the hessian of the logistic regression likelihood f
    implicit none
    integer n,m
    !f2py intent(hide),depend(x) :: n=shape(x,0), m = shape(x,1)
    real*8, intent(in):: x(n,m), beta(m)
    real*8, intent(out):: h(m,m), s(n), xbeta(n)
    call dgemm('n', 'n', n, 1, m, 1d0, x, n, beta, m, 0d0, xbeta, n)
    call logistic(xbeta, n, s)
    call dgemm('t', 'n', m, m, n, 1d0, x, n, spread(s*(1d0-s)/float(n), 2, m)*x, n, 0d0, h, m)
end

subroutine newtons_logreg(x, y, beta_i, max_iter, tol, n, m) 
    ! compute the maximum likelihood estimate of logistic regression w/ design
    !   matrix x and response variable y
    implicit none
    integer n,m
    !f2py intent(hide),depend(x) :: n=shape(x,0), m = shape(x,1)
    integer, intent(in):: max_iter
    integer o, it
    real*8, intent(in):: y(n), x(n,m), tol
    real*8, intent(inout):: beta_i(m)
    real*8 h(m,m), g(m)
    real*8 s(n), xbeta(n)
    it=0
    do while ( it .lt. max_iter)
        call hessf(x, beta_i, n, m, h, xbeta, s)
        call gradf(x, y, n, m, g, s)
        call dposv('u', m, 1, h, m, g, m, o)
        beta_i = beta_i - g
        if (norm2(g) .lt. tol) return
        it = it + 1
    enddo
end    
        
! end of file fortran_version.f

Answer 1

Despite my initial conviction that I was including all relevant compiler flags, in fact I wasn't. I forgot the -fopenmp flag on the C++ compiler, which enables multithreading and a fairer comparison between implementations. I also found that -march=native makes the code three times faster than the more general -march=nocona which I was originally using, and also slightly faster than the Fortran version. In conclusion, neither implementation is faster than the other, because it's the compiler options that make all of the difference.

For future reference, the following stack overflow question's top answer has a full description of all the relevant flags and the speedup they provide when compiling C++ written with Eigen.

https://stackoverflow.com/questions/51656818/benchmarking-matrix-multiplication-performance-c-eigen-is-much-slower-than