title: Neural Net with softmax output failing to converge
body: <p>I've been working on Stamford's Deep Learning Tutorial and I'm having an issue with one of the exercises, the neural network with the softmax output layer. Here is my implementation in R:</p>
  
  <pre><code>train &lt;- function(training.set, labels, costFunc, activationFunc, outputActivationFunc, activationDerivative, hidden.unit.count = 7, learningRate = 0.3, decayRate=0.02, momentumRate=0.02, samples.count, batch.size, verbose=F, debug=F){
  
    #initialize weights and biases
    w1 &lt;- matrix( rnorm(hidden.unit.count * input.unit.count, sd=0.5), nrow=hidden.unit.count, ncol=input.unit.count)
    b1 &lt;- matrix(-1, nrow=hidden.unit.count, ncol=1)
    w2 &lt;- matrix(rnorm(output.unit.count * hidden.unit.count, sd=0.5), nrow=output.unit.count, ncol=hidden.unit.count)
    b2 &lt;- matrix(-1, nrow=output.unit.count, ncol=1)
  
    cost.list&lt;- matrix(rep(seq(1:floor(samples.count / batch.size)), each=2), byrow=T, ncol=2)
    cost.list[, 2] &lt;- 0
  
    i &lt;- 1
    while(i &lt; samples.count){
      z2 &lt;- w1 %*% training.set[, i: (i + batch.size - 1)] + matrix(rep(b1, each=batch.size), ncol=batch.size,byrow=T)
      a2 &lt;- activationFunc(z2)
  
      z3 &lt;- w2 %*% a2 + matrix(rep(b2, each=batch.size), ncol=batch.size,byrow=T)
      h  &lt;- outputActivationFunc(z3)
  
      #calculate error
      output.error &lt;- (h - labels[, i: (i + batch.size - 1)]) 
      hidden.error &lt;- (t(w2) %*% output.error) * sigmoidPrime(z2)
  
      # calculate gradients for both layers
      gradW2 &lt;- hidden.error %*% t(training.set[ ,i: (i + batch.size - 1)]) - momentumRate * gradW2.prev - decayRate * w1
      gradw2 &lt;- output.error %*% t(a2) - momentumRate * gradw2.prev - decayRate * w2
  
      gradW2.prev &lt;- gradW2
      gradw2.prev &lt;- gradw2
  
      #update weights and biases
      w1 &lt;- w1 - learningRate * gradW2 / batch.size
      w2 &lt;- w2 - learningRate * gradW3 / batch.size
  
      b1 &lt;- b1 - learningRate * rowSums(gradW2) / batch.size
      b2 &lt;- b2 - learningRate * rowSums(gradW3) / batch.size
  
      i &lt;- i + batch.size
    }
  
    return (list(w1, w2, b1, b2, cost.list))
  }
  </code></pre>
  
  <p>Here is the softmax function I use on the output layer and also the cost function I use with softmax:</p>
  
  <pre><code>softmax &lt;- function(a){
    a &lt;- a - apply(a, 1, function(row){ 
        return (max(row))
    })
  
    a &lt;- exp(a)
  
    return (sweep(a, 2, colSums(a), FUN='/'))
  }
  
  softmaxCost &lt;- function(w, b, x, y, decayRate, batch.size){
    a &lt;- w %*% x + matrix(rep(b, each=dim(x)[2]), byrow = T, ncol=dim(x)[2])
  
    h &lt;- softmax(a)
  
    cost &lt;- -1/batch.size * (sum(y * log(h))) + decayRate/2 * sum((w * w))
  
    return (cost)
  }
  </code></pre>
  
  <p>I've checked the gradients computed by my program against numerical gradients and they are different. However, I can't find the source of the incorrect gradient calculation.</p>
  
  <p>Also, I've successfully used this network using sigmoid activation at the output layer on the MNIST whereas the using the softmax layer simply doesn't work(11% accuracy). This leads me to believe that the issue lies in my softmax implementation.</p>