# Artificial Neural Networks (ANNs) and Backpropagation

Classify images into ten digits by using an artificial neural network. Let $x \in \Re^{784\times N}$ be a collection of N input images (stacked in columns) and $y \in [0, 1]^{10\times N}$ be the collection of our class predictions made for all images and encoded by so-called one-hot codes meant to approximate class probabilities.  

A neural network makes the prediction y from the collection x of images by forward propagation which is defined recursively as follows
$$
\boldsymbol{y}=\boldsymbol{h}_L \quad with \quad \boldsymbol{h}_k=g_k(\boldsymbol{a}_k), \quad \boldsymbol{a}_k=\boldsymbol{W}_k \boldsymbol{h}_{k-1} + \boldsymbol{b}_k \quad and \quad \boldsymbol{h}_0=\boldsymbol{x}
$$
with $L$ the number of layers, $k$ the layer index, $g_k$ a so-called activation function, $\boldsymbol{h}_k$ an array of $N$ hidden feature vectors (also stacked in columns), $\boldsymbol{a}_k$ an array of $N$ vectors of activations (also called potentials), $\boldsymbol{W}_k$ a matrix of synaptic weights and $\boldsymbol{b}_k$ a vector of biases.  

We will use a shallow neural network which means $L = 2$. The hidden layer will use ReLU activation function and the output layer will use Softmax.  

The matrices $\boldsymbol{W}_k$ and vectors $\boldsymbol{b}_k$ will be learned/estimated from the MNIST training set by using logistic regression. This consists of minimizing the cross-entropy loss, between the collection d of desired one-hot codes ($d_{ij} = 1$ if the j-th image represents the digit i, 0 otherwise) and the prediction $y$, which is defined as follows
$$
E=-\Sigma_{j=1}^{N} \Sigma_{i=1}^{10} d_{ij} \log{y_{ij}}
$$  

The optimization will be performed by gradient descent with backpropagation that can be implemented
iteratively, for $t$ = 0, 1, ..., and some initializations $\boldsymbol{W}_k^0$ and $\boldsymbol{b}_k^0$ , as:


```python
import numpy as np
from MNISTtools import load, show
from matplotlib import pyplot as plt
```

## 1. Load dataset


```python
xtrain, ltrain = load(dataset='training', path='dataset/')
```


```python
print(f'xtrain.shape = {xtrain.shape}')
print(f'ltrain.shape = {ltrain.shape}')
print(f'size of training set = {xtrain.shape[1]}')
print(f'feature dimension = {xtrain.shape[0]}')
```

    xtrain.shape = (784, 60000)
    ltrain.shape = (60000,)
    size of training set = 60000
    feature dimension = 784


## 2. Display the image of index 42


```python
print(f'label of index 42 is {ltrain[42]}')
show(xtrain[:,42].reshape((28,28)))
```

## 3. The range and type of `xtrain`


```python
print(f'range of xtrain = from {xtrain.min()} to {xtrain.max()}')
print(f'type of xtrain = {type(xtrain)}')
print(f'dtype of xtrain = {xtrain.dtype}')
```

    range of xtrain = from 0 to 255
    type of xtrain = <class 'numpy.ndarray'>
    dtype of xtrain = uint8


## 4. Normalize `xtrain`
Normalize `xtrain` s.t. it is in the range `[-1,1]` of type `float32`.


```python
def normalize_MNIST_images(x):
    '''
    Inputs:
        x: data
    '''
    x_norm = x.astype(np.float32)
    return x_norm*2/255-1
```


```python
xtrain = normalize_MNIST_images(xtrain)
print(f'min normalized xtrain = {np.min(xtrain)}')
print(f'max normalized xtrain = {np.max(xtrain)}')
```

    min normalized xtrain = -1.0
    max normalized xtrain = 1.0


## 5. One-hot encoding


```python
def label2onehot(lbl):
    '''
    Convert label (n,) to one-hot form (d,n).
    '''
    d = np.zeros((lbl.max() + 1, lbl.size))
    d[lbl, np.arange(lbl.size)] = 1
    return d
```


