我尝试移植 NN 代码 here给 Julia ，希望能提高网络训练的速度。在我的桌面上，事实证明是这样。

然而，在我的 MacBook 上，Python + numpy 远远击败了 Julia。
使用相同的参数进行训练，Python 的速度是 Julia 的两倍多(4.4s vs 10.6s for one epoch)。考虑到 Julia 在我的桌面上比 Python 快(大约 2 秒)，似乎 Python/numpy 在 mac 上使用了一些 Julia 没有使用的资源。即使并行化代码也只能让我下降到 ~6.6s(尽管这可能是因为我在编写并行代码方面没有那么丰富的经验)。我认为问题可能是 Julia 的 BLAS 比 mac 中本地使用的 vecLib 库慢，但尝试不同的构建似乎并没有让我更接近。我尝试使用 USE_SYSTEM_BLAS = 1 进行构建，并使用 MKL 进行构建，其中 MKL 给出了更快的结果(上面发布的时间)。

我将在下面发布笔记本电脑的版本信息以及我的 Julia 实现以供引用。我此时无法访问桌面，但我在 Windows 上运行相同版本的 Julia，使用 openBLAS，与同样使用 openBLAS 的全新安装的 Python 2.7 相比。

我在这里缺少什么吗？

编辑:我知道我的 Julia 代码在优化方面还有很多不足之处，我真的很感谢任何让它更快的提示。然而，这并不是 Julia 在我的笔记本电脑上变慢的情况，而是 Python 快得多的情况。在我的桌面上，Python 在大约 13 秒内运行一个 epoch，在笔记本电脑上只需要大约 4.4 秒。我最感兴趣的是这种差异从何而来。我意识到这个问题可能有些表述不当。

笔记本电脑版本:

julia> versioninfo()
Julia Version 0.6.2
Commit d386e40c17 (2017-12-13 18:08 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin17.4.0)
  CPU: Intel(R) Core(TM) i5-7360U CPU @ 2.30GHz
  WORD_SIZE: 64
  BLAS: libmkl_rt
  LAPACK: libmkl_rt
  LIBM: libopenlibm
  LLVM: libLLVM-3.9.1 (ORCJIT, broadwell)

