Solutions to exercises and problems from chapter 2 of Neural networks and deep learning

Book by Michael Nielsen

Solutions by Peter Bessman

Problem 1: Alternate presentation of the equations of backpropagation

Part 1

Given \(\boldsymbol{\delta}^L = \nabla_{\mathbf{a}} C \odot \sigma'(\mathbf{z}^L)\), show that \(\boldsymbol{\delta}^L = \mathbf{\Sigma}'(\mathbf{z}^L)\nabla_{\mathbf{a}} C\).

Note that Nielsen uses column vectors for gradients, so \(\nabla_{\mathbf{a}} C\) is an \(N \times 1\) column vector:

\[\nabla_{\mathbf{a}} C = \left[ \frac{\partial C}{\partial a^L_1}, \frac{\partial C}{\partial a^L_2}, \dots, \frac{\partial C}{\partial a^L_N} \right]^T\]

Therefore \(\sigma'(\mathbf{z}^L)\) must also be an \(N \times 1\) column vector.

Note also that Nielsen uses \(\mathbf{\Sigma}'(\mathbf{z}^L)\) to represent an \(N \times N\) diagonal matrix populated by the elements of \(\sigma'(\mathbf{z}^L)\):

\[\mathbf{\Sigma}'(\mathbf{z}^L) = \text{diag}(\sigma'(\mathbf{z}^L)) = \begin{bmatrix} \sigma'(z^L_1) & & 0 \\ & \ddots & \\ 0 & & \sigma'(z^L_N) \end{bmatrix}\]

The Hadamard product denoted by \(\odot\) executes element-wise multiplication between two vectors of identical dimensions.

We will achieve our proof by expanding the two forms of \(\boldsymbol{\delta}^L\) and demonstrating equivalence:

\[ \begin{align*} \\ \nabla_{\mathbf{a}} C \odot \sigma'(\mathbf{z}^L) &= \begin{bmatrix} \frac{\partial C}{\partial a^L_1} \\ \vdots \\ \frac{\partial C}{\partial a^L_N} \end{bmatrix} \odot \begin{bmatrix} \sigma'(z^L_1) \\ \vdots \\ \sigma'(z^L_N) \end{bmatrix} &= \begin{bmatrix} \frac{\partial C}{\partial a^L_1} \sigma'(z^L_1) \\ \vdots \\ \frac{\partial C}{\partial a^L_N} \sigma'(z^L_N) \end{bmatrix} \\ \\ \mathbf{\Sigma}'(\mathbf{z}^L)\nabla_{\mathbf{a}} C &= \begin{bmatrix} \sigma'(z^L_1) & & 0 \\ & \ddots & \\ 0 & & \sigma'(z^L_N) \end{bmatrix} \begin{bmatrix} \frac{\partial C}{\partial a^L_1} \\ \vdots \\ \frac{\partial C}{\partial a^L_N} \end{bmatrix} &= \begin{bmatrix} \frac{\partial C}{\partial a^L_1} \sigma'(z^L_1) \\ \vdots \\ \frac{\partial C}{\partial a^L_N} \sigma'(z^L_N) \end{bmatrix} \\ \end{align*} \]

Part 2

As with part 1, given that \(\boldsymbol{\delta}^l = ((\mathbf{W}^{l+1})^T \boldsymbol{\delta}^{l+1}) \odot \sigma'(\mathbf{z}^l)\) we must show that \(\boldsymbol{\delta}^l = \mathbf{\Sigma}'(\mathbf{z}^l) ((\mathbf{W}^{l+1})^T \boldsymbol{\delta}^{l+1})\).

We begin by expanding the given equation for \(\boldsymbol{\delta}^l\); assuming \(\mathbf{W}^{l+1}\) is an \(N \times M\) matrix this results in:

