I tried porting the NN code presented here to Julia, hoping for a speed increase in training the network. On my desktop, this proved to be the case.
However, on my MacBook, Python + numpy beats Julia by miles.
Training with the same parameters, Python is more than twice as fast as Julia (4.4s vs 10.6s for one epoch). Considering that Julia is faster than Python (by ~2s) on my desktop, it seems like there's some resource that Python/numpy is utilizing on the mac that Julia isn't. Even parallelizing the code only gets me down to ~6.6s (although this might be due to me not being that experienced in writing parallel code). I thought the problem might be that Julia's BLAS was slower than the vecLib library used natively in mac, but experimenting with different builds didn't seem to get me much closer. I tried building both with USE_SYSTEM_BLAS = 1, and building with MKL, of which MKL gave the faster result (the times posted above).
I'll post my version info for the laptop as well as my Julia implementation below for reference. I don't have access to the desktop at this time, but I was running the same version of Julia on Windows, using openBLAS, comparing with a clean installation of Python 2.7 also using openBLAS.
Is there something I'm missing here?
EDIT: I know that my Julia code leaves a lot to be desired in terms of optimization, I really appreciate any tips to make it faster. However, this is not a case of Julia being slower on my laptop but rather Python being much faster. On my desktop, Python runs one epoch in ~13 seconds, on the laptop it only takes ~4.4s. What I'm interested in the most is where this difference comes from. I realize the question may be somewhat poorly formulated.
Versions on laptop:
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 code (sequential):
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
z = w*activation .+ b
that kind of code is creating arrays when they aren't required. Why not just have a cache array to make this non-allocating? Then your arrays of arrays should be StaticVectors of arrays (that's a small difference). Also, are you timing with@btime
or including Julia's startup + JIT time? If you're using amain
function and calling it from the command line, likely half of the time is just booting up Julia+LLVM and not actually running your script. This is not a recommended way to run Julia. – Chris Rackauckas Apr 8 '18 at 15:19a .= σ.(w*a .+ b)
instead ofa = σ.(w*a + b)
. Alsoactivation .= σ.(z)
and thenactivations[i+1] .= activation
. Remember, slicing without a view produces a copy. – Chris Rackauckas Apr 8 '18 at 16:23