Multiclass Neural Networks

The focus of this article will be on the math behind simple neural networks and implementing the code in python from scratch. Also, the differences between binary and multiclass models will be highlighted.

This article tries to be as detailed as necessary in covering each derivation of the formulas, so that every line of code can be comprehended. This might not be a guide for complete beginners, as some basic knowledge in linear algebra is required.

primer

A neural network, simply put, a function $\phi$ that maps input to output. The dimensions of the input/output and its inner topology all are variable. The function itself is not explicitly defined beforehand, but inferred from the dataset.

$\begin{aligned} \R^{d_1} \; &\;\to \quad \R^{d_2} \\ x = \pmatrix{x_1 \\ x_2} &\overset \phi \longmapsto \pmatrix{a_1 \\ a_2 \\ a_3} = a\\ \end{aligned}$

Example of a neural network $\phi$ mapping input $x$ to output $a$.

The inputs and outputs will be column vectors, as I find that more intuitive when working with matrices. This is an arbitrary convention and can is easily changed by transposition.
We can also run a batch of samples through the neural network.

$\begin{aligned} \R^{d_1 \x N} \quad &\to \qquad \R^{d_2 \x N}\\ X = \bmatrix{ | & & |\\ x^{(1)} & \cdots & x^{(N)}\\ | & & |\\ } &\overset \phi \longmapsto \bmatrix{ | & & |\\ a^{(1)} & \cdots & a^{(N)}\\ | & & |\\ } = A \end{aligned}$

Different datapoints will be noted as upper indices in parentheses, the individual entries as lower indices.

The propagation from one layer to the next consists of an affine linear transformation, $x \mapsto Wx + b$, followed by a non-linear activation function $\s$, applied element-wise.

$\begin{aligned} \varphi(x) &= \sigma (Wx + b) = \s(z) \\ &= \veceval{\sigma(z_{\i})}{d_2} = \vecfill{a}{d_2} = a \in \R^{d_2} \end{aligned}$

$W$ is a matrix $\in \mathbb{R}^{d_2\times d_1}$ (here, a linear Transformation from $\R^{d_1}$ to $\R^{d_2}$) and the bias $b$ is a scalar $\in \R$.

A neural network may be composed of multiple layers. Each layer acting on the previous layer’s output or “neural activation”.

$\varphi _ l(a^{[l-1]}) = \s^{[l]}(W^{[l-1]}x + b^{[l]}) = a^{[l]}$

$\begin{aligned} \phi (x) &= (\varphi _ L \circ \cdots \circ \varphi _2 \circ \varphi _1) (x) \\ &= \s^{[L]}(W^{[L]}\s^{[L-1]}(\cdots\s^{[1]}(W^{[1]}x + b^{[1]})) + b^{[L]}) \\ &= a^{[L]} \\ \end{aligned}$

Upper indices in brackets denote the layer from $1$ to $L$.

It would seem more logical to start with the general multiclass model and then to study the binary model as a special case. Nonetheless, I will begin with the binary model, since calculating derivatives and predicting is simpler in this case.

binary classification

The regular multiclass model assumes two output dimensions, each representing the probability of the sample belonging to the respective class $0$ or $1$. Yet, as the probability of belonging to 'class $1$’ is $1$ minus the probability of belonging to 'class $0$': $a_2 = 1 - a_1$, and vice-versa, one output variable becomes redundant. Only one output ‘neuron’ $a$ is needed.
The estimated probability of whether the input belongs to each class, given our model $\phi$ then is

$\begin{aligned} \P(x \text{ of class } &1 \enspace|\; \phi)= a \\ \P(x \text{ of class } &0 \enspace | \;\phi)= 1 - a \\ \end{aligned}$

, where $a \in (0,1)$.

imports

Numpy will be used for matrix calculations and matplotlib for plots. Later, scikit will help in creating a toy dataset.

import numpy as np

initialization

The function nn_init() initializes a neural network with random weights and zero bias. It takes a parameter layer_dims, an array specifying the number of neurons for each layer. layer_dims = [2, 8, 1] for example. The first entry is the input dimension and reflects the dimension of the data, the last is the network’s output dimension and is equal to the number of classes in the general case.
The function returns a dictionary param containing the network’s information.
The weights connecting the input to the first layer are stored in param['W1'], the bias in param['b1'], and so on. The number of layers $L$ is len(layer_dims) - 1 (not counting the input layer).

def nn_init(layer_dims):

