计算科学 - Eigen C++ 库比 Fortran 慢 4 倍以上 - 吾爱随笔录

我正在尝试从 Python 调用编译程序的方法。我的主要兴趣是迭代优化方法，所以我正在测试牛顿方法的实现来解决逻辑回归问题。

在我的测试中，使用 f2py 的 Fortran 实现比使用 C++、Eigen 和 Pybind11 的实现快 4 倍以上。我确保包含常见的编译器标志，如 -DNDEBUG、-mtune=My_Architecture 等。我没有包含 ffast-math，因为我的理解是它不安全。

实现之间的这种性能差距比我预期的要大。这仅仅是因为 Fortran 在科学计算方面的优势，还是我在 C++ 实现中遗漏了一些东西？

更新1：

我尝试了评论中的一些建议。这个问题仍然悬而未决，所以我感谢其他建议。

根据 Spencer 的建议，我在 C++ 代码中实现了第二个计时器，以确保问题与包装无关。它不是。C++ 中的时间比 Python 调用报告的时间少千分之一秒。
根据查理的建议，我对代码做了一些小的修改。代码没有明显更快。
我仔细检查了两个实现是否为相同的输入数据执行相同数量的牛顿迭代。他们是这样。
我强迫 eigen 调用 lapack/blas 而不是使用它自己的线性代数实现。令人惊讶的是，该代码仍然比 Fortran 慢得多。

更新 2：

有罪线在 Hessian 计算中。线

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

每次运行需要半秒以上。Charlie S 从更改H为的评论中提出的建议H.no_alias()并没有改善运行时间，但他的直觉是正确的。Fortran 中的可比行是

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

其中spread函数替换了 eigen 对的使用asDiagonal()。我尝试了以下两种计算 H 的替代方法，但都没有提高速度。

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

目前尚不清楚如何使用 eigen 的语法来解决这个瓶颈。

原始代码和运行时间。

以下是 10 次运行的运行时间，设计矩阵大小为 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

C++ 实现：

#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);
}

和 Fortran 实现：

! 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