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

Question

我试图优化一个简单的动力系统的仿真，其中网络的响应以及它的参数(权值)都是根据简单的线性方程组来演化的。模拟需要运行数以千万计的时间步骤，但网络规模通常很小。因此，性能不受矩阵向量乘积的约束，而是由临时数组、绑定检查和其他较不明显的因素所限制。由于我是新的朱莉娅，我希望任何提示，以进一步优化性能。

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的建议，我去掉了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！一级更新。请注意，如果使用多个线程(OpenBLAS问题？)，则这两个调用都相当慢，所以最好设置

blas_set_num_threads(1)

如果网络大小为20，则时间上的增益为15%，而规模为50的网络的时间增益为50%。没有更多的内存分配，新的时间安排是

1.638287 seconds (11 allocations: 1.266 KB)

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

更新4

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

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

这与原来的版本大致相同(到目前为止，有50%是剃掉的)。

Answer 1

虽然您使用的是Devectorize.jl包，但我建议您将所有这些矢量化操作显式地写出为简单循环。我希望这将给你一个显着的性能提升。

Devectorize包无疑是一个很大的贡献，但是要查看它为您完成脏工作的循环，您可以这样做(来自包自述的一个例子)：

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的这种用途。