forward pass

A nonlinear activation function is applied to each layer after the affine linear transformation. The hidden layers will all be equipped with the rectified linear unit, or relu-activation function.

$relu(z) = z_{\{z> 0\}} = \begin{cases} z & \text{, if $z > 0$} \\ 0 & \text{otherwise} \end{cases}$

def relu(Z):

The last layer will use the $sigmoid$ activation function. This way, the final output will be squished to values between $0$ and $1$ which then can be used for predicting the classes.

$sigmoid(z) = \frac{1}{1 + e^{-z}}$

def sigmoid(Z):

Each layer’s activation $A$ is calculated as follows:

$A = \s( WA_{prev} + b )$

The forward pass begins with the input $X$ (a single sample or a batch) and iterates over the layers.

def nn_forward(X, param):

All activations $A^{[l]}$ ($l \in 1, .. L$) are stored in the dictionary named cache, and the input $X$ in cache['A0'].

loss function

For our network to be able to learn anything, we need to formulate a loss function. For classification tasks, the cross-entropy loss function is a common choice. Its formulation for one sample is:

$\mathcal{H}(p,q) = -\sum_{i \in C} p_i \log(q_i)$

, where the sum is over all possible classes $C$.

$p_i = \P(x \text{ of class } i)$

is the true probability of the sample belonging to class $i$, and

$q_i = \P(x \text{ of class } i \enspace|\; \phi)$

is the predicted probability given the neural network model $\phi$.
For the binary case, where the class labels $y$ are either $0$ or $1$, the sum expands to

$\mathcal{H}(p,q) = -p_1 \log(q_1) - p_0 \log(q_0)$

The true probability of a sample belonging to ‘class 1’ is reflected by its label.

$p_1 = y$

And the true probability of not belonging to ‘class 1’, but ‘class 0’ is

$p_0 = 1 - y$

Similarly, $q_i$ are the network’s predictions, the last layer’s output $a$.

$\begin{aligned} q_1 &= a \\ q_0 &= 1 - a\end{aligned}$

The binary loss function then is simplified to:

$\mathcal{L}(a,y) = -y\log a - (1-y)\log (1-a)$

For a batch of samples $X = [x^{(1)}, .., x^{(N)}]$, the loss can be averaged.

$\mathcal{L}(A,Y) = - \frac{1}{N} \sum_{n=0}^N \left(y^{(n)}\log (a^{(n)}) + (1-y^{(n)})\log(1-a^{(n)})\right)$

As a reminder:
For now, $A = (a^{(1)}, .. a^{(N)})$ and $Y = (y^{(1)}, .. y^{(N)})$ are row-vectors, and $a^{(n)}$ and $y^{(n)}$ scalars, since the output has only one dimension. Later, in the multivariate case, these will be matrices and vectors respectively.

In python:

def loss(A, Y):

backward pass

Gradient descent is like a blind man trying to reach the bottom of a hill. He is only able to feel the hill’s slope with his feet and take small steps in the hill’s downward direction.

A derivative of a function similarly reflects the rate of change (or slope when graphed) of a parameter. It essentially says how much the output will change by a small tweak in input. In the multivariable case, the gradient points towards the direction of biggest (positive) change in output. In order to drive the network to a lower loss, all parameters are moved in the opposite direction of their gradients, in the ‘direction of negative slope’, by a small amount. Although, in general, this does not guarantee a reduce in cost.

The network’s parameters are its weights and biases. Thus, the gradients $dW := \pd{\L}{W}^\top$ and $db := \pd{\L}{b}^\top$ need to be evaluated for each layer. For this, the derivatives of the activation functions, the loss and the affine linear functions need to be calculated.

In the case of $\mathcal L(A(Z))$ (the loss function), where $\L$ (a scalar) is dependent on $A$ (a matrix), which, in turn is dependent on $Z$ (a matrix), evaluating the derivative w.r.t. $Z$ can be done for all indices using the following chain rule:

$\pd{\L}{z_{ij}} = \sum _{kl} \pd{\L}{a_{kl}}\pd{a_{kl}}{z_{ij}}$

In matrix form this is often written as

$\pd{\L}{Z} = \pd{\L}{A}\pd{A}{Z}$

