The Math inside Neural Networks

10 min read

This article is intended for those who have an undergraduate level understanding of linear algebra, calculus, as well as some introductory knowledge of vanilla neural networks. I’ve tried my best to keep things simple without turning this into an ‘Introduction to Neural Networks’ article.

I. Introduction

There are three primary processes that go on inside a vanilla multilayered neural network in a supervised learning setup:

  1. Feedforward - The network guesses an answer that is almost always incorrect.
  2. Error Calculation - The network calculates the amount of error between its answer and the correct answer.
  3. Backpropagation - The network updates its parameters (or weights) to minimise the error.

We will analyse the mathematics of all three of these processes in order to reach an understanding of the overall system.

II. Feedforward

A neural network converts a vector X in an input space Rx to a vector Y in an output space Ry. Since a vector is being transformed from one space to another, we know a matrix is involved. This matrix is actually the weights of that particular layer. Thus, each layer after the input layer represents a vector transformation with a corresponding weight matrix.

We take a closer look by analysing the set of operations inside the first hidden layer. Let’s say the input layer has a total ‘N’ number of neurons and the first hidden layer has a total ‘M’ number of neurons.

W1,1W1,2W2,1W3,3W3,mWn,4Wn,mX₁X₂X₃···XₙH₁H₂H₃H₄···Hₘvector Xvector Hmatrix WN elementsM elementsN x MconnectionsInput Layerof Neural NetworkFirst Hidden Layerof Neural Network

Mathematically, this is just standard matrix multiplication:

X⃗=(x1x2⋯xn),W⃗=(w1,1w1,2⋯w1,mw2,1w2,2⋯w2,m⋮⋮⋱⋮wn,1wn,2⋯wn,m),H⃗=(h1h2⋯hm)\vec{X} = \begin{pmatrix} x_1 & x_2 & \cdots & x_n \end{pmatrix}, \quad \vec{W} = \begin{pmatrix} w_{1,1} & w_{1,2} & \cdots & w_{1,m} \\ w_{2,1} & w_{2,2} & \cdots & w_{2,m} \\ \vdots & \vdots & \ddots & \vdots \\ w_{n,1} & w_{n,2} & \cdots & w_{n,m} \end{pmatrix}, \quad \vec{H} = \begin{pmatrix} h_1 & h_2 & \cdots & h_m \end{pmatrix}
X⃗[1×N]×W⃗[N×M]=H⃗[1×M]\vec{X}_{[1 \times N]} \times \vec{W}_{[N \times M]} = \vec{H}_{[1 \times M]}

Note: real networks also add a bias vector to each layer’s weighted sum (H = X·W + b). We omit bias terms throughout because b has the same dimensions as H

We now have a system that converts vectors from one vector space to another. If this was the only operation in each layer, the entire neural network could be represented as a single weight matrix by post-multiplying the weight matrices of each layer. Of course, that is not what happens!

The reason we have multiple layers in the network is because each layer has an activation function. An activation function squishes the outputs of a layer into a certain range. Because there are no restrictions on how large the values inside X or W can be, after a few multiplications the magnitude of values inside H can increase to the point where computers run out of space to store those values. Besides, we don’t need very large values for tasks like classification or probabilistic prediction.

A⃗[1×M]=f(H⃗[1×M])=(f(h1)f(h2)⋯f(hm))\vec{A}_{[1 \times M]} = f(\vec{H}_{[1 \times M]}) = \begin{pmatrix} f(h_1) & f(h_2) & \cdots & f(h_m) \end{pmatrix}

We thus obtain the final activation values of the first layer after passing H through an activation function. Generally, we want the following characteristics in an activation function:

(i) It should be monotonic increasing (in other words, an increase in h should result in an increase in a).

(ii) It should be bounded to a certain interval. This could either be a lower bound, upper bound, or both.

(iii) It should be differentiable. We will discuss why this is important during backpropagation.

Two functions that fit this requirement are the sigmoid function and the tanh function. Of course, there are more activation functions for other use cases but let’s consider these two for now. The graph below shows us how these functions constrain the input. The sigmoid function is in blue and the tanh function is in green.

