# Neural Networks

<a id='Table_of_Content'></a>

**[1. Neural Networks](#1.Neural_Networks)**
 * [1.1. Perceptron](#1.1.Perceptron)
 * [1.2. Sigmoid](#1.2.Sigmoid)
 
**[2. Neural Networks Architecture](#2.Neural_Networks_Architecture)**

**[3. Training Neural Network](#3.Training_Neural_Network)**
 * [3.1. Forward Propagation](#3.1.Forward_Propagation)
 * [3.2. Compute Error](#3.2.Compute_Error)
 * [3.3. Back Propagation](#3.3.Back_Propagation)
 * [3.4. Gradient Descent](#3.4.Gradient_Descent)
 * [3.5. Computational Graph](#3.5.Computational_Graph)
 * [3.6. Gradient_Checking](#3.6.Gradient_Checking)
 * [3.7. Parameter Update](#3.7.Parameter_Update)
 * [3.8. Learning Rate](#3.8.Learning_Rate)



<a id='1.Neural_Networks'></a>

# 1. Neural Networks



Neural networks (NN) are a broad family of algorithms that have formed the basis for the recent resurgence in the computational field called deep learning. Early work on neural networks actually began in the 1950s and 60s. And just recently, neural network has experienced a resurgence of interest, as deep learning has achieved impressive state-of-the-art results. 

Neural network is basically a mathematical model built from simple functions with changing parameters. Just like a biological neuron has dendrites to receive signals, a cell body to process them, and an axon to send signals out to other neurons, an artificial neuron has a number of input channels, a processing stage, and one output that can branch out to multiple other artificial neurons. Neurons are interconnected and pass message to each other.

To understand neural networks, Let's get started with **Perceptron**.

<center></center>

<a id='1.1.Perceptron'></a>

## 1.1. Perceptron

A perceptron takes several binary inputs, $x_1,x_2,…,x_n$, and produces a single binary output.

<center></center>

The example above shows a perceptron taking three inputs $x_1, x_2, x_3$. Each input is given a $weight$ $W \in \mathbb{R}$ and it serves to express the importance of its corresponding input in the computation of output for that perceptron. The perceptron output, 0 or 1, is determined by the weighted sum $\sum_i w_ix_i$ with respect to a $threshold$ value as follows:

\begin{equation}
    output = \left\{
            \begin{array}{rl}
            0 & \text{if } \sum_iw_ix_i \leq \text{threshold}\\
            1 & \text{if } \sum_iw_ix_i > \text{threshold}
            \end{array} \right.
\end{equation}

The weighted sum can be categorically defined as a dot product between $w$ and $x$ as follows: 

$$\sum_i w_ix_i \equiv w \cdot x$$

where $w$ and $x$ are vectors corresponding to weights and inputs respectively. Introducing a bias term $b \equiv -threshold$ results in

\begin{equation}
    output = \left\{
            \begin{array}{rl}
            0 & \text{if } w \cdot x + b \leq 0\\
            1 & \text{if } w \cdot x + b > 0
            \end{array} \right.
\end{equation}

You can think of the $bias$ as a measure of how easy it is to get the perceptron to output 1. For a perceptron with a high positive $bias$, it is extremely easy for the perceptron to output 1. In constrast, if the $bias$ is relatively  a negative value, it is difficult for the perceptron to output 1.

<center></center>

A way to think about the perceptron is that it is a device that makes **decisions** by weighing up evidence. 



```python
import numpy as np
X = np.array([0, 1, 1])
W = np.array([5, 1, -3])
b=5

def perceptron_neuron(X, W, b):
    return int(X.dot(W)+b > 0)

perceptron_neuron(X,W,b)
```




    1



<a id='1.2.Sigmoid'></a>
## 1.2. Sigmoid 

Small changes to $weights$ and $bias$ of any perceptron in a network can cause the output to flip from 0 to 1 or 1 to 0. This flip can cause the behaviour of the rest of the network to change in a complicated way.


```python
x = np.array([100])
b = np.array([9])
w1 = np.array([-0.08])
w2 = np.array([-0.09])

print(perceptron_neuron(x,w1,b))
print(perceptron_neuron(x,w2,b))

```

    1
    0


The problem above can be overcome by using a Sigmoid neuron. It functions similarly to a Perceptron but modified such that small changes in $weights$ and $bias$ cause only a small change in the output.

As with a Perceptron, a Sigmoid neuron also computes $w \cdot x + b $, but now with the Sigmoid function being incorporated as follows:

\begin{equation}
    z = w \cdot x + b \\
    \sigma(z) = \frac{1}{1+e^{-z}}
\end{equation}

<center></center>

A Sigmoid function produces output between 0 and 1, and the figure below shows  the function. If $z$ is large and positive, the output of a sigmoid neuron approximates to 1, just as it would for a perceptron. Alternatively if $z$ is highly negative, the output approxiates to 0.

<center></center>


```python
def sigmoid(x):
    return 1/(1 + np.exp(-x))

def sigmoid_neuron(X, W, b):
    z = X.dot(W)+b
    return sigmoid(z)

print(sigmoid_neuron(x,w1,b))
print(sigmoid_neuron(x,w2,b))
```

    [0.73105858]
    [0.5]


Click here to go back [Table of Content](#Table_of_Content).

<a id='2.Neural_Networks_Architecture'></a>
                         
# 2. Neural Networks Architecture

A neural network can take many forms. A typical architecture consists of an input layer (leftmost), an output layer (rightmost), and a middle layer (hidden layer). Each layer can have multiple neurons while the number of neurons in the output layer is dependent on the number of classes.

<center></center>


Click here to go back [Table of Content](#Table_of_Content).


```python
# Create dataset
from sklearn.datasets import make_moons, make_circles
import matplotlib.pyplot as plt
import pandas as pd

seed = 123

np.random.seed(seed)
X, y = make_circles(n_samples=1000, factor=.5, noise=.1, random_state=seed)
```


```python
colors = {0:'red', 1:'blue'}
df = pd.DataFrame(dict(x=X[:,0], y=X[:,1], label=y))

fig, ax = plt.subplots()
grouped = df.groupby('label')
for key, group in grouped:
    group.plot(ax=ax, kind='scatter', x='x', y='y', label=key, color=colors[key])
plt.show()
```

<a id='3.Training_Neural_Network'></a>

## 3. Training Neural Network

Previously we have learnt that $weights$ express the importance of variables, and $bias$ is a threshold to control the behaviour of neurons. So, how can we determine these $weights$ and $bias$?

Consider these steps: 
1. Since we do not know the ideal $weights$ and $bias$, we initialize them using random numbers **(parameters initialization).**
2. Let the data flow through the network with these initialized $weights$ and $bias$ to get a predicted output. This process is known as **forward propagation**.
3. Compare the predicted output with the actual output. An error is computed if there is a difference between them. A high error thus indicates that current $weights$ and $bias$ do not give an accurate prediction. **(compute error)**
4. To fix these $weights$ and $bias$, a backward computation is carried out by finding the partial derivative of error with respect to each $weight$ and $bias$ and then updating their values accordingly. This process is known as **backpropagation**.
5. Repeat steps (2) to (4) until the error is below a pre-defined threshold to obtain the optimized $weights$ and $bias$.


<a id='3.1.Forward_Propagation'></a>

### 3.1. Forward Propagation

<center></center>

The model above has two neurons in the input layer, four neurons (also known as **activation units**) in the hidden layer, and one neuron in the output layer. 

$X = [x_1, x_2]$ is the input matrix

$W^{j+1} =$ matrix of $weights$ controlling function mapping from layer $j$ to layer $j + 1$ 

\begin{align}
     W^{j+1} \equiv 
\begin{bmatrix}
w_{11} & w_{12} \\
w_{21} & w_{22} \\
w_{31} & w_{32} \\
w_{41} & w_{42} \\
\end{bmatrix}^{j+1}
\end{align}
$W^{j+1}_{kl}$, $k$ is node in layer $j+1$, $l$ is node in layer $j$

$B^{j+1} = $ matrix of $bias$ controlling function mapping from layer $j$ to layer $j + 1$

\begin{align}
     B^{j+1} \equiv 
\begin{bmatrix}
b_{1} & b_{2} & b_{3} & b_{4} 
\end{bmatrix}^{j+1}
\end{align}

$B^{j+1}_k$, $k$ is node in layer $j + 1$

the activation units can be label as $a_i^j =$ "activation" of unit $i$ in layer $j$. 

if $j=0,\quad a_i^j, $ is equivalent to input layer.  


Finally the activation function in layer 1 can be denode as, 

\begin{align}
     a_1^1 = \sigma(W_{11}^1x_1+W_{12}^1x_2+B^1_{1}) \\
     a_2^1 = \sigma(W_{21}^1x_1+W_{22}^1x_2+B^1_{2}) \\
     a_3^1 = \sigma(W_{31}^1x_1+W_{32}^1x_2+B^1_{3}) \\
     a_4^1 = \sigma(W_{41}^1x_1+W_{42}^1x_2+B^1_{4})
\end{align}

Simplified using vectorization,

\begin{align}
     a^1 = \sigma(X \cdot W^{1T}+B^1) \\
\end{align}

\begin{align}
     output = \sigma( a^1 \cdot W^{2T}+B^2) \\
\end{align}

Implement forward propagation and feed in the generated data

Tips:
- Using numpy function *random.randn()* to generate a Gaussian distribution with mean 0, and variance 1.


```python
#step 1: parameters initialization
def initialize_params():
    params = {
        'W1': np.random.randn(4,2),
        'B1': np.random.randn(1,4),
        'W2': np.random.randn(1,4),
        'B2': np.random.randn(1,1),
    }
    return params
```


```python
#step 2: forward propagation
x_ = np.array([X[0]])
y_ = np.array([y[0]])

np.random.seed(0)
params = initialize_params()

def forward(X, params):
    a1 = sigmoid(X.dot(params['W1'].T)+params['B1'])
    output = sigmoid(a1.dot(params['W2'].T)+params['B2'])
    cache={'a1':a1, 'params': params}
    return output, cache

output, cache = forward(x_, params)
print('Actual output: ', y_)
print('Predicted output: ',output)
```

    Actual output:  [0]
    Predicted output:  [[0.91620003]]


<a id='3.2.Compute_Error'></a>

### 3.2. Compute Error

Error is also known as Loss or Cost.

**Loss function** is usually a function defined on a data point, prediction and label, and measures the penalty. 

**Cost function** is usually more general. It might be a sum of loss functions over your training set plus some model complexity penalty (regularization).

To compute the error, we should first define a $cost$ $function$. For simplicity, we will use **One Half Mean Squared Error** as our cost function. The equation is listed below:

\begin{equation}
    MSE = \frac{1}{2n} \sum (\hat y - y)^2
\end{equation}

where $n$ is the number of training samples, $\hat y$ is the predicted output, and $y$ is the actual output. A low cost results is returned if the predicted output is close to the actual output, which indicates a good measure of accuracy. 


```python
#step 3: cost function
def mse(yhat, y):
    n = yhat.shape[0]
    return (1/(2*n)) * np.sum(np.square(yhat-y))

mse(output, y_)
```




    0.4197112460136526



<a id='3.3.Back_Propagation'></a>

### 3.3. Back Propagation

Now we know that to get a good prediction, the cost should be as low/small as possible. To minimize the cost, we have to tune the $weights$ and $bias$, but how can we do that? Do we go with random trial and error or is there a better way to do it? Fortunately, there is a better way and it is called **Gradient Descent**.



<a id='3.4.Gradient_Descent'></a>

### 3.4. Gradient Descent

Gradient descent is an optimization algorithm that iteratively looks for optimal $weights$ and $bias$ so that the cost gets smaller and eventually equals zero.

In the interative process, the gradient (of the cost function with respect to $weights$ and $bias$) is computed. The gradient is the change in cost when $weights$ and $bias$ are changed. This helps us update $weights$ and $bias$ in the direction in which the cost is minimized.

Let's recall the forward propagation equation:

\begin{align}
     a^1 = \sigma(X \cdot W^{1T}+B^1) \\
     output = \sigma( a^1 \cdot W^{2T}+B^2) \\
     cost = \frac{1}{2n} \sum (output - y)^2
\end{align}

Arrange them into a single equation, and $cost$, $L$ can be defined as follows:

\begin{align}
     L = \frac{1}{2n} \sum (\sigma( \sigma(X \cdot W^{1T}+B^1) \cdot W^{2T}+B^2)  - y)^2
\end{align}

From the equation, we want to find the gradient or derivative of $L$ with respect to $W^1, W^2, B^1, B^2$.
\begin{align}
      \frac{\partial L}{\partial W^1}, \frac{\partial L}{\partial W^2}, \frac{\partial L}{\partial B^1}, \frac{\partial L}{\partial B^2}
\end{align}

Computation of partial derivatives of $L$ with respect to the $weights$ and $bias$ can become very complex if the layer of the network grows. To make it simple, we can actually break the equation into smaller compenents and use **chain rule** to derive the partial derivative.

<a id='3.5.Computational_Graph'></a>

### 3.5. Computational Graph
Eventually, we can think of the forward propagation equation as a computational graph.

<center></center>
<center></center>
<center></center>





#### Scalar Example
\begin{align}
     a_1^1 = \sigma(W_{11}^1x_1+W_{12}^1x_2+B^1_{1}) \quad \equiv \quad \frac{1}{1+\exp^{-(W_1x_1+W_2x_2+B_1)}}
\end{align}

<center></center>

#### Note:
\begin{align}
   L \quad &\rightarrow \quad \frac{\partial L}{\partial L} = 1 \\
   L = \frac{1}{2n} \sum (output - y)^2 \quad &\rightarrow \quad \frac{\partial L}{\partial output} = \frac{1}{n} (output - y) \\
     f(x) = e^x \quad &\rightarrow \quad  \frac{\partial f}{\partial x} = e^x \\
     f(x) = xy \quad &\rightarrow \quad  \frac{\partial f}{\partial x} = y, \quad \frac{\partial f}{\partial y} = x \\
     f(x) = 1/x \quad &\rightarrow \quad  \frac{\partial f}{\partial x} = -1/x^2 \\
     f(x) = x+c \quad &\rightarrow \quad  \frac{\partial f}{\partial x} = 1 \\
     \sigma(x) = \frac{1}{1+e^{-x}} \quad &\rightarrow \quad \frac{\partial f}{\partial \sigma} = \sigma(1-\sigma)
\end{align}


```python
print('X: ', X[0])
print('W11: ', params['W1'][0])
print('B11: ', params['B1'][0][0])

# how is the graph look like in our case?
# calculate forward and backward flows
```

    X:  [-0.08769568  1.08597835]
    W11:  [1.76405235 0.40015721]
    B11:  -0.10321885179355784




#### A Vectorized Example

\begin{align}
     a^1 = \sigma(X \cdot W^{1T}+B^1) \\
\end{align}

<center></center>

#### Note:
\begin{align}
    q = X\cdot(W^T) \quad &\rightarrow \quad \frac{\partial f}{\partial X} = \frac{\partial f}{\partial q} \cdot W , \quad \frac{\partial f}{\partial W} = X^T \cdot \frac{\partial f}{\partial q} \\
    l = q+B \quad &\rightarrow \quad \frac{\partial f}{\partial B} = \begin{bmatrix}1 & 1 \end{bmatrix} \cdot \frac{\partial f}{\partial l} , \quad \frac{\partial f}{\partial q} = \frac{\partial f}{\partial l}
\end{align}



```python
def dmse(output, y):
    return (output - y)/output.shape[0]
    
def backward(X, output, y, cache):
    grads={}
    a1=cache['a1']
    params=cache['params']
    
    dloss = dmse(output, y)
    
    doutput = output*(1-output)*dloss
    
    #compute gradient of B2 and W2
    dW2 = a1.T.dot(doutput)
    dB2 = np.sum(doutput, axis=0, keepdims=True)
    
    dX2 = doutput.dot(params['W2'])
    da1 = a1*(1-a1)*dX2
    
    #compute gradient of B1 and W1
    dW1 = X.T.dot(da1)
    dB1 = np.sum(da1, axis=0, keepdims=True)
    
    grads['W1'] = dW1.T
    grads['W2'] = dW2.T
    grads['B1'] = dB1
    grads['B2'] = dB2
    
    return grads
    
```


```python
X_ = X[:3]
Y_ = y[:3].reshape(-1,1)

def step(X,y,params):

    output, cache = forward(X, params)

    cost = mse(output, y)

    grads = backward(X, output, y, cache)

    return (cost, grads)

np.random.seed(0)
params = initialize_params()

cost, grads = step(X_, Y_, params)

print(cost)
print(grads)
```

    0.4184302966572367
    {'W1': array([[-0.00208634,  0.0076568 ],
           [-0.00053028,  0.00068688],
           [-0.0004186 ,  0.00346495],
           [-0.00170289,  0.00300797]]), 'W2': array([[0.03212358, 0.05777044, 0.02222057, 0.05217427]]), 'B1': array([[0.01009638, 0.00114439, 0.00470392, 0.00425563]]), 'B2': array([[0.07033491]])}


<a id='3.6.Gradient_Checking'></a>

### 3.6. Gradient Checking
We can use Numerical Gradient to evaluate the gradient (Analytical gradient) that we have calculated.

<center></center>

Consider the image above, where the red line is our function, the blue line is the gradient derived from the point $x$, the green line is the approximated gradient from the point of $x$, and $h$ is the step size. It can then be shown that:


$$ \frac{\partial f}{\partial x} \approx \frac{Y_C-Y_B}{X_C-X_B} \quad = \quad \frac{f(x+h) - f(x-h)}{(x+h)-(x-h)} \quad = \quad \frac{f(x+h) - f(x-h)}{2h} $$



##### EXAMPLE



```python
w1 = 3; x1 = 1; w2 = 2; x2 = -2; b1 = 2

h = 1e-4

def f(w1, x1, w2, x2, b1):
    linear = (w1*x1)+(w2*x2)+b1
    return 1/(1+np.exp(-linear))


num_grad_x1 = (f(w1, x1+h, w2, x2, b1) - f(w1, x1-h, w2, x2, b1))/(2*h)
print(num_grad_x1)

```

    0.5898357981348745



```python
# vectorized gradient checking
def gradient_check(f, x, h=0.00001):
    grad = np.zeros_like(x)
  # iterate over all indexes in x
    it = np.nditer(x, flags=['multi_index'])
    while not it.finished:
    # evaluate function at x+h
        ix = it.multi_index
        oldval = x[ix]
        x[ix] = oldval + h # increment by h
        fxph = f(x) # evalute f(x + h)
        x[ix] = oldval - h
        fxnh = f(x) # evaluate f(x - h)
        x[ix] = oldval # restore

        # compute the partial derivative with centered formula
        grad[ix] = ((fxph - fxnh) / (2 * h)).sum() # the slope
        it.iternext() # step to next dimension

    return grad

def rel_error(x, y):
    return np.max(np.abs(x - y) / (np.maximum(1e-8, np.abs(x) + np.abs(y))))
```


```python
for param_name in grads:
    f = lambda W: step(X_, Y_, params)[0]
    
    param_grad_num = gradient_check(f, params[param_name])
    print('%s max relative error: %e' % (param_name, rel_error(param_grad_num, grads[param_name])))

```

    W1 max relative error: 3.437088e-09
    W2 max relative error: 5.328225e-11
    B1 max relative error: 2.182159e-09
    B2 max relative error: 3.440171e-11


<a id='3.7.Parameter_Update'></a>

### 3.7. Parameter Update
After getting the gradient, parameters are updated depend on the gradients; if it is positive, then the updated parameters reduces in value, and if it is negative, then the updated parameter increases in value. Regardless of the gradient, the main goal is to reach the global minimum.

<a id='3.8.Learning_Rate'></a>

### 3.8. Learning Rate
Large or small updates are controlled by the learning rate known as $\alpha$. Hence, gradient descent equations are as follows:

\begin{align}
    W^1 &= W^1 - \alpha * \frac{\partial L}{\partial W^1} \\
    B^1 &= B^1 - \alpha * \frac{\partial L}{\partial B^1} \\
    W^2 &= W^2 - \alpha * \frac{\partial L}{\partial W^2} \\
    B^2 &= B^2 - \alpha * \frac{\partial L}{\partial B^2} \\
\end{align}



```python
def update_parameter(params, grads, learning_rate):
    params['W1'] += -learning_rate * grads['W1']
    params['B1'] += -learning_rate * grads['B1']
    params['W2'] += -learning_rate * grads['W2']
    params['B2'] += -learning_rate * grads['B2']
    
```


```python
params = initialize_params()

def train(X,y,learning_rate=0.1,num_iters=30000,batch_size=256):
    num_train = X.shape[0]
    costs = []
    for it in range(num_iters):
        random_indices = np.random.choice(num_train, batch_size)
        X_batch = X[random_indices]
        y_batch = y[random_indices]
        
        cost, grads = step(X_batch, y_batch, params)
        costs.append(cost)
        
        # update parameters 
        update_parameter(params, grads, learning_rate)
            
    return costs

costs = train(X,y.reshape(-1,1))
plt.plot(costs)
```


```python
def predict(X, params):
    W1 = params['W1']
    B1 = params['B1']
    W2 = params['W2']
    B2 = params['B2']
    
    output, _ = forward(X, params)
    return output
```


```python
# test on training samples
y_pred = []
for i in range(len(X)):
    pred = np.squeeze(predict(X[i],params)).round()
    y_pred.append(pred)
    
plt.scatter(X[:,0], X[:,1], c=y_pred, linewidths=0, s=20);
```


```python
# test on new samples
X_new, _ = make_circles(n_samples=1000, factor=.5, noise=.1)
y_pred = []

for i in range(len(X_new)):
    pred = np.squeeze(predict(X_new[i],params)).round()
    y_pred.append(pred)
    
plt.scatter(X_new[:,0], X_new[:,1], c=y_pred, linewidths=0, s=20);
```

# References:

- http://neuralnetworksanddeeplearning.com/

- http://cs231n.github.io/optimization-2/

- http://kineticmaths.com/index.php?title=Numerical_Differentiation

- https://google-developers.appspot.com/machine-learning/crash-course/backprop-scroll/

Click here to go back [Table of Content](#Table_of_Content).

cbf67ea71443059e40d8aa49aedc89b3082dd38d
ipynb
Jupyter Notebook
01-neural-networks/neural-network-training.ipynb
williamardianto/neural-network-from-scratch
9e01bcd2d6663a074e8152a83ffe671992b14a89
MIT
01-neural-networks/neural-network-training.ipynb
williamardianto/neural-network-from-scratch
9e01bcd2d6663a074e8152a83ffe671992b14a89
MIT
01-neural-networks/neural-network-training.ipynb
williamardianto/neural-network-from-scratch
9e01bcd2d6663a074e8152a83ffe671992b14a89
MIT
Qwen/Qwen-72B
1. YES
2. YES