运行 julia 函数的巨大内存分配？

Question

我尝试在 julia 命令中运行以下函数，但是在对函数计时时，我看到太多的内存分配，我不知道为什么。

function pdpf(L::Int64, iters::Int64)
snr_dB = -10 
snr = 10^(snr_dB/10) 
Pf = 0.01:0.01:1      
thresh = rand(100)
Pd     = rand(100)

for m = 1:length(Pf)
    i = 0
    for k = 1:iters 
        n = randn(L) 
        s = sqrt(snr) * randn(L) 
        y = s + n 

        energy_fin = (y'*y) / L
        @inbounds thresh[m] = erfcinv(2Pf[m]) * sqrt(2/L) + 1

        if energy_fin[1] >= thresh[m]
            i += 1   
        end
    end
    @inbounds Pd[m] = i/iters
end

#thresh = erfcinv(2Pf) * sqrt(2/L) + 1
#Pd_the = 0.5 * erfc(((thresh - (snr + 1)) * sqrt(L)) / (2*(snr + 1)))

end

在我笔记本电脑上的 julia 命令中运行该函数，我得到以下令人震惊的数字：

julia> @time pdpf(1000, 10000)
 17.621551 seconds (9.00 M allocations: 30.294 GB, 7.10% gc time)

我的代码有什么问题？任何帮助表示赞赏。

Answer 1

我认为这种内存分配并不令人惊讶。例如，考虑执行内部循环的所有时间：

for m = 1:length(Pf) 这会给你 100 次处决

for k = 1:iters 这会根据您提供给函数的参数为​​您提供 10,000 次执行。

randn(L) 这会根据您提供给函数的参数为​​您提供一个长度为 1,000 的随机向量。

因此，仅考虑这些，您就会生成 100*10,000*1000 = 10 亿个 Float64 随机数。它们中的每一个都占用 64 位 = 8 字节。 IE。 8GB就在那里。而且，您已经收到了两次对 randn(L) 的调用，这意味着您已经分配了 16GB。

然后您有y = s + n，这意味着另外 8GB 的​​分配，最多为 24GB。我没有详细查看剩余的代码来让你从 24GB 分配到 30GB，但这应该告诉你，GB 分配开始在你的代码中加起来并不难。

如果您正在寻找需要改进的地方，我会提示您可以通过使用正常随机变量的属性来改进这些行：

    n = randn(L) 
    s = sqrt(snr) * randn(L) 
    y = s + n 

通过这种方式，您应该可以轻松地将此处的分配从 24GB 减少到 8GB。请注意，y 在此处将是您定义的正态随机变量，并想出一种方法来生成与 y 现在具有相同分布的正态随机变量。

另一件小事，snr 是函数内部的常量。然而，你继续使用它的sqrt 100 万次。在某些情况下，“检查您的工作”可能会有所帮助，但我认为您可以确信计算机会在第一时间正确处理，因此您无需让它继续重新计算； ）。还有其他类似的地方你可以改进你的代码以避免重复计算，我会留给你去定位。

Answer 2

aireties 给出了一个很好的答案来解释为什么你有这么多的分配。您可以采取更多措施来减少分配次数。使用this property，我们知道y = s+n 确实是y = sqrt(snr) * randn(L) + randn(L)，因此我们可以改为使用y = rvvar*randn(L)，其中rvvar= sqrt(1+sqrt(snr)^2) 在循环外定义（感谢修复！）。这将使所需的随机变量数量减半。

在循环之外，您可以保存 sqrt(2/L) 以减少一点时间。

我不认为转置是特殊情况，所以尝试使用dot(y,y) 而不是y'*y。我知道dot 肯定只是一个无需转置的循环，而另一个可能会根据 Julia 的版本转置。

有助于提高性能（但不是分配）的方法是使用一个大的 randn(L,iters) 并循环遍历它。原因是如果你一次生成所有的随机数，它会更快，因为它可以使用 SIMD 和一堆其他好东西。如果你想隐式地做到这一点而不改变你的代码，你可以使用ChunkedArrays.jl，在那里你可以使用rands = ChunkedArray(randn,L)来初始化它，然后每次你想要randn(L)，你可以使用next(rands)。在 ChunkedArray 内部，它实际上会生成更大的向量并根据需要对其进行补充，但这样您就可以获取您的 randn(L) 而无需跟踪所有这些。

编辑：

ChunkedArrays 可能仅在 L 较小时才能节省时间。这给出了代码：

function pdpf(L::Int64, iters::Int64)
snr_dB = -10
snr = 10^(snr_dB/10)
Pf = 0.01:0.01:1
thresh = rand(100)
Pd     = rand(100)
rvvar= sqrt(1+sqrt(snr)^2)

for m = 1:length(Pf)
    i = 0
    for k = 1:iters
      y = rvvar*randn(L)
      energy_fin = (y'*y) / L
      @inbounds thresh[m] = erfcinv(2Pf[m]) * sqrt(2/L) + 1

      if energy_fin[1] >= thresh[m]
          i += 1
      end
    end
    @inbounds Pd[m] = i/iters
end
end

它的运行时间是使用两个 randn 调用的一半。事实上，我们从 ProfileViewer 中得到：

@profile pdpf(1000, 10000)
using ProfileView
ProfileView.view()

我圈出了y = rvvar*randn(L)这行的两部分，所以绝大多数时间都是随机数生成。上次我检查您仍然可以通过更改为VSL.jl library 来获得不错的随机数生成速度，但是您需要将 MKL 链接到您的 Julia 构建。请注意from the Google Summer of Code page，您可以看到有一个项目可以使用更快的 psudo-rngs 创建一个 repo RNG.jl。看起来它已经实施了一些新的。您可能想查看它们，看看它们是否提供加速（或帮助该项目！）