```python
dtrain = label2onehot(ltrain)
print(f'dtrain shape = {dtrain.shape}')
print(f'{np.where(dtrain[:,42] == 1)}')
print(f'{np.where(dtrain[:,42] == 1) == ltrain[42]}')
```

    dtrain shape = (10, 60000)
    (array([7]),)
    [[ True]]


## 6. One-hot decoding


```python
def onehot2label(d):
    '''
    Inputs:
        d: one-hot encoding labels
    '''
    lbl = d.argmax(axis=0)
    return lbl
```


```python
all(ltrain == onehot2label(dtrain))
```




    True



## 7. The softmax function
$$
y_i=g(\boldsymbol{a})_i=\frac{e^{a_i}}{\Sigma_{j=1}^{10} e^{a_j}} \quad where \quad i\in[1,10]
$$


```python
def softmax(a):
    '''
    Inputs: 
        a: data
    '''
    exp_a = np.exp(a - np.max(a, axis=0))
    return exp_a / np.sum(exp_a, axis=0)
```

## 8. $$ 
\frac{\partial g(\boldsymbol{a})_i}{\partial a_i} = {g(\boldsymbol{a})}_i (1-{g(\boldsymbol{a})_i})
$$

Proof:

$$
\begin{align}
\frac{\partial g(\boldsymbol{a})_i}{\partial a_i} & = 
  \frac{\partial}{\partial a_i}\frac{e^{a_i}}{\Sigma_{k=1}^{N} e^{a_k}} \\
& = \frac{e^{a_i}\Sigma_{k=1}^{N}e^{a_k}-e^{a_i}e^{a_i}}{(\Sigma_{k=1}^{N}e^{a_k})^2} \\
& = \frac{e^{a_i}}{\Sigma_{k=1}^{N}e^{a_k}} \frac{\Sigma_{k=1}^{N}e^{a_k}-e^{a_i}}{\Sigma_{k=1}^{N}e^{a_k}} \\
& = {g(\boldsymbol{a})}_i (1-{g(\boldsymbol{a})_i})
\end{align}
$$

## Q9$$
\frac{\partial{g(\boldsymbol{a})_i}}{\partial a_j}=-g(\boldsymbol{a})_i g(\boldsymbol{a})_j \quad for \, j \neq i
$$

Proof:

$$
\begin{align}
\frac{\partial g(\boldsymbol{a})_i}{\partial a_j} & = 
  \frac{\partial}{\partial a_j} \frac{e^{a_i}}{\Sigma_{k=1}^{N} e^{a_k}} \\
& = \frac{0-e^{a_i}e^{a_j}}{(\Sigma_{k=1}^{N}e^{a_k})^2} \\
& = -\frac{e^{a_i}}{\Sigma_{k=1}^{N}e^{a_k}} \frac{e^{a_j}}{\Sigma_{k=1}^{N}e^{a_k}} \\
& = -{g(\boldsymbol{a})}_i g(\boldsymbol{a})_j
\end{align}
$$

## 10. Compute the $\delta$
$$
\delta = g(a)\otimes e-<g(a),e>g(a)
$$


```python
def softmaxp(a, e):
    '''
    Compute delta for the backpropagation.
    
    Inputs:
        a: predictions before softmax
        e: onehot labels
        
    Returns:
        delta
    '''
    softmax_a = softmax(a)
    
    # element-wise product
    e_prod = softmax_a*e
    
    # dot product
    dot_prod = np.sum(e_prod, axis=0)
    
    return e_prod - dot_prod*softmax_a
```

# 11. Numerical approximations for $\delta\$
$$
\delta=\frac{\partial g(a)}{\partial a}\times e = \lim_{\epsilon \to 0}\frac{g(a+\epsilon e)-g(a)}{\epsilon}
$$


```python
eps = 1e-6 # finite difference step
a = np.random.randn(10, 200) # random inputs
e = np.random.randn(10, 200) # random directions

diff = softmaxp(a, e)
diff_approx = (softmax(a + eps*e) - softmax(a)) / eps
rel_error = np.abs(diff-diff_approx).mean() / np.abs(diff_approx).mean()
print(rel_error, 'should be smaller than 1e-6')
```

    4.888515969765918e-07 should be smaller than 1e-6