sig(x)=ex1+ex,tanh(x)=ex−e−xex+e−xsig(x) = \frac{e^x}{1 + e^x}, \quad tanh(x) = \frac{e^x - e^{-x}}{e^x + e^{-x}}
-2-1012-2-112

This process of weight multiplication followed by activation is repeated for each layer, till we reach the final layer. If we want to express this in form of pseudocode, we can do it as follows:

def forwardPass(input):
    current_activation = input
    for layer in layers[1:]:
        weight_multiplication = current_activation.matrixMultiply(layer.weightMatrix)
        current_activation = weight_multiplication.activate()
    return current_activation

The layers include the output layer of the network. This means at the end of the forward pass we have the output of the entire network. The number of elements in the output will be equal to the number of neurons in the output layer. At this point our network has guessed an answer, which is most likely wrong since it is an untrained network. Up next, we will describe how the network quantifies its error.

III. Error Calculation

In supervised learning we have labelled outputs for our inputs during training. We call these labels ‘targets’. To calculate the error in the network’s output Y we need to compare it to the corresponding target T. But how exactly do we determine the error. Do we simply subtract Y from T?

Whatever the error is, we ideally want to minimise it. Error is also referred to as the ‘Loss’ of the network. Whatever this loss function is, it should be as low as possible for maximum accuracy of the network. Let us plot the target T and the output Y on a graph.

Output YTarget T

The red line on the graph denotes y = t, which is the optimal condition. The blue point denotes any point (y0,t0) for any input x0. What we really want is to minimise the distance between the line (y = t) and the point (y0,t0). This gives us a good starting point for defining the loss function.

The distance between any point (y0,t0) and the line (y=t) can be described using the formula:

d=(y0−t0)22d = \sqrt{\frac{(y_0 - t_0)^2}{2}}

To minimise d, we minimise |y0 - t0|. Assuming there are K elements in the output layer, we define the loss function E as:

E⃗[1×K]=((y1−t1)2(y2−t2)2⋯(yk−tk)2)\vec{E}_{[1 \times K]} = \begin{pmatrix} (y_1 - t_1)^2 & (y_2 - t_2)^2 & \cdots & (y_k - t_k)^2 \end{pmatrix}

We have described E mathematically. That’s nice, but we took all this trouble to answer the question - How much should the weights of the network be adjusted in a way that reduces E? Introducing Backpropagation.

IV. Backpropagation

Before we adjust the weights of the entire network, we need to determine the error for each layer. The loss function E denotes the error for the network, that’s true. But more specifically, it is the error of the output layer.

Backpropagation is a popular approach to solve this problem. Starting from the output layer, the error propagates backwards until the weights of all layers have been updated. Let us first look at how we update the weights in the output layer.

To do that we first express E as a function of the weights of the output layer. If Wyz denotes the weight matrix of the output layer, and A denotes the activations of the second last layer:

A note on notation: we name the last three layers x, y and z. The output layer therefore has Z neurons, which means the K from Section III equals Z from here on.

W1,1W1,2W2,1W4,3W4,zWy,4Wy,zAy₁Ay₂Ay₃Ay₄···AyyZ₁Z₂Z₃Z₄···Zzvector Ayvector Zmatrix WY elementsZ elementsY x ZconnectionsSecond Last Layerof the NetworkFinal Layer ofthe Network
Y⃗=f(A⃗y×W⃗yz)=f(Z⃗)=(y1y2⋯yz)\vec{Y} = f(\vec{A}_y \times \vec{W}_{yz}) = f(\vec{Z}) = \begin{pmatrix} y_1 & y_2 & \cdots & y_z \end{pmatrix}
A⃗y=(a1a2⋯ay),W⃗yz=(w1,1w1,2⋯w1,zw2,1w2,2⋯w2,z⋮⋮⋱⋮wy,1wy,2⋯wy,z)\vec{A}_y = \begin{pmatrix} a_1 & a_2 & \cdots & a_y \end{pmatrix}, \quad \vec{W}_{yz} = \begin{pmatrix} w_{1,1} & w_{1,2} & \cdots & w_{1,z} \\ w_{2,1} & w_{2,2} & \cdots & w_{2,z} \\ \vdots & \vdots & \ddots & \vdots \\ w_{y,1} & w_{y,2} & \cdots & w_{y,z} \end{pmatrix}

