时间戳

首页 > 解决方案 > 如何在 Fastor 或 Xtensor 中编写快速的 c++ 惰性评估代码？

c++ - 如何在 Fastor 或 Xtensor 中编写快速的 c++ 惰性评估代码？

问题描述

我是 C++ 的新手，听说eigen、blaze、Fastor和Xtensor 等具有惰性求值和 simd 的库对于矢量化操作来说很快。

我通过以下函数测量了一些进行基本数值运算的时间：

（快速）

using namespace Fastor;

template<typename T, size_t num>
T func2(Tensor<T,num> &u) {

    Tensor<T,num> z;
    for (auto k=0; k<100; ++k){
        z = u * u;
        z /= exp(u+u);
        z *= 1.;
        z *= sin(u) * cos(z);
    }
    return z(last);
}

(Xtensor)

template<typename T, size_t num>
T func2(xt::xtensor_fixed<T, xt::xshape<num>> &u) {

    xt::xtensor_fixed<T, xt::xshape<num>> z;

    for (auto k=0; k<100; ++k){
        z = u * u;
        z /= xt::exp(u+u);
        z *= 1.;
        z *= xt::sin(u) * xt::cos(z);
    }
    return z(0);
}

编译标志：

（快速）

-std=c++14 -O3 -march=native -funroll-loops -DNDEBUG -mllvm -inline-threshold=10000000 -ffp-contract=fast  -mfma -I/Path/to/Fastor -DFASTOR_NO_ALIAS -DFASTOR_DISPATCH_DIV_TO_MUL_EXPR

(Xtensor)

 -std=c++14 -O3 -march=native -funroll-loops -DNDEBUG -mllvm -inline-threshold=10000000 -ffp-contract=fast  -mfma -I/Path/to/xsimd/include/ -I/Path/to/xtl/include/ -I/Path/to/xtensor/include/ -I/Path/to/xtensor-blas/include/ -DXTENSOR_USE_XSIMD -lblas -llapack -DHAVE_CBLAS=1

编译器：Apple LLVM version 10.0.0 (clang-1000.11.45.5)

处理器：2.6 GHz Intel Core i5

为了比较，我还测量了用 python 编写的函数，它优化了numba.vectorize

@numba.vectorize(['float64(float64)'],nopython=True)
def func(x):
    for k in range(100):
        z = x * x
        z /= np.exp(x + x)
        z *= 1.0
        z *= np.sin(x) * np.cos(x)
    return z

结果（以 usec 为单位）表明

---------------------------------------
num     |  Fastor  |  Xtensor | numba
---------------------------------------
100     |  286     |  201     | 13
1000    |  2789    |  1202    | 65
10000   |  29288   |  20468   | 658
100000  |  328033  |  165263  | 3166
---------------------------------------

我做错什么了吗？Fastor 和 Xtensor 怎么会慢 50 倍。

如何通过使用auto关键字来使用表达式模板和惰性求值？

谢谢你的帮助！

@Jérôme Richard 感谢您的帮助！

有趣的是，Fastor 和 Xtensor 无法忽略冗余的 for 循环。无论如何，我对每个数字运算进行了更公平的比较。

SIMD 的因子 2 也很有意义。

（快速）

template<typename T, size_t num>
T func_exp(Tensor<T,num> &u) {
    Tensor<T,num> z=u;
    for (auto k=0; k<100; ++k){
        z += exp( u );
    }
    return z(0);
}
template<typename T, size_t num>
T func_sin(Tensor<T,num> &u) {
    Tensor<T,num> z=u;
    for (auto k=0; k<100; ++k){
        z += sin( u );
    }
    return z(0);
}
template<typename T, size_t num>
T func_cos(Tensor<T,num> &u) {
    Tensor<T,num> z=u;
    for (auto k=0; k<100; ++k){
        z += cos( u );
    }
    return z(0);
}
template<typename T, size_t num>
T func_add(Tensor<T,num> &u) {
    Tensor<T,num> z=u;
    for (auto k=0; k<100; ++k){
        z += u;
    }
    return z(0);
}
template<typename T, size_t num>
T func_mul(Tensor<T,num> &u) {
    Tensor<T,num> z=u;
    for (auto k=0; k<100; ++k){
        z *= u;
    }
    return z(0);
}
template<typename T, size_t num>
T func_div(Tensor<T,num> &u) {
    Tensor<T,num> z=u;
    for (auto k=0; k<100; ++k){
        z /= u;
    }
    return z(0);
}