## 12. ReLU and its directional derivative


```python
def relu(a):
    return np.maximum(a, 0)

def relup(a, e):
    '''
    The directional derivative of the ReLU function.
    '''
    eps = 1e-6
    return (relu(a + eps * e) - relu(a)) / eps
```


```python
eps = 1e-6 # finite difference step
a = np.random.randn(10, 200) # random inputs
e = np.random.randn(10, 200) # random directions

diff = relup(a, e)
diff_approx = (relu(a + eps*e) - relu(a)) / eps
rel_error = np.abs(diff-diff_approx).mean() / np.abs(diff_approx).mean()
print(rel_error, 'should be smaller than 1e-6')
```

    0.0 should be smaller than 1e-6


## 13. Initialization of the net
Utilize `He and Xavier initializations`.


```python
def init_shallow(Ni, Nh, No):
    '''
    Inputs:
        Ni: dimension of input layer
        Nh: dimension of output layer
        No: number of unit of output layer
    
    Returns:
        parameters of the net
    '''
    b1 = np.random.randn(Nh, 1) / np.sqrt((Ni+1.)/2.)
    W1 = np.random.randn(Nh, Ni) / np.sqrt((Ni+1.)/2.)
    b2 = np.random.randn(No, 1) / np.sqrt((Nh+1.))
    W2 = np.random.randn(No, Nh) / np.sqrt((Nh+1.))
    return W1, b1, W2, b2
```


```python
Ni = xtrain.shape[0] # 784
Nh = 64
No = dtrain.shape[0] # 10
netinit = init_shallow(Ni, Nh, No)
```

## 14. Foward propagation


```python
def forwardprop_shallow(x, net):
    '''
    Inputs:
        x: data
        net: parameters of the net
    
    Returns:
        prediction    
    '''
    W1 = net[0] # (64, 784)
    b1 = net[1] # (64, 1)
    W2 = net[2] # (10, 64)
    b2 = net[3] # (10, 1)
    
    a1 = W1.dot(x) + b1 # (64, 60000)
    h1 = relu(a1) # (64, 60000)
    
    a2 = W2.dot(h1) + b2 # (10, 60000)
    y = softmax(a2) # (10, 60000)
    return y
```


```python
yinit = forwardprop_shallow(xtrain, netinit)
```

## 15. Loss


```python
def eval_loss(y, d):
    '''
    Inputs:
        y (10, 60000): prediction
        d (10, 60000): ground truth
    
    Returns:
        Average cross-entropy loss (1, 1)
    '''
    return -np.sum(d*np.log(y))/y.shape[0]/y.shape[1]
```


```python
print(eval_loss(yinit, dtrain), 'should be around .26')
```

    0.26297468338637725 should be around .26


## 16. Error rate


```python
def eval_perfs(y, lbl):
    '''
    Compute the percentage of misclassified samples.
    
    Inputs:
        y (10, 60000): prediction
        lbl (10, 60000): ground truth
    
    Returns:
        Error rate
    '''
    return np.sum(onehot2label(y)!=lbl)/lbl.size
```


```python
print(eval_perfs(yinit, ltrain))
```

    0.8715


## 17. Update the parameters 
Perform one backpropagation update for the network.

$$
E=-\Sigma_{j=1}^{N} \Sigma_{i=1}^{10} d_{ij} \log{y_{ij}} \\
(\nabla_y E)_i=-\frac{d_i}{y_i}
$$


