Answers.Body: <p>I'm not familiar with python I'm afraid, but I'll give a go at presenting some of the maths...</p>
  
  <p>In order to perform the <strong>back propagation</strong> using matrices, matrix <a href="https://en.wikipedia.org/wiki/Transpose" rel="nofollow noreferrer">transpose</a> is used for different sized layers. Also note that when using matrices for this, we need to distinguish different types of matrix multiplication, namely matrix <a href="https://en.wikipedia.org/wiki/Matrix_multiplication" rel="nofollow noreferrer">product</a>, and the <a href="https://en.wikipedia.org/wiki/Hadamard_product_(matrices)" rel="nofollow noreferrer">Hadamard product</a> which operates differently (the latter of which is denoted by a circle with a dot in the centre).</p>
  
  <p>Note the <strong>back propagation</strong> formulas:</p>
  
  <p>(EQ1)
  <span class="math-container">\begin{equation*}\delta ^{l} = (w^{l+1})^{T} \delta ^{l+1} \odot \sigma {}' (z ^{l})\end{equation*}</span></p>
  
  <p>(EQ2)
  <span class="math-container">\begin{equation*}\frac{\partial E}{\partial w}=a^{l-1}\delta ^{l}\end{equation*}</span></p>
  
  <p>As you can see the transpose of the weight matrix is used to multiply against the delta of the layer. </p>
  
  <p>As an example, consider a simple network with an input layer of 3 neurons, and an output layer of 2 neurons.</p>
  
  <p>The <em>unactivated</em> <strong>feed forward</strong> output is given by...</p>
  
  <p><span class="math-container">\begin{equation*}z_{1} = h_{1} w_{1} + h_{2} w_{3} + h_{3} w_{5}\end{equation*}</span>
  <span class="math-container">\begin{equation*}z_{2} = h_{1} w_{2} + h_{2} w_{4} + h_{3} w_{6}\end{equation*}</span></p>
  
  <p>Which can be represented in matrix form as...</p>
  
  <p><span class="math-container">\begin{equation*}
  z
  =
  \begin{bmatrix}z_{1}\\z_{2}\end{bmatrix}
  =
  \begin{bmatrix}h_{1} w_{1} + h_{2} w_{3} + h_{3} w_{5}\\h_{1} w_{2} + h_{2} w_{4} + h_{3} w_{6}\end{bmatrix}
  \end{equation*}</span></p>
  
  <p>(EQ3)
  <span class="math-container">\begin{equation*}
  =
  \begin{bmatrix}w_{1}&amp;w_{3}&amp;w_{5}\\w_{2}&amp;w_{4}&amp;w_{6}\end{bmatrix}
  \begin{bmatrix}h_{1}\\h_{2}\\h_{3}\end{bmatrix}
  \end{equation*}</span></p>
  
  <p>Which allows us to <strong>forward propagate</strong> from a layer of 3 neurons, to a layer of 2 neurons (via matrix multiplication).</p>
  
  <p>For the <strong>back propagation</strong>, we need the <em>error</em> of each neruon...</p>
  
  <p>N.B <em>Cost</em> is a metric used to denote the <em>error</em> of the entire network and depends on your cost function. It is usually the mean of the sum of the <em>errors</em> of each nuerons, but of course depends on your cost function used.</p>
  
  <p>e.g for <a href="https://en.wikipedia.org/wiki/Mean_squared_error" rel="nofollow noreferrer">MSE</a> is...
  <span class="math-container">\begin{equation*}
  C_{MSE}=\frac{1}{N}\sum (o_{n}-t_{n})^{2}
  \end{equation*}</span></p>
  
  <p>We are interested in the derivative of the <em>error</em> for each neuron (not the <em>cost</em>), which by the chain rule is...</p>
  
  <p><span class="math-container">\begin{equation*}
  \frac{\partial E}{\partial z}
  =
  \frac{\partial E}{\partial o}
  \frac{\partial o}{\partial z}
  \end{equation*}</span></p>
  
  <p>Expressed in matrix form...</p>
  
  <p><span class="math-container">\begin{equation*}
  \frac{\partial E}{\partial z}
  =
  \begin{bmatrix}
  \frac{\partial E}{\partial z_{1}}\\
  \frac{\partial E}{\partial z_{2}}
  \end{bmatrix}
  =
  \begin{bmatrix}
  \frac{\partial E}{\partial o_{1}}
  \frac{\partial o_{1}}{\partial z_{1}}\\
  \frac{\partial E}{\partial o_{2}}
  \frac{\partial o_{2}}{\partial z_{2}}
  \end{bmatrix}
  \end{equation*}</span></p>
  
  <p>(EQ4)
  <span class="math-container">\begin{equation*}
  =
  \begin{bmatrix}
  \frac{\partial E}{\partial o_{1}}
  \\
  \frac{\partial E}{\partial o_{2}}
  \end{bmatrix}
  \odot 
  \begin{bmatrix}
  \frac{\partial o_{1}}{\partial z_{1}}
  \\
  \frac{\partial o_{2}}{\partial z_{2}}
  \end{bmatrix}
  \end{equation*}</span></p>
  
  <p>Note the use of the hadamard product here. In fact, given these are simply vectors, vector dot product would work, but hadamard is used as this becomes important later when using this expression with the matrix equations as we want to distinguish the use of hadamard product as opposed to matrix product.</p>
  
  <p>We start off with our first delta error, which for the first layer <strong>back propagation</strong> is...</p>
  
  <p><span class="math-container">\begin{equation*}
  \delta ^{L} = \frac{\partial E}{\partial o}
  \end{equation*}</span></p>
  
  <p>And then we want to calculate the next delta error using the formula (EQ1)...</p>
  
  <p><span class="math-container">\begin{equation*}
  \delta ^{l} = (w^{l+1})^{T} \delta ^{l+1} \odot \sigma {}' (z ^{l}) = \frac{\partial E}{\partial h}
  \end{equation*}</span></p>
  
  <p>Explicitly, the equations are...</p>
  
  <p><span class="math-container">\begin{equation*}\frac{\partial E}{\partial h_{1}}=\frac{\partial E}{\partial z_{1}} w_{1} + \frac{\partial E}{\partial z_{2}} w_{2}\end{equation*}</span>
  <span class="math-container">\begin{equation*}\frac{\partial E}{\partial h_{2}}=\frac{\partial E}{\partial z_{1}} w_{3} + \frac{\partial E}{\partial z_{2}} w_{4}\end{equation*}</span>
  <span class="math-container">\begin{equation*}\frac{\partial E}{\partial h_{3}}=\frac{\partial E}{\partial z_{1}} w_{5} + \frac{\partial E}{\partial z_{2}} w_{6}\end{equation*}</span></p>
  
  <p>Also given the transpose of the weight matrix which allows us to <strong>back propagate</strong> from a layer of 2 neurons to a layer of 3 neurons (via matrix multiplication) ...</p>
  
  <p><span class="math-container">\begin{equation*}
  \left (w  \right )^{T}
  =
  \left ( \begin{bmatrix}
  w_{1} &amp; w_{3} &amp; w_{5}\\ 
  w_{2} &amp; w_{4} &amp; w_{6}
  \end{bmatrix}
   \right )^{T}
  =
  \begin{bmatrix}
  w_{1} &amp; w_{2}\\ 
  w_{3} &amp; w_{4}\\ 
  w_{5} &amp; w_{6}
  \end{bmatrix}
  \end{equation*}</span></p>
  
  <p>So in a similar way to how we represented the forward pass (EQ3), this can be represented in matrix form...</p>
  
  <p><span class="math-container">\begin{equation*}
  \frac{\partial E}{\partial h}
  =
  \begin{bmatrix}
  \frac{\partial E}{\partial h_{1}}\\ 
  \frac{\partial E}{\partial h_{2}}\\ 
  \frac{\partial E}{\partial h_{3}}
  \end{bmatrix}
  =
  \begin{bmatrix}
  \frac{\partial E}{\partial z_{1}} w_{1} + \frac{\partial E}{\partial z_{2}} w_{2}
  \\ 
  \frac{\partial E}{\partial z_{1}} w_{3} + \frac{\partial E}{\partial z_{2}} w_{4}
  \\
  \frac{\partial E}{\partial z_{1}} w_{5} + \frac{\partial E}{\partial z_{2}} w_{6}
  \end{bmatrix}
  \end{equation*}</span></p>
  
  <p>(EQ5)
  <span class="math-container">\begin{equation*}
  =
  \begin{bmatrix}
  w_{1} &amp; w_{2}\\ 
  w_{3} &amp; w_{4}\\ 
  w_{5} &amp; w_{6}
  \end{bmatrix}
  \begin{bmatrix}
  \frac{\partial E}{\partial z_{1}}
  \\ 
  \frac{\partial E}{\partial z_{2}}
  \end{bmatrix}
  \end{equation*}</span></p>
  
  <p>And then plugging the hadamard version of the delta (EQ4) into this, we get...</p>
  
  <p><span class="math-container">\begin{equation*}
  \begin{bmatrix}
  w_{1} &amp; w_{2}\\ 
  w_{3} &amp; w_{4}\\ 
  w_{5} &amp; w_{6}
  \end{bmatrix}
  \begin{bmatrix}
  \frac{\partial E}{\partial o_{1}}
  \\
  \frac{\partial E}{\partial o_{2}}
  \end{bmatrix}
  \odot 
  \begin{bmatrix}
  \frac{\partial o_{1}}{\partial z_{1}}
  \\
  \frac{\partial o_{2}}{\partial z_{2}}
  \end{bmatrix}
  \end{equation*}</span></p>
  
  <p>aka (EQ1) ...</p>
  
  <p><span class="math-container">\begin{equation*}(w^{l+1})^{T} \delta ^{l+1} \odot \sigma {}' (z ^{l})\end{equation*}</span></p>
  
  <p>And thus we have <strong>back propagated</strong> from a layer of 2 neurons, to a layer of 3 neurons (via matrix multiplication) thanks to transpose.</p>
  
  <p>For completeness... the other aspect of the back propagation to uses matrices, is the delta weight matrix....</p>
  
  <p><span class="math-container">\begin{equation*}
  \frac{\partial E}{\partial w}
  =
  \begin{bmatrix}
  \frac{\partial E}{\partial w_{1}} &amp; \frac{\partial E}{\partial w_{3}} &amp; \frac{\partial E}{\partial w_{5}}\\ 
  \frac{\partial E}{\partial w_{2}} &amp; \frac{\partial E}{\partial w_{4}} &amp; \frac{\partial E}{\partial w_{6}}
  \end{bmatrix}
  \end{equation*}</span></p>
  
  <p>As mentioned before, you need to cache weight matrix and the activated layer output of a forward pass of the network.</p>
  
  <p>In a similar vein to (EQ3)...</p>
  
  <p>We explcitily have the equations...</p>
  
  <p><span class="math-container">\begin{equation*}\frac{\partial E}{\partial w_{1}}=\frac{\partial E}{\partial z_{1}}h_{1}\end{equation*}</span>
  <span class="math-container">\begin{equation*}\frac{\partial E}{\partial w_{2}}=\frac{\partial E}{\partial z_{2}}h_{1}\end{equation*}</span>
  <span class="math-container">\begin{equation*}\frac{\partial E}{\partial w_{3}}=\frac{\partial E}{\partial z_{1}}h_{2}\end{equation*}</span>
  <span class="math-container">\begin{equation*}\frac{\partial E}{\partial w_{4}}=\frac{\partial E}{\partial z_{2}}h_{2}\end{equation*}</span>
  <span class="math-container">\begin{equation*}\frac{\partial E}{\partial w_{5}}=\frac{\partial E}{\partial z_{1}}h_{3}\end{equation*}</span>
  <span class="math-container">\begin{equation*}\frac{\partial E}{\partial w_{6}}=\frac{\partial E}{\partial z_{2}}h_{3}\end{equation*}</span></p>
  
  <p>Note the use of h1, h2 and h3, which are the activated outputs of the previous layer.. (or in the case of our example these are the inputs).</p>
  
  <p>Which we represent in matrix form...</p>
  
  <p><span class="math-container">\begin{equation*}
  \frac{\partial E}{\partial w}
  =
  \begin{bmatrix}
  \frac{\partial E}{\partial z_{1}}h_{1} &amp; \frac{\partial E}{\partial z_{1}}h_{2} &amp; \frac{\partial E}{\partial z_{1}}h_{3}\\ 
  \frac{\partial E}{\partial z_{2}}h_{1} &amp; \frac{\partial E}{\partial z_{2}}h_{2} &amp; \frac{\partial E}{\partial z_{2}}h_{3}
  \end{bmatrix}
  \end{equation*}</span>
  (EQ6)
  <span class="math-container">\begin{equation*}
  =
  \begin{bmatrix}
  \frac{\partial E}{\partial z_{1}} \\ \frac{\partial E}{\partial z_{2}}
  \end{bmatrix}
  \begin{bmatrix}
  h_{1} &amp;
  h_{2} &amp; 
  h_{3}
  \end{bmatrix}
  \end{equation*}</span></p>
  
  <p>Which just so happens to be (EQ2) ...</p>
  
  <p><span class="math-container">\begin{equation*}\frac{\partial E}{\partial w}=\delta ^{l}a^{l-1}\end{equation*}</span></p>
  
  <p>:)</p>
  
  <p>Since the delta weight matrix and the original weight matrix are the same dimensions, it is trivia to apply the learning rate...</p>
  
  <p><span class="math-container">\begin{equation*}
  w =
  \begin{bmatrix}
  w_{1} - \alpha \frac{\partial E}{\partial w_{1}} &amp; w_{3} - \alpha \frac{\partial E}{\partial w_{3}} &amp; w_{5} - \alpha \frac{\partial E}{\partial w_{5}}\\ 
  w_{2} - \alpha \frac{\partial E}{\partial w_{2}} &amp; w_{4} - \alpha \frac{\partial E}{\partial w_{4}} &amp; w_{6} - \alpha \frac{\partial E}{\partial w_{6}}
  \end{bmatrix}
  \end{equation*}</span></p>
  
  <p><span class="math-container">\begin{equation*}
  =
  \begin{bmatrix}w_{1}&amp;w_{3}&amp;w_{5}\\w_{2}&amp;w_{4}&amp;w_{6}\end{bmatrix}
  -
  \alpha
  \begin{bmatrix}
  \frac{\partial E}{\partial w_{1}} &amp; \frac{\partial E}{\partial w_{3}} &amp; \frac{\partial E}{\partial w_{5}}\\ 
  \frac{\partial E}{\partial w_{2}} &amp; \frac{\partial E}{\partial w_{4}} &amp; \frac{\partial E}{\partial w_{6}}
  \end{bmatrix}
  \end{equation*}</span></p>
  
  <p><span class="math-container">\begin{equation*} = w - \alpha \frac{\partial E}{\partial w}\end{equation*}</span></p>
  
  <p>I'll omit the equations for the bias, as these are easy to put in and don't require any other equations other than the ones presented here. 
  So in summary, the equations you want are (EQ3), (EQ5) and (EQ6) to use matrices, and the ability to change move between layers of differing number of neurons is the matrix transpose. Hope this helps, and let me know if you want me to expand on anything. </p>
  
  <p>It might be also important to note that you appear to be using MSE as the cost. This is perhaps not optimal for the MNIST data set which is a clasification. A better cost function to use would be <a href="https://www.google.com/search?client=firefox-b-d&amp;q=cross%20entropy" rel="nofollow noreferrer">cross entropy</a>.</p>
  
  <p>Happy propagating!</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-04-08T23:53:35.957Z