Next we calculate the gradient of E with respect to Wyz using the chain rule:

∂E⃗∂W⃗yz=∂E⃗∂Y⃗×∂Y⃗∂Z⃗×∂Z⃗∂W⃗yz\frac{\partial \vec{E}}{\partial \vec{W}_{yz}} = \frac{\partial \vec{E}}{\partial \vec{Y}} \times \frac{\partial \vec{Y}}{\partial \vec{Z}} \times \frac{\partial \vec{Z}}{\partial \vec{W}_{yz}}

We can resolve these three terms individually very easily. Let’s start with dE/dY:

∂E⃗∂Y⃗=(∂e1∂y1∂e2∂y2⋯∂ez∂yz)=(∂(y1−t1)2∂y1∂(y2−t2)2∂y2⋯∂(yz−tz)2∂yz)\frac{\partial \vec{E}}{\partial \vec{Y}} = \begin{pmatrix} \frac{\partial e_1}{\partial y_1} & \frac{\partial e_2}{\partial y_2} & \cdots & \frac{\partial e_z}{\partial y_z} \end{pmatrix} = \begin{pmatrix} \frac{\partial (y_1 - t_1)^2}{\partial y_1} & \frac{\partial (y_2 - t_2)^2}{\partial y_2} & \cdots & \frac{\partial (y_z - t_z)^2}{\partial y_z} \end{pmatrix}
=(2(y1−t1)2(y2−t2)⋯2(yz−tz))=2(Y⃗−T⃗)= \begin{pmatrix} 2(y_1 - t_1) & 2(y_2 - t_2) & \cdots & 2(y_z - t_z) \end{pmatrix} = 2(\vec{Y} - \vec{T})

To calculate dY/dZ, recall that Y=f(Z) where f is the activation function. This is why we insist on having a differentiable activation function.

∂Y⃗∂Z⃗=∂f(Z⃗)∂Z⃗=f′(Z⃗)\frac{\partial \vec{Y}}{\partial \vec{Z}} = \frac{\partial f(\vec{Z})}{\partial \vec{Z}} = f'(\vec{Z})

We have a trick to quickly compute the derivatives. The derivatives of the sigmoid and tanh functions can be represented as follows:

∂ sig(x)∂x=sig(x)×(1−sig(x))\frac{\partial \, sig(x)}{\partial x} = sig(x) \times (1 - sig(x))
∂ tanh(x)∂x=1−tanh(x)2\frac{\partial \, tanh(x)}{\partial x} = 1 - tanh(x)^2

dZ/dWyz is actually just Ay

∂Z⃗∂W⃗yz=A⃗y\frac{\partial \vec{Z}}{\partial \vec{W}_{yz}} = \vec{A}_y

We can rewrite the dE/dWyz as:

∂E⃗∂W⃗yz=2(Y⃗−T⃗)×f′(Z⃗)×A⃗y\frac{\partial \vec{E}}{\partial \vec{W}_{yz}} = 2(\vec{Y} - \vec{T}) \times f'(\vec{Z}) \times \vec{A}_y

Treat the multiplications above as element-wise for now. The dimensions do not strictly line up as matrix products yet. We will make this equation dimensionally rigorous shortly.

We can also determine the gradient for each connection Wij in the matrix Wyz. Also note that we can safely ignore the ‘2’ in the equation because we are not concerned with the coefficients.

∂E⃗∂wij=2(yj−tj)×f′(zj)×ai=δj×ai\frac{\partial \vec{E}}{\partial w_{ij}} = 2(y_j - t_j) \times f'(z_j) \times a_i = \delta_j \times a_i