\[((\mathbf{W}^{l+1})^T\boldsymbol{\delta}^{l+1}) \odot \sigma'(\mathbf{z}^l) = \begin{bmatrix} ((\mathbf{W}^{l+1})^T)_1 \boldsymbol{\delta}^{l+1} \\ \vdots \\ ((\mathbf{W}^{l+1})^T)_M \boldsymbol{\delta}^{l+1} \end{bmatrix} \odot \begin{bmatrix} \sigma'(z^l_1) \\ \vdots \\ \sigma'(z^l_M) \end{bmatrix} = \begin{bmatrix} ((\mathbf{W}^{l+1})^T)_1 \boldsymbol{\delta}^{l+1} \sigma'(z^l_1) \\ \vdots \\ ((\mathbf{W}^{l+1})^T)_M \boldsymbol{\delta}^{l+1} \sigma'(z^l_M) \end{bmatrix}\]

Factoring out the scalar \(\sigma'(z^l_j)\) components into an \(M \times M\) diagonal matrix demonstrates equivalence.

\[\boldsymbol{\delta}^l = \begin{bmatrix} \sigma'(z^l_1) & & 0 \\ & \ddots & \\ 0 & & \sigma'(z^l_M) \end{bmatrix} \begin{bmatrix} ((\mathbf{W}^{l+1})^T)_1 \boldsymbol{\delta}^{l+1} \\ \vdots \\ ((\mathbf{W}^{l+1})^T)_M \boldsymbol{\delta}^{l+1} \end{bmatrix} = \mathbf{\Sigma}'(\mathbf{z}^l) ((\mathbf{W}^{l+1})^T \boldsymbol{\delta}^{l+1})\]

Part 3

We are tasked with showing that:

\[ \begin{align*} \boldsymbol{\delta}^l = \mathbf{\Sigma}'(\mathbf{z}^l) (\mathbf{W}^{l+1})^T \ldots \mathbf{\Sigma}'(\mathbf{z}^{L-1}) (\mathbf{W}^L)^T \mathbf{\Sigma}'(\mathbf{z}^L) \nabla_{\mathbf{a}} C \tag{PART3} \end{align*} \]

This is simply achieved by unrolling the last two equations we proved:

\[ \begin{align*} \boldsymbol{\delta}^L &= \mathbf{\Sigma}'(\mathbf{z}^L)\nabla_{\mathbf{a}} C \tag{PART1} \\ \boldsymbol{\delta}^l &= \mathbf{\Sigma}'(\mathbf{z}^l) ((\mathbf{W}^{l+1})^T \boldsymbol{\delta}^{l+1}) \tag{PART2} \end{align*} \]

We already have a base case from \(\text{(PART1)}\). Next we can use the recurrence relation from \(\text{(PART2)}\) to evaluate one step back at \(l = L - 1\):

\[ \begin{align*} \boldsymbol{\delta}^{L-1} &= \mathbf{\Sigma}'(\mathbf{z}^{L-1}) ((\mathbf{W}^{(L-1)+1})^T \boldsymbol{\delta}^{(L-1)+1}) \\ &= \mathbf{\Sigma}'(\mathbf{z}^{L-1}) ((\mathbf{W}^{L})^T \boldsymbol{\delta}^{L}) \\ \therefore \boldsymbol{\delta}^{L-1} &= \mathbf{\Sigma}'(\mathbf{z}^{L-1}) (\mathbf{W}^{L})^T \mathbf{\Sigma}'(\mathbf{z}^L)\nabla_{\mathbf{a}} C \end{align*} \]

Taking two steps back to \(l = L - 2\) we obtain a general form:

\[ \begin{align*} \boldsymbol{\delta}^{L-2} &= \mathbf{\Sigma}'(\mathbf{z}^{L-2}) ((\mathbf{W}^{(L-2)+1})^T \boldsymbol{\delta}^{(L-2)+1}) \\ &= \mathbf{\Sigma}'(\mathbf{z}^{L-2}) ((\mathbf{W}^{L-1})^T \boldsymbol{\delta}^{L-1}) \\ &= \mathbf{\Sigma}'(\mathbf{z}^{L-2}) (\mathbf{W}^{L-1})^T \mathbf{\Sigma}'(\mathbf{z}^{L-1}) (\mathbf{W}^{L})^T \mathbf{\Sigma}'(\mathbf{z}^L)\nabla_{\mathbf{a}} C \\ \therefore \boldsymbol{\delta}^{L-n} &= \mathbf{\Sigma}'(\mathbf{z}^{L-n}) (\mathbf{W}^{L-n+1})^T \ldots \mathbf{\Sigma}'(\mathbf{z}^{L-1}) (\mathbf{W}^{L})^T \mathbf{\Sigma}'(\mathbf{z}^L)\nabla_{\mathbf{a}} C \end{align*} \]

If we let \(l = L - n\) we arrive at the form in \(\text{(PART3)}\) above:

\[ \boldsymbol{\delta}^{L-n} = \boldsymbol{\delta}^{l} = \mathbf{\Sigma}'(\mathbf{z}^{l}) (\mathbf{W}^{l+1})^T \ldots \mathbf{\Sigma}'(\mathbf{z}^{L-1}) (\mathbf{W}^{L})^T \mathbf{\Sigma}'(\mathbf{z}^L)\nabla_{\mathbf{a}} C \]

Exercise 1: Proving the backpropagation equations

We will prove the following four equations:

\[ \begin{align*} \\ \boldsymbol{\delta}^L &= \nabla_{\mathbf{a}} C \odot \sigma'(\mathbf{z}^L) \tag{BP1} \\ \\ \boldsymbol{\delta}^l &= ((\mathbf{W}^{l+1})^T \boldsymbol{\delta}^{l+1}) \odot \sigma'(\mathbf{z}^l) \tag{BP2} \\ \\ \frac{\partial C}{\partial b^l_j} &= \delta^l_j \tag{BP3} \\ \\ \frac{\partial C}{\partial w^l_{jk}} &= a^{l-1}_k \delta^l_j \tag{BP4} \\ \end{align*} \]

Throughout, we will make use of the following definitions provided in the text:

\[ \begin{align*} \\ \delta^l_j &\equiv \frac{\partial C}{\partial z^l_j} \tag{ERR} \\ \\ \mathbf{z}^l &\equiv \mathbf{W}^l \mathbf{a}^{l-1} + \mathbf{b}^l \tag{INPUT} \\ \\ \mathbf{a}^l &\equiv \sigma(\mathbf{z}^l) \tag{ACT} \\ \end{align*} \]

Nielsen proves the first two, but I offer my own take below in addition to proofs for the other equations.

And while we don’t care about the definiton of the cost function \(C\), we do need to consider its key characteristics:

  1. It can be written as an average \(C = \frac{1}{n} \sum_x C_x\) over cost functions \(C_x\) for individual training examples, \(x\).
  2. it can be written as a function of the output activations \(\mathbf{a}^L\).

Proving BP1

To simplify the proof, let’s first decompose \(\text{(BP1)}\) into its \(j\text{-th}\) component function:

\[ \begin{align*} \\ \boldsymbol{\delta}^L &= \nabla_{\mathbf{a}} C \odot \sigma'(\mathbf{z}^L) \tag{BP1} \\ \\ &= \left[ \frac{\partial C}{\partial a^L_1}, \frac{\partial C}{\partial a^L_2}, \dots \right]^T \odot \left[ \sigma'(z^L_1), \sigma'(z^L_2), \dots \right]^T \\ \\ \therefore \delta^L_j &= \frac{\partial C}{\partial a^L_j} \sigma'(z^L_j) \tag{BP1'} \\ \end{align*} \]

We already have a definition for \(\delta^L_j\) from \(\text{(ERR)}\), and if we expand that with the chain rule we see that it looks very similar to \(\text{(BP1')}\):

\[ \begin{align*} \\ \delta^L_j \equiv \frac{\partial C}{\partial z^L_j} = \sum_{j'} \frac{\partial C}{\partial a^L_{j'}} \frac{\partial a^L_{j'}}{\partial z^L_j} \tag{ERR2} \\ \end{align*} \]

Hence our strategy will be to show that \(\text{(ERR)}\) is equivalent to \(\text{(BP1')}\) when \(l = L\), which by extension would prove that \(\text{(BP1)}\) is correct.

Note: \({C}\) is defined as a function of \({\mathbf{a}^L}\), which in turn is a function of \({\mathbf{z}^L}\). This is why we can state \(\delta^L_j\) in terms of \({\mathbf{a}^L}\) using the chain rule.

Also, we use \({\sum_{j'}}\) because a change in any \({a^L_{j'}}\) can potentially contribute to a change in \({C}\); so by the law of total derivatives we must sum up each partial derivative of \({C}\) with respect to \({a^L_{j'}}\) to determine the total contribution of changes in \({\mathbf{a}^L}\) to \({C}\). For more information, see The Matrix Calculus You Need For Deep Learning by Terrence Parr and Jeremy Howard.

Finally, we use \({j'}\) to indicate that although \({a^L_{j'}}\) represents the same kind of subscript as \({z^L_{j}}\), they are nevertheless different variables, with \({j'}\) being an index of summation while \({j}\) is a variable of the \({\delta^L_j}\) identity.

We first note that we can eliminate the sum operator from \(\text{(ERR2)}\). To see why, first decompose \(\text{(ACT)}\) down to \(a^l_{j} = \sigma(z^l_{j})\), then evaluate \(\partial a^L_{j'} / \partial z^L_{j}\):

\[ \begin{align*} \\ \frac{\partial a^L_{j'}}{\partial z^L_{j}} = \frac{\partial}{\partial z^L_{j}} \sigma(z^L_{j'}) = \begin{cases} \sigma'(z^L_j) & \text{if } j' = j \\ 0 & \text{if } j' \neq j \end{cases} \\ \end{align*} \]

Now expand \(\text{(ERR2)}\) and evaluate:

\[ \begin{align*} \\ \delta^L_j &= \sum_{j'} \frac{\partial C}{\partial a^L_{j'}} \frac{\partial a^L_{j'}}{\partial z^L_j} \tag{ERR2} \\ \\ &= \frac{\partial C}{\partial a^L_{1}} \frac{\partial a^L_{1}}{\partial z^L_j} + \frac{\partial C}{\partial a^L_{2}} \frac{\partial a^L_{2}}{\partial z^L_j} + \dots + \frac{\partial C}{\partial a^L_{j}} \frac{\partial a^L_{j}}{\partial z^L_j} + \dots \\ \\ &= \frac{\partial C}{\partial a^L_{1}} (0) + \frac{\partial C}{\partial a^L_{2}} (0) + \dots + \frac{\partial C}{\partial a^L_{j}} \sigma'(z^L_j) + \dots \\ \\ \delta^L_j &= \frac{\partial C}{\partial a^L_{j}} \sigma'(z^L_j) \tag{ERR3} \\ \end{align*} \]

This is equivalent to \(\text{(BP1')}\) that we worked out previously, so by extension \(\text{(BP1)}\) is proved correct.

Proving BP2

As with our proof for \(\text{(BP1)}\), we first decompose \(\text{(BP2)}\) into its \(j\text{-th}\) component function:

\[ \begin{align*} \boldsymbol{\delta}^l &= \left( (\mathbf{W}^{l+1})^\top \boldsymbol{\delta}^{l+1} \right) \odot \sigma'(\mathbf{z}^l) \tag{BP2} \\ &= \left[ [(\mathbf{W}^{l+1})^\top]_1 \boldsymbol{\delta}^{l+1}, [(\mathbf{W}^{l+1})^\top]_2 \boldsymbol{\delta}^{l+1}, \dots \right]^\top \odot \left[ \sigma'(z^l_1), \sigma'(z^l_2), \dots \right]^\top \\ \therefore \delta^l_j &= [(\mathbf{W}^{l+1})^\top]_j \boldsymbol{\delta}^{l+1} \sigma'(z^l_j) \tag{BP2'} \end{align*} \]

We worked out previously that \(\partial a^L_j / \partial z^L_j = \sigma'(z^L_j)\), which implies that \(\sigma'(z^l_j) = \partial a^l_j / \partial z^l_j\). Note also that we can expand \([(\mathbf{W}^{l+1})^\top]_j\) and \(\boldsymbol{\delta}^{l+1}\):

\[ \begin{align*} [(\mathbf{W}^{l+1})^\top]_j &= \left[ w^{l+1}_{1j}, w^{l+1}_{2j}, \dots \right] \\ \boldsymbol{\delta}^{l+1} &= \left[ \delta^{l+1}_1, \delta^{l+1}_2, \dots \right]^\top \\ \end{align*} \]

Therefore \(\text{(BP2')}\) becomes:

\[ \begin{align*} \delta^l_j &= [(\mathbf{W}^{l+1})^\top]_j \boldsymbol{\delta}^{l+1} \frac{\partial a^l_j}{\partial z^l_j} \\ &= \left[ w^{l+1}_{1j}, w^{l+1}_{2j}, \dots \right] \left[ \delta^{l+1}_1, \delta^{l+1}_2, \dots \right]^\top \frac{\partial a^l_j}{\partial z^l_j} \\ &= \left( w^{l+1}_{1j}\delta^{l+1}_1 + w^{l+1}_{2j}\delta^{l+1}_2 + \dots \right) \frac{\partial a^l_j}{\partial z^l_j} \\ \delta^l_j &= \sum_k w^{l+1}_{kj}\delta^{l+1}_k \frac{\partial a^l_j}{\partial z^l_j} \tag{BP2''} \end{align*} \]

We also previously established that:

\[ \begin{align*} \delta^L_j &= \frac{\partial C}{\partial a^L_j} \sigma'(z^L_j) \\ \therefore \delta^l_j &= \frac{\partial C}{\partial a^l_j} \sigma'(z^l_j) = \frac{\partial C}{\partial a^l_j} \frac{\partial a^l_j}{\partial z^l_j} \end{align*} \]

This implies that if \(\text{(BP2'')}\) is true then \(\sum_k w^{l+1}_{kj}\delta^{l+1}_k\) must resolve to \(\partial C / \partial a^l_j\):

\[ \begin{align*} \sum_k w^{l+1}_{kj}\delta^{l+1}_k = \frac{\partial C}{\partial a^l_j} \tag{PROP1} \end{align*} \]

If we prove \(\text{PROP1}\), we will have proven \(\text{(BP2'')}\) and therefore \(\text{(BP2')}\) and \(\text{(BP2)}\). Its form suggests using the chain rule again, so it may be useful to expand the \(\delta^{l+1}_k\) term:

\[ \begin{align*} \sum_k w^{l+1}_{kj}\delta^{l+1}_k = \sum_k w^{l+1}_{kj} \frac{\partial C}{\partial z^{l+1}_k} \tag{PROP1'} \end{align*} \]

This gives us a very strong hint about how the sum term must resolve:

\[ \begin{align*} \sum_k w^{l+1}_{kj} \frac{\partial C}{\partial z^{l+1}_k} &= \frac{\partial C}{\partial a^l_j} \\ \therefore w^{l+1}_{kj} &\stackrel{?}\rightsquigarrow \frac{\partial z^{l+1}_k}{\partial a^l_j} \tag{PROP2} \end{align*} \]

So let’s investigate \(\partial z^{l+1}_k / \partial a^l_j\). First recall the definition of \(\text{(INPUT)}\), and use that to express its \(k\text{-th}\) component function at layer \(l+1\):

\[ \begin{align*} \mathbf{z}^l &= \mathbf{W}^l \mathbf{a}^{l-1} + \mathbf{b}^l \tag{INPUT} \\ \therefore z^{l+1}_k &= \mathbf{w}^{l+1}_k \mathbf{a}^l + b^{l+1}_k \\ &= \sum_{j'} w^{l+1}_{kj'} a^l_{j'} + b^{l+1}_k \\ \end{align*} \]

Note: \(j'\) relates to \(j\) as in the proof of \(\text{(BP1)}\), see the note there for details.

We can use this to expand \(\partial z^{l+1}_k / \partial a^l_j\) and evaluate it:

\[ \begin{align*} \frac{\partial z^{l+1}_k}{\partial a^l_j} &= \frac{\partial}{\partial a^l_j} \left( \sum_{j'} w^{l+1}_{kj'} a^l_{j'} + b^{l+1}_k \right) \\ &= \frac{\partial}{\partial a^l_j} \left( \left( w^{l+1}_{k1} a^l_1 + w^{l+1}_{k2} a^l_2 + \dots + w^{l+1}_{kj} a^l_j + \dots \right) + b^{l+1}_k \right) \\ &= \left( \frac{\partial}{\partial a^l_j} (w^{l+1}_{k1} a^l_1) + \frac{\partial}{\partial a^l_j} (w^{l+1}_{k2} a^l_2) + \dots + \frac{\partial}{\partial a^l_j} (w^{l+1}_{kj} a^l_j) + \dots \right) + \frac{\partial}{\partial a^l_j} (b^{l+1}_k) \\ &= \left( 0 + 0 + \dots + w^{l+1}_{kj}(1) + \dots \right) + 0 \\ \therefore \frac{\partial z^{l+1}_k}{\partial a^l_j} &= w^{l+1}_{kj} \quad \checkmark \tag{PROP2} \end{align*} \]

Having verified \(\text{PROP2}\), we can combine it with \(\text{PROP1'}\) to express \(\text{PROP1}\) entirely as a sum of chained partial derivatives: \[ \begin{align*} \sum_k w^{l+1}_{kj}\delta^{l+1}_k &= \frac{\partial C}{\partial a^l_j} \tag{PROP1} \\ \sum_k w^{l+1}_{kj}\delta^{l+1}_k &= \sum_k w^{l+1}_{kj} \frac{\partial C}{\partial z^{l+1}_k} \tag{PROP1'} \\ w^{l+1}_{kj} &= \frac{\partial z^{l+1}_k}{\partial a^l_j} \tag{PROP2} \\ \therefore \sum_k w^{l+1}_{kj}\delta^{l+1}_k &= \sum_k \frac{\partial z^{l+1}_k}{\partial a^l_j} \frac{\partial C}{\partial z^{l+1}_k} = \frac{\partial C}{\partial a^l_j} \tag{PROP1''} \end{align*} \]

We can see that \(\text{(PROP1'')}\) is trivially correct as it’s simply an expansion of \(\partial C / \partial a^l_j\) using the definition of the total derivative. This proves \(\text{(PROP1')}\) and \(\text{(PROP1)}\), and therefore \(\text{(BP2'')}\) , \(\text{(BP2')}\) and ultimately \(\text{(BP2)}\)

Proving BP3

For this proof it will once again be helpful to recall the definition for \(\text{(ERR)}\):

\[ \begin{align*} \frac{\partial C}{\partial b_j^l} &= \delta_j^l \tag{BP3} \\ \delta_j^l &\equiv \frac{\partial C}{\partial z_j^l} \tag{ERR} \end{align*} \]

This suggests a proof through the now-familiar total derivative:

\[ \begin{align*} \frac{\partial C}{\partial b_j^l} = \sum_{j'} \frac{\partial C}{\partial z_{j'}^l} \frac{\partial z_{j'}^l}{\partial b_j^l} \stackrel{?}\rightsquigarrow \frac{\partial C}{\partial z_j^l} \tag{PROP3} \end{align*} \]

Note: \(j'\) relates to \(j\) as in the proof of \(\text{(BP1)}\), see the note there for details.

If we prove \(\text{(PROP3)}\) we will have proven \(\text{(BP3)}\). We begin by using \(\text{(INPUT)}\) to derive \(z_{j'}^l\) from \(\mathbf{z}\):

\[ \begin{align*} \mathbf{z}^l &= \mathbf{W}^l \mathbf{a}^{l-1} + \mathbf{b}^l \tag{INPUT} \\ \therefore z_{j'}^l &= \mathbf{w}^{l}_{j'} \mathbf{a}^{l-1} + b^{l}_{j'} \\ &= \sum_{k} (w_{j'k}^l a_k^{l-1}) + b_{j'}^l \\ \end{align*} \]

We can then substitute this into \(\partial z_{j'}^l / \partial b_j^l\):

\[\frac{\partial z_{j'}^l}{\partial b_j^l} = \frac{\partial}{\partial b_j^l} \left( \sum_{k} (w_{j'k}^l a_k^{l-1}) + b_{j'}^l \right) \]

This will evaluate to \({0}\) except when \({j' = j}\), in which case our result is simply \({1}\). We can use this to expand \(\text{(PROP3)}\) and prove it correct, thereby proving \(\text{(BP3)}\):

\[ \begin{align*} \\ \frac{\partial C}{\partial b_j^l} &= \sum_{j'} \frac{\partial C}{\partial z_{j'}^l} \frac{\partial z_{j'}^l}{\partial b_j^l} \stackrel{?}\rightsquigarrow \frac{\partial C}{\partial z_j^l} \tag{PROP3} \\ &= \frac{\partial C}{\partial z_1^l} \frac{\partial z_1^l}{\partial b_j^l} + \frac{\partial C}{\partial z_2^l} \frac{\partial z_2^l}{\partial b_j^l} + \dots + \frac{\partial C}{\partial z_j^l} \frac{\partial z_j^l}{\partial b_j^l} + \dots \\ &= \frac{\partial C}{\partial z_1^l} (0) + \frac{\partial C}{\partial z_2^l} (0) + \dots + \frac{\partial C}{\partial z_j^l} (1) + \dots \\ \therefore \frac{\partial C}{\partial b_j^l} &= \frac{\partial C}{\partial z_j^l} (1) = \frac{\partial C}{\partial z_j^l} \quad \checkmark \\ \end{align*} \]

Proving BP4

The setup for this proof is very similar to the proof for \(\text{BP3}\) in that we will try to prove equivalence with the definition \(\text{(ERR)}\) using the total derivative:

\[ \begin{align*} \\ \frac{\partial C}{\partial w_{jk}^l} &= a_k^{l-1} \delta_j^l \tag{BP4} \\ \delta_j^l &\equiv \frac{\partial C}{\partial z_j^l} \tag{ERR} \\ \frac{\partial C}{\partial w_{jk}^l} &= \sum_{j'} \frac{\partial C}{\partial z_{j'}^l} \frac{\partial z_{j'}^l}{\partial w_{jk}^l} \stackrel{?}\rightsquigarrow a_k^{l-1} \frac{\partial C}{\partial z_j^l} \tag{PROP4} \end{align*} \]

Note: \(j'\) relates to \(j\) as in the proof of \(\text{(BP1)}\), see the note there for details.

If we can prove \(\text{(PROP4)}\), we will have proven \(\text{(BP4)}\). Once again, we begin by using \(\text{(INPUT)}\) to derive \(z_{j'}^l\) from \(\mathbf{z}\):

\[ \begin{align*} \mathbf{z}^l &= \mathbf{W}^l \mathbf{a}^{l-1} + \mathbf{b}^l \tag{INPUT} \\ \therefore z_{j'}^l &= \mathbf{w}^{l}_{j'} \mathbf{a}^{l-1} + b^{l}_{j'} \\ &= \sum_{k'} w_{j'k'}^l a_{k'}^{l-1} + b_{j'}^l \\ \end{align*} \]

Note: \(k'\) relates to \(k\) in the same way that \(j'\) relates to \(j\). We have to do this trick for not just \(j\) but also \(k\) because they are both named parameters of \(\text{(BP4)}\).

We can substitute this into \(\text{(PROP4)}\) to expand it and — with a bit of patience and brute force — prove it correct, thereby proving \(\text{(BP4)}\):

\[ \begin{align*} \\ \frac{\partial C}{\partial w_{jk}^l} &= \sum_{j'} \frac{\partial C}{\partial z_{j'}^l} \frac{\partial}{\partial w_{jk}^l} \left( \left( \sum_{k'} w_{j'k'}^l a_{k'}^{l-1} \right) + b_{j'}^l \right) \stackrel{?}\rightsquigarrow a_k^{l-1} \frac{\partial C}{\partial z_j^l} \tag{PROP4} \\ \\ &= \sum_{j'} \frac{\partial C}{\partial z_{j'}^l} \frac{\partial}{\partial w_{jk}^l} \left( (w_{j'1}^l a_1^{l-1} + w_{j'2}^l a_2^{l-1} + \dots + w_{j'k}^l a_k^{l-1} + \dots) + b_{j'}^l \right) \\ \\ &= \sum_{j'} \frac{\partial C}{\partial z_{j'}^l} \left(\left( 0 + 0 + \dots + \frac{\partial}{\partial w_{jk}^l} \left( w_{j'k}^l a_k^{l-1} \right) + \dots \right) + 0 \right) \\ \\ &= \sum_{j'} \frac{\partial C}{\partial z_{j'}^l} \frac{\partial}{\partial w_{jk}^l} \left( w_{j'k}^l a_k^{l-1} \right) \\ \\ &= \frac{\partial C}{\partial z_{1}^l} \frac{\partial}{\partial w_{jk}^l} \left( w_{1k}^l a_k^{l-1} \right) + \frac{\partial C}{\partial z_{2}^l} \frac{\partial}{\partial w_{jk}^l} \left( w_{2k}^l a_k^{l-1} \right) + \dots + \frac{\partial C}{\partial z_{j}^l} \frac{\partial}{\partial w_{jk}^l} \left( w_{jk}^l a_k^{l-1} \right) + \dots \\ \\ &= \frac{\partial C}{\partial z_{1}^l} (0) + \frac{\partial C}{\partial z_{2}^l} (0) + \dots + \frac{\partial C}{\partial z_{j}^l} \left( (1) a_k^{l-1} \right) + \dots \\ \\ \frac{\partial C}{\partial w_{jk}^l} &= \frac{\partial C}{\partial z_{j}^l} a_k^{l-1} \quad \checkmark \\ \\ \end{align*} \]

Exercise 2: Backpropagation with a single modified neuron

Suppose we modify a single neuron in a feedforward network so that the output from the neuron is given by \(f(\sum_j w_j x_j + b)\) , where \(f\) is some function other than the sigmoid. How should we modify the backpropagation algorithm in this case?

The general point is that backpropagation does not much depend on the nature of the activation function, which has been the sigmoid in the discussion thus far. The algorithm will “work” for any continuously differentiable function (though it may become less effective). Thus replacing the sigmoid with \(f\) only requires a simple modification to the way we calculate a neuron’s activation in the feedforward pass, as well as the first two fundamental backpropagation equations used in the backward pass:

\[ \begin{align*} \mathbf{a}^{l} &= \sigma(\mathbf{z}^{l}) &&\longrightarrow f(\mathbf{z}^{l}) \\ \boldsymbol{\delta}^L &= \nabla_{\mathbf{a}} C \odot \sigma'(\mathbf{z}^L) &&\longrightarrow \nabla_{\mathbf{a}} C \odot f'(\mathbf{z}^L) \tag{BP1} \\ \boldsymbol{\delta}^l &= ((\mathbf{W}^{l+1})^T \boldsymbol{\delta}^{l+1}) \odot \sigma'(\mathbf{z}^l) &&\longrightarrow ((\mathbf{W}^{l+1})^T \boldsymbol{\delta}^{l+1}) \odot f'(\mathbf{z}^l) \tag{BP2} \\ \end{align*} \]

The specific point is that using a different activation function for a single neuron will only affect the output of that single neuron along the lines described above:

\[ \begin{align*} a^{l}_j &= \sigma(z^{l}_j) &&\longrightarrow f(z^{l}_j) \\ \delta^L_j &= \frac{\partial C}{\partial a^L_j} \sigma'(z^L_j) &&\longrightarrow \frac{\partial C}{\partial a^L_j} f'(z^L_j) \tag{BP1} \\ \delta^l_j &= ((\mathbf{W}^{l+1})^T \boldsymbol{\delta}^{l+1})_j \sigma'(z^l_j) &&\longrightarrow ((\mathbf{W}^{l+1})^T \boldsymbol{\delta}^{l+1})_j f'(z^l_j) \tag{BP2} \\ \end{align*} \]

Exercise 3: Backpropagation with linear neurons

Suppose we replace the usual non-linear \(\sigma\) function with \(\sigma(z) = z\) throughout the network. Rewrite the backpropagation algorithm for this case.

As per the above answer, this is a matter of straightforward substitution, made even simpler by the fact that \(\frac{d}{dz}(z) = 1\):

  1. Input x: Set the corresponding activation \(\mathbf{a}^{1}\) for the input layer.

  2. Feedforward: For each \(l = 2, 3, \ldots, L\) compute \(\mathbf{z}^{l} = \mathbf{W}^l \mathbf{a}^{l-1}+\mathbf{b}^l\) and \(\mathbf{a}^{l} = \mathbf{z}^{l}\).

  3. Output error \(\boldsymbol{\delta}^L\): Compute the vector \(\boldsymbol{\delta}^{L} = \nabla_{\mathbf{a}} C\).

  4. Backpropagate the error: For each \(l = L-1, L-2, \ldots, 2\) compute \(\boldsymbol{\delta}^{l} = ((\mathbf{W}^{l+1})^T \boldsymbol{\delta}^{l+1})\).

  5. Output: The gradient of the cost function is given by \(\frac{\partial C}{\partial w^l_{jk}} = a^{l-1}_k \delta^l_j\) and \(\frac{\partial C}{\partial b^l_j} = \delta^l_j\).

Problem 2: Fully matrix-based approach to backpropagation over a mini-batch

Theory

We will follow the outline that Nielsen provides for generic mini-batch stochastic gradient descent, along the way noting the changes we make to adapt this procedure to a fully matrix-based approach.

Note: After spending so much time reading and writing equations where \(\mathbf{a}^x\) means “the vector \(\mathbf{a}\) at layer \(x\),” I find it very jarring to have it now mean “the vector \(\mathbf{a}\) with mini-batch index \(x\).” So while Nielsen moves the layer index into the outer superscript, I keep it at the inner superscript and use the outer superscript for the mini-batch index.

Likewise, Nielsen defines layer \(l = 1\) as the input data or “pseudo-neurons,” with \(l = 2\) being the first layer of “proper” neurons. This feels odd to me, and while I could work in a codebase like this I have the luxury of writing my own code here and will therefore use what feels more natural to me: \(l = 0\) represents the input data, and \(l = 1\) is the first layer of real neurons.

I want it to be easy to cross-reference the following against Nielsen’s outline. However I’m willing to compromise on that a bit in order to preserve cohesiveness here. I have therefore taken the compromise position of modifying quotes from Nielsen to account for my notation.

Our procedure begins with the following:

  1. Input a set of training examples
  2. For each training example x: Set the corresponding input activation \(a^{0,x}\), and perform the following steps:

Let’s first define a few symbols, after which we will show how to package multiple training example vectors into a single matrix:

Symbol Shape Description
\(l\) Layer index (\(0\) = input)
\(N^l\) Number of neurons in layer \(l\)
\(M\) Mini-batch size
\(\mathbf{x}^{l,m}\) \({N^{l-1} \times 1}\) Inputs for layer \(l\) at minibatch position \(m\)

This is the general form for inputs into a layer:

\[ \underbrace{\begin{bmatrix} \mathbf{x}^{l,1} , \mathbf{x}^{l,2} , \dots , \mathbf{x}^{l,M} \end{bmatrix}}_{N^{l-1} \times M} \]

Our training data provides the inputs into layer \(l = 1\):

\[ \underbrace{\begin{bmatrix} \mathbf{x}^{1,1} , \mathbf{x}^{1,2} , \dots , \mathbf{x}^{1,M} \end{bmatrix}}_{N^{0} \times M} \]

Now it is time to:

Feedforward: For each \(l = 1, 2, \ldots, L\) compute \(z^{x,l} = w^l a^{x,l-1}+b^l\) and \(a^{x,l} = \sigma(z^{x,l})\).

We know \(\mathbf{a}\) takes \(\mathbf{z}\) as input, so let’s start there. Studying the definition for \(\text{(INPUT)}\) reveals that it will naturally parallelize:

\[ \begin{align*} \\ \mathbf{z}^l &\equiv \mathbf{W}^l \mathbf{a}^{l-1} + \mathbf{b}^l \tag{INPUT} \\ \\ \end{align*} \]

The input to our first layer of actual neurons is just our selected training examples for the batch, i.e. when \(l=1\) then \(\mathbf{a}^{l-1} = \mathbf{a}^{0}\) corresponds to the input training data.

Multiplying our matrix of inputs against the weight matrix is, by definition, the same as multiplying each vector of input against the weight matrix. We can then add the same vector of bias terms to every vector in the resulting matrix. This will give us a matrix of pre-activations where every column is the vector of pre-activations for the corresponding column in the input matrix.

Symbol Shape Description
\(\mathbf{W}^{l}\) \(N^l \times N^{l-1}\) Weights in layer \(l\)
\(\mathbf{w}^{l}_n\) \(1 \times N^{l-1}\) Weights for neuron \(n\) in layer \(l\)
\(\mathbf{b}^l\) \(N^l \times 1\) Biases for layer \(l\)
\(\mathbf{z}^{l,m}\) \(N^l \times 1\) Pre-activations for layer \(l\) at minibatch position \(m\)

\[ \underbrace{\begin{bmatrix} \mathbf{z}^{l,1}, \mathbf{z}^{l,2}, \dots, \mathbf{z}^{l,M} \end{bmatrix}}_{N^l \times M} = \underbrace{\begin{bmatrix} \mathbf{w}^{l}_1 \mathbf{w}^{l}_2 \dots \mathbf{w}^{l}_{N^l} \end{bmatrix}}_{N^{l} \times N^{l-1}} \underbrace{\begin{bmatrix} \mathbf{x}^{l,1}, \mathbf{x}^{l,2}, \dots, \mathbf{x}^{l,M} \end{bmatrix}}_{N^{l-1} \times M} + \underbrace{\begin{bmatrix} \mathbf{b}^l, \mathbf{b}^l, \dots, \mathbf{b}^l \end{bmatrix}}_{N^l \times M} \]

Now recall our definition for \(\text{(ACT)}\):

\[ \begin{align*} \\ \mathbf{a}^l &\equiv \sigma(\mathbf{z}^l) \tag{ACT} \\ \end{align*} \]

The \(\sigma\) is applied elementwise, so we can easily simplify our preceding notation into a fully matrix-based form and use that to clearly indicate how we will arrive at our activations.

Symbol Shape Description
\(\mathbf{B}^l\) \(N^l \times M\) Mini-batch biases for layer \(l\)
\(\mathbf{X}^l\) \(N^{l-1} \times M\) Mini-batch inputs for layer \(l\)
\(\mathbf{Z}^l\) \(N^l \times M\) Mini-batch pre-activations for layer \(l\)
\(\mathbf{A}^l\) \(N^l \times M\) Mini-batch activations for layer \(l\)

\[ \begin{align*} \mathbf{Z}^l &= \mathbf{W}^{l} \mathbf{X}^l + \mathbf{B}^l \\ \mathbf{A}^l &= \sigma (\mathbf{Z}^l) \\ \mathbf{X}^{l+1} &= \mathbf{A}^l \end{align*} \]

A brief detour to discuss the cost function: While the fundamental equations of backpropagation are largely agnostic to the specifics of the cost function, in order to actually train a neural network we will need to know how to calculate the gradient of our particular cost function with respect to the output layer activations. In other words: we need to know what \(\nabla_{\mathbf{a}} C\) really is.

We are using the quadratic cost function:

\[C = \frac{1}{2} \| \mathbf{y} - \mathbf{a}^L \|^2 = \frac{1}{2} \sum_j (y_j - a_j^L)^2\]

And we can therefore find that \(\nabla_{\mathbf{a}} C\) is:

\[ \begin{align*} \\ \nabla_{\mathbf{a}} C &= \begin{bmatrix} \partial C / \partial a_1^L \\ \partial C / \partial a_2^L \\ \vdots \end{bmatrix} \\ \frac{\partial C}{\partial a_k^L} &= \frac{\partial}{\partial a_k^L} \left[ \frac{1}{2} \sum_j (y_j - a_j^L)^2 \right] \\ &= \frac{\partial}{\partial a_k^L} \left[ (y_1 - a_1^L)^2 + \cdots + (y_k - a_k^L)^2 + \cdots \right] \frac{1}{2} \\ &= \frac{1}{2} \left[ 0 + \cdots + \frac{\partial}{\partial a_k^L} (y_k - a_k^L)^2 + \cdots \right] = \frac{\partial}{\partial a_k^L} \frac{1}{2} (y_k - a_k^L)^2 \\ &= (2) \frac{1}{2} (y_k - a_k^L) \frac{\partial}{\partial a_k^L} (y_k - a_k^L) = (y_k - a_k^L)(-1) = a_k^L - y_k \\ \nabla_{\mathbf{a}} C &= \begin{bmatrix} a_1^L - y_1 \\ a_2^L - y_2 \\ \vdots \\ a_k^L - y_k \end{bmatrix} \\ \nabla_{\mathbf{a}} C &= \mathbf{a}^L - \mathbf{y} \\ \end{align*} \]

We now have everything we need for the next step:

Output error \(\delta^{L,x}\): Compute the vector \(\delta^{L,x} = \nabla_a C_x \odot \sigma'(z^{L,x})\).

Having determined that \(\nabla_{\mathbf{a}} C = \mathbf{a}^L - \mathbf{y}\), and with the knowledge that both subtraction and the Hadamard product are elementwise operations in linear algebra, we can simply continue our approach of packing mini-batch vectors into matrices to achieve a fully matrix-based approach.

Symbol Shape Description
\(L\) Output layer index
\(\mathbf{a}^{l,m}\) \(N^l \times 1\) Activations for layer \(l\), mini-batch sample \(m\)
\(\mathbf{y}^m\) \(N^L \times 1\) Targets for mini-batch sample \(m\)
\(\boldsymbol{\delta}^{l,m}\) \(N^l \times 1\) Errors for layer \(l\), mini-batch sample \(m\)

\[ \underbrace{\begin{bmatrix} \boldsymbol{\delta}^{L,1}, \boldsymbol{\delta}^{L,2}, \dots, \boldsymbol{\delta}^{L,M} \end{bmatrix}}_{N^L \times M} =\\ \left( \underbrace{\begin{bmatrix} \mathbf{a}^{L,1}, \mathbf{a}^{L,2}, \dots, \mathbf{a}^{L,M} \end{bmatrix}}_{N^L \times M} -\\ \underbrace{\begin{bmatrix} \mathbf{y}^1, \mathbf{y}^2, \dots, \mathbf{y}^M \end{bmatrix}}_{N^L \times M} \right) \]

Once again, we will tidy this up into fully matrix-based notation.

Symbol Shape Description
\(\mathbf{\Delta}^l\) \(N^l \times M\) Mini-batch errors for layer \(l\)
\(\mathbf{Y}\) \(N^L \times M\) Mini-batch targets

\[ \mathbf{\Delta}^L = (\mathbf{A}^L - \mathbf{Y}) \odot \sigma'(\mathbf{Z}^L) \]

We have now crossed the inflection point and are moving into the backward pass proper:

Backpropagate the error: For each \(l = L-1, L-2, \ldots, 1\) compute \(\delta^{l,x} = ((w^{l+1})^T \delta^{l+1,x}) \odot \sigma'(z^{l,x})\).

This can be straightforwardly expressed in a fully matrix-based form without defining any new symbols:

\[ \underbrace{\begin{bmatrix} \boldsymbol{\delta}^{l,1}, \boldsymbol{\delta}^{l,2}, \dots, \boldsymbol{\delta}^{l,M} \end{bmatrix}}_{N^l \times M} = \left( \underbrace{(\mathbf{W}^{l+1})^T}_{N^l \times N^{l+1}} \underbrace{\begin{bmatrix} \boldsymbol{\delta}^{l+1,1}, \boldsymbol{\delta}^{l+1,2}, \dots, \boldsymbol{\delta}^{l+1,M} \end{bmatrix}}_{N^{l+1} \times M} \right) \odot \sigma' \left( \underbrace{\begin{bmatrix} \mathbf{z}^{l,1}, \mathbf{z}^{l,2}, \dots, \mathbf{z}^{l,M} \end{bmatrix}}_{N^l \times M} \right) \]

\[ \mathbf{\Delta}^l = \left( (\mathbf{W}^{l+1})^T \mathbf{\Delta}^{l+1} \right) \odot \sigma'(\mathbf{Z}^l) \]

We are on the home stretch now. The time has come for:

Gradient descent: For each \(l = L, L-1, \ldots, 1\) update the weights according to the rule \(w^l \rightarrow w^l-\frac{\eta}{m} \sum_x \delta^{l,x} (a^{l-1,x})^T\), and the biases according to the rule \(b^l \rightarrow b^l-\frac{\eta}{m} \sum_x \delta^{l,x}\).

In essence, we are re-arranging \(\text{(BP4)}\) as a vector outer product, which produces a matrix of weight errors for a given mini-batch at layer \(l\). We then sum up all of these weight error matrices and perform elementwise division by \(M\), which produces an average weight error matrix. We scale this matrix by the learning rate \(\eta\) and use that to update the layer weights.

We can represent this as a very elegant matrix multiplication, although the reason why that works is not so obvious as in the prior steps. The basic idea is to pack the layer error vectors for the mini-batch into a matrix, and multiply that by another matrix consisting of the transposed layer input vectors for the mini-batch stacked vertically (which we can do because when transposed they form row vectors). I spell this all out in somewhat excruciating detail below.

Symbol Shape Description
\(\eta\) Learning rate
\((\mathbf{X}^l)^T\) \(M \times N^{l-1}\) Transposed mini-batch input matrix (\(\mathbf{X}^l = \mathbf{A}^{l-1}\))

\[ \begin{align*} \\ \underbrace{ \sum_{m=1}^M \boldsymbol{\delta}^{l,m} (\mathbf{x}^{l,m})^T }_{N^l \times N^{l-1}} &= \begin{bmatrix} \delta^{l,1}_1 x^{l,1}_1 & \delta^{l,1}_1 x^{l,1}_2 & \dots & \delta^{l,1}_1 x^{l,1}_{N^{l-1}} \\ \delta^{l,1}_2 x^{l,1}_1 \\ \dots & & \ddots \\ \delta^{l,1}_{N^l} x^{l,1}_1 & & & \delta^{l,1}_{N^l} x^{l,1}_{N^{l-1}} \end{bmatrix} + \begin{bmatrix} \delta^{l,2}_1 x^{l,2}_1 & \delta^{l,2}_1 x^{l,2}_2 & \dots & \delta^{l,2}_1 x^{l,2}_{N^{l-1}} \\ \delta^{l,2}_2 x^{l,2}_1 \\ \dots & & \ddots \\ \delta^{l,2}_{N^l} x^{l,2}_1 & & & \delta^{l,2}_{N^l} x^{l,2}_{N^{l-1}} \end{bmatrix} + \dots + \begin{bmatrix} \delta^{l,M}_1 x^{l,M}_1 & \delta^{l,M}_1 x^{l,M}_2 & \dots & \delta^{l,M}_1 x^{l,M}_{N^{l-1}} \\ \delta^{l,M}_2 x^{l,M}_1 \\ \dots & & \ddots \\ \delta^{l,M}_{N^l} x^{l,M}_1 & & & \delta^{l,M}_{N^l} x^{l,M}_{N^{l-1}} \end{bmatrix} \\ &= \sum_{m=1}^M \begin{bmatrix} \delta^{l,m}_1 x^{l,m}_1 & \delta^{l,m}_1 x^{l,m}_2 & \dots & \delta^{l,m}_1 x^{l,m}_{N^{l-1}} \\ \delta^{l,m}_2 x^{l,m}_1 \\ \dots & & \ddots \\ \delta^{l,m}_{N^l} x^{l,m}_1 & & & \delta^{l,m}_{N^l} x^{l,m}_{N^{l-1}} \end{bmatrix} \\ &= \underbrace{\begin{bmatrix} \delta^{l,1}_1 & \delta^{l,2}_1 & \dots & \delta^{l,M}_1 \\ \delta^{l,1}_2 \\ \dots & & \ddots \\ \delta^{l,1}_{N^l} & & & \delta^{l,M}_{N^l} \end{bmatrix}}_{N^l \times M} \underbrace{\begin{bmatrix} x^{l,1}_1 & x^{l,1}_2 & \dots & x^{l,1}_{N^{l-1}} \\ x^{l,2}_1 \\ \dots & & \ddots \\ x^{l,M}_{1} & & & x^{l,M}_{N^{l-1}} \end{bmatrix}}_{M \times N^{l-1}} \\ &= \underbrace{\begin{bmatrix} \boldsymbol{\delta}^{l,1} , \boldsymbol{\delta}^{l,2} , \dots , \boldsymbol{\delta}^{l,M} \end{bmatrix}}_{N^l \times M} \underbrace{\begin{bmatrix} (\mathbf{x}^{l,1})^T \\ (\mathbf{x}^{l,2})^T \\ \vdots \\ (\mathbf{x}^{l,M})^T \end{bmatrix}}_{M \times N^{l-1}} \\ \frac{\partial C}{\partial \mathbf{W}^l} &= \mathbf{\Delta}^l (\mathbf{X}^l)^T \end{align*} \]

Expressing our update step in fully matrix-based form we obtain:

\[ \begin{align*} \mathbf{W}^l &\rightarrow \mathbf{W}^l - \frac{\eta}{M} \frac{\partial C}{\partial \mathbf{W}^l} \\ \therefore \mathbf{W}^l &\rightarrow \mathbf{W}^l - \frac{\eta}{M} \mathbf{\Delta}^l (\mathbf{X}^l)^T \end{align*} \]

The same basic idea holds for the biases, only simpler: note what happens if we multiply the matrix of mini-batch layer error vectors with a vector of \(1\text{'s}\).

Symbol Shape Description
\(\mathbf{1}_M\) \(M \times 1\) Column vector of ones

\[ \begin{align*} \sum_{m=1}^M \boldsymbol{\delta}^{l,m} &= \underbrace{\begin{bmatrix} \delta^{l,1}_1 + \delta^{l,2}_1 + \dots + \delta^{l,M}_1 \\ \delta^{l,1}_2 + \delta^{l,2}_2 + \dots + \delta^{l,M}_2 \\ \dots \\ \delta^{l,1}_{N^l} + \delta^{l,2}_{N^l} + \dots + \delta^{l,M}_{N^l} \end{bmatrix}}_{N^l \times 1} \\ &= \underbrace{\begin{bmatrix} \boldsymbol{\delta}^{l,1} , \boldsymbol{\delta}^{l,2} , \dots , \boldsymbol{\delta}^{l,M} \end{bmatrix}}_{N^l \times M} \underbrace{\begin{bmatrix} 1 \\ 1 \\ \vdots \\ 1 \end{bmatrix}}_{M \times 1} \\ \frac{\partial C}{\partial \mathbf{b}^l} &= \mathbf{\Delta}^l \mathbf{1}_M \end{align*} \]

Now imagine if we instead multiplied by an \(M \times M\) matrix of \(1\text{'s}\).

Symbol Shape Description
\(\mathbf{1}_{M \times M}\) \(M \times M\) Square matrix of ones

\[ \begin{align*} & \underbrace{\begin{bmatrix} \boldsymbol{\delta}^{l,1} , \boldsymbol{\delta}^{l,2} , \dots , \boldsymbol{\delta}^{l,M} \end{bmatrix}}_{N^l \times M} \underbrace{\begin{bmatrix} \mathbf{1}_M, \mathbf{1}_M, \dots, \mathbf{1}_M \end{bmatrix}}_{M \times M} \\ & \underbrace{\begin{bmatrix} \mathbf{\Delta}^l \mathbf{1}_M , \mathbf{\Delta}^l \mathbf{1}_M , \dots , \mathbf{\Delta}^l \mathbf{1}_M \end{bmatrix}}_{N^l \times M} \\ & \underbrace{\begin{bmatrix} \frac{\partial C}{\partial \mathbf{b}^l} , \frac{\partial C}{\partial \mathbf{b}^l} , \dots , \frac{\partial C}{\partial \mathbf{b}^l} \end{bmatrix}}_{N^l \times M} = \frac{\partial C}{\partial \mathbf{B}^l} \\ \therefore \frac{\partial C}{\partial \mathbf{B}^l} &= \mathbf{\Delta}^l \mathbf{1}_{M \times M} \end{align*} \]

Finally, we can express our bias update step in fully matrix-based form:

\[ \begin{align*} \mathbf{b}^l &\rightarrow \mathbf{b}^l - \frac{\eta}{M}\frac{\partial C}{\partial \mathbf{b}^l} \\ \therefore \mathbf{B}^l &\rightarrow \mathbf{B}^l - \frac{\eta}{M} \mathbf{\Delta}^l \mathbf{1}_{M \times M} \end{align*} \]

Implementation

Parameters and architecture

I tried to write code that was primarily legible, and secondarily performant. Rather than worry about reuse and generality (i.e. proper functions and data structures), I just gathered the key configuration up here. You can tweak these settings and then use Run Selected Cell and All Below to see your new results at the bottom.

N = [784, 30, 10] # Neurons per layer
M = 128 # Minibatch size
NUM_EPOCHS = 30
ETA = 3.0

Loading the data

import numpy as np
import gzip

# mnist_modern.pkl.gz should work with modern Python, and is adapted from
# mnist.pkl.gz, which comes from:
#   https://github.com/MichalDanielDobrzanski/DeepLearningPython
# and is distributed under the MIT license.
mnist_file = gzip.open('mnist_modern.pkl.gz', 'rb')
mnist_data = np.load(mnist_file, allow_pickle=True, encoding='latin1')

# Each entry contains an array of 784 float
# values on [0.0, 1.0] representing the 28x28 digit bitmap.
training_images = mnist_data[0][0]
validation_images = mnist_data[1][0]
test_images = mnist_data[2][0]

# Integer value 0-9 of the corresponding image.
training_labels = mnist_data[0][1]
validation_labels = mnist_data[1][1]
test_labels = mnist_data[2][1]

# Setup 1-hot encoding:
# Each entry contains a size 10 array of floats with a 1.0 in the
# position corresponding to the digit and a 0.0 everywhere else.
# This represents what perfect output from our neural network
# should look like for this digit.
training_targets = np.eye(10)[training_labels]
validation_targets = np.eye(10)[validation_labels]
test_targets = np.eye(10)[test_labels]

mnist_file.close()

Building the network

Again, I wanted code that was clear and performant, and to my mind, that calls for a struct of arrays design. First, when in doubt, it’s a good default for performance in the modern “parallel everything” era. Two, it naturally aligns closely with mathematical notation, e.g \(A^l\) becomes A[l]. Third: I just like it.

So in reading the following, cross-reference it with the math worked out above, and all should become clear.

L = len(N)-1 # Output layer
FORWARD_LAYERS = [x for x in range(1, L+1)] # Everything except the input
BACKWARD_LAYERS = [x for x in reversed(range(1, L))] # Neither output nor input

rng = np.random.default_rng()

# The first layer is blank except for A, which is where input arrives.
W = [""]
B = [""]
X = [""]
Z = [""]
A = [np.empty((N[0], M))]
D = [""]

for l in FORWARD_LAYERS:
    # Note that data "goes in the columns and comes out the rows"
    # so all matrices have to be sized accordingly: columns in, rows out.
    W.append(rng.standard_normal((N[l], N[l-1])) * 0.1)
    b = rng.standard_normal((N[l], 1))
    B.append(np.repeat(b, M, axis=1))
    X.append(np.empty((N[l-1], M)))
    Z.append(np.empty((N[l], M)))
    A.append(np.empty((N[l], M)))
    D.append(np.empty((N[l], M)))

Sanity check

The next couple of blocks ensure that we have all the parameters we ought to have in each layer, and also that their dimensions correctly relate to each other and across layers.

print("W:", len(W))
print("B:", len(B))
print("X:", len(X))
print("Z:", len(Z))
print("A:", len(A))
print("D:", len(D))
assert(len(W) == len(B))
assert(len(B) == len(X))
assert(len(X) == len(Z))
assert(len(Z) == len(A))
assert(len(A) == len(D))
W: 3
B: 3
X: 3
Z: 3
A: 3
D: 3
print(f"A[{0}]: {A[0].shape} -> copies into...")
assert(A[0].shape == X[1].shape)

for l in FORWARD_LAYERS:
    print(f"X[{l}]: {X[l].shape} -> multiplies against...")
    print(f"W[{l}]: {W[l].shape} -> added with...")
    print(f"B[{l}]: {B[l].shape} -> result stored in...")
    print(f"Z[{l}]: {Z[l].shape} -> activated and stored in...")
    print(f"A[{l}]: {A[l].shape} -> copies into...")
    assert(X[l].shape[0] == W[l].shape[1])
    assert(X[l].shape[1] == B[l].shape[1] and W[l].shape[0] == B[l].shape[0])
    assert(B[l].shape == Z[l].shape)
    assert(Z[l].shape == A[l].shape)
    assert(A[l].shape == D[l].shape)
A[0]: (784, 128) -> copies into...
X[1]: (784, 128) -> multiplies against...
W[1]: (30, 784) -> added with...
B[1]: (30, 128) -> result stored in...
Z[1]: (30, 128) -> activated and stored in...
A[1]: (30, 128) -> copies into...
X[2]: (30, 128) -> multiplies against...
W[2]: (10, 30) -> added with...
B[2]: (10, 128) -> result stored in...
Z[2]: (10, 128) -> activated and stored in...
A[2]: (10, 128) -> copies into...

Forward pass

We will need to use the forward pass logic both to train the network and to validate it, so I broke it into a helper function. I also defined helper functions for the activation just to make the code more readable and better aligned with the math.

def sigmoid(z):
    return 1.0/(1.0+np.exp(-z))

def sigmoid_prime(z):
    return sigmoid(z)*(1-sigmoid(z))
def forward_pass():
    for l in FORWARD_LAYERS:
        X[l][:] = A[l-1]

        # The @ does matrix multiplication.
        Z[l][:] = W[l] @ X[l] + B[l]
        A[l][:] = sigmoid(Z[l])

Training and validation

With all the housekeeping sorted, we are ready for the main event. The general idea is to create a shuffled list of training sample indices, then step through them one mini-batch at a time until we have completed an epoch. Each epoch shuffles, completes the forward pass, calculates the output error, completes the backward pass, and then updates all the weights and biases.

Validation proceeds similarly, except that we use the validation data with its original indices and stop after the forward pass. We conclude by printing out performance statistics for the epoch.

Note: Why are we using validation data here, and not the test data like in chapter 1? Because that is actually the best practice, as Nielsen discusses in chapter 3.

RATE = ETA / M

# Ensure we never have a minibatch that is only partially full
# by defining our number of samples as an integer multiple of our
# minibatch size (hence our use of the // "floor division" operator).
NUM_TRAINING_SAMPLES = (len(training_images) // M) * M
NUM_VALIDATION_SAMPLES = (len(validation_images) // M) * M

# We only need one matrix to store minibatch targets; that is, this
# doesn't need to exist per layer.
Y = np.empty((N[L], M))

# This square matrix of ones equal to the minibatch size is used
# for updating the biases efficiently.
MxM = np.ones((M, M))

for epoch in range(NUM_EPOCHS):
    shuffled_indices = rng.permutation(NUM_TRAINING_SAMPLES)

    #### TRAINING
    for m in range(0, NUM_TRAINING_SAMPLES, M):
        indices = shuffled_indices[m : m + M]
        images = training_images[indices]

        # Our data is "one image per row," for our forward pass we need
        # "one image per column," so we transpose and copy it over.
        # Note that the [:] below is Python for "copy the data into
        # previously allocated memory," without it we would be reassigning
        # pointers.
        A[0][:] = images.T

        forward_pass()

        #### OUTPUT ERROR
        # We do the same transpose trick with the targets
        # that we did with the images.
        targets = training_targets[indices]
        Y[:] = targets.T

        # The * does elementwise multiplication (Hadamard product).
        D[L][:] = (A[L] - Y) * sigmoid_prime(Z[L])

        #### BACKWARD PASS
        for l in BACKWARD_LAYERS:
            D[l] = (W[l+1].T @ D[l+1]) * sigmoid_prime(Z[l])

        #### GRADIENT DESCENT
        for l in FORWARD_LAYERS:
            W[l] = W[l] - RATE * (D[l] @ X[l].T)
            B[l] = B[l] - RATE * (D[l] @ MxM)

    #### VALIDATION
    total_correct = 0
    for m in range(0, NUM_VALIDATION_SAMPLES, M):
        images = validation_images[m : m + M]
        targets = validation_targets[m : m + M]
        
        A[0][:] = images.T
        Y[:] = targets.T

        forward_pass()

        # We use argmax() columnwise (axis=0) to get the 1-hot encoded
        # result of our feedforward pass. We compare that against our
        # targets and store a list equal to the minibatch size, where
        # each entry has a 1 for a match and a 0 for miss. We then 
        # total up that list and add it to our running total for the epoch.
        expected = np.argmax(Y, axis=0)
        actual = np.argmax(A[L], axis=0)
        correct = (expected == actual) * 1
        subtotal = np.sum(correct)
        total_correct += subtotal

    percent_correct = total_correct / NUM_VALIDATION_SAMPLES
    print("{0}:\t{1} / {2}\t({3:.2%})".format(
        epoch, total_correct, NUM_VALIDATION_SAMPLES, percent_correct))
0:  9056 / 9984 (90.71%)
1:  9224 / 9984 (92.39%)
2:  9310 / 9984 (93.25%)
3:  9358 / 9984 (93.73%)
4:  9400 / 9984 (94.15%)
5:  9450 / 9984 (94.65%)
6:  9465 / 9984 (94.80%)
7:  9476 / 9984 (94.91%)
8:  9506 / 9984 (95.21%)
9:  9516 / 9984 (95.31%)
10: 9542 / 9984 (95.57%)
11: 9533 / 9984 (95.48%)
12: 9546 / 9984 (95.61%)
13: 9550 / 9984 (95.65%)
14: 9569 / 9984 (95.84%)
15: 9567 / 9984 (95.82%)
16: 9582 / 9984 (95.97%)
17: 9590 / 9984 (96.05%)
18: 9594 / 9984 (96.09%)
19: 9591 / 9984 (96.06%)
20: 9608 / 9984 (96.23%)
21: 9604 / 9984 (96.19%)
22: 9611 / 9984 (96.26%)
23: 9616 / 9984 (96.31%)
24: 9625 / 9984 (96.40%)
25: 9603 / 9984 (96.18%)
26: 9618 / 9984 (96.33%)
27: 9628 / 9984 (96.43%)
28: 9636 / 9984 (96.51%)
29: 9630 / 9984 (96.45%)