Julia 的广播速度是 Matlab 的两倍

Posted 2023-03-23

【中文标题】Julia 的广播速度是 Matlab 的两倍

【英文标题】：Julia broadcasting twice as slow as Matlab

【发布时间】：2021-07-02 07:35:16

【问题描述】：

我正在努力让自己熟悉 Julia，以便从 Matlab 迁移，到目前为止一切顺利，直到我开始使用广播来移植一个特定函数，该函数的执行速度或多或少是 Matlab 的两倍。

function features(X::VectorFloat64,M::Int,hyper::Float64,mid::Float64)
    X = X.-mid
    H = 4.0.*hyper.+maximum(abs.(X))
    X = (X.+H)./(2.0.*H)
    w = transpose(1:M)
    S = (sqrt.(2.0.*pi).*hyper).*exp.((-0.5.*hyper.^2).*((pi.*w./(2.0.*H)).^2))
    f = H.^(-0.5).*sin.(pi.*X.*w).*sqrt.(S)
end

任何帮助将不胜感激！

【问题讨论】：

到底是什么问题？我从未使用过 Julia，但语言之间的速度提高 2 倍是完全正常的

问题是什么？如果您想获得有用的答案，请具体明确！

如果您能提供实际的时间和制作方式，那就太好了。一开始很常见，测量编译时间或基准全局变量性能（这是出了名的慢）

有些广播实际上是不需要的，比如第二行，2.0 .* H，第 5 行的大部分，等等...但我确实尝试过，删除它们只节省 4%。您是否使用@btime features($X,$M,$hyper,$mid) 进行基准测试？否则，您包括编译时间或在基准测试中将参数视为全局变量。如果您在这里没有得到答案，请尝试official discourse forum。像这样的问题很快就会得到解答！

在最后一行中，您通过在列向量 (X) 和行向量 (S) 上广播来创建矩阵。那是你要的吗？如果您能提供示例输入和相应的所需输出，那就太好了。

【参考方案1】：

首先，您对广播的使用不是最佳的。你用的太多了，还不够；）

其次，几乎所有运行时 (99.9%) 都发生在广播的 sin 表达式中，因此应该集中精力。

第三，在这种情况下，您不应该真的期望 Julia 能胜过 Matlab。这正是 Matlab 的优化目标：直接按元素调用优化的 C/Fortran 例程。此外，Matlab 默认是多线程的，隐式地并行运行元素调用，而 Julia 要求您明确多线程。

目前，2 倍的差异似乎并非不合理。

不过，还是努力吧。这里先介绍几个cmets：

X = X .- mid

您错过了就地操作，使用

X .= X .- mid

相反。这节省了中间数组的分配。

H = 4.0.*hyper.+maximum(abs.(X))

通过标量广播 (hyper) 是徒劳的，最坏的情况是浪费。并且abs.(X) 创建了一个不必要的临时数组。而是使用带有函数输入的maximum 的版本，这样效率更高：

H = 4 * hyper + maximum(abs, X)

这里还有一些不必要的点：

S = (sqrt.(2.0.*pi).*hyper).*exp.((-0.5.*hyper.^2).*((pi.*w./(2.0.*H)).^2))

避免再次通过标量广播并在大多数地方使用整数而不是浮点数：

S = (sqrt(2pi) * hyper) .* exp.((-0.5 * hyper^2 * (pi/2H)^2) .* w.^2)

注意x^(-0.5) 比1/sqrt(x) 慢很多，所以

f = H.^(-0.5).*sin.(pi.*X.*w).*sqrt.(S)

应该是

f = sin.(pi .* X .* w') .* (sqrt.(S)' ./ sqrt(H))

让我们把它放在一起：

function features2(X::VectorFloat64,M::Int,hyper::Float64,mid::Float64)
    X .= X .- mid
    H = 4 * hyper + maximum(abs, X)
    X .= (X .+ H) ./ (2 * H)
    w = 1:M
    S = (sqrt(2pi) * hyper) .* exp.((-0.5 * hyper^2 * (pi/2H)^2) .* w.^2)
    f = sin.(pi .* X .* w') .* (sqrt.(S)' ./ sqrt(H))
    return f
end

基准测试：

jl> X = rand(10000);

jl> M = 100;

jl> hyper = rand();

jl> mid = 0.4;

jl> @btime features($X, $M, $hyper, $mid);
  17.339 ms (9 allocations: 7.86 MiB)