(Xtensor)

template<typename T, size_t nn>
T func_exp(xt::xtensor_fixed<T, xt::xshape<nn>> &u) {
    xt::xtensor_fixed<T, xt::xshape<nn>> z=u;
    for (auto k=0; k<100; ++k){
        z += xt::exp( u );
    }
    return z(0);
}
template<typename T, size_t nn>
T func_sin(xt::xtensor_fixed<T, xt::xshape<nn>> &u) {
    xt::xtensor_fixed<T, xt::xshape<nn>> z=u;
    for (auto k=0; k<100; ++k){
        z += xt::sin( u );
    }
    return z(0);
}
template<typename T, size_t nn>
T func_cos(xt::xtensor_fixed<T, xt::xshape<nn>> &u) {
    xt::xtensor_fixed<T, xt::xshape<nn>> z=u;
    for (auto k=0; k<100; ++k){
        z += xt::sin( u );
    }
    return z(0);
}
template<typename T, size_t nn>
T func_add(xt::xtensor_fixed<T, xt::xshape<nn>> &u) {
    xt::xtensor_fixed<T, xt::xshape<nn>> z=u;
    for (auto k=0; k<100; ++k){
        z += u;
    }
    return z(0);
}
template<typename T, size_t nn>
T func_mul(xt::xtensor_fixed<T, xt::xshape<nn>> &u) {
    xt::xtensor_fixed<T, xt::xshape<nn>> z=u;
    for (auto k=0; k<100; ++k){
        z *= u;
    }
    return z(0);
}
template<typename T, size_t nn>
T func_div(xt::xtensor_fixed<T, xt::xshape<nn>> &u) {
    xt::xtensor_fixed<T, xt::xshape<nn>> z=u;
    for (auto k=0; k<100; ++k){
        z /= u;
    }
    return z(0);
}

（努巴）

@numba.vectorize(['float64(float64)'],nopython=True)
def func_exp(u):
    z = u
    for k in range(100):
        z += exp(u)
    return z
@numba.vectorize(['float64(float64)'],nopython=True)
def func_sin(u):
    z = u
    for k in range(100):
        z += sin(u)
    return z
@numba.vectorize(['float64(float64)'],nopython=True)
def func_cos(u):
    z = u
    for k in range(100):
        z += cos(u)
    return z
@numba.vectorize(['float64(float64)'],nopython=True)
def func_add(u):
    z = u
    for k in range(100):
        z += u
    return z
@numba.vectorize(['float64(float64)'],nopython=True)
def func_mul(u):
    z = u
    for k in range(100):
        z *= u
    return z
@numba.vectorize(['float64(float64)'],nopython=True)
def func_div(u):
    z = u
    for k in range(100):
        z *= u
    return z

结果显示

----------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
unit [1E-6 sec] |          exp              |         sin               |           cos             |         add           |           mul         |          div          |
----------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
                |   F     |   X     |   N   |   F     |   X     |   N   |   F     |   X     |   N   |   F   |   X   |   N   |   F   |   X   |   N   |   F   |   X   |   N   |
----------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
n=100           | 135/135 | 38/38   | 10    | 162/162 | 65/32   | 9     | 111/110 | 34/58   | 9     | 0.07  | 0.06  | 6.2   | 0.06  | 0.05  | 9.6   | 0.06  | 0.05  | 9.6   |
n=1000          | 850/858 | 501/399 | 110   | 1004/961| 522/491 | 94    | 917/1021| 486/450 | 92    | 20    | 43    | 57    | 22    | 40    | 91    | 279   | 275   | 91    |
n=10000         | 8113    | 4160    | 830   | 10670   | 4052    | 888   | 10094   | 3436    | 1063  | 411   | 890   | 645   | 396   | 922   | 1011  | 2493  | 2735  | 914   |
n=100000        | 84032   | 46173   | 8743  | 104808  | 48203   | 8745  | 102868  | 53948   | 8958  | 6138  | 18803 | 5672  | 6039  | 13851 | 9204  | 23404 | 33485 | 9149  |
----------------------------------------------------------------------------------------------------------------------------------------------------------------------------|

像这样的格式135/135表示结果。without/with-ffast-math

事实证明

  1. Fastor/Xtensor 在 , , 中的表现非常糟糕exp，sin这cos令人惊讶。
  2. +=Fastor/Xtensor 在, *=,中的缩放比 Numba 差/=。

这是 Fastor/Xtensor 的本质吗？

