xtensor 的“operator/”比 numpy 的“/”慢

问题描述

我正在尝试将我之前用 python 编写的一些代码转换为 C++，我目前正在测试 xtensor，看看它是否可以比 numpy 更快地完成我需要的工作。

我的一个函数接受一个方阵 d 和一个标量 alpha，并执行元素操作 alpha/(alpha+d)。背景：此函数用于测试 alpha 的哪个值是“最佳”，因此处于循环中，其中 d 始终相同，但 alpha 不同。

以下所有时间尺度均为运行该函数的 100 个实例的平均值。

在numpy中，大约需要0.27秒，代码如下：

def kfun(d,alpha):
    k = alpha /(d+alpha)
    return k

但是xtensor大约需要0.36秒，代码如下：

xt::xtensor<double,2> xk(xt::xtensor<double,2> d,double alpha){
    return alpha/(alpha+d);
}

我还使用 std::vector 尝试了以下版本，但从长远来看，我不想使用这个版本，即使只用了 0.22 秒。

std::vector<std::vector<double>> kloops(std::vector<std::vector<double>> d,double alpha,int d_size){
    for (int i = 0; i<d_size; i++){
        for (int j = 0; j<d_size; j++){
            d[i][j] = alpha/(alpha + d[i][j]);
        }
    }
    return d;
}

我注意到 xtensor 中的 operator/ 使用“延迟广播”，有没有办法让它立即生效？

编辑：

在Python中，函数调用如下，并使用“time”包计时

t0 = time.time()
for i in range(100):
    kk = k(dsquared,alpha_squared)
print(time.time()-t0)

在 C++ 中，我调用的函数如下，并使用 chronos 计时：

//d is saved as a 1D npy file,an artefact from old code
auto sd2 = xt::load_npy<double>("/path/to/d.npy");

shape = {7084,7084};
    xt::xtensor<double,2> xd2(shape);
    for (int i = 0; i<7084;i++){
        for (int j=0; j<7084;j++){
            xd2(i,j) = (sd2(i*7084+j));
        }
    }

auto start = std::chrono::steady_clock::now();
for (int i = 0;i<10;i++){
    matrix<double> kk = kfun(xd2,4000*4000,7084);
}
auto end = std::chrono::steady_clock::now();
std::chrono::duration<double> elapsed_seconds = end-start;
std::cout << "k takes: " << elapsed_seconds.count() << "\n";

如果您想运行此代码，我建议使用 xd2 作为对称的 7084x7084 随机矩阵，对角线上的零点。

函数的输出，一个名为 k 的矩阵，然后继续在其他函数中使用，但我仍然需要 d 保持不变，因为它将被重用。

结束编辑

要运行我的 C++ 代码，我在终端中使用以下行：

cd "/path/to/src/" && g++ -mavx2 -ffast-math -DXTENSOR_USE_XSIMD -O3 ccode.cpp -o ccode -I/path/to/xtensorinclude && "/path/to/src/"ccode

提前致谢！

解决方法

C++ 实现的一个问题可能是它创建了一个或什至可能是两个可以避免的临时副本。第一个副本来自没有通过引用（或完美转发）传递参数。如果不查看其余代码，很难判断这是否对性能有影响。如果保证在方法 d 之后不被使用，编译器可能会将 xk() 移动到方法中，但更有可能将数据复制到 d 中。

引用传递，方法可以改为

xt::xtensor<double,2> xk(const xt::xtensor<double,2>& d,double alpha){
    return alpha/(alpha+d);
}

要使用完美转发（并启用其他 xtensor 容器，如 xt::xarray 或 xt::xtensor_fixed），可以将该方法更改为

template<typename T>
xt::xtensor<double,2> xk(T&& d,double alpha){
    return alpha/(alpha+d);
}

此外，您可以避免为返回值保留内存。同样，如果没有看到其余的代码，很难判断。但是如果在循环内部使用该方法，并且返回值始终具有相同的形状，那么在循环外部创建返回值并通过引用返回可能会有所帮助。为此，可以将方法更改为：

template<typename T,typename U>
void xk(T& r,U&& d,double alpha){
    r = alpha/(alpha+d);
}

如果保证 d 和 r 不指向同一个内存，可以进一步将 r 包裹在 xt::noalias() 中以避免在分配之前临时复制结果。如果不按引用返回，函数的返回值也是如此。

祝你好运，编码愉快！