如何在Fastor或Xtensor中编写快速的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ômeRichard）感谢您的帮助！

有趣的是，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> 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> 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> 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> 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> 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>> 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,size_t nn>
T func_add(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>> 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>> z=u;
    for (auto k=0; k<100; ++k){
        z /= u;
    }
    return z(0);
}

（Numba）

@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> z=u + 100. * exp(u);;
    return z;
}

template<typename T,size_t nn>
auto func_exp2(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无法通过惰性求值向朴素的100*exp(u)表达for循环？
为什么随着张量大小的增加，Fastor / Xtensor变得更快？

解决方法

Numpy实现之所以快得多，是因为它与其他两个实现计算的东西不同。

实际上，python版本不读取表达式z中的np.sin(x) * np.cos(x)。结果，Numba JIT足够聪明，可以仅执行一次循环，从而证明Fastor和Numba之间的系数为100。您可以通过将range(100)替换为range(10000000000)并遵守相同的时间来进行检查。

最后，在本基准测试中，XTensor的速度比Fastor快，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 ++编译器要好，因为它处理的是可识别SIMM的高级代码，而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与the floating-point addition is not associative并不严格等于100 * u。尽管-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;（前者强制对表达式求值）。