Answers.LastActivityDate: 2020-04-09T11:36:09.210Z
Answers.LastEditDate: 2020-04-09T11:36:09.210Z
Answers.UserDisplayName: lfgtm
Body: <p>I'm developing my first neural network, using the well known MNIST database of handwritten digit. I want the NN to be able to classify a number from 0 to 9 given an image.</p>
  
  <p>My neural network consists of three layers: the input layer (784 neurons, each one for every pixel of the digit), a hidden layer of 30 neurons (it could also be 100 or 50, but I'm not too worried about hyperparameter tuning yet), and the output layer, 10 neurons, each one representing the activation for every digit. That gives to me two weight matrices: one of 30x724 and a second one of 10x30.</p>
  
  <p>I know and understand the theory behind back propagation, optimization and the mathematical formulas behind that, that's not a problem as such. I can optimize the weights for the second matrix of weights, and the cost is indeed being reduced over time. But I'm not able to keep propagating that back because of the matrix structure.</p>
  
  <p>Knowing that I have find the derivative of the cost w.r.t. the weights:</p>
  
  <pre><code>d(cost) / d(w) = d(cost) / d(f(z)) * d(f(z)) / d(z) * d(z) / d(w)
  </code></pre>
  
  <p>(Being <code>f</code> the activation function and <code>z</code> the dot product plus the bias of a neuron)</p>
  
  <p>So I'm in the rightmost layer, with an output array of 10 elements. <code>d(cost) / d(f(z))</code> is the subtraction of the observed an predicted values. I can multiply that by <code>d(f(z)) / d(z)</code>, which is just <code>f'(z)</code> of the rightmost layer, also an unidimensional vector of 10 elements, having now d(cost) / d(z) calculated. Then, <code>d(z)/d(w)</code> is just the input to that layer, i.e. the output of the previous one, which is a vector of 30 elements. I figured that I can transpose d(cost) / <code>d(z)</code> so that <code>T( d(cost) / d(z) ) * d(z) / d(w)</code> gives me a matrix of (10, 30), which makes sense because it matches with the dimension of the rightmost weight matrix.</p>
  
  <p>But then I get stuck. The dimension of <code>d(cost) / d(f(z))</code> is (1, 10), for <code>d(f(z)) / d(z)</code> is (1, 30) and for <code>d(z) / d(w)</code> is (1, 784). I don't know how to come up with a result for this.</p>
  
  <p>This is what I've coded so far. The incomplete part is the <code>_propagate_back</code> method. I'm not caring about the biases yet because I'm just stuck with the weights and first I want to figure this out. </p>
  
  
  
  <pre><code>import random
  from typing import List, Tuple
  
  import numpy as np
  from matplotlib import pyplot as plt
  
  import mnist_loader
  
  np.random.seed(42)
  
  NETWORK_LAYER_SIZES = [784, 30, 10]
  LEARNING_RATE = 0.05
  BATCH_SIZE = 20
  NUMBER_OF_EPOCHS = 5000
  
  
  def sigmoid(x):
      return 1 / (1 + np.exp(-x))
  
  
  def sigmoid_der(x):
      return sigmoid(x) * (1 - sigmoid(x))
  
  
  class Layer:
  
      def __init__(self, input_size: int, output_size: int):
          self.weights = np.random.uniform(-1, 1, [output_size, input_size])
          self.biases = np.random.uniform(-1, 1, [output_size])
          self.z = np.zeros(output_size)
          self.a = np.zeros(output_size)
          self.dz = np.zeros(output_size)
  
      def feed_forward(self, input_data: np.ndarray):
          input_data_t = np.atleast_2d(input_data).T
          dot_product = self.weights.dot(input_data_t).T[0]
          self.z = dot_product + self.biases
          self.a = sigmoid(self.z)
          self.dz = sigmoid_der(self.z)
  
  
  class Network:
  
      def __init__(self, layer_sizes: List[int], X_train: np.ndarray, y_train: np.ndarray):
          self.layers = [
              Layer(input_size, output_size)
              for input_size, output_size
              in zip(layer_sizes[0:], layer_sizes[1:])
          ]
          self.X_train = X_train
          self.y_train = y_train
  
      @property
      def predicted(self) -&gt; np.ndarray:
          return self.layers[-1].a
  
      def _normalize_y(self, y: int) -&gt; np.ndarray:
          output_layer_size = len(self.predicted)
          normalized_y = np.zeros(output_layer_size)
          normalized_y[y] = 1.
  
          return normalized_y
  
      def _calculate_cost(self, y_observed: np.ndarray) -&gt; int:
          y_observed = self._normalize_y(y_observed)
          y_predicted = self.layers[-1].a
  
          squared_difference = (y_predicted - y_observed) ** 2
  
          return np.sum(squared_difference)
  
      def _get_training_batches(self, X_train: np.ndarray, y_train: np.ndarray) -&gt; Tuple[np.ndarray, np.ndarray]:
          train_batch_indexes = random.sample(range(len(X_train)), BATCH_SIZE)
  
          return X_train[train_batch_indexes], y_train[train_batch_indexes]
  
      def _feed_forward(self, input_data: np.ndarray):
          for layer in self.layers:
              layer.feed_forward(input_data)
              input_data = layer.a
  
      def _propagate_back(self, X: np.ndarray, y_observed: int):
          """
          der(cost) / der(weight) = der(cost) / der(predicted) * der(predicted) / der(z) * der(z) / der(weight)
          """
          y_observed = self._normalize_y(y_observed)
          d_cost_d_pred = self.predicted - y_observed
  
          hidden_layer = self.layers[0]
          output_layer = self.layers[1]
  
          # Output layer weights
          d_pred_d_z = output_layer.dz
          d_z_d_weight = hidden_layer.a  # Input to the current layer, i.e. the output from the previous one
  
          d_cost_d_z = d_cost_d_pred * d_pred_d_z
          d_cost_d_weight = np.atleast_2d(d_cost_d_z).T * np.atleast_2d(d_z_d_weight)
  
          output_layer.weights -= LEARNING_RATE * d_cost_d_weight
  
          # Hidden layer weights
          d_pred_d_z = hidden_layer.dz
          d_z_d_weight = X
  
          # ...
  
      def train(self, X_train: np.ndarray, y_train: np.ndarray):
          X_train_batch, y_train_batch = self._get_training_batches(X_train, y_train)
          cost_over_epoch = []
  
          for epoch_number in range(NUMBER_OF_EPOCHS):
              X_train_batch, y_train_batch = self._get_training_batches(X_train, y_train)
  
              cost = 0
              for X_sample, y_observed in zip(X_train_batch, y_train_batch):
                  self._feed_forward(X_sample)
                  cost += self._calculate_cost(y_observed)
                  self._propagate_back(X_sample, y_observed)
  
              cost_over_epoch.append(cost / BATCH_SIZE)
  
          plt.plot(cost_over_epoch)
          plt.ylabel('Cost')
          plt.xlabel('Epoch')
          plt.savefig('cost_over_epoch.png')
  
  
  training_data, validation_data, test_data = mnist_loader.load_data()
  X_train, y_train = training_data[0], training_data[1]
  
  network = Network(NETWORK_LAYER_SIZES, training_data[0], training_data[1])
  network.train(X_train, y_train)
  </code></pre>
  
  <p>This is the code for <code>mnist_loader</code>, in case someone wanted to reproduce the example:</p>
  
  <pre><code>import pickle
  import gzip
  
  
  def load_data():
      f = gzip.open('data/mnist.pkl.gz', 'rb')
      training_data, validation_data, test_data = pickle.load(f, encoding='latin-1')
      f.close()
  
      return training_data, validation_data, test_data
  
  </code></pre>
  
Tags: <neural-networks><backpropagation>
Title: How to perform back propagation with different sized layers?
UserDisplayName: Julen
ThreadText: # How to perform back propagation with different sized layers?
  
  (Asked by: username_0 on 2020-04-06)
  
  I'm developing my first neural network, using the well known MNIST database of handwritten digit. I want the NN to be able to classify a number from 0 to 9 given an image.
  My neural network consists of three layers: the input layer (784 neurons, each one for every pixel of the digit), a hidden layer of 30 neurons (it could also be 100 or 50, but I'm not too worried about hyperparameter tuning yet), and the output layer, 10 neurons, each one representing the activation for every digit. That gives to me two weight matrices: one of 30x724 and a second one of 10x30.
  I know and understand the theory behind back propagation, optimization and the mathematical formulas behind that, that's not a problem as such. I can optimize the weights for the second matrix of weights, and the cost is indeed being reduced over time. But I'm not able to keep propagating that back because of the matrix structure.
  Knowing that I have find the derivative of the cost w.r.t. the weights:
  d(cost) / d(w) = d(cost) / d(f(z)) * d(f(z)) / d(z) * d(z) / d(w)
  
  (Being f the activation function and z the dot product plus the bias of a neuron)
  So I'm in the rightmost layer, with an output array of 10 elements. d(cost) / d(f(z)) is the subtraction of the observed an predicted values. I can multiply that by d(f(z)) / d(z), which is just f'(z) of the rightmost layer, also an unidimensional vector of 10 elements, having now d(cost) / d(z) calculated. Then, d(z)/d(w) is just the input to that layer, i.e. the output of the previous one, which is a vector of 30 elements. I figured that I can transpose d(cost) / d(z) so that T( d(cost) / d(z) ) * d(z) / d(w) gives me a matrix of (10, 30), which makes sense because it matches with the dimension of the rightmost weight matrix.
  But then I get stuck. The dimension of d(cost) / d(f(z)) is (1, 10), for d(f(z)) / d(z) is (1, 30) and for d(z) / d(w) is (1, 784). I don't know how to come up with a result for this.
  This is what I've coded so far. The incomplete part is the _propagate_back method. I'm not caring about the biases yet because I'm just stuck with the weights and first I want to figure this out. 
  import random
  from typing import List, Tuple
  
  import numpy as np
  from matplotlib import pyplot as plt
  
  import mnist_loader
  
  np.random.seed(42)
  
  NETWORK_LAYER_SIZES = [784, 30, 10]
  LEARNING_RATE = 0.05
  BATCH_SIZE = 20
  NUMBER_OF_EPOCHS = 5000
  
  
  def sigmoid(x):
      return 1 / (1 + np.exp(-x))
  
  
  def sigmoid_der(x):
      return sigmoid(x) * (1 - sigmoid(x))
  
  
  class Layer:
  
      def __init__(self, input_size: int, output_size: int):
          self.weights = np.random.uniform(-1, 1, [output_size, input_size])
          self.biases = np.random.uniform(-1, 1, [output_size])
          self.z = np.zeros(output_size)
          self.a = np.zeros(output_size)
          self.dz = np.zeros(output_size)
  
      def feed_forward(self, input_data: np.ndarray):
          input_data_t = np.atleast_2d(input_data).T
          dot_product = self.weights.dot(input_data_t).T[0]
          self.z = dot_product + self.biases
          self.a = sigmoid(self.z)
          self.dz = sigmoid_der(self.z)
  
  
  class Network:
  
      def __init__(self, layer_sizes: List[int], X_train: np.ndarray, y_train: np.ndarray):
          self.layers = [
              Layer(input_size, output_size)
              for input_size, output_size
              in zip(layer_sizes[0:], layer_sizes[1:])
          ]
          self.X_train = X_train
          self.y_train = y_train
  
      @property
      def predicted(self) -> np.ndarray:
          return self.layers[-1].a
  
      def _normalize_y(self, y: int) -> np.ndarray:
          output_layer_size = len(self.predicted)
          normalized_y = np.zeros(output_layer_size)
          normalized_y[y] = 1.
  
          return normalized_y
  
      def _calculate_cost(self, y_observed: np.ndarray) -> int:
          y_observed = self._normalize_y(y_observed)
          y_predicted = self.layers[-1].a
  
          squared_difference = (y_predicted - y_observed) ** 2
  
          return np.sum(squared_difference)
  
      def _get_training_batches(self, X_train: np.ndarray, y_train: np.ndarray) -> Tuple[np.ndarray, np.ndarray]:
          train_batch_indexes = random.sample(range(len(X_train)), BATCH_SIZE)
  
          return X_train[train_batch_indexes], y_train[train_batch_indexes]
  
      def _feed_forward(self, input_data: np.ndarray):
          for layer in self.layers:
              layer.feed_forward(input_data)
              input_data = layer.a
  
      def _propagate_back(self, X: np.ndarray, y_observed: int):
          """
          der(cost) / der(weight) = der(cost) / der(predicted) * der(predicted) / der(z) * der(z) / der(weight)
          """
          y_observed = self._normalize_y(y_observed)
          d_cost_d_pred = self.predicted - y_observed
  
          hidden_layer = self.layers[0]
          output_layer = self.layers[1]
  
          # Output layer weights
          d_pred_d_z = output_layer.dz
          d_z_d_weight = hidden_layer.a  # Input to the current layer, i.e. the output from the previous one
  
          d_cost_d_z = d_cost_d_pred * d_pred_d_z
          d_cost_d_weight = np.atleast_2d(d_cost_d_z).T * np.atleast_2d(d_z_d_weight)
  
          output_layer.weights -= LEARNING_RATE * d_cost_d_weight
  
          # Hidden layer weights
          d_pred_d_z = hidden_layer.dz
          d_z_d_weight = X
  
          # ...
  
      def train(self, X_train: np.ndarray, y_train: np.ndarray):
          X_train_batch, y_train_batch = self._get_training_batches(X_train, y_train)
          cost_over_epoch = []
  
          for epoch_number in range(NUMBER_OF_EPOCHS):
              X_train_batch, y_train_batch = self._get_training_batches(X_train, y_train)
  
              cost = 0
              for X_sample, y_observed in zip(X_train_batch, y_train_batch):
                  self._feed_forward(X_sample)
                  cost += self._calculate_cost(y_observed)
                  self._propagate_back(X_sample, y_observed)
  
              cost_over_epoch.append(cost / BATCH_SIZE)
  
          plt.plot(cost_over_epoch)
          plt.ylabel('Cost')
          plt.xlabel('Epoch')
          plt.savefig('cost_over_epoch.png')
  
  
  training_data, validation_data, test_data = mnist_loader.load_data()
  X_train, y_train = training_data[0], training_data[1]
  
  network = Network(NETWORK_LAYER_SIZES, training_data[0], training_data[1])
  network.train(X_train, y_train)
  
  This is the code for mnist_loader, in case someone wanted to reproduce the example:
  import pickle
  import gzip
  
  
  def load_data():
      f = gzip.open('data/mnist.pkl.gz', 'rb')
      training_data, validation_data, test_data = pickle.load(f, encoding='latin-1')
      f.close()
  
      return training_data, validation_data, test_data
  
  **Tags:** neural-networks, backpropagation
  
  ---
  
  ## 1 answers
  
  ### Answer by username_1 on 2020-04-08. Score: 2 | Accepted answer
  
  I'm not familiar with python I'm afraid, but I'll give a go at presenting some of the maths...
  
  In order to perform the **back propagation** using matrices, matrix [transpose](https://en.wikipedia.org/wiki/Transpose) is used for different sized layers. Also note that when using matrices for this, we need to distinguish different types of matrix multiplication, namely matrix [product](https://en.wikipedia.org/wiki/Matrix_multiplication), and the [Hadamard product](https://en.wikipedia.org/wiki/Hadamard_product_(matrices)) which operates differently (the latter of which is denoted by a circle with a dot in the centre).
  
  Note the **back propagation** formulas:
  
  (EQ1)
  \begin{equation\*}\delta ^{l} = (w^{l+1})^{T} \delta ^{l+1} \odot \sigma {}' (z ^{l})\end{equation\*}
  
  (EQ2)
  \begin{equation\*}\frac{\partial E}{\partial w}=a^{l-1}\delta ^{l}\end{equation\*}
  
  As you can see the transpose of the weight matrix is used to multiply against the delta of the layer.
  
  As an example, consider a simple network with an input layer of 3 neurons, and an output layer of 2 neurons.
  
  The *unactivated* **feed forward** output is given by...
  
  \begin{equation\*}z\_{1} = h\_{1} w\_{1} + h\_{2} w\_{3} + h\_{3} w\_{5}\end{equation\*}
  \begin{equation\*}z\_{2} = h\_{1} w\_{2} + h\_{2} w\_{4} + h\_{3} w\_{6}\end{equation\*}
  
  Which can be represented in matrix form as...
  
  \begin{equation\*}
  z
  =
  \begin{bmatrix}z\_{1}\\z\_{2}\end{bmatrix}
  =
  \begin{bmatrix}h\_{1} w\_{1} + h\_{2} w\_{3} + h\_{3} w\_{5}\\h\_{1} w\_{2} + h\_{2} w\_{4} + h\_{3} w\_{6}\end{bmatrix}
  \end{equation\*}
  
  (EQ3)
  \begin{equation\*}
  =
  \begin{bmatrix}w\_{1}&w\_{3}&w\_{5}\\w\_{2}&w\_{4}&w\_{6}\end{bmatrix}
  \begin{bmatrix}h\_{1}\\h\_{2}\\h\_{3}\end{bmatrix}
  \end{equation\*}
  
  Which allows us to **forward propagate** from a layer of 3 neurons, to a layer of 2 neurons (via matrix multiplication).
  
  For the **back propagation**, we need the *error* of each neruon...
  
  N.B *Cost* is a metric used to denote the *error* of the entire network and depends on your cost function. It is usually the mean of the sum of the *errors* of each nuerons, but of course depends on your cost function used.
  
  e.g for [MSE](https://en.wikipedia.org/wiki/Mean_squared_error) is...
  \begin{equation\*}
  C\_{MSE}=\frac{1}{N}\sum (o\_{n}-t\_{n})^{2}
  \end{equation\*}
  
  We are interested in the derivative of the *error* for each neuron (not the *cost*), which by the chain rule is...
  
  \begin{equation\*}
  \frac{\partial E}{\partial z}
  =
  \frac{\partial E}{\partial o}
  \frac{\partial o}{\partial z}
  \end{equation\*}
  
  Expressed in matrix form...
  
  \begin{equation\*}
  \frac{\partial E}{\partial z}
  =
  \begin{bmatrix}
  \frac{\partial E}{\partial z\_{1}}\\
  \frac{\partial E}{\partial z\_{2}}
  \end{bmatrix}
  =
  \begin{bmatrix}
  \frac{\partial E}{\partial o\_{1}}
  \frac{\partial o\_{1}}{\partial z\_{1}}\\
  \frac{\partial E}{\partial o\_{2}}
  \frac{\partial o\_{2}}{\partial z\_{2}}
  \end{bmatrix}
  \end{equation\*}
  
  (EQ4)
  \begin{equation\*}
  =
  \begin{bmatrix}
  \frac{\partial E}{\partial o\_{1}}
  \\
  \frac{\partial E}{\partial o\_{2}}
  \end{bmatrix}
  \odot
  \begin{bmatrix}
  \frac{\partial o\_{1}}{\partial z\_{1}}
  \\
  \frac{\partial o\_{2}}{\partial z\_{2}}
  \end{bmatrix}
  \end{equation\*}
  
  Note the use of the hadamard product here. In fact, given these are simply vectors, vector dot product would work, but hadamard is used as this becomes important later when using this expression with the matrix equations as we want to distinguish the use of hadamard product as opposed to matrix product.
  
  We start off with our first delta error, which for the first layer **back propagation** is...
  
  \begin{equation\*}
  \delta ^{L} = \frac{\partial E}{\partial o}
  \end{equation\*}
  
  And then we want to calculate the next delta error using the formula (EQ1)...
  
  \begin{equation\*}
  \delta ^{l} = (w^{l+1})^{T} \delta ^{l+1} \odot \sigma {}' (z ^{l}) = \frac{\partial E}{\partial h}
  \end{equation\*}
  
  Explicitly, the equations are...
  
  \begin{equation\*}\frac{\partial E}{\partial h\_{1}}=\frac{\partial E}{\partial z\_{1}} w\_{1} + \frac{\partial E}{\partial z\_{2}} w\_{2}\end{equation\*}
  \begin{equation\*}\frac{\partial E}{\partial h\_{2}}=\frac{\partial E}{\partial z\_{1}} w\_{3} + \frac{\partial E}{\partial z\_{2}} w\_{4}\end{equation\*}
  \begin{equation\*}\frac{\partial E}{\partial h\_{3}}=\frac{\partial E}{\partial z\_{1}} w\_{5} + \frac{\partial E}{\partial z\_{2}} w\_{6}\end{equation\*}
  
  Also given the transpose of the weight matrix which allows us to **back propagate** from a layer of 2 neurons to a layer of 3 neurons (via matrix multiplication) ...
  
  \begin{equation\*}
  \left (w \right )^{T}
  =
  \left ( \begin{bmatrix}
  w\_{1} & w\_{3} & w\_{5}\\
  w\_{2} & w\_{4} & w\_{6}
  \end{bmatrix}
  \right )^{T}
  =
  \begin{bmatrix}
  w\_{1} & w\_{2}\\
  w\_{3} & w\_{4}\\
  w\_{5} & w\_{6}
  \end{bmatrix}
  \end{equation\*}
  
  So in a similar way to how we represented the forward pass (EQ3), this can be represented in matrix form...
  
  \begin{equation\*}
  \frac{\partial E}{\partial h}
  =
  \begin{bmatrix}
  \frac{\partial E}{\partial h\_{1}}\\
  \frac{\partial E}{\partial h\_{2}}\\
  \frac{\partial E}{\partial h\_{3}}
  \end{bmatrix}
  =
  \begin{bmatrix}
  \frac{\partial E}{\partial z\_{1}} w\_{1} + \frac{\partial E}{\partial z\_{2}} w\_{2}
  \\
  \frac{\partial E}{\partial z\_{1}} w\_{3} + \frac{\partial E}{\partial z\_{2}} w\_{4}
  \\
  \frac{\partial E}{\partial z\_{1}} w\_{5} + \frac{\partial E}{\partial z\_{2}} w\_{6}
  \end{bmatrix}
  \end{equation\*}
  
  (EQ5)
  \begin{equation\*}
  =
  \begin{bmatrix}
  w\_{1} & w\_{2}\\
  w\_{3} & w\_{4}\\
  w\_{5} & w\_{6}
  \end{bmatrix}
  \begin{bmatrix}
  \frac{\partial E}{\partial z\_{1}}
  \\
  \frac{\partial E}{\partial z\_{2}}
  \end{bmatrix}
  \end{equation\*}
  
  And then plugging the hadamard version of the delta (EQ4) into this, we get...
  
  \begin{equation\*}
  \begin{bmatrix}
  w\_{1} & w\_{2}\\
  w\_{3} & w\_{4}\\
  w\_{5} & w\_{6}
  \end{bmatrix}
  \begin{bmatrix}
  \frac{\partial E}{\partial o\_{1}}
  \\
  \frac{\partial E}{\partial o\_{2}}
  \end{bmatrix}
  \odot
  \begin{bmatrix}
  \frac{\partial o\_{1}}{\partial z\_{1}}
  \\
  \frac{\partial o\_{2}}{\partial z\_{2}}
  \end{bmatrix}
  \end{equation\*}
  
  aka (EQ1) ...
  
  \begin{equation\*}(w^{l+1})^{T} \delta ^{l+1} \odot \sigma {}' (z ^{l})\end{equation\*}
  
  And thus we have **back propagated** from a layer of 2 neurons, to a layer of 3 neurons (via matrix multiplication) thanks to transpose.
  
  For completeness... the other aspect of the back propagation to uses matrices, is the delta weight matrix....
  
  \begin{equation\*}
  \frac{\partial E}{\partial w}
  =
  \begin{bmatrix}
  \frac{\partial E}{\partial w\_{1}} & \frac{\partial E}{\partial w\_{3}} & \frac{\partial E}{\partial w\_{5}}\\
  \frac{\partial E}{\partial w\_{2}} & \frac{\partial E}{\partial w\_{4}} & \frac{\partial E}{\partial w\_{6}}
  \end{bmatrix}
  \end{equation\*}
  
  As mentioned before, you need to cache weight matrix and the activated layer output of a forward pass of the network.
  
  In a similar vein to (EQ3)...
  
  We explcitily have the equations...
  
  \begin{equation\*}\frac{\partial E}{\partial w\_{1}}=\frac{\partial E}{\partial z\_{1}}h\_{1}\end{equation\*}
  \begin{equation\*}\frac{\partial E}{\partial w\_{2}}=\frac{\partial E}{\partial z\_{2}}h\_{1}\end{equation\*}
  \begin{equation\*}\frac{\partial E}{\partial w\_{3}}=\frac{\partial E}{\partial z\_{1}}h\_{2}\end{equation\*}
  \begin{equation\*}\frac{\partial E}{\partial w\_{4}}=\frac{\partial E}{\partial z\_{2}}h\_{2}\end{equation\*}
  \begin{equation\*}\frac{\partial E}{\partial w\_{5}}=\frac{\partial E}{\partial z\_{1}}h\_{3}\end{equation\*}
  \begin{equation\*}\frac{\partial E}{\partial w\_{6}}=\frac{\partial E}{\partial z\_{2}}h\_{3}\end{equation\*}
  
  Note the use of h1, h2 and h3, which are the activated outputs of the previous layer.. (or in the case of our example these are the inputs).
  
  Which we represent in matrix form...
  
  \begin{equation\*}
  \frac{\partial E}{\partial w}
  =
  \begin{bmatrix}
  \frac{\partial E}{\partial z\_{1}}h\_{1} & \frac{\partial E}{\partial z\_{1}}h\_{2} & \frac{\partial E}{\partial z\_{1}}h\_{3}\\
  \frac{\partial E}{\partial z\_{2}}h\_{1} & \frac{\partial E}{\partial z\_{2}}h\_{2} & \frac{\partial E}{\partial z\_{2}}h\_{3}
  \end{bmatrix}
  \end{equation\*}
  (EQ6)
  \begin{equation\*}
  =
  \begin{bmatrix}
  \frac{\partial E}{\partial z\_{1}} \\ \frac{\partial E}{\partial z\_{2}}
  \end{bmatrix}
  \begin{bmatrix}
  h\_{1} &
  h\_{2} &
  h\_{3}
  \end{bmatrix}
  \end{equation\*}
  
  Which just so happens to be (EQ2) ...
  
  \begin{equation\*}\frac{\partial E}{\partial w}=\delta ^{l}a^{l-1}\end{equation\*}
  
  :)
  
  Since the delta weight matrix and the original weight matrix are the same dimensions, it is trivia to apply the learning rate...
  
  \begin{equation\*}
  w =
  \begin{bmatrix}
  w\_{1} - \alpha \frac{\partial E}{\partial w\_{1}} & w\_{3} - \alpha \frac{\partial E}{\partial w\_{3}} & w\_{5} - \alpha \frac{\partial E}{\partial w\_{5}}\\
  w\_{2} - \alpha \frac{\partial E}{\partial w\_{2}} & w\_{4} - \alpha \frac{\partial E}{\partial w\_{4}} & w\_{6} - \alpha \frac{\partial E}{\partial w\_{6}}
  \end{bmatrix}
  \end{equation\*}
  
  \begin{equation\*}
  =
  \begin{bmatrix}w\_{1}&w\_{3}&w\_{5}\\w\_{2}&w\_{4}&w\_{6}\end{bmatrix}
  -
  \alpha
  \begin{bmatrix}
  \frac{\partial E}{\partial w\_{1}} & \frac{\partial E}{\partial w\_{3}} & \frac{\partial E}{\partial w\_{5}}\\
  \frac{\partial E}{\partial w\_{2}} & \frac{\partial E}{\partial w\_{4}} & \frac{\partial E}{\partial w\_{6}}
  \end{bmatrix}
  \end{equation\*}
  
  \begin{equation\*} = w - \alpha \frac{\partial E}{\partial w}\end{equation\*}
  
  I'll omit the equations for the bias, as these are easy to put in and don't require any other equations other than the ones presented here.
  So in summary, the equations you want are (EQ3), (EQ5) and (EQ6) to use matrices, and the ability to change move between layers of differing number of neurons is the matrix transpose. Hope this helps, and let me know if you want me to expand on anything.
  
  It might be also important to note that you appear to be using MSE as the cost. This is perhaps not optimal for the MNIST data set which is a clasification. A better cost function to use would be [cross entropy](https://www.google.com/search?client=firefox-b-d&q=cross%20entropy).
  
  Happy propagating!
  
  
  ---
  
  