A word of warning. It is easy to make mistakes using the matrix version of the chain rule without knowing what is going on under the hood. Note that for example, $\pd A Z$ is a 4-tensor and the derivative $\pd \L A$ is a linear transformation from $\R^{d\x N}$ to $\R$, which cannot even be represented as a matrix using conventional notation. Even if we, nonetheless, decide to write the derivative as a matrix, keeping consistency with the numerator layout would make it a matrix of size $N\x d$. That is why I choose to write $\pdt \L A$ when referring to the gradient which has the same dimensions as $A$ $\in \R^{d\x N}$. I further clarify this .

Rewriting the above equation yields

$\begin{aligned} \pd{\L}{Z}^\top &= \pd{A}{Z}^\top\pd{\L}{A}^\top \\ \end{aligned}$

or

$\begin{aligned} dZ &= \pd{A}{Z}^\top dA \end{aligned}$

$\begin{aligned} &\R^{d\x N}& &\to& &\R^{d\x N}& &\to& \R \\ &Z& &\mapsto& &A& &\mapsto& \L \\ &dZ& &\underset{\pd{A}{Z}^\top}{\leftmapsto}& &dA& && \end{aligned}$

And it can be seen how, in order to calculate the desired derivative $dZ$, the gradient $dA$ is ‘backpropagated’.

cross-entropy loss

The cross-entropy loss function

$\mathcal{L}(A,Y) = - \frac{1}{N} \sum_{n=0}^N \left ( y^{(n)}\log (a^{(n)}) + (1-y^{(n)})\log(1-a^{(n)}) \right )$

has the following gradient:

$dA := \pdt \L A = \frac{1}{N}(-Y/A + (1-Y)/(1-A))$

Where the expression $Y/A$ above is meant to be performed element-wise.
This can be seen by calculating the derivative w.r.t. one sample. Then $a^{(n)}$ and $y^{(n)}$ become scalars:

$\frac{\partial \mathcal L}{\partial a^{(n)}} = \frac{1}{N}(-\frac{y^{(n)}}{a^{(n)}} + \frac{1-y^{(n)}}{1-a^{(n)}})$

As the summands are independent of each other, the gradient can be written in the above, vectorized form.

activation functions

relu:

$relu'(z) = \mathbb{1}_{\{z\ge 0\}} = \begin{cases} 1 & \text{, if $z \ge 0$} \\ 0 & \text{otherwise} \end{cases}$

sigmoid:

$\sigma '(z) = \sigma(z)(1 - \sigma(z)) = a(1 - a)$

$\begin{aligned} \sigma '(z) &= \left (\frac{1}{1+e^{-z}}\right )' = -1 (1+e^{-z})^{-2}e^{-z}(-1) \\ &= \frac{1}{1+e^{-z}}\frac{e^{-z}}{1+e^{-z}} = \frac{1}{1+e^{-z}}\frac{1 + e^{-z} - 1}{1+e^{-z}} \\ &= \frac{1}{1+e^{-z}}(1 - \frac{1}{1+e^{-z}}) = \sigma(z)(1 - \sigma(z)) \end{aligned}$

Where the chain-rule was applied in the first equation.

affine linear functions

Given an affine linear function $Z = WX + b$, the derivatives are $\pdt Z W = \cdot X^\top$ and $\pdt Z b = \cdot \vec\1$.
A derivation can be found . Many matrix derivatives can also be looked up in the matrix cookbook.

piecing together

Before implementing the code in python, I believe it is helpful to run through a small 2-layer example, from which the general algorithm can be derived. The aim is to evaluate $dW$ and $db$ for each layer.

$forward$		$backward$
$Z^{[1]} = W^{[1]} X + b^{[1]} \\ \downarrow$		$dW^{[1]} = dZ^{[1]} X^T \\ db^{[1]} = dZ^{[1]}\vec\1 \\ \uparrow$
$A^{[1]} = relu(Z^{[1]}) \\ \downarrow$		$dZ^{[1]} = dA^{[1]} * (Z^{[1]} > 0) \\ \uparrow$
$Z^{[2]} = W^{[2]} A^{[1]} + b^{[2]} \\ \downarrow$		$dW^{[2]} = dZ^{[2]} {A^{[1]}}^T \\ db^{[2]} = dZ^{[2]}\vec\1 \\ dA^{[1]} = {W^{[2]}}^T dZ^{[2]} \\ \uparrow$
$A^{[2]} = sigmoid(Z^{[2]}) \\ \downarrow$		$dZ^{[2]} = dA * A^{[2]} *(1 - A^{[2]}) \\ \uparrow$
$\L = - \frac{1}{N} \sum [Y\log (A)$ $+ (1-Y)\log(1-A)]$	$\longrightarrow$	$dA^{[2]} = \frac{1}{N}(-Y/A^{[2]}$ $+ (1-Y)/(1-A^{[2]}) )$

