如何在 Fastor 或 Xtensor 中编写快速的 C++ 惰性评估代码?

How to write fast c++ lazy evaluation code in Fastor or Xtensor?

我是 c++ 的新手,听说像 eigenblazeFastor[=102 这样的库=] 和 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 英特尔酷睿 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 在 expsincos 中的表现非常糟糕,这令人惊讶。
  2. Fastor/Xtensor 在 +=*=/=.
  3. 方面比 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 无法通过惰性求值将 for 循环表达为天真的 100*exp(u)
  2. 为什么 Fastor/Xtensor 随着张量大小的增加而变得更快?

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

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

最后,XTensor 在此基准测试中比 Fastor 更快,因为它似乎 use its own fast SIMD implementation of exp/sin/cos while Fastor seems to use a scalar implementation from libm 证明 XTensor 和 Fastor 之间的因子 2 是合理的。


回复更新:

Fastor/Xtensor performs really bad in exp, sin, cos, which was surprising.

没有。我们不能从基准中得出结论。 你比较的是编译器优化代码的能力。在这种情况下,Numba 优于普通 C++ 编译器,因为它处理 high-level SIMD-aware 代码,而 C++ 编译器必须处理巨大的 low-level template-based 代码 来自 Fastor/Xtensor 库。 从理论上讲,我认为 C++ 编译器应该可以应用与 Numba 相同的 high-level 优化,但它更难。此外,请注意 Numpy 倾向于 create/allocate 临时数组,而 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,因为 the floating-point addition is not associative。虽然 -ffast-math 有助于克服此问题,但编译器可能仍会由于优化传递冲突而无法执行此类优化。例如,太多的迭代会阻止循环展开,从而阻止表达式的因式分解。

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

Fastor/Xtensor scales worse than Numba in +=, *=, /=.

在这种情况下,

Numba 可以用乘法代替常量除法(即 1/u 可能是 预计算的 )。除此之外,请注意 Fastor 和 Numba 彼此相对接近。

why Fastor/Xtensor is not able to express the for-loop to a naive 100*exp(u) by lazy-evaluation?

我认为lazy-evaluation并不意味着表达式会自动因式分解/优化。相反,它意味着只应在需要时才计算结果。然而,表达式分解可能是一个很好的特性,可以在未来的 Fastor/Xtensor 版本中添加(显然还没有)。

why Fastor/Xtensor gets faster as tensor size increases?

我认为它们一样快,而不是更快(时间变化可能是噪音)。因此,我猜这些表达式实际上并没有计算出来。这可能是由于 lazy-evaluation 因为 z 从未被读取过。尝试使用 return z(0); 而不是 return z; (前者强制计算表达式)。