title: Extremely small or NaN values appear in training neural network
body: <p>I'm trying to implement a neural network architecture in Haskell, and use it on MNIST.</p>
  <p>I'm using the <code>hmatrix</code> package for linear algebra.
  My training framework is built using the <code>pipes</code> package.</p>
  <p>My code compiles and doesn't crash. But the problem is, certain combinations of layer size (say, 1000), minibatch size, and learning rate give rise to <code>NaN</code> values in the computations. After some inspection, I see that extremely small values (order of <code>1e-100</code>) eventually appear in the activations. But, even when that doesn't happen, the training still doesn't work. There's no improvement over its loss or accuracy.</p>
  <p>I checked and rechecked my code, and I'm at a loss as to what the root of the problem could be.</p>
  <p>Here's the backpropagation training, which computes the deltas for each layer:</p>
  <pre><code>backward lf n (out,tar) das = do
      let δout = tr (derivate lf (tar, out)) -- dE/dy
          deltas = scanr (\(l, a') δ -&gt;
                           let w = weights l
                           in (tr a') * (w &lt;&gt; δ)) δout (zip (tail $ toList n) das)
      return (deltas)
  </code></pre>
  <p><code>lf</code> is the loss function, <code>n</code> is the network (<code>weight</code> matrix and <code>bias</code> vector for each layer), <code>out</code> and <code>tar</code> are the actual output of the network and the <code>target</code> (desired) output, and <code>das</code> are the activation derivatives of each layer.</p>
  <p>In batch mode, <code>out</code>, <code>tar</code> are matrices (rows are output vectors), and <code>das</code> is a list of the matrices.</p>
  <p>Here's the actual gradient computation:</p>
  <pre><code>  grad lf (n, (i,t)) = do
      -- Forward propagation: compute layers outputs and activation derivatives
      let (as, as') = unzip $ runLayers n i
          (out) = last as
      (ds) &lt;- backward lf n (out, t) (init as') -- Compute deltas with backpropagation
      let r  = fromIntegral $ rows i -- Size of minibatch
      let gs = zipWith (\δ a -&gt; tr (δ &lt;&gt; a)) ds (i:init as) -- Gradients for weights
      return $ GradBatch ((recip r .*) &lt;$&gt; gs, (recip r .*) &lt;$&gt; squeeze &lt;$&gt; ds)
  </code></pre>
  <p>Here, <code>lf</code> and <code>n</code> are the same as above, <code>i</code> is the input, and <code>t</code> is the target output (both in batch form, as matrices).</p>
  <p><code>squeeze</code> transforms a matrix into a vector by summing over each row. That is, <code>ds</code> is a list of matrices of deltas, where each column corresponds to the deltas for a row of the minibatch. So, the gradients for the biases are the average of the deltas over all the minibatch. The same thing for <code>gs</code>, which corresponds to the gradients for the weights.</p>
  <p>Here's the actual update code:</p>
  <pre><code>move lr (n, (i,t)) (GradBatch (gs, ds)) = do
      -- Update function
      let update = (\(FC w b af) g δ -&gt; FC (w + (lr).*g) (b + (lr).*δ) af)
          n' = Network.fromList $ zipWith3 update (Network.toList n) gs ds
      return (n', (i,t))
  </code></pre>
  <p><code>lr</code> is the learning rate. <code>FC</code> is the layer constructor, and <code>af</code> is the activation function for that layer.</p>
  <p>The gradient descent algorithm makes sure to pass in a negative value for the learning rate. The actual code for the gradient descent is simply a loop around a composition of <code>grad</code> and <code>move</code>, with a parameterized stop condition.</p>
  <p>Finally, here's the code for a mean square error loss function:</p>
  <pre><code>mse :: (Floating a) =&gt; LossFunction a a
  mse = let f (y,y') = let gamma = y'-y in gamma**2 / 2
            f' (y,y') = (y'-y)
        in  Evaluator f f'
  </code></pre>
  <p><code>Evaluator</code> just bundles a loss function and its derivative (for calculating the delta of the output layer).</p>
  <p>The rest of the code is up on GitHub: <a href="https://github.com/DrPyser/NeuralNetwork" rel="noreferrer">NeuralNetwork</a>.</p>
  <p>So, if anyone has an insight into the problem, or even just a sanity check that I'm correctly implementing the algorithm, I'd be grateful.</p>