```python
def update_shallow(x, d, net, gamma=.05):
    '''
    Inputs:
        x: data
        d: ground truth
        net: parameters of the net
        gamma: learning rate
    
    Returns:
        parameters of the net
    '''
    W1 = net[0]
    b1 = net[1]
    W2 = net[2]
    b2 = net[3]
    Ni = W1.shape[1]
    Nh = W1.shape[0]
    No = W2.shape[0]
    gamma = gamma / x.shape[1] # normalized by the training dataset size
    
    a1 = W1.dot(x) + b1
    h1 = relu(a1)
    a2 = W2.dot(h1) + b2
    y = softmax(a2)
    
    # derivative of cross-entropy loss
    e = -d/y + (1-d)/(1-y)
    
    # backprop
    delta2 = softmaxp(a2, e)
    delta1 = relup(a1, W2.T.dot(delta2))
    
    W2 = W2 - gamma * delta2.dot(h1.T)
    W1 = W1 - gamma * delta1.dot(x.T)
    b2 = b2 - gamma * delta2.sum(axis=1, keepdims=True)
    b1 = b1 - gamma * delta1.sum(axis=1, keepdims=True)
    
    return W1, b1, W2, b2
```

## 18. Backpropagation


```python
def backprop_shallow(x, d, net, T, gamma=.05):
    '''
    Inputs:
        x: data
        d: ground truth
        net: parameters of the net
        T: number of updates
        gamma: learning rate
    
    Returns:
        parameters of the net
    '''
    lbl = onehot2label(d)
    for t in range(T):
        # UPDATE NET
        y = forwardprop_shallow(x, net)
        # DISPLAY LOSS AND PERFS
        net = update_shallow(x, d, net, gamma)
        print(f'\niter = {t+1}')
        print(f'loss = {eval_loss(y, d):.4f}')
        print(f'error rate = {eval_perfs(y, lbl):.4f}')
    return net
```


