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

Question

我正在尝试将我之前用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++中，我调用的函数如下所示，并使用计时器计时：

//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

提前感谢！

Answer 1

C++实现的一个问题可能是它创建了一个甚至可能是两个可以避免的临时副本。第一个副本来自于没有通过引用传递参数(或完美转发)。如果不查看代码的其余部分，很难判断这是否会对性能产生影响。如果在方法xk()之后保证不使用d，编译器可能会将它移动到方法中，但它更有可能将数据复制到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不会指向相同的内存，那么可以在xt::noalias()中进一步包装r，以避免在分配结果之前进行临时复制。如果不通过引用返回，则函数的返回值也是如此。

祝你好运，祝你编码愉快！