如何在Fastor或Xtensor中编写快速的c++惰性计算代码？

rslzwgfq  于 2023-01-06  发布在  其他
关注(0)|答案(2)|浏览(110)

我是一个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英特尔酷睿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

结果（以微秒为单位）表明

---------------------------------------
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的格式表示-ffast-math的结果without/with。
原来

  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无法通过惰性求值将for循环表示为原始100*exp(u)？
1.为什么Fastor/Xtensor会随着Tensor大小的增加而变快？

c++

2条答案

按热度按时间
inkz8wg9

inkz8wg91#

Numpy实现要快得多的原因是它计算的东西与其他两个不同。
实际上，python版本没有在表达式np.sin(x) * np.cos(x)中读取z，因此，Numba JIT足够聪明，只执行一次循环证明Fastor和Numba之间的系数为100，您可以通过将range(100)替换为range(10000000000)并观察相同的计时来检查。
最后，在此基准测试中，XTensor比Fastor更快，因为XTensor看起来像use its own fast SIMD implementation of exp/sin/cos，而Fastor看起来像use a scalar implementation from 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（如the floating-point addition is not associative）。虽然-ffast-math有助于克服此问题，但编译器仍可能由于优化传递冲突而无法执行此类优化。例如，过多的迭代可能会阻止循环展开，从而阻止表达式的因式分解。
我强烈建议您执行更现实的基准测试。
Fastor/Xtensor在+=、=、/=方面的缩放比Numba差。
在这种情况下，Numba可以用乘法代替常数除法（即1/u可能是预先计算的）。除此之外，请注意Fastor和Numba彼此相对接近。
为什么Fastor/Xtensor不能通过延迟求值将for循环表示为一个简单的100exp（u）？
我认为lazy-evaluation并不意味着表达式会被自动分解/优化，而是意味着结果应该只在需要的时候计算。然而，表达式分解可能是未来Fastor/Xtensor版本中添加的一个很好的特性（显然还没有）。
为什么Fastor/Xtensor会随着Tensor大小的增加而变快？
我认为它们一样快，而不是更快（时间变化可能是噪声）。因此，我猜表达式实际上没有计算。这可能是由于lazy-evaluation，因为z从未被读取。尝试使用return z(0);而不是return z;（前者强制表达式被求值）。

赞(0）分享  回复(0）举报  2023-01-06
ymdaylpp

ymdaylpp2#

我想你误解了惰性求值的工作原理。C++是一种强类型语言，而Python不是。当你执行一个生成表达式模板的操作时，它会生成一个新类型。
此代码：

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);
}

不会产生您所认为的结果。z = u * u会产生一个表示u * u的表达式模板，立即调用它并将其赋值给z，因为z的类型为Tensor〈T，num〉。要使表达式模板继续循环，z的类型必须随着每次迭代而改变!这在Python中是可能的，因为Python是一种动态类型语言。Fastor和Xtensor假设您尝试在每一步计算表达式，从而破坏了它们执行约简的任何机会（许多库例如Blaze、Eigen、Fastor、xtensor、Fastor的文档甚至指出，它将尝试在可能的情况下使用einsum符号和约简来自动约简表达式。它的实现有一个相对复杂的成本模型来说明库实际上有多简单）。
要在C++中实现这一点，需要展开循环，并且在准备求值之前不要给z赋值，可以使用std::make_index_sequence：

template<std::size_t ... Is, typename T>
constexpr auto unroll(std::index_sequence<Is...>, T&& expr) noexcept -> decltype(auto)
{
    return ((Is, expr), ...);
}

template<std::size_t ... Is, typename T>
constexpr auto unroll_add(std::index_sequence<Is...>, T&& expr) noexcept -> decltype(auto)
{
    return ((Is, expr) + ...);
}

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

    Tensor<T,num> z = unroll(std::make_index_sequence<100>{},
        u * u
        / exp(u+u)
        * 1.
        * sin(u) * cos(z)
    );
    return z(last);
}

template<typename T, size_t num>
T func_exp(Tensor<T,num> &u) {
    Tensor<T,num> z = u + unroll_add(std::make_index_sequence<100>{},
        exp( u );
    );
    return z(0);
}

使用表达式模板，可以执行如下多步操作，而无需立即调用表达式：

auto&& a = u + u; // add<T&, T&>
auto&& b = a * u; // mul<add<T&, T&>, T&>
Tensor c = b;.  // Tensor is evaluated here
赞(0）分享  回复(0）举报  2023-01-06
我来回答

相关问题

    查看更多

    热门标签

    更多
    JavaquerypythonNode开发语言requestUtil数据库Table后端算法LoggerMessageElementParser

    最新问答

    更多
    • 蜀ICP备13028337号-1 大数据知识库 https://www.saoniuhuo.com © All rights reserved
    • 本站内容来源互联网,如果侵犯您的权益请联系我们删除, 联系方式：448109455@qq.com