疯狂的空间分配求解随机微分方程组

Question

第一次在这里提问。我以前使用一个简单的 MATLAB 脚本对 90 个 Hopf 振荡器进行建模，这些振荡器通过矩阵耦合，具有 randn 噪声，并具有简单的欧拉阶跃积分。我想升级这个，所以我进入了 Julia，似乎有很多令人兴奋的特性。

这是我正在求解的方程组

我有点迷失了。我开始使用 Differentialequations.jl (stochastic solver) ，找到了一个解决方案，并发现自己有一个基准，告诉我解决 200 秒占用 4 Gb ！（2.5 Gb with alg_hints=[:stiff]）（我没有修复 dt，以前我使用 dt=0.1）

function Shopf(du,u,p,t)

  du[1:90,1]=(p[1:90,1]-u[1:90,1].^2.0-u[1:90,2].^2.0).*u[1:90,1]-p[1:90,2].*u[1:90,2] + 0.5*(-p[: , end].*u[:,1]+p[:,4:end-1] *u[:,1])
  du[1:90,2]=(p[1:90,1]-u[1:90,1].^2.0-u[1:90,2].^2.0).*u[1:90,1]+p[1:90,2].*u[1:90,1] + 0.5*(-p[: , end].*u[:,2]+p[:,4:end-1] *u[:,2])

end


function ?_Shopf(du,u,p,t)

du[1:90,1]=0.04*ones(90,1)
du[1:90,2]=0.04*ones(90,1)


end


#initial condition
u0=-0.1*ones(90,2);
#initial time
t0=0.0;
#final time
tend=200.0;
#setting parameter matrix
p0=[0.1 , 2*pi*0.04]
push!(p0,-p0[2])
p=p0'.*ones(90,3);
SC=SC;
p=[p SC]
p=[p sum(SC,dims=2)]
#


#col 1 :alpha
#col 2-3 : [w0 -w0]

#col 3-93 : coupling matrix
#col 94: col-wise sum of coupling matrix


@benchmark solve(prob_sde_Shopf,nlsolver=Rosenbrock23(),alg_hints=[:stiff])

BenchmarkTools.Trial：内存估计：2.30 GiB

分配估计：722769

---

最短时间：859.224 毫秒（13.24% GC）

中位时间：942.707 毫秒（13.10% GC）

平均时间：975.430 毫秒（12.99% GC）

最长时间：1.223 秒（13.00% GC）

---

样品：6

评估/样本：1

有什么想法吗？我正在检查几种解决方案，但它们都没有将内存量减少到合理的数量。提前致谢。

Answer 1

您正在创建数量惊人的临时数组。每个切片都会创建一个临时的。你在这里和那里放一个点，但你必须把所有的点都加进去才能得到融合广播。相反，您可以使用@.宏来为您完成。此外，使用@views将确保切片不会复制：

function Shopf(du, u, p, t)
    @. du[1:90, 1] = @views (p[1:90, 1] - u[1:90, 1]^2 - u[1:90, 2]^2) * u[1:90, 1] -
        p[1:90, 2] * u[1:90,2] + 0.5 * (-p[:, end] * u[:, 1] + p[:, 4:end-1] * u[:,1])
    @. du[1:90, 2] = @views (p[1:90, 1] - u[1:90, 1]^2 - u[1:90, 2]^2) * u[1:90, 1] + 
        p[1:90, 2] * u[1:90,1] + 0.5 * (-p[:, end] * u[:, 2] + p[:, 4:end-1] * u[:,2])  
end

还有，不要写x^2.0，用x^2，前者是慢浮力，后者是快x * x。事实上，在乘法、加法等任何可能的地方尝试使用整数。

这是另一件事

function ?_Shopf(du,u,p,t)

du[1:90,1]=0.04*ones(90,1)
du[1:90,2]=0.04*ones(90,1)


end

无需在赋值右侧创建两个临时数组。只写这个：

function ?_Shopf(du, u, p, t)
    du[1:90, 1:2] .= 0.04
end

更快更简单。请注意，我尚未对此进行测试，因此请更正任何错别字。

（最后，请使用缩进并在运算符周围放置空格，它使您的代码更易于阅读。）

更新：我真的不知道你的代码应该做什么，奇怪的索引是什么，但这里有一个可能的改进，它只使用循环（我认为它实际上更干净，并且会让你做进一步的优化）：

生产操作A是一个矩阵乘积，所以你不能避免在那里分配，除非你可以传入一个缓存数组来处理，使用mul!. 除此之外，您在下面不应有任何分配。

function shopf!(du, u, p, t)
    A = @view p[:, 4:end-1] * u
    # mul!(A, view(p, 4:end-1), u)  # in-place matrix product
    for i in axes(u, 1)
        val = (p[i, 1] - u[i, 1]^2 - u[i, 2]^2) * u[i, 1]  # don't calculate this twice
        du[i, 1] = val - (p[i, 2] * u[i, 2]) - (0.5 * p[i, end] * u[i, 1]) + 
            (0.5 * A[i, 1])
        du[i, 2] = val + (p[i, 2] * u[i, 1]) - (0.5 * p[i, end] * u[i, 2]) + 
            (0.5 * A[i, 2])
    end
end

在此之后，您可以添加各种优化，@inbounds如果您确定数组大小、多线程，@simd甚至@avx来自 LoopVectorization 实验包。