Thus the new value of any connection Wij in the matrix Wyz is written as follows. The symbol gamma (γ) represents the learning rate and absorbs all constant coefficients. This is also why we discarded the ‘2’ previously.

wijnew=wijold−γ∂E⃗∂wij=wijold−γδjaiw_{ij_{new}} = w_{ij_{old}} - \gamma \frac{\partial \vec{E}}{\partial w_{ij}} = w_{ij_{old}} - \gamma \delta_j a_i

We saw how a given weight Wij in the output layer can be updated. But to understand backpropagation for hidden layers we need to understand the matrices at play. In order to satisfy the dimensionality of the vectors, we take the transpose on the RHS and rewrite the matrix equation as follows:

∂E⃗∂W⃗yz=A⃗y T⋅(Y⃗−T⃗)⋅Dz\frac{\partial \vec{E}}{\partial \vec{W}_{yz}} = \vec{A}_y^{\,T} \cdot (\vec{Y} - \vec{T}) \cdot D_z
f′(Z⃗)=Dz=(f′(z1)0⋯00f′(z2)⋯0⋮⋮⋱⋮00⋯f′(zz))f'(\vec{Z}) = D_z = \begin{pmatrix} f'(z_1) & 0 & \cdots & 0 \\ 0 & f'(z_2) & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & f'(z_z) \end{pmatrix}
(Y⃗−T⃗)=e⃗z=((y1−t1)(y2−t2)⋯(yz−tz))(\vec{Y} - \vec{T}) = \vec{e}_z = \begin{pmatrix} (y_1 - t_1) & (y_2 - t_2) & \cdots & (y_z - t_z) \end{pmatrix}

Here Dz is a ZxZ diagonal matrix containing the derivatives of the output activations.

∂E⃗∂W⃗yz=A⃗y T⋅e⃗z⋅Dz=A⃗y T⋅δ⃗z\frac{\partial \vec{E}}{\partial \vec{W}_{yz}} = \vec{A}_y^{\,T} \cdot \vec{e}_z \cdot D_z = \vec{A}_y^{\,T} \cdot \vec{\delta}_z
W⃗yznew=W⃗yzold−γ A⃗y T⋅δ⃗z\vec{W}_{yz_{new}} = \vec{W}_{yz_{old}} - \gamma \, \vec{A}_y^{\,T} \cdot \vec{\delta}_z

This is how computers calculate gradients. We finally calculate the gradients for the hidden layer with the weight matrix Wxy.

δ⃗y=δ⃗z⋅W⃗yz T⋅Dy\vec{\delta}_y = \vec{\delta}_z \cdot \vec{W}_{yz}^{\,T} \cdot D_y
W⃗xynew=W⃗xyold−γ A⃗x T⋅δ⃗y\vec{W}_{xy_{new}} = \vec{W}_{xy_{old}} - \gamma \, \vec{A}_x^{\,T} \cdot \vec{\delta}_y

The above process is repeated for the previous layers until we finally update the weights of the first hidden layer. Backpropagation can be summarised by pseudocode as follows:

def backPropagate(error):
    current_layer = layers[-1]
    previous_layer = layers[-2]
    differential_matrix = current_layer.getDifferential()
    delta = error.matrixMultiply(differential_matrix)
    current_layer.weightMatrix -= learning_rate * (previous_layer.activations.transpose().matrixMultiply(delta))
    for i in range(1, len(layers) - 1):
        next_layer = layers[-i]
        current_layer = layers[-1 - i]
        previous_layer = layers[-2 - i]
        differential_matrix = current_layer.getDifferential()
        delta = delta.matrixMultiply(next_layer.weightMatrix.transpose()).matrixMultiply(differential_matrix)
        current_layer.weightMatrix -= learning_rate * (previous_layer.activations.transpose().matrixMultiply(delta))

V. Conclusion

We discussed how neural networks are carefully coordinated matrix operations. If you found this article useful and want to study the concepts behind Neural Networks in detail, you may read ‘Neural Networks – A Systematic Introduction’ by Raul Rojas.