我将表达式修改为

template<typename T, size_t num>
auto func_exp2(Tensor<T,num> &u) {
    Tensor<T,num> z=u + 100. * exp(u);;
    return z;
}

template<typename T, size_t nn>
auto func_exp2(xt::xtensor_fixed<T, xt::xshape<nn>> &u) {
    xt::xtensor_fixed<T, xt::xshape<nn>> z=u + 100.*xt::exp(u);
    return z;
}

@numba.vectorize(['float64(float64)'],nopython=True)
def func_exp2(u):
    z = u + 100 * exp(u)
    return z

它给了

-----------------------------------------------------------------
unit [1E-6 sec] |     Fastor    |     Xtensor   |      Numba    |
-----------------------------------------------------------------
n=100           |     0.100     |     0.066     |       1.8     |
n=1000          |     0.073     |     0.057     |       3.6     |
n=10000         |     0.086     |     0.089     |      26.7     |
n=100000        |     0.056     |     0.065     |     275.7     |
-----------------------------------------------------------------

发生了什么？

  1. 为什么 Fastor/Xtensor 无法100*exp(u)通过惰性评估将 for 循环表达为天真的？
  2. 为什么随着张量大小的增加 Fastor/Xtensor 变得更快？

标签： c++optimizationlazy-evaluationexpression-templatesxtensor

解决方案

Numpy 实现更快的原因是它计算的东西与其他两个不同。

事实上，python 版本并没有z在表达式中读取np.sin(x) * np.cos(x)。因此，Numba JIT 足够聪明，只需执行一次循环即可证明 Fastor 和 Numba 之间的系数为 100。您可以通过替换range(100)并range(10000000000)观察相同的时间来检查。

最后，在这个基准测试中，XTensor 比 Fastor 更快，因为它似乎使用了自己的 exp/sin/cos的快速SIMD 实现，而 Fastor 似乎使用了 libm 的标量实现，证明 XTensor 和 Fastor 之间的因子为 2 是合理的。

回复更新：

Fastor/Xtensor 在 exp、sin、cos 方面的表现非常糟糕，这令人惊讶。

不，我们不能从基准中得出结论。您所比较的是编译器优化您的代码的能力。在这种情况下，Numba 比普通的 C++ 编译器要好，因为它处理高级 SIMD 感知代码，而 C++ 编译器必须处理来自 Fastor/Xtensor 库的大量基于模板的低级代码。从理论上讲，我认为 C++ 编译器应该可以应用与 Numba 相同的高级优化，但它更难。此外，请注意 Numpy 倾向于创建/分配临时数组，而 Fastor/Xtensor 不应该。

在实践中，Numba 更快，因为u是一个常数exp(u)，sin(u)和 也是cos(u)。因此，Numba 预先计算表达式（仅计算一次）并仍然在循环中执行求和。以下代码给出了相同的时间：

@numba.vectorize(['float64(float64)'],nopython=True)
def func_exp(u):
    z = u
    tmp = exp(u)
    for k in range(100):
        z += tmp
    return z

我猜 C++ 实现由于惰性评估而不会执行此优化。在两个 github 项目上报告这个优化问题可能是个好主意。

此外，请注意u + u + ... + u不严格等于，100 * u因为浮点加法不是 associative。虽然-ffast-math有助于克服这个问题，但编译器仍可能由于优化传递冲突而无法执行这种优化。例如，太多的迭代会阻止循环展开，从而阻止表达式的分解。

我强烈建议您执行更现实的基准测试。

Fastor/Xtensor 在 +=、*=、/= 中的扩展性比 Numba 差。

在这种情况下，Numba 可以用乘法代替常数除法（即1/u可以预先计算）。除此之外，请注意 Fastor 和 Numba 彼此相对较近。

为什么 Fastor/Xtensor 无法通过惰性求值将 for 循环表达为天真的 100*exp(u)？

我认为惰性求值并不意味着表达式会自动分解/优化。这意味着只在需要时才计算结果。然而，表达式分解可能是一个很好的特性，可以在未来的 Fastor/Xtensor 版本中添加（显然还没有）。

为什么随着张量大小的增加 Fastor/Xtensor 变得更快？

我认为它们同样快，而不是更快（时间变化可能是噪音）。因此，我猜这些表达式实际上并没有计算出来。这可能是由于从未读取过的惰性评估。z尝试使用return z(0);而不是return z;（前者强制评估表达式）。

推荐阅读