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

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

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)

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)
    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
! end of file fortran_version.f
  • 4
    $\begingroup$ One naturally wonders what part of the difference is associated with the different wrappers. Did you compile and call them natively and compare? $\endgroup$ Nov 18 '21 at 5:13
  • 2
    $\begingroup$ @Spencer. I agree this is a natural thing to look into. My current tests are in Python, but I will implement tests in C++ and Fortran and update the question with the results. $\endgroup$ Nov 18 '21 at 5:26
  • 2
    $\begingroup$ Maybe you can report start and end times from the c++/Fortran code. These should be one-liners and you can stil conveniently use the wrappers to call everything. $\endgroup$
    – MPIchael
    Nov 18 '21 at 7:47
  • 2
    $\begingroup$ One thing I've noticed: you are using LDLT in the eigen version, while dposv performs an LLT decomposition, which is faster (though probably not 4x faster). To make things more comparable, substitute G = H.ldlt().solve(G); for LLT<Ref<MatrixXd>> H_LLT(H); H_LLT.solveInPlace(G); This will avoid a copy of H and G. $\endgroup$
    – Charlie S
    Nov 18 '21 at 12:16
  • 2
    $\begingroup$ @SpencerBryngelson the noalias() renders the matrix-matrix product closer to that of ?gemm. Likewise with in-place factorization vs copying the entire matrix and then factorizing it. The fortran-style BLAS and LAPACK calls are intended to minimize allocations. Eigen will not avoid this by default, which may widen the performance gap considering this occurs each iteration. $\endgroup$
    – Charlie S
    Nov 18 '21 at 17:59

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.



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.