Project intention

Why rebuild an MLP without deep-learning libraries?

The intention of this project is simple: for friends who have just started learning neural networks, it is often difficult to see what forward propagation and backward propagation are really doing. A framework can train a model in a few lines, but it can also hide the chain rule, the gradients, the shapes of matrices, and the reason parameters change.

So this project implements a multilayer perceptron with explicit derivation and explicit Python code. The network uses Fashion-MNIST as the dataset, and contains two derivation paths: sigmoid with mean squared error, and softmax with cross entropy.

Build the network once by hand, and the equations stop being symbols on a page. They become the code path of learning itself.

1 · Introduction / 2.1 · Forward propagation

Forward propagation turns input information into a prediction.

This is an introduction of the theory of MLP, illustrated by derivation with both MSE using sigmoid and cross entropy using softmax. Forward propagation can be understood as the process of network prediction: it transforms input information into classification results.

Suppose the input layer information is \([x_1, x_2, x_3]\). For layer \(l\), \(L_l\) represents all neurons of the layer, the output is \(y_l\), the output of the \(j\)-th node is \(y_l^{(j)}\), the input of that node is \(u_l^{(j)}\), the weight matrix connecting layer \(l\) and layer \(l-1\) is \(W_l\), and the weight from the \(i\)-th node of the previous layer to the \(j\)-th node of layer \(l\) is \(w_l^{ji}\).

Equation 1 \[ \begin{cases} y_2^{(1)} = f(u_2^{(1)}) = f\left(\sum_{i=1}^{n} w_2^{1i}x_i + b_2^{(1)}\right) \\ y_2^{(2)} = f(u_2^{(2)}) = f\left(\sum_{i=1}^{n} w_2^{2i}x_i + b_2^{(2)}\right) \end{cases} \]

To make this easy to code, the first equation can be written as a matrix expression:

Equation 2 \[ y_2 = \begin{bmatrix} y_2^{(1)} \\ y_2^{(2)} \end{bmatrix} = f\left( \begin{bmatrix} w_2^{11} & w_2^{12} & w_2^{13} \\ w_2^{21} & w_2^{22} & w_2^{23} \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} + \begin{bmatrix} b_2^{(1)} \\ b_2^{(2)} \end{bmatrix} \right) = f(W_2X + b_2) \]

Extending the forward propagation calculation process of the second layer to any layer in the network gives:

Equation 3 \[ \begin{cases} y_l^{(j)} = f(u_l^{(j)}) \\ u_l^{(j)} = \sum_{i\in L_{l-1}} w_l^{ji}y_{l-1}^{(i)} + b_l^{(j)} \\ y_l = f(u_l) = f(W_ly_{l-1}+b_l) \end{cases} \]

Here, \(f(\cdot)\) is the activation function, and \(b_l^{(j)}\) is the bias of the \(j\)-th node in layer \(l\).

2.2 · Backward propagation

Backward propagation asks how each parameter changed the loss.

After the basic model is built, the process of training is to update the model parameters. Due to the multi-layer network structure, it is not possible to directly update the parameters of the hidden layer by the loss function. Instead, the loss is propagated from the top layer to the bottom layer to estimate the parameters.

Suppose there are multiple neurons in the output layer of the multilayer perceptron, and each neuron corresponds to a label. The input sample is \(x=[x_1,x_2,\ldots,x_n]\), and the label is \(t\). For the output layer, this note includes two loss functions:

Equation 4 · MSE \[ E_{MSE}=\frac{1}{2}\sum_{j\in L_k}\left(t^j-y_k^{(j)}\right)^2 \]
Equation 5 · Cross entropy \[ E_{Cross\ Entropy}=-\sum_j t^j\cdot \log(y_k^{(j)}) \]

The loss function determines how the output layer parameters are updated. In order to improve training accuracy, the loss function with sigmoid as the output activation is defined as MSE, while the loss function with softmax as the output activation is defined as cross entropy. To minimize the loss, gradient descent is used.