‘*’ and ‘/’ on matrices are used to indicate element-wise actions, keeping close to python’s syntax.

Making sure that the gradients match the dimensions of the variables also can prove helpful. For example one should be able to verify that A2.shape == dA2.shape.

The backpropagation through the last layer can be further simplified, as it is only necessary to calculate $dZ^{[2]}$.

$\begin{aligned} dZ^{[2]} &= dA * A^{[2]} *(1 - A^{[2]}) \\ &= -\frac{1}{N}(Y*(1-A^{[2]}) - (1-Y) * A^{[2]}) \\ &= \frac{1}{N}(A^{[2]} * (1-Y) - Y * (1-A^{[2]})) \end{aligned}$

This step also prevents numerical errors that could arise when $A^{[2]}$ is close to $0$ or $1$.

Implementation:

def relu_backward(Z):

parameter update

The parameters are updated using the following rule:

$W = W - \alpha \pdt \L W = W - \alpha dW$

, given the learning rate $\alpha$.

def update_param(param, grad, learning_rate=0.01):

training

Training the neural network is achieved by iterating over the following steps:

1. forward pass (evaluating $A^{[l]}$ for each layer)
2. backward pass (computing the gradients $dW^{[l]}$ and $db^{[l]}$)
3. parameter update (tuning the weights and biases)

A list of the computed costs can be stored for logging/graphing the network’s progress.