Python 2.7.14 (default, Mar 22 2018, 14:43:05) 
[GCC 4.2.1 Compatible Apple LLVM 9.0.0 (clang-900.0.39.2)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> import numpy
>>> numpy.show_config()
lapack_opt_info:
    extra_link_args = ['-Wl,-framework', '-Wl,Accelerate']
    extra_compile_args = ['-msse3']
    define_macros = [('NO_ATLAS_INFO', 3), ('HAVE_CBLAS', None)]
openblas_lapack_info:
  NOT AVAILABLE
atlas_3_10_blas_threads_info:
  NOT AVAILABLE
atlas_threads_info:
  NOT AVAILABLE
openblas_clapack_info:
  NOT AVAILABLE
atlas_3_10_threads_info:
  NOT AVAILABLE
atlas_blas_info:
  NOT AVAILABLE
atlas_3_10_blas_info:
  NOT AVAILABLE
atlas_blas_threads_info:
  NOT AVAILABLE
openblas_info:
  NOT AVAILABLE
blas_mkl_info:
  NOT AVAILABLE
blas_opt_info:
    extra_link_args = ['-Wl,-framework', '-Wl,Accelerate']
    extra_compile_args = ['-msse3', '-I/System/Library/Frameworks/vecLib.framework/Headers']
    define_macros = [('NO_ATLAS_INFO', 3), ('HAVE_CBLAS', None)]
blis_info:
  NOT AVAILABLE
atlas_info:
  NOT AVAILABLE
atlas_3_10_info:
  NOT AVAILABLE
lapack_mkl_info:
  NOT AVAILABLE

Julia 代码(顺序):

using MLDatasets

mutable struct network
    num_layers::Int64
    sizearr::Array{Int64,1}
    biases::Array{Array{Float64,1},1}
    weights::Array{Array{Float64,2},1}
end

function network(sizes)
    num_layers = length(sizes)
    sizearr = sizes
    biases = [randn(y) for y in sizes[2:end]]
    weights = [randn(y, x) for (x, y) in zip(sizes[1:end-1], sizes[2:end])]

    network(num_layers, sizearr, biases, weights)
end

σ(z) = 1/(1+e^(-z))
σ_prime(z) = σ(z)*(1-σ(z))

function (net::network)(a)
    for (w, b) in zip(net.weights, net.biases)
        a = σ.(w*a + b)
    end
    return a
end

function SGDtrain(net::network, training_data, epochs, mini_batch_size, η, test_data=nothing)
    n_test = test_data != nothing ? length(test_data):nothing
    n = length(training_data)

    for j in 1:epochs
        training_data = shuffle(training_data)
        mini_batches = [training_data[k:k+mini_batch_size-1] for k in 1:mini_batch_size:n]

        @time for batch in mini_batches
            update_batch(net, batch, η)
        end

        if test_data != nothing
            println("Epoch ", j,": ", evaluate(net, test_data), "/", n_test)
        else
            println("Epoch ", j," complete.")
        end
    end
end

function update_batch(net::network, batch, η)
    ∇_b = net.biases .- net.biases
    ∇_w = net.weights .- net.weights

    for (x, y) in batch
        δ_∇_b, δ_∇_w = backprop(net, x, y)
        ∇_b += δ_∇_b
        ∇_w += δ_∇_w
    end

    net.biases -= (η/length(batch))∇_b
    net.weights -= (η/length(batch))∇_w
end

function backprop(net::network, x, y)
    ∇_b = copy(net.biases)
    ∇_w = copy(net.weights)

    len = length(net.sizearr)
    activation = x
    activations = Array{Array{Float64,1}}(len)
    activations[1] = x
    zs = copy(net.biases)

    for i in 1:len-1
        b = net.biases[i]; w = net.weights[i]
        z = w*activation .+ b
        zs[i] = z
        activation = σ.(z)
        activations[i+1] = activation[:]
    end

    δ = (activations[end] - y) .* σ_prime.(zs[end])
    ∇_b[end] = δ[:]
    ∇_w[end] = δ*activations[end-1]'

    for l in 1:net.num_layers-2
        z = zs[end-l]
        δ = net.weights[end-l+1]'δ .* σ_prime.(z)
        ∇_b[end-l] = δ[:]
        ∇_w[end-l] = δ*activations[end-l-1]'
    end
    return (∇_b, ∇_w)
end

function evaluate(net::network, test_data)
    test_results = [(findmax(net(x))[2] - 1, y) for (x, y) in test_data]
    return sum(Int(x == y) for (x, y) in test_results)
end

function loaddata(rng = 1:50000)
    train_x, train_y = MNIST.traindata(Float64, Vector(rng))
    train_x = [train_x[:,:,x][:] for x in 1:size(train_x, 3)]
    train_y = [vectorize(x) for x in train_y]
    traindata = [(x, y) for (x, y) in zip(train_x, train_y)]

    test_x, test_y = MNIST.testdata(Float64)
    test_x = [test_x[:,:,x][:] for x in 1:size(test_x, 3)]
    testdata = [(x, y) for (x, y) in zip(test_x, test_y)]
    return traindata, testdata
end

function vectorize(n)
    ev = zeros(10,1)
    ev[n+1] = 1
    return ev
end

function main()
    net = network([784, 30, 10])
    traindata, testdata = loaddata()
    SGDtrain(net, traindata, 10, 10, 1.25, testdata)
end

最佳答案

我首先运行您的代码:

7.110379 seconds (1.37 M allocations: 20.570 GiB, 19.81%gc time)
Epoch 1: 7960/10000
6.147297 seconds (1.27 M allocations: 20.566 GiB, 18.33%gc time)

哎呀，每个时期分配 21GiB？那是你的问题。它对垃圾收集产生了很大的影响，你的计算机内存越少，它就必须处理越多。因此，让我们来解决这个问题。

主要思想是预先分配缓冲区，然后修改数组而不是创建新数组。在您的代码中，您开始使用 backprop:

∇_b = copy(net.biases)
∇_w = copy(net.weights)

len = length(net.sizearr)
activation = x
activations = Array{Array{Float64,1}}(len)
activations[1] = x
zs = copy(net.biases)

您正在使用 copy 的事实意味着您可能应该预先分配一些东西!因此，让我们从 zs 和 activations 开始。我扩展了您的网络以容纳这些缓存阵列:

mutable struct network
    num_layers::Int64
    sizearr::Array{Int64,1}
    biases::Array{Array{Float64,1},1}
    weights::Array{Array{Float64,2},1}
    zs::Array{Array{Float64,1},1}
    activations::Array{Array{Float64,1},1}
end

function network(sizes)
    num_layers = length(sizes)
    sizearr = sizes
    biases = [randn(y) for y in sizes[2:end]]
    weights = [randn(y, x) for (x, y) in zip(sizes[1:end-1], sizes[2:end])]
    zs = [randn(y) for y in sizes[2:end]]
    activations = [randn(y) for y in sizes[1:end]]
    network(num_layers, sizearr, biases, weights, zs, activations)
end

然后我更改了您的 backprop 以使用这些缓存:

function backprop(net::network, x, y)
    ∇_b = copy(net.biases)
    ∇_w = copy(net.weights)

    len = length(net.sizearr)
    activations = net.activations
    activations[1] .= x
    zs = net.zs

    for i in 1:len-1
        b = net.biases[i]; w = net.weights[i];
        z = zs[i]; activation = activations[i+1]
        z .= w*activations[i] .+ b
        activation .= σ.(z)
    end

    δ = (activations[end] - y) .* σ_prime.(zs[end])
    ∇_b[end] = δ[:]
    ∇_w[end] = δ*activations[end-1]'

    for l in 1:net.num_layers-2
        z = zs[end-l]
        δ = net.weights[end-l+1]'δ .* σ_prime.(z)
        ∇_b[end-l] = δ[:]
        ∇_w[end-l] = δ*activations[end-l-1]'
    end
    return (∇_b, ∇_w)
end

这导致分配的内存大幅减少。但是还有很多事情要做。首先让我们将 * 更改为 A_mul_B!。这个函数是一个矩阵乘法，它写入数组 C (A_mul_B!(C,A,B)) 而不是创建一个新矩阵，这可以大大减少你的内存分配。所以我做了:

for l in 1:net.num_layers-2
    z = zs[end-l]
    δ = net.weights[end-l+1]'δ .* σ_prime.(z)
    ∇_b[end-l] .= vec(δ)
    atransp = activations[end-l-1]'
    A_mul_B!(∇_w[end-l],δ,atransp)
end

但是我没有使用分配的 '，而是使用 reshape 因为我只想要一个 View :

for l in 1:net.num_layers-2
    z = zs[end-l]
    δ = net.weights[end-l+1]'δ .* σ_prime.(z)
    ∇_b[end-l] .= vec(δ)
    atransp = reshape(activations[end-l-1],1,length(activations[end-l-1]))
    A_mul_B!(∇_w[end-l],δ,atransp)
end

(它还实现了更快的 OpenBLAS 调度。尽管 MKL 可能有所不同)。但是你还在复制

    ∇_b = copy(net.biases)
    ∇_w = copy(net.weights)

并且您在每一步都分配了一堆 δs，所以我所做的下一个更改是预先分配这些 δs 并就地完成(它看起来就像之前的更改一样)。

然后我做了一些分析。在 Juno 中，这只是:

@profile main()
Juno.profiler()

或者如果您不使用 Juno，您可以将第二部分替换为 ProfileView.jl .我得到了:

所以大部分时间都花在了BLAS上，但是有一个问题。看到像 ∇_w += δ_∇_w 这样的操作正在创建一堆矩阵!相反，我们希望遍历每个矩阵并通过其更改矩阵就地更新每个矩阵。这扩展为:

function update_batch(net::network, batch, η)
    ∇_b = net.∇_b
    ∇_w = net.∇_w

    for i in 1:length(∇_b)
      fill!(∇_b[i],0.0)
    end

    for i in 1:length(∇_w)
      fill!(∇_w[i],0.0)
    end

    for (x, y) in batch
        δ_∇_b, δ_∇_w = backprop(net, x, y)
        ∇_b .+= δ_∇_b
        for i in 1:length(∇_w)
          ∇_w[i] .+= δ_∇_w[i]
        end
    end

    for i in 1:length(∇_b)
      net.biases[i] .-= (η/length(batch)).*∇_b[i]
    end

    for i in 1:length(∇_w)
      net.weights[i] .-= (η/length(batch)).*∇_w[i]
    end
end

我沿着相同的路线做了一些更多的更改，我的最终代码如下:

mutable struct network
    num_layers::Int64
    sizearr::Array{Int64,1}
    biases::Array{Array{Float64,1},1}
    weights::Array{Array{Float64,2},1}
    weights_transp::Array{Array{Float64,2},1}
    zs::Array{Array{Float64,1},1}
    activations::Array{Array{Float64,1},1}
    ∇_b::Array{Array{Float64,1},1}
    ∇_w::Array{Array{Float64,2},1}
    δ_∇_b::Array{Array{Float64,1},1}
    δ_∇_w::Array{Array{Float64,2},1}
    δs::Array{Array{Float64,2},1}
end

function network(sizes)
    num_layers = length(sizes)
    sizearr = sizes
    biases = [randn(y) for y in sizes[2:end]]
    weights = [randn(y, x) for (x, y) in zip(sizes[1:end-1], sizes[2:end])]
    weights_transp = [randn(x, y) for (x, y) in zip(sizes[1:end-1], sizes[2:end])]
    zs = [randn(y) for y in sizes[2:end]]
    activations = [randn(y) for y in sizes[1:end]]
    ∇_b = [zeros(y) for y in sizes[2:end]]
    ∇_w = [zeros(y, x) for (x, y) in zip(sizes[1:end-1], sizes[2:end])]
    δ_∇_b = [zeros(y) for y in sizes[2:end]]
    δ_∇_w = [zeros(y, x) for (x, y) in zip(sizes[1:end-1], sizes[2:end])]
    δs = [zeros(y,1) for y in sizes[2:end]]
    network(num_layers, sizearr, biases, weights, weights_transp, zs, activations,∇_b,∇_w,δ_∇_b,δ_∇_w,δs)
end

function update_batch(net::network, batch, η)
    ∇_b = net.∇_b
    ∇_w = net.∇_w

    for i in 1:length(∇_b)
      ∇_b[i] .= 0.0
    end

    for i in 1:length(∇_w)
      ∇_w[i] .= 0.0
    end

    δ_∇_b = net.δ_∇_b
    δ_∇_w = net.δ_∇_w

    for (x, y) in batch
        backprop!(net, x, y)
        for i in 1:length(∇_b)
          ∇_b[i] .+= δ_∇_b[i]
        end
        for i in 1:length(∇_w)
          ∇_w[i] .+= δ_∇_w[i]
        end
    end

    for i in 1:length(∇_b)
      net.biases[i] .-= (η/length(batch)).*∇_b[i]
    end

    for i in 1:length(∇_w)
      net.weights[i] .-= (η/length(batch)).*∇_w[i]
    end
end

function backprop!(net::network, x, y)
    ∇_b = net.δ_∇_b
    ∇_w = net.δ_∇_w

    len = length(net.sizearr)
    activations = net.activations
    activations[1] .= x
    zs = net.zs
    δs = net.δs

    for i in 1:len-1
        b = net.biases[i]; w = net.weights[i];
        z = zs[i]; activation = activations[i+1]
        A_mul_B!(z,w,activations[i])
        z .+= b
        activation .= σ.(z)
    end

    δ = δs[end]
    δ .= (activations[end] .- y) .* σ_prime.(zs[end])
    ∇_b[end] .= vec(δ)
    atransp = reshape(activations[end-1],1,length(activations[end-1]))
    A_mul_B!(∇_w[end],δ,atransp)

    for l in 1:net.num_layers-2
        z = zs[end-l]
        transpose!(net.weights_transp[end-l+1],net.weights[end-l+1])
        A_mul_B!(δs[end-l],net.weights_transp[end-l+1],δ)
        δ = δs[end-l]
        δ .*= σ_prime.(z)
        ∇_b[end-l] .= vec(δ)
        atransp = reshape(activations[end-l-1],1,length(activations[end-l-1]))
        A_mul_B!(∇_w[end-l],δ,atransp)
    end
    return nothing
end

其他一切保持不变。为确保完成，我将 @time 添加到 backprop 调用并获取:

0.000070 seconds (8 allocations: 352 bytes)
0.000066 seconds (8 allocations: 352 bytes)
0.000090 seconds (8 allocations: 352 bytes)

所以这是非分配的。我将 @time 添加到 for (x, y) in batch 循环并获取

0.000636 秒(80 次分配:3.438 KiB) 0.000610 秒(80 次分配:3.438 KiB) 0.000624 秒(80 次分配:3.438 KiB)

所以这告诉我基本上所有剩余的分配都来自迭代器(这可以改进，但可能不会改进时间)。所以最后的时间是:

Epoch 2: 8428/10000
  4.005540 seconds (586.87 k allocations: 23.925 MiB)
Epoch 1: 8858/10000
  3.488674 seconds (414.49 k allocations: 17.082 MiB)
Epoch 2: 9104/10000

这在我的机器上几乎快了 2 倍，但每个循环的内存分配减少了 1200 倍。这意味着在 RAM 较慢和较小的机器上，这种方法应该做得更好(我的桌面有相当多的内存，所以它真的不太关心!)。

最终的配置文件显示大部分时间都在 A_mul_B! 调用中，因此几乎所有内容现在都受到我的 OpenBLAS 速度的限制，所以我完成了。我可以做的一些额外的事情是多线程一些其他循环，但是给出分析的返回会很小所以我会把它留给你(基本上只是把 Threads.@threads 放在像 ∇_w[i] .+= δ_∇_w[i]).

希望这不仅能改进您的代码，还能教会您如何分析、预分配、使用就地操作以及考虑性能。

关于python - macOS Python 和 numpy 在训练神经网络方面比 Julia 更快，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/49719076/