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

Question

我是 c++ 的新手，听说像 eigen、blaze、Fastor[=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.

原来是

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

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

发生了什么事？

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

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

Answer 1

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; （前者强制计算表达式）。

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

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

c++

optimization

lazy-evaluation

expression-templates

xtensor