存储微分方程的解耗时太长

Question

我正在解一个微分方程，其解的大小为(2000,200000)，例如,大小(sol) =(2000, 2000000)。现在我想把前1000个变量的时间演变存储到一个数组中。例如，

x=sol[1:1000,:];

，从这里我需要计算一个值

r1= abs.(mean(exp.(x*1im),dims=1)[1,:])[1:end]

现在将解存储到x，然后计算r1花费太多时间，甚至导致系统崩溃。有什么解决办法吗?

我使用Atom来运行代码，我的系统的RAM是48 GB。

Answer 1

我将假设数组sol是Matrix{Float64}，尽管这将是在您的问题中包含的重要内容。

然后我将指出您在哪些地方创建了不必要的分配:

x = sol[1:1000,:]
r1 = abs.(mean(exp.(x*1im),dims=1)[1,:])[1:end]

从最里面开始:

1. x*1im创建一个与x相同大小的全新数组。
2. exp.()创建另一个与x相同大小的完整数组。
3. mean创建另一个数组，尽管这是一个合理的期望。
4. [1,:]创建一个与x相同大小的新数组
5. [1:end]还创建了输出矢量的副本，这似乎根本没有任何目的。
6. x = sol[1:1000, :]分配一个大数组

我们想去掉这些:

1。您可以在x后面写exp.(x .* 1im)，但使用cis函数更有效，它与exp(1im * x)做同样的事情

3. 广播很好，但由于您只想沿列平均，因此不需要临时数组，只需写入mean(cis, x; dims=1)即可，这不会创建临时数组。
4. 你可以使用vec函数从矩阵中强制一个矢量而不复制数据。
5. 尝试制作view以避免复制。

结果:

# original
function absmean(sol)
x = sol[1:1000,:]
r1 = abs.(mean(exp.(x*1im),dims=1)[1,:])[1:end]
return r1
end
# new version
function absmean2(sol)
x = @view sol[1:1000, :]
r1 = abs.(mean(cis, x; dims=1)) |> vec
return r1
end

基准:

# making a smaller version to reduce benchmarking time
julia> sol = randn(2000, 2000);
julia> @btime absmean($sol);
72.293 ms (18 allocations: 76.39 MiB)
julia> @btime absmean2($sol);
29.985 ms (15 allocations: 47.88 KiB)

不是，快得多，但小于1/1000的内存使用。你可以做很多事情来提高性能，但这些都是"容易摘的果实"。

最后的话。似乎您的解决方案数据是沿着sol的行而不是沿着列定向的。这有点奇怪，因为Julia数组是列为主的。在这个特定的计算中，这无关紧要，但在其他情况下，它可能对性能有害。

编辑:我将添加一个版本来展示如何获得更显著的加速。这使用了LoopVectorization。Jl具有多线程(6核，12线程)。为了让LoopVectorization接受代码，需要重新编写一些代码:

using LoopVectorization
function absmean4(sol)
x = @view sol[1:1000, :]
r = similar(x, size(x, 2))
# use @turbo instead of @tturbo for single-threaded
@tturbo for j in axes(x, 2)  
s = zero(eltype(x))
c = zero(eltype(x))
for i in axes(x, 1)
(a, b) = sincos(x[i, j])
s += a
c += b
end
r[j] = sqrt(s^2 + c^2) / size(x, 1)
end
return r
end
julia> @btime absmean4($sol);
1.830 ms (1 allocation: 15.75 KiB)

这比原来的速度提高了40倍，内存使用减少了5000倍。

如果你想完全全力以赴，你可以预先分配输出向量r并将其传递给你的函数。那么你将有零分配。