Neural Network in Julia

Question

I am currently trying to implement a Neural Net in Julia with the goal of eventually implementing a stacked autoencoder. My code seems to work but I would appreciate any constructive criticism. If there exists a style guide for Julia, I am not concerned with that. However, any other comments would be very much welcome. I would also like to be able to write an implementation that can be extended to more complicated architectures without making significant alterations to the basics of the code. This is not that in any way but ideas on how to do this would be very helpful.

type ANN2


#
# Neural Network type...
#


# define vars
weights::Dict
bias::Dict
As::Dict
Ns::Dict
Fs::Dict
Ss::Dict
weightdelta::Dict
biasdelta::Dict
shape::Array{Int64,1}
numlayers::Int64
averror::Float64

# define methods
forward::Function
calcuate_deltas::Function
init::Function
setshape::Function
sgm::Function
updateone::Function
updateepoch::Function
calculate_error::Function

# Constructer
function ANN2()
    this = new ()

    this.weights = Dict{Int64,Any}()
    this.bias = Dict{Int64,Any}()
    this.As = Dict{Int64,Any}()
    this.Ns = Dict{Int64,Any}()
    this.Fs = Dict{Int64,Any}()
    this.weightdelta = Dict{Int64,Any}()
    this.biasdelta = Dict{Int64,Any}()
    this.Ss = Dict{Int64,Any}()
    this.numlayers = 0

    # Set the shape of the network
    this.setshape = function(shape)
        this.shape = shape
        this.numlayers = size(this.shape)[1] - 1
        return nothing
    end

    # initialise weights and bias
    this.init = function()
        for (ind,(a,b)) in enumerate(zip(this.shape[1:end-1],this.shape[2:end]))
            this.weights[ind] = rand(b,a)
            this.bias[ind] = rand(b)
        end
        return nothing
    end

    # Calculate output of network given one input
    this.forward = function (input::Array{Float64,1})
        this.As[0] = input 
        for i = 1:this.numlayers                
            this.Ns[i] = net.weights[i]*this.As[i-1] + net.bias[i]
            this.As[i] = this.sgm(this.Ns[i])
            this.Fs[i] = this.As[i].*(1-this.As[i])
        end
        return this.As[this.numlayers]
    end

    # calculate weight and bias updates
    # if avg is true then updates are accumulated 
    # if avg is false then updates are overwritten
    this.calcuate_deltas = function (input::Array{Float64,1},target::Array{Float64,1},rate::Float64,avg::Bool)
        this.forward(input)
        for i in reverse(1:this.numlayers)
            if i == this.numlayers
                this.Ss[i] = this.Fs[i].*(this.As[i] - target)
                if avg
                    this.weightdelta[i] = this.weightdelta[i]+rate.*(this.Ss[i]*this.As[i-1]')
                    this.biasdelta[i] = this.biasdelta[i]+rate.*this.Ss[i]
                else
                    this.weightdelta[i] = rate.*(this.Ss[i]*this.As[i-1]')
                    this.biasdelta[i] = rate.*this.Ss[i]
                end
            else
                this.Ss[i] = this.Fs[i].*(this.weights[i+1]'*this.Ss[i+1])
                if avg
                    this.weightdelta[i] = this.weightdelta[i]+rate.*(this.Ss[i]*this.As[i-1]')
                    this.biasdelta[i] = this.biasdelta[i]+rate.*this.Ss[i]
                else
                    this.weightdelta[i] = rate.*(this.Ss[i]*this.As[i-1]')
                    this.biasdelta[i] = rate.*this.Ss[i]
                end
            end
        end
        return nothing
    end

    # calculate new weights and bias from one input target pair
    this.updateone = function(input::Array{Float64,1},target::Array{Float64,1},rate::Float64)
        this.calcuate_deltas(input,target,rate,false)
        for i in 1:this.numlayers
            this.weights[i] = this.weights[i] - this.weightdelta[i]
            this.bias[i] = this.bias[i] - this.biasdelta[i]
        end
        return nothing
    end

    # calculate new weights and bias from training set
    # randomly sample from training set n (cases) input target pairs
    # update weights and bias by averaging updates for each pair
    this.updateepoch = function(cases::Int64,inputs::Dict,targets::Dict,rate::Float64)
        this.updateone(inputs[1],targets[1],rate)
        for i in 1:cases
            ind = rand(1:length(inputs))
            input = inputs[ind]
            target = targets[ind]
            this.calcuate_deltas(input,target,rate,true)
        end

        for i in 1:this.numlayers
            this.weightdelta[i] = (1/cases).*this.weightdelta[i]
            this.weights[i] = this.weights[i] - this.weightdelta[i]
            this.biasdelta[i] = (1/cases).*this.biasdelta[i]
            this.bias[i] = this.bias[i] - this.biasdelta[i]
        end
    end

    # sigmoid function
    this.sgm = function(x::Array{Float64,1})
            return 1./(1+exp(-x))
    end

    # calculate current error for one input target pair
    this.calculate_error = function(input::Array{Float64,1},target::Array{Float64,1})
        this.forward(input)
        return (this.As[this.numlayers] - target)'*(this.As[this.numlayers] - target)
    end

    return this
end

Answer 1

Coding Style

• "Methods", as in Functions which are part of a type are not Julianic.
  • Instead use Functions with typed parameters
  • Rather than ann.forward(input::Array{Float64,1})
  • use forward(ann::ANN2, input::Array{Float64,1})
• Functions which mutate there inputs should end with a bang (!)
  • So infact: forward!(ann::ANN2, input::Array{Float64,1})
• It Vector{T} is a type alias for Array{T,1} it is cleaner to read
  • So forward!(ann::ANN2, input::Vector{Float64})

• Don't over type specify. Input doesn't really have to be a Float64, it could be any kind of number. So instead use

  • So forward!{N<:Number}(ann::ANN2, input::Vector{N})
  • The most important thing this will allow you to do is Gradient Checking with ForwardDiff.jl
    • Gradient Checking your neural networks is very important as the back propagation algorithm is very hard, and fiddly, and it kinda works even if you have it wrong -- making the mistake hard to notice. See: here. People don't appreciate just how important doing a gradient check is.

• You do not need return statements at the end of functions, functions implictly return the last statement. So rather than return 1./(1+exp(-x)), you just write 1./(1+exp(-x)) as the last line of the sigmoid function

• Returning nothing is not a normal practice in Julia.

  • Functions which Mutate (modify) one of their inputs normally return the modified version of that input for fluid programming eg sort!(xs) will sort xs inplace and then return (the now sorted) xs.
  • Functions which don't modify their inputs, normally have some output to return (otherwise why were they called?). The exception is logging functions and that kinda thing. But they can just be allowed to return the return value of the last function they called. Which will likely be nothing

• use clearer names: As, Ns, Fs, Ss what are these? I don't know.

  • Ws and bs can be understood to be Weights and Baises as that is a very common notation -- but you don't use those they are named in full.
  • If you are using naming from a particular paper, you should add a link to that paper in a comment and have a comment signing what each name is for.
  • I can guess As is the activations of each layer.
  • on reading closer, looks like Fs[i] is the derivivitive of sigmoid(As[i]). Why are you calculating that doing feedforward? It is part of the back propagation step?
  • Ss is the error signal, I think
  • Ns is the input to each sigmoid
  • Working that out took me over 10 minutes
• Int keyed Dicts are a bit of a code smell. Particularly when they are consecutively indexed.
  • I see why you are using them -- you want zero indexes lists
  • I would rethink that, and if you still decide that it is the best way at least leave comment so your future self knows that are actually lists.

• It would be clearer (and might be marginally faster) if you gave those Dicts (by which I mean lists) a fixed return type rather than Any. They do have a constant return type (I think they are all Vector{Float64}).

  • This would also allow the type checker to catch some logic errors

• calcuate_deltas is spelt wrong. correct is calculate_deltas

• rather than having if statements inside a for-loop, checking if this is the first index and do almost entirely different code, you could instead just loop over the later indexes. So rather than:

  for i in reverse(1:this.numlayers)
      if i == this.numlayers
          this.Ss[i] = this.Fs[i].*(this.As[i] - target)
          if avg
              this.weightdelta[i] = this.weightdelta[i]+rate.*(this.Ss[i]*this.As[i-1]')
              this.biasdelta[i] = this.biasdelta[i]+rate.*this.Ss[i]
          else
              this.weightdelta[i] = rate.*(this.Ss[i]*this.As[i-1]')
              this.biasdelta[i] = rate.*this.Ss[i]
          end
      else
          this.Ss[i] = this.Fs[i].*(this.weights[i+1]'*this.Ss[i+1])
          if avg
              this.weightdelta[i] = this.weightdelta[i]+rate.*(this.Ss[i]*this.As[i-1]')
              this.biasdelta[i] = this.biasdelta[i]+rate.*this.Ss[i]
          else
              this.weightdelta[i] = rate.*(this.Ss[i]*this.As[i-1]')
              this.biasdelta[i] = rate.*this.Ss[i]
          end
      end
  end
  

  do

  this.Ss[this.numlayers] = this.Fs[this.numlayers].*(this.As[this.numlayers] - target)
  if avg
      this.weightdelta[this.numlayers] = this.weightdelta[this.numlayers]+rate.*(this.Ss[this.numlayers]*this.As[this.numlayers-1]')
      this.biasdelta[this.numlayers] = this.biasdelta[this.numlayers]+rate.*this.Ss[this.numlayers]
  else
      this.weightdelta[this.numlayers] = rate.*(this.Ss[this.numlayers]*this.As[this.numlayers-1]')
      this.biasdelta[this.numlayers] = rate.*this.Ss[this.numlayers]
  end
  
  for i in this.numlayers-1 :-1: 1
         this.Ss[i] = this.Fs[i].*(this.weights[i+1]'*this.Ss[i+1])
         if avg
              this.weightdelta[i] = this.weightdelta[i]+rate.*(this.Ss[i]*this.As[i-1]')
              this.biasdelta[i] = this.biasdelta[i]+rate.*this.Ss[i]
         else
              this.weightdelta[i] = rate.*(this.Ss[i]*this.As[i-1]')
              this.biasdelta[i] = rate.*this.Ss[i]
         end
  end
  

Neural Network/Architecture

• You are initialising your weights and biases with uniformly distributed values between 0.0 and 1.0. That is not going to go well for training. The spread is too large so the network will saturate. You also have no negative values.

  • Yann Lecun provides some solid and motivated guidelines in his very well known paper on Efficient BackProp. There is a lot more than just choising inital weights there.
  • as a very rough rule of thumb, rather than rand I suggest 0.1*randn so you are sampling from the normal distribution.

• I have not checked your backpropergation logic. A gradient check will do a better job than I could do anyway. You are looking for accuracy of your calculated gradient to match the automatic gradient to within <10^-15. Anything else is an indication of a mistake in your logic. (And I know how sucky that is, spent several days stuck on trying to match my gradient check)

• You say avg which is short for average but you mean momentum. That is what it is called in literature.

• The Neural Network Type should not be storing information directly tied to its transient state during execution. (particularly not for a nonrecoccurant network).

  • Mutating the state like that means you can be running multiple feedforwards at once. ie no parallelism.
  • As, Ns, Fs, Ss, weightdelta, biasdelta all have to go.
  • Restructure forward, calcuate_deltas and updateone/updateepoch into
    • forward: takes input as input still, returns As, Ns. Does not mutate any state variables.
    • get_gradient: Takes input and expected output as inputs, calls forward, calculates Fs, Ss and also weightdelta, biasdelta Does not mess with momentum (avg) or a learning rate. Does not mutate any state variables.
    • updateepoch, takes a list of training cases (inputs and expected outputs), as well as learning rate, and wither or not to do momentum. Performs a map (or a parallel map) over them through get_gradient, to have a list of gradients for all Weights/biases. It mutliplies that by a learning rate, and applies the changes to the weights/biases, optionally using the momentum to keep value. It could be extended to take a "minibatch" size, and then pocess one minibatch at a time. It could also be extended to do L1/L2 regularisation. Futhermore it could be replaced with a call to a strong nonlinear optimisation library like Optim. Though rearranging your data so thing works is a little fiddly. LeCunn talks a little about doing this in Effient BackProp.
    • updateone! simply an overload for updateepoch! that has a list of training cases contain just one.