def train(X, Y, param, epochs=50000, learning_rate=0.001,

predicting and evaluation

When making a class prediction, the class with the higher probability is chosen. Thus, if the estimated probability for ‘class 1’ is greater than 0.5, the prediction is set to ‘class 1’, otherwise it will be ‘class 0’.

The networks empirical correct classification rate is

$\frac{\text{\# correct classifications}}{\text{\# samples}}$

def predict(X, param):

toy data

Sci-kit sklearn has built in functions to easily create datasets.

def load_dataset():

The dataset is shuffled and then split into two parts, a training and a testing dataset. The test set is used for evaluating the network’s ability to extrapolate on data it has never seen before.

def split_dataset(X, Y, ratio=0.8):

X_train shape: (2, 240)
Y_train shape: (240,)
X_test shape: (2, 60)
Y_test shape: (60,)

binary model

Now that all necessary building components are finished, it is time to assemble the network. The number of input dimensions can be read directly from the dataset (in this case the data is 2-dimensional) and the output uses one neuron in the binary case. The dimension of the hidden layer is somewhat arbitrarily set to 8. Passing the specified layer dimensions to the nn_init() function, the parameters of a randomly initialized neural network are returned. The network is then trained on the training dataset.

d_in = X_train.shape[0]

epoch 0, loss: 0.8943569024225951
epoch 100000, loss: 0.47308653964775604
epoch 200000, loss: 0.46855221922894785
epoch 300000, loss: 0.40372674060427743
epoch 400000, loss: 0.36625403094480086
epoch 500000, loss: 0.35062643425360634

The next part is about visualizing the performance of our network. It is not essential for understanding neural networks. The code of plot_decision_boundary is mostly adapted from scikit.

code for plot_costs, plot_decision_boundary

def plot_costs(costs, plot_every=1000, learning_rate=.001):

plot_costs(costs)

accuracy on train set: 0.8666666666666667
accuracy on test set: 0.8

multiclass neural network

The current network architecture only has one output neuron. Given, say, 10 classes, the network’s i-th output neuron then is to reflect the probability of the input belonging to ‘class i’.

$a = \vecfill{a}{c} = \veceval{\P(\text{x of class }\i \;|\; \phi)}{c}$

As an example

$a = \begin{pmatrix} 0.2 \\ \vdots \\ 0.6 \\ \vdots \\ 0.1\end{pmatrix} \text{ \} class i probability}$

Changing the current code-framework to be able to handle multiple classes involves the following:

• implementing the softmax layer, to be used as the last layer’s activation function
• revisiting the original cross-entropy loss formulation
• calculating the backward pass for the new last layer
• modifying the prediction function

softmax layer

After propagating the the input through an affine linear layer $Z = WX + b$, $Z$ can take on all sorts of values, positive or negative. If the 10 outputs are to reflect a probability distribution (i.e. have values between 0 and 1 that all sum to 1), the following can be done.

Each value is mapped to the exponential function, to ensure that it will be positive, and then divided by the sum of all outputs. The result can then be seen as a probability distribution.
This is exactly what softmax does:

$softmax(z) = \frac{e^z}{\sum_i e^{z_i}} = a$

$\R^d \to \R^d$

As a side note: in the binary case it was redundant to work with two outputs, since the probability of not belonging to ‘class 1’ is simply

$a_1 = 1 - a_1$

In that case

$\begin{aligned} a_1 &= \frac{e^{z_1}}{e^{z_1} + e^{z_2}} = \frac{1}{1 + e^{-z_1}e^{z_2}} \\ &= \frac{1}{1+e^{-(z_1-z_2)}} \\ &= sigmoid(z_1-z_2) \end{aligned}$

which can be substituted as

$sigmoid(z)$

Again, we see the binary case emerge as a special case of the general case.

A numerical trick that is often employed in order to prevent the exponential term from blowing up (which may result in unstable computations) is the following:

$softmax(z) = \frac{e^{-c}}{e^{-c}}\frac{e^z}{\sum_i e^{z_i-c}} = \frac{e^{z-c}}{\sum_i e^{z_i-c}}$

$c$ is set to $\underset{i}{max}\{ z_i\}$, the biggest value, to ensure numerical stability.

def softmax(Z):

cross-entropy loss

As a reminder, the loss for one sample is

$\L(a, \hat y) = -\sum_{i \in C} \hat y_i \log(a_i)$

, where $\hat y$ is the one hot encoding of y.

The ‘one hat’ above the $y$ is supposed to remind of the one hot encoding.

Before, a scalar $y$ was used as a label for sample $x$. $y$ = 3, as an example.
Now, with one hot encoding, it is assumed that $\hat y$ is a vector of, say, length 10, with a ‘1’ at the index corresponding to the sample’s class, and ‘0’ everywhere else.

$x$ of ‘class 3’:

$y = 3 \quad \Longrightarrow \quad \hat y = \begin{pmatrix} 0 \\ \vdots \\ 1 \\ \vdots \\ 0 \end{pmatrix}\leftarrow \text{third index}$

The class label data shape is the same as the neural network’s output.

To reduce overhead we can avoid explicitly converting our labels to one hot encodings by using numpy’s indexing.

$\mathcal L = -\sum_{i \in C} \hat y_i \log(a_i) = -\log(a_y)$

Using numpy’s notation to access elements in arrays this would look like -log(a[y]).

For a batch of samples:

$\L = -\frac 1 N\sum_{n = 1}^N \log(a_{y^{(n)}}^{(n)})$

Using numpy’s advanced indexing, accessing these elements is done via A[Y, range(N)]. This expands to A[[Y(1), Y(2), ..], [1, 2, ..]], which then can be summed with np.sum.

def loss(A, Y):

In deriving the theory, it is useful to stick to the one-hot formulation:

$\begin{aligned} \L(A, \hat Y) &= -\frac 1 N\sum_{n = 1}^N \sum_{i \in C} \hat y_i^{(n)} \log(a_i^{(n)}) \\&= -\frac 1 N\sum \hat Y \odot \log(A) \end{aligned}$

where the sum $\sum$ is over all free indices and ‘$\odot$’ is the element-wise multiplication

From now on, to make it clearer, ‘$\odot$’ and ‘$\oslash$’ will be used for element-wise multiplication and division.

backward pass

The goal is to calculate $dZ = \pdt \L Z = \pdt A Z \pdt \L A = \pdt A Z dA$, the derivative of the last layer.

cross-entropy loss

The derivative of the cross-entropy loss

$\L(A, \hat Y) = - \frac 1 N \sum \hat Y \odot log(A)$

remains as before

$dA = \pdt \L A = -\frac{1}{N} \hat Y \oslash A$

, since all operations are still performed element-wise.

softmax layer

The derivative of the softmax layer $\pdt A Z$ is a bit more involved. The derivative of a matrix-function w.r.t. another matrix is a rank-4 tensor. I cover its full derivation via the chain rule .

Yet, we really only are interested in calculating $dZ = \pdt \L Z$ (a matrix). And since the sum in the loss function involves independent terms and the softmax function acts on each sample individually, it is possible to calculate the derivative $dz = \pdt \L z = \pdt a z \pdt \L a$ w.r.t. one input and then piece together the full derivative $dZ$. The full derivative then becomes $\pdt \L Z = \bmatrix{\pdt \L {z^{(1)}} & \cdots & \pdt \L {z^{(N)}}}$. And now, the term $\pdt a z$ is a (more manageable) matrix.

With respect to one sample, the derivative of the loss function is

$da = \pdt \L a = -\frac{1}{N} \hat y \oslash a = -\frac{1}{N} \pmatrix{ 0 \\ \vdots \\\frac{1}{a_y} \\ \vdots \\ 0 }$

And the derivative $\pdt a z$ of the softmax layer

$A = softmax(Z) = \frac{e^Z}{\sum_i e^{Z_i}}$

is calculated as follows:

$\begin{aligned} i = j: \quad \pd {a_i} {z_i} &= \frac{\partial}{\partial z_i}\left (\frac{e^{z_i}}{\sum_k e^{z_k}}\right ) \\ &= \frac{1}{\sum_i{e^{z_i}}} e^{z_i} - \frac{e^{z_i}}{(\sum_k e^{z_k})^2}e^{z_i} \\ &= a_i(1-a_i) \\ i \neq j: \quad \pd {a_i} {z_j} &= \frac{\partial}{\partial z_j}\left (\frac{e^{z_i}}{\sum_k e^{z_k}}\right ) \\ &= - \frac{e^{z_i}}{(\sum_k e^{z_k})^2}e^{z_j} \\ &= -a_ia_j \end{aligned}$

The product rule is used in after the first equation.

In full:

$\begin{aligned} \pdt a z &= - \begin{bmatrix} a_1(a_1-1) & a_1a_2 & \cdots \\ a_2a_1 & a_2(a_2-1) & & \\ \vdots & & \ddots \end{bmatrix} \\ & = - (aa^\top - diag(a)) \end{aligned}$

Combining both derivatives to obtain $dz$ gives:

$\begin{aligned} dz &= \pdt \L z = \pdt a z \pdt \L a \\ &= \frac{1}{N} (aa^\top - diag(a))(\hat y \oslash a) \\ &= \frac 1 N (a\mathbb (a \oslash a) ^\top \hat y - diag(a\oslash a)\hat y) \\ &= \frac 1 N (a \mathbb 1 ^\top \hat y - I\hat y) = \frac 1 N (a - \hat y) \end{aligned}$

This is easier seen when expanding the matrices.

$\begin{aligned} dz &= \pdt \L z = \pdt a z \pdt \L a \\ &= \begin{bmatrix} a_1(a_1-1) & a_1a_2 & \cdots \\ a_2a_1 & a_2(a_2-1) & & \\ \vdots & & \ddots \end{bmatrix} \pmatrix{ 0 \\ \vdots \\\frac{1}{a_y} \\ \vdots \\ 0 }\\ &= \frac{1}{N} \begin{pmatrix} a_1 \\ a_2 \\ \vdots \\ (a_y-1) \\ \vdots \\ a_d \end{pmatrix} = \frac 1 N (a - \hat y) \end{aligned}$

$a - \hat y$ results in subtracting $1$ from the element of $a$ with the index corresponding to its true class.

The full derivative now is

$dZ = \frac 1 N ( A - \hat Y)$

In python, again, using advanced indexing:

def last_layer_backward(A, Y):

predicting

The output layer contains the predicted class probabilities.

$a = \vecfill{a}{n} = \veceval{\P(\text{x of class }\i \;|\; \phi)}{n}$

The final prediction is done by selecting the class with the highest estimated probability. This corresponds to the index of the highest value.
For a batch of samples this is done for each sample, i.e. each column of the matrix (axis 0 in numpy).

def predict(X, param):

toy data

A toy dataset can be created by reusing the toy dataset for the binary case and adding multivariate normal distributed samples as a third class.

def load_dataset_multi():

multiclass model

Finally, the multiclass neural network is initialized by specifying the layer dimensions and then is trained on the dataset, just as in the binary case.

d_in = X_train.shape[0]

epoch 0, loss: 0.8492644177201928
epoch 10000, loss: 0.584044408860389
epoch 20000, loss: 0.5478712675155261
epoch 30000, loss: 0.5135834238863085
epoch 40000, loss: 0.500427787036695
epoch 50000, loss: 0.4898009397411907

accuracy on train set: 0.8125
accuracy on test set: 0.75

The full python code can be found on my Github. Alternatively, the code can be run directly in a Google Colab notebook here.