```python
print('\niter_max = 100\n')
nettrain = backprop_shallow(xtrain, dtrain, netinit, 100)
```

    
    iter_max = 100
    
    
    iter = 1
    loss = 0.2630
    error rate = 0.8715
    
    iter = 2
    loss = 0.2224
    error rate = 0.8073
    
    iter = 3
    loss = 0.2094
    error rate = 0.7125
    
    iter = 4
    loss = 0.2003
    error rate = 0.6544
    
    iter = 5
    loss = 0.1919
    error rate = 0.5986
    
    iter = 6
    loss = 0.1840
    error rate = 0.5502
    
    iter = 7
    loss = 0.1765
    error rate = 0.5087
    
    iter = 8
    loss = 0.1695
    error rate = 0.4803
    
    iter = 9
    loss = 0.1634
    error rate = 0.4573
    
    iter = 10
    loss = 0.1589
    error rate = 0.4576
    
    iter = 11
    loss = 0.1629
    error rate = 0.4871
    
    iter = 12
    loss = 0.1730
    error rate = 0.5871
    
    iter = 13
    loss = 0.1893
    error rate = 0.5526
    
    iter = 14
    loss = 0.1635
    error rate = 0.5360
    
    iter = 15
    loss = 0.1357
    error rate = 0.3711
    
    iter = 16
    loss = 0.1290
    error rate = 0.3540
    
    iter = 17
    loss = 0.1291
    error rate = 0.3766
    
    iter = 18
    loss = 0.1290
    error rate = 0.3786
    
    iter = 19
    loss = 0.1385
    error rate = 0.4335
    
    iter = 20
    loss = 0.1280
    error rate = 0.3814
    
    iter = 21
    loss = 0.1240
    error rate = 0.3857
    
    iter = 22
    loss = 0.1174
    error rate = 0.3435
    
    iter = 23
    loss = 0.1181
    error rate = 0.3681
    
    iter = 24
    loss = 0.1080
    error rate = 0.3079
    
    iter = 25
    loss = 0.1064
    error rate = 0.3271
    
    iter = 26
    loss = 0.1007
    error rate = 0.2932
    
    iter = 27
    loss = 0.1031
    error rate = 0.3241
    
    iter = 28
    loss = 0.0951
    error rate = 0.2821
    
    iter = 29
    loss = 0.0975
    error rate = 0.3120
    
    iter = 30
    loss = 0.0922
    error rate = 0.2815
    
    iter = 31
    loss = 0.0948
    error rate = 0.3071
    
    iter = 32
    loss = 0.0892
    error rate = 0.2742
    
    iter = 33
    loss = 0.0898
    error rate = 0.2942
    
    iter = 34
    loss = 0.0861
    error rate = 0.2661
    
    iter = 35
    loss = 0.0865
    error rate = 0.2833
    
    iter = 36
    loss = 0.0826
    error rate = 0.2547
    
    iter = 37
    loss = 0.0824
    error rate = 0.2675
    
    iter = 38
    loss = 0.0795
    error rate = 0.2461
    
    iter = 39
    loss = 0.0792
    error rate = 0.2559
    
    iter = 40
    loss = 0.0768
    error rate = 0.2387
    
    iter = 41
    loss = 0.0761
    error rate = 0.2450
    
    iter = 42
    loss = 0.0747
    error rate = 0.2339
    
    iter = 43
    loss = 0.0735
    error rate = 0.2349
    
    iter = 44
    loss = 0.0730
    error rate = 0.2298
    
    iter = 45
    loss = 0.0710
    error rate = 0.2251
    
    iter = 46
    loss = 0.0714
    error rate = 0.2264
    
    iter = 47
    loss = 0.0688
    error rate = 0.2159
    
    iter = 48
    loss = 0.0700
    error rate = 0.2228
    
    iter = 49
    loss = 0.0670
    error rate = 0.2084
    
    iter = 50
    loss = 0.0687
    error rate = 0.2192
    
    iter = 51
    loss = 0.0654
    error rate = 0.2009
    
    iter = 52
    loss = 0.0673
    error rate = 0.2145
    
    iter = 53
    loss = 0.0643
    error rate = 0.1947
    
    iter = 54
    loss = 0.0660
    error rate = 0.2091
    
    iter = 55
    loss = 0.0633
    error rate = 0.1913
    
    iter = 56
    loss = 0.0646
    error rate = 0.2040
    
    iter = 57
    loss = 0.0624
    error rate = 0.1885
    
    iter = 58
    loss = 0.0632
    error rate = 0.1981
    
    iter = 59
    loss = 0.0614
    error rate = 0.1844
    
    iter = 60
    loss = 0.0617
    error rate = 0.1912
    
    iter = 61
    loss = 0.0602
    error rate = 0.1805
    
    iter = 62
    loss = 0.0601
    error rate = 0.1847
    
    iter = 63
    loss = 0.0587
    error rate = 0.1757
    
    iter = 64
    loss = 0.0585
    error rate = 0.1787
    
    iter = 65
    loss = 0.0572
    error rate = 0.1693
    
    iter = 66
    loss = 0.0570
    error rate = 0.1725
    
    iter = 67
    loss = 0.0557
    error rate = 0.1630
    
    iter = 68
    loss = 0.0555
    error rate = 0.1662
    
    iter = 69
    loss = 0.0542
    error rate = 0.1572
    
    iter = 70
    loss = 0.0541
    error rate = 0.1605
    
    iter = 71
    loss = 0.0529
    error rate = 0.1530
    
    iter = 72
    loss = 0.0531
    error rate = 0.1573
    
    iter = 73
    loss = 0.0519
    error rate = 0.1515
    
    iter = 74
    loss = 0.0529
    error rate = 0.1585
    
    iter = 75
    loss = 0.0519
    error rate = 0.1555
    
    iter = 76
    loss = 0.0541
    error rate = 0.1708
    
    iter = 77
    loss = 0.0540
    error rate = 0.1713
    
    iter = 78
    loss = 0.0566
    error rate = 0.1869
    
    iter = 79
    loss = 0.0578
    error rate = 0.1914
    
    iter = 80
    loss = 0.0568
    error rate = 0.1855
    
    iter = 81
    loss = 0.0580
    error rate = 0.1931
    
    iter = 82
    loss = 0.0551
    error rate = 0.1737
    
    iter = 83
    loss = 0.0554
    error rate = 0.1807
    
    iter = 84
    loss = 0.0527
    error rate = 0.1605
    
    iter = 85
    loss = 0.0525
    error rate = 0.1659
    
    iter = 86
    loss = 0.0506
    error rate = 0.1497
    
    iter = 87
    loss = 0.0503
    error rate = 0.1547
    
    iter = 88
    loss = 0.0491
    error rate = 0.1431
    
    iter = 89
    loss = 0.0487
    error rate = 0.1477
    
    iter = 90
    loss = 0.0479
    error rate = 0.1386
    
    iter = 91
    loss = 0.0476
    error rate = 0.1432
    
    iter = 92
    loss = 0.0471
    error rate = 0.1365
    
    iter = 93
    loss = 0.0469
    error rate = 0.1402
    
    iter = 94
    loss = 0.0465
    error rate = 0.1349
    
    iter = 95
    loss = 0.0463
    error rate = 0.1384
    
    iter = 96
    loss = 0.0461
    error rate = 0.1338
    
    iter = 97
    loss = 0.0459
    error rate = 0.1369
    
    iter = 98
    loss = 0.0458
    error rate = 0.1333
    
    iter = 99
    loss = 0.0455
    error rate = 0.1362
    
    iter = 100
    loss = 0.0455
    error rate = 0.1331