jl> @btime features2($X, $M, $hyper, $mid);
  17.173 ms (4 allocations: 7.63 MiB)

这并不是什么加速。不过，分配较少。问题是运行时很大程度上受sin 广播的支配。

让我们尝试多线程。我有 8 个内核，所以我使用 8 个线程：

function features3(X::VectorFloat64,M::Int,hyper::Float64,mid::Float64)
    X .= X .- mid
    H = 4 * hyper + maximum(abs, X)
    X .= (X .+ H) ./ (2 * H)
    w = transpose(1:M)
    S = (sqrt(2pi) * hyper) .* exp.((-0.5 * hyper^2 * (pi/2H)^2) .* w.^2)
    f = similar(X, length(X), M)
    temp = sqrt.(S) ./ sqrt(H)
    Threads.@threads for j in axes(f, 2)
        wj = w[j]
        tempj = temp[j]
        for i in axes(f, 1)
            @inbounds f[i, j] = tempj * sin(pi * X[i] * w[j])
        end
    end
    return f
end

基准测试：

jl> @btime features3($X, $M, $hyper, $mid);
  1.919 ms (45 allocations: 7.63 MiB)

这好多了，循环和显式线程快 9 倍。

但仍有一些选项可供选择：例如 LoopVectorization.jl。您可以安装这个惊人的软件包，但您需要一个新版本，可能存在一些安装问题，具体取决于您拥有的其他软件包。 LoopVectorization 有两个特别感兴趣的宏，@avx 和 @avxt，前者做了很多工作来矢量化（在模拟意义上）你的代码，单线程，而后者做同样的事情，但是多线程。

using LoopVectorization

function features4(X::VectorFloat64,M::Int,hyper::Float64,mid::Float64)
    X .= X .- mid
    H = 4 * hyper + maximum(abs, X)
    X .= (X .+ H) ./ (2 * H)
    w = collect(1:M)  # I have to use collect here due to some issue with LoopVectorization
    S = (sqrt(2pi) * hyper) .* exp.((-0.5 * hyper^2 * (pi/2H)^2) .* w.^2)
    f = @avx sin.(pi .* X .* w') .* (sqrt.(S)' ./ sqrt(H))
    return f
end

function features4t(X::VectorFloat64,M::Int,hyper::Float64,mid::Float64)
    X .= X .- mid
    H = 4 * hyper + maximum(abs, X)
    X .= (X .+ H) ./ (2 * H)
    w = collect(1:M)  # I have to use collect here due to some issue with LoopVectorization
    S = (sqrt(2pi) * hyper) .* exp.((-0.5 * hyper^2 * (pi/2H)^2) .* w.^2)
    f = @avxt sin.(pi .* X .* w') .* (sqrt.(S)' ./ sqrt(H))
    return f
end

这些函数之间的唯一区别是 @avx 与 @avxt。

基准测试：

jl> @btime features4($X, $M, $hyper, $mid);
  2.695 ms (5 allocations: 7.63 MiB)

单线程情况下的一个非常好的加速。

jl> @btime features4t($X, $M, $hyper, $mid);
  431.700 μs (5 allocations: 7.63 MiB)

多线程 avx 代码的速度是我笔记本电脑上原始代码的 40 倍。还不错吧？

【讨论】：

也值得一提sinpi

我的理解是sinpi更准确，但速度较慢。将 sin(pi*..) 替换为 sinpi(..) 时，我发现速度降低了 70%。

太棒了，非常感谢！但是有一个问题：当乘以标量常量时，为什么将它们保留为整数而不是浮点数（附加 .0）更好？如果我运行简单的测试，让它们作为整数比让它们作为浮点数表现更差

性能应该是一样的，所以我怀疑这种情况下的基准有问题。使用整数的原因有两个：它更具可读性和简洁性。并且它不会无意中更改类型（尽可能多）。例如，如果 x 是整数，2*x 是整数，如果 x 是单精度浮点数，等等。而2.0*x 在这些情况下将结果转换为双精度浮点数。整数的破坏性较小。

@DNF 目前sinpi 速度较慢，但​​它应该一样快，因为它具有更简单的减少。目前它使用sin_kernel 和cos_kernel 来保存代码，但如果它有自己的内核，我敢肯定它会更快。我希望很快能有一个 PR 来解决这个问题。