Equation 6 · Chain rule for each parameter \[ \begin{cases} \dfrac{\partial E}{\partial w_l^{(ji)}} = \dfrac{\partial E}{\partial y_l^{(j)}}\dfrac{\partial y_l^{(j)}}{\partial w_l^{(ji)}} = \dfrac{\partial E}{\partial y_l^{(j)}}\dfrac{\partial y_l^{(j)}}{\partial u_l^{(j)}}\dfrac{\partial u_l^{(j)}}{\partial w_l^{(ji)}} \\ \dfrac{\partial E}{\partial b_l^{(j)}} = \dfrac{\partial E}{\partial y_l^{(j)}}\dfrac{\partial y_l^{(j)}}{\partial b_l^{(j)}} = \dfrac{\partial E}{\partial y_l^{(j)}}\dfrac{\partial y_l^{(j)}}{\partial u_l^{(j)}}\dfrac{\partial u_l^{(j)}}{\partial b_l^{(j)}} \end{cases} \]

The following three derivatives are easy to derive from the definition of the forward pass:

Equation 7 \[ \begin{cases} \dfrac{\partial y_l^{(j)}}{\partial u_l^{(j)}} = f'(u_l^{(j)}) \\ \dfrac{\partial u_l^{(j)}}{\partial w_l^{(ji)}} = y_{l-1}^{(i)} \\ \dfrac{\partial u_l^{(j)}}{\partial b_l^{(j)}} = 1 \end{cases} \]

Therefore, the parameter gradients become:

Equation 8 \[ \begin{cases} \dfrac{\partial E}{\partial w_l^{(ji)}} = \dfrac{\partial E}{\partial y_l^{(j)}}\dfrac{\partial y_l^{(j)}}{\partial u_l^{(j)}}\dfrac{\partial u_l^{(j)}}{\partial w_l^{(ji)}} = \dfrac{\partial E}{\partial y_l^{(j)}}f'(u_l^{(j)})y_{l-1}^{(i)} \\ \dfrac{\partial E}{\partial b_l^{(j)}} = \dfrac{\partial E}{\partial y_l^{(j)}}\dfrac{\partial y_l^{(j)}}{\partial u_l^{(j)}}\dfrac{\partial u_l^{(j)}}{\partial b_l^{(j)}} = \dfrac{\partial E}{\partial y_l^{(j)}}f'(u_l^{(j)}) \end{cases} \]

Each node in the next layer is related to all nodes from the former layer. The loss function can therefore be considered as a function of the input of each node in the next layer:

Equation 9 \[ \begin{aligned} \dfrac{\partial E}{\partial y_l^{(j)}} &= \dfrac{\partial E\left(u_{l+1}^{(1)},u_{l+1}^{(2)},\ldots,u_{l+1}^{(k)},\ldots,u_{l+1}^{(K)}\right)}{\partial y_l^{(j)}} \\ &= \sum_{k\in L_{l+1}} \dfrac{\partial E}{\partial u_{l+1}^{(k)}}\dfrac{\partial u_{l+1}^{(k)}}{\partial y_l^{(j)}} \\ &= \sum_{k\in L_{l+1}} \dfrac{\partial E}{\partial y_{l+1}^{(k)}}\dfrac{\partial y_{l+1}^{(k)}}{\partial u_{l+1}^{(k)}}\dfrac{\partial u_{l+1}^{(k)}}{\partial y_l^{(j)}} \\ &= \sum_{k\in L_{l+1}} \dfrac{\partial E}{\partial y_{l+1}^{(k)}}\dfrac{\partial y_{l+1}^{(k)}}{\partial u_{l+1}^{(k)}}w_{l+1}^{kj} \end{aligned} \]

Node sensitivity

The gradient becomes easier once the node sensitivity is named.

For better understanding of calculation, define \(\delta=\partial E/\partial u\) as the rate of change of error to input. This is node sensitivity. Therefore, the node sensitivity of the \(j\)-th node from layer \(l\) is:

Equation 10 \[ \delta_l^{(j)} = \dfrac{\partial E}{\partial u_l^{(j)}} = \dfrac{\partial E}{\partial y_l^{(j)}}\dfrac{\partial y_l^{(j)}}{\partial u_l^{(j)}} = \dfrac{\partial E}{\partial y_l^{(j)}}f'(u_l^{(j)}) \]

Applying node sensitivity to the derivation of the loss function simplifies the upstream error term:

Equation 11 \[ \begin{aligned} \dfrac{\partial E}{\partial y_l^{(j)}} &= \sum_{k\in L_{l+1}} \dfrac{\partial E}{\partial y_{l+1}^{(k)}}\dfrac{\partial y_{l+1}^{(k)}}{\partial u_{l+1}^{(k)}}w_{l+1}^{kj} \\ &= \sum_{k\in L_{l+1}}\delta_{l+1}^{(k)}w_{l+1}^{kj} \end{aligned} \]

For the output layer, there is no next layer to consider. The output of the final layer is directly related to the error. Therefore, after multiplying by \(f'(u_l^{(j)})\), the sensitivity can be written as:

Equation 12 \[ \delta_l^{(j)} = \dfrac{\partial E}{\partial y_l^{(j)}}f'(u_l^{(j)}) = \begin{cases} f'(u_l^{(j)})\displaystyle\sum_{k\in L_{l+1}}\delta_{l+1}^{(k)}w_{l+1}^{kj}, & \text{if } l \text{ is hidden layer} \\ f'(u_l^{(j)})\dfrac{\partial E}{\partial y_l^{(j)}}, & \text{if } l \text{ is output layer} \end{cases} \]

The gradient of the loss function with respect to each parameter is:

Equation 13 \[ \begin{cases} \dfrac{\partial E}{\partial w_l^{(ji)}} = \dfrac{\partial E}{\partial u_l^{(j)}}\dfrac{\partial u_l^{(j)}}{\partial w_l^{(ji)}} = \delta_l^{(j)}y_{l-1}^{(i)} \\ \dfrac{\partial E}{\partial b_l^{(j)}} = \dfrac{\partial E}{\partial u_l^{(j)}}\dfrac{\partial u_l^{(j)}}{\partial b_l^{(j)}} = \delta_l^{(j)} \end{cases} \]

To make it suitable for every node and easy to code, the matrix form is \(\partial E/\partial W_l = \delta_l y_{l-1}^T\) and \(\partial E/\partial b_l = \delta_l\), where:

Equation 14 \[ \delta_l = \begin{cases} (W_{l+1}^{T}\delta_{l+1})\circ f'(u_l), & \text{if } l \text{ is hidden layer} \\ \dfrac{\partial E}{\partial y_l}\circ f'(u_l), & \text{if } l \text{ is output layer} \end{cases} \]

The symbol \(\circ\) represents multiplication of corresponding elements in a matrix or vector. The update equations of weights from each layer are:

Equation 15 \[ \begin{cases} W_l = W_l - \eta\dfrac{\partial E}{\partial W_l} = W_l - \eta\delta_l y_{l-1}^{T} \\ b_l = b_l - \eta\dfrac{\partial E}{\partial b_l} = b_l - \eta\delta_l \end{cases} \]

2.3 · Derivation of sigmoid and MSE

Sigmoid turns the derivative into a reusable local term.

The sigmoid function can be simplified as:

Equation 16 \[ y_l = \dfrac{1}{1+e^{-x}} = \dfrac{e^x}{e^x+1} = 1-(e^x+1)^{-1} \]

Therefore, the derivation of the sigmoid function is:

Equation 17 \[ (y_l)' = (-1)(-1)(e^x+1)^{-2}e^x = (e^x+1)^{-2}e^{-2x}e^x = (e^x+1)^{-1}\dfrac{1}{1+e^{-x}} = y_l(1-y_l) \]

Referring to the \(E_{MSE}\) equation above, the derivative of the MSE loss function is:

Equation 18 \[ \dfrac{\partial E}{\partial y_l^{(j)}} = y_l^{(j)} - t^{(j)} \]

In the sigmoid implementation of this project, the output layer sensitivity becomes the MSE derivative multiplied by the sigmoid derivative:

Sigmoid output sensitivity in code 
sigma_out = (layer_2 - target) * (layer_2 * (1 - layer_2))

2.3.1 · Derivation of softmax and cross entropy

Softmax plus cross entropy collapses into a clean error term.

The derivative of softmax and cross entropy can be considered as:

Setup \[ \dfrac{\partial E}{\partial u_i} = \sum_j\left(\dfrac{\partial E_j}{\partial y_j}\dfrac{\partial y_j}{\partial u_i}\right) \]

For the first derivative \(\partial E_j/\partial y_j\):

Equation 19 \[ \dfrac{\partial E_j}{\partial y_j} = \dfrac{\partial(-t_j\ln y_j)}{\partial y_j} = -t_j\dfrac{1}{y_j} \]

As for the softmax part, there are two conditions. When \(i=j\):

Equation 20 \[ \dfrac{\partial y_j}{\partial u_i} = \dfrac{\partial\left(\dfrac{e^{z_i}}{\sum_k e^{z_k}}\right)}{\partial z_i} = \dfrac{\sum_k e^{z_k}e^{z_i}-(e^{z_i})^2}{(\sum_k e^{z_k})^2} = \left(\dfrac{e^{z_i}}{\sum_k e^{z_k}}\right)\left(1-\dfrac{e^{z_i}}{\sum_k e^{z_k}}\right) = y_j(1-y_j) \]

However, if \(i\neq j\):

Equation 21 \[ \dfrac{\partial y_j}{\partial u_i} = \dfrac{\partial\left(\dfrac{e^{z_j}}{\sum_k e^{z_k}}\right)}{\partial z_i} = -e^{z_j}\left(\dfrac{1}{\sum_k e^{z_k}}\right)^2e^{z_i} = -y_i y_j \]

Combining the equations above, the final derivation of softmax and cross entropy is:

Equation 22 \[ \begin{aligned} \dfrac{\partial E}{\partial u_i} &= \sum_j\left(\dfrac{\partial E_j}{\partial y_j}\dfrac{\partial y_j}{\partial u_i}\right) \\ &= \sum_{j\neq i}\left(\dfrac{\partial E_j}{\partial y_j}\dfrac{\partial y_j}{\partial u_i}\right) + \sum_{j=i}\left(\dfrac{\partial E_j}{\partial y_j}\dfrac{\partial y_j}{\partial u_i}\right) \\ &= \left(\sum_{j\neq i}-t_j\dfrac{1}{y_j}(-y_iy_j)\right) + \left(-t_i\dfrac{1}{y_i}\right)(y_i(1-y_i)) \\ &= \left(\sum_{j\neq i}y_it_j\right)+(-t_i(1-y_i)) \\ &= \left(\sum_{j\neq i}y_it_j\right)+y_it_i-t_i \\ &= y_i\sum_j t_j-t_i \end{aligned} \]

There is only one target value \(t_i\) equal to 1, therefore \(\sum_j t_j=1\). In the end, the final derivative is:

Softmax · cross entropy result \[ \dfrac{\partial E}{\partial u_i}=y_i-t_i \]

In the softmax implementation of this project, this appears directly as the output sensitivity:

Softmax output sensitivity in code 
sigma_out = layer_2 - target

2.3.2 · Problems of coding sigmoid and softmax

The final step is numerical stability.

The program overflowed when designing sigmoid and softmax methods in this project. To solve these problems, the properties of the two activation functions were studied and the equations were rewritten into numerically safer forms.

The original sigmoid form is \(\sigma(z)=1/(1+e^{-z})\). When the input is negative with a large absolute value, computing \(e^{-z}\) directly can overflow. To solve this, the sigmoid code is modified as follows:

Stable sigmoid implementation 
def sigmoid_function(input):
    fz = []
    input = input.tolist()
    for num in input:
        if num >= 0:
            fz.append(1.0 / (1 + math.exp(-num)))
        else:
            fz.append(math.exp(num) / (1 + math.exp(num)))
    output = np.array(fz)
    return output

Similarly, the original form of softmax is:

Softmax original form \[ y_k = \dfrac{e^{a_k}}{\sum_{i=1}^{n}e^{a_i}} \]

If the input signal is too large, the program can overflow. To avoid this, the formula can be modified without changing the result of the activation function:

Softmax stable form \[ y_k = \dfrac{e^{a_k-C}}{\sum_{i=1}^{n}e^{a_i-C}} \]

Here, \(C\) is the maximum number in the input array. The corresponding code is:

Stable softmax implementation 
def softmax_function(input):
    max_input = np.max(input)
    input = np.exp(input - max_input)
    sum_input = np.sum(input)
    output = input / sum_input
    return output
View project repository Back to home