performance - Julia :优化简单动力学系统的仿真

我正在尝试优化简单动态系统的仿真，在该系统中，网络的响应及其参数(权重)根据简单的线性方程式发展。仿真需要运行数千万个时间步，但是网络规模通常很小。因此，性能受矩阵向量乘积的约束较小，而是受临时数组，边界检查和其他不那么明显的因素约束。由于我是 Julia (Julia)的新手，所以我希望能有任何进一步优化性能的提示。

function train_network(A, T, Of, cs, dt)
    N, I = size(T)
    z    = zeros(I)
    r    = zeros(N)

    @inbounds for t in 1:size(cs, 1)
        # precompute
        Az  = A*z
        Ofr = Of*r

        # compute training signal
        @devec z += dt.*(Az + cs[t] - 0.5.*z)
        I_teach   = T*(Az + cs[t])
        Tz        = T*z

        # rate updates
        @devec r += dt.*(I_teach - Ofr - 0.1.*r)

        # weight updates
        for i in 1:I
            @devec T[:, i] += dt.*1e-3.*(z[i].*r - T[:, i])
        end

        for n in 1:N
            @devec Of[:, n] += dt.*1e-3.*(Tz.*r[n] - Of[:, n])
        end
    end
end

# init parameters
N, I = 20, 2
dt  = 1e-3

# init weights
T = rand(N, I)*N
A = rand(I, I)
Of = rand(N, N)/N

# simulation time & input
sim_T = 2000
ts = 0:dt:sim_T
cs = randn(size(ts, 1), I)

对网络计时(2.000.000步)

@time train_network(A, T, Of, cs, dt)

产生时间

3.420486 seconds (26.12 M allocations: 2.299 GB, 6.65% gc time)

更新1

遵循David Sanders的建议，我摆脱了devec宏并写下了循环。这确实减少了阵列分配，并将性能提高了约25％，以下是新数字:

2.648113 seconds (18.00 M allocations: 1.669 GB, 5.60% gc time)

网络规模越小，提升幅度越大。可以在here中找到更新后的仿真代码的要点。

更新2

许多内存分配是由于矩阵向量乘积。因此，为了摆脱这些缺陷，我通过就地BLAS操作BLAS.genv!替换了这些产品，该操作将时序再减少了25％，并将内存分配减少了90％，

1.990031 seconds (2.00 M allocations: 152.589 MB, 0.69% gc time)

更新了代码here。

更新3

最大的1级更新也可以由对本地BLAS函数的两个调用BLAS.scal!代替。用于缩放和BLAS.ger!进行1级更新。需要注意的是，如果使用了一个以上的线程(OpenBLAS是否有问题？)，这两个调用都将非常缓慢。

blas_set_num_threads(1)

网络规模为20时，时序增加15％，网络规模为50时，时序增加50％。没有更多的内存分配，新的时序为

1.638287 seconds (11 allocations: 1.266 KB)

同样，可以在here中找到更新的代码。

更新4

我写了一个基本的Cython script来比较到目前为止的结果。主要区别在于，我不使用任何对BLAS的调用，而是有循环:在Cython中，注入(inject)低级BLAS调用是很痛苦的，而对于较小的网络，调用numpy dot会产生太多开销(我尝试过...)。时间是

CPU times: user 3.46 s, sys: 6 ms, total: 3.47 s, Wall time: 3.47 s

与原始版本大致相同(到目前为止，已从原始版本中删除了50％)。

最佳答案

尽管您使用的是Devectorize.jl包，但我建议您仅将所有这些矢量化操作显式地写出为简单循环。我希望这将为您带来显着的性能提升。
Devectorize软件包无疑是一个很大的贡献，但是要查看它跳到麻烦的圆环上来为您做些肮脏的工作，您可以执行以下操作(来自README包的示例):

using Devectorize

a = rand(2,2);
b = rand(2,2);
c = rand(2,2);

julia> macroexpand(:(@devec r = exp(a + b) .* sum(c)))

在这里，macroexpand是一个函数，用于告诉您@devec宏将其自变量扩展为的代码(该行其余部分的代码)。
我不会在这里显示输出来破坏惊喜，但是您将要手工编写的不仅仅是简单的for循环。

此外，您拥有大量分配的事实表明并非所有矢量操作都得到了正确处理。

顺便说一句，不要忘了先进行少量运行，这样就不会安排编译步骤的时间。

[切记:这里，exp是将通常的指数函数应用于矩阵的每个元素的函数，等效于map(exp, a+b)。 expm给出矩阵的指数。有人谈论过不赞成使用exp。]

关于performance - Julia :优化简单动力学系统的仿真，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/32845996/