## 19.  Evaluate the performance on the testing dataset


```python
xtest, ltest = load(dataset='testing', path='dataset/')
print(f'xtest.shape = {xtest.shape}')
print(f'ltest.shape = {ltest.shape}')
```

    xtest.shape = (784, 10000)
    ltest.shape = (10000,)



```python
# use trained network to evaluate the performance on testing data
xtest = normalize_MNIST_images(xtest)
pred_test = forwardprop_shallow(xtest, nettrain)
print(f'error rate = {eval_perfs(pred_test, ltest):.4f}')
```

    error rate = 0.1286


## 20. Backpropagation based on minibatch gradient descent


```python
def backprop_minibatch_shallow(x, d, net, T, B=100, gamma=.05):
    '''
    Inputs:
        x: data
        d: ground truth
        net: parameters of net
        T: number of epoch
        B: minibatches
        gamma: learning rate
    
    Returns:
        parameters of net
    '''
    N = x.shape[1]
    NB = int((N+B-1)/B)
    lbl = onehot2label(d)
    for t in range(T): # epoch
        shuffled_indices = np.random.permutation(range(N))
        for l in range(NB):
            minibatch_indices = shuffled_indices[B*l:min(B*(l+1), N)]
            # UPDATE NET
            net = update_shallow(x[:,minibatch_indices], d[:,minibatch_indices], net, gamma)
        y = forwardprop_shallow(x, net)
        # DISPLAY LOSS AND PERFS
        print(f'\nepoch = {t+1}')
        print(f'loss = {eval_loss(y, d):.4f}')
        print(f'error rate = {eval_perfs(y, lbl):.4f}')
    return net
```

## 21. Evaluate the performance on the testing dataset based on network with minibatch


```python
netminibatch = backprop_minibatch_shallow(xtrain, dtrain, netinit, 5, B=100)
```

    
    epoch = 1
    loss = 0.0267
    error rate = 0.0790
    
    epoch = 2
    loss = 0.0205
    error rate = 0.0584
    
    epoch = 3
    loss = 0.0191
    error rate = 0.0574
    
    epoch = 4
    loss = 0.0148
    error rate = 0.0432
    
    epoch = 5
    loss = 0.0124
    error rate = 0.0361



```python
# use trained network to evaluate the performance on testing data
pred_test_mini = forwardprop_shallow(xtest, netminibatch)
print(f'error rate = {eval_perfs(pred_test_mini, ltest):.4f}')
```

    error rate = 0.0399


25c59eb0a95f3e87212c9bdbf82412e4edad2204
ipynb
Jupyter Notebook
src/Artificial-Neural-Networks.ipynb
lychengr3x/Artificial-Neuarl-Networks
f2d02b55e55361339940ffdf28e9c745f515160f
MIT
src/Artificial-Neural-Networks.ipynb
lychengr3x/Artificial-Neuarl-Networks
f2d02b55e55361339940ffdf28e9c745f515160f
MIT
src/Artificial-Neural-Networks.ipynb
lychengr3x/Artificial-Neuarl-Networks
f2d02b55e55361339940ffdf28e9c745f515160f
MIT
Qwen/Qwen-72B
1. YES
2. YES