# Neural Network

The following example follows Andrew Trask’s old blog post, which is nice because it tries to demonstrate a neural net in very few lines of code, much like this document’s goal.

The data setup is very simple (only 4 observations!), and I keep the Python code essentially identical outside of very slight cosmetic (mostly name/space) changes. For more detail I suggest following the original posts, but I’ll add some context here and there. In addition, see the logistic regression chapter and gradient descent chapter.

- https://iamtrask.github.io/2015/07/12/basic-python-network/
- https://iamtrask.github.io/2015/07/27/python-network-part2/

## Example 1

### Python

In this initial example, while it can serve as instructive starting point for backpropagation, we’re not really using what most would call a neural net, but rather just an alternative way to estimate a logistic regression. `layer_1`

in this case is just the linear predictor after a nonlinear transformation (sigmoid).

Note that in this particular example however, that the first column is perfectly correlated with the target `y`

, which would cause a problem if no regularization were applied (a.k.a. separation).

Description of the necessary objects from the blog. These will be consistent throughout the demo.

XInput dataset matrix where each row is a training exampleyOutput dataset matrix where each row is a training examplelayer_0First Layer of the Network, specified by the input datalayer_1Second Layer of the Network, otherwise known as the hidden layersynapse_0First layer of weights, Synapse 0, connecting layer_0 to layer_1.*Elementwise multiplication, so two vectors of equal size are multiplying corresponding values 1-to-1 to generate a final vector of identical size.-Elementwise subtraction, so two vectors of equal size are subtracting corresponding values 1-to-1 to generate a final vector of identical size.x.dot(y)If x and y are vectors, this is a dot product. If both are matrices, it’s a matrix-matrix multiplication. If only one is a matrix, then it’s vector matrix multiplication.

```
import numpy as np
# sigmoid function
def nonlin(x, deriv = False):
if(deriv == True):
return x*(1 - x)
return 1/(1 + np.exp(-x))
# input dataset
= np.array([
X 0, 0, 1],
[0, 1, 1],
[1, 0, 1],
[1, 1, 1]
[
])
# output dataset
= np.array([[0, 0, 1, 1]]).T
y
# seed random numbers to make calculation
# deterministic (just a good practice)
1)
np.random.seed(
# initialize weights randomly with mean 0 (or just use np.random.uniform)
= 2*np.random.random((3, 1)) - 1
synapse_0
for iter in np.arange(10000):
# forward propagation
= X
layer_0 = nonlin(np.dot(layer_0, synapse_0))
layer_1
# how much did we miss?
= y - layer_1
layer_1_error
# multiply how much we missed by the
# slope of the sigmoid at the values in layer_1
= layer_1_error * nonlin(layer_1, True)
l1_delta
# update weights
+= np.dot(layer_0.T, l1_delta)
synapse_0
print("Output After Training:")
```

`Output After Training:`

`print(np.append(layer_1, y, axis = 1))`

```
[[0.00966449 0. ]
[0.00786506 0. ]
[0.99358898 1. ]
[0.99211957 1. ]]
```

### R

For R I make a couple changes, but it should be easy to follow from the original Python, and I keep the original comments. I convert the code into a function so that settings can be altered more easily.

```
= matrix(
X c(0, 0, 1,
0, 1, 1,
1, 0, 1,
1, 1, 1),
nrow = 4,
ncol = 3,
byrow = TRUE
)
# output dataset
= c(0, 0, 1, 1)
y
# seed random numbers to make calculation
# deterministic (just a good practice)
set.seed(1)
# initialize weights randomly with mean 0
= matrix(runif(3, min = -1, max = 1), 3, 1)
synapse_0
# sigmoid function
<- function(x, deriv = FALSE) {
nonlin if (deriv)
* (1 - x)
x else
plogis(x)
}
<- function(X, y, synapse_0, maxiter = 10000) {
nn_1
for (iter in 1:maxiter) {
# forward propagation
= X
layer_0 = nonlin(layer_0 %*% synapse_0)
layer_1
# how much did we miss?
= y - layer_1
layer_1_error
# multiply how much we missed by the
# slope of the sigmoid at the values in layer_1
= layer_1_error * nonlin(layer_1, deriv = TRUE)
l1_delta
# update weights
= synapse_0 + crossprod(layer_0, l1_delta)
synapse_0
}
list(layer_1 = layer_1, layer_1_error = layer_1_error, synapse_0 = synapse_0)
}
= nn_1(X, y, synapse_0)
fit_nn
message("Output After Training: \n",
paste0(capture.output(cbind(fit_nn$layer_1, y)), collapse = '\n'))
```

```
Output After Training:
y
[1,] 0.009670417 0
[2,] 0.007864211 0
[3,] 0.993590571 1
[4,] 0.992115835 1
```

A key takeaway from this demonstration regards the update step. Using the derivative, we are getting the slope of the sigmoid function at the point of interest, e.g. the three points below. If a predicted probability (`layer 1`

) is close to zero or 1, this suggests the prediction ‘confidence’ is high, the slopes are shallow, and the result is that there is less need for updating. For the others, e.g. close to `x = 0`

, there will be relatively more updating, as the error will be greater. The last step is to compute the weight updates for each weight for each observation, sum them, and update the weights accordingly. So our update is a function of the error and the slope. If both are small, then there will be little update, but if there is larger error and/or less confident prediction, the result will ultimately be a larger update.

## Example 2

In the following example, we have similarly simple data, but technically a bit harder problem to estimate. In this case, `y = 1`

only if we apply the XOR function to columns 1 and 2. This makes the relationship between inputs and output a nonlinear one.

To account for the nonlinearity we will have two ‘hidden’ layers, the first of four ‘nodes’, and the second a single node, which relates the transformed information to the target variable. In the hidden layer, we can think of each node as a linear combination of the input, with weights represented by paths (in the original post these weights/connections are collectively called ‘synapses’). This is followed by a nonlinear transformation for each node (e.g. sigmoid). The final predictions are a linear combination of the hidden layer, followed by another nonlinear transformation (sigmoid again). The last transformation needs to put the result between 0 and 1, but you could try other transformations for the first hidden layer. In fact sigmoid is rarely used except for the final layer.

### Python

For the following, I remove defining `layer_0`

, as it is just the input matrix `X`

. But otherwise, only minor cleanup and more explicit names as before.

```
import numpy as np
def nonlin(x, deriv = False):
if(deriv == True):
return x*(1 - x)
return 1/(1 + np.exp(-x))
= np.array([
X 0, 0, 1],
[0, 1, 1],
[1, 0, 1],
[1, 1, 1]
[
])
= np.array([
y 0],
[1],
[1],
[0]
[
])
1)
np.random.seed(
# randomly initialize our weights with mean 0
= 2*np.random.random((3, 4)) - 1
synapse_0 = 2*np.random.random((4, 1)) - 1
synapse_1
for j in np.arange(30000):
# Feed forward through layers 0, 1, and 2
= nonlin(np.dot(X, synapse_0))
layer_1 = nonlin(np.dot(layer_1, synapse_1))
layer_2
# how much did we miss the target value?
= y - layer_2
layer_2_error
if (j% 10000) == 0:
print("Error:" + str(np.mean(np.abs(layer_2_error))))
# in what direction is the target value?
# were we really sure? if so, don't change too much.
= layer_2_error*nonlin(layer_2, deriv=True)
layer_2_delta
# how much did each layer_1 value contribute to the layer_2 error (according to the weights)?
= layer_2_delta.dot(synapse_1.T)
layer_1_error
# in what direction is the target layer_1?
# were we really sure? if so, don't change too much.
= layer_1_error * nonlin(layer_1, deriv=True)
layer_1_delta
+= layer_1.T.dot(layer_2_delta)
synapse_1 += X.T.dot(layer_1_delta) synapse_0
```

```
Error:0.4964100319027255
Error:0.008584525653247153
Error:0.005789459862507812
```

`print('Final error: ' + str(np.round(np.mean(np.abs(layer_2_error)), 5)))`

`Final error: 0.00463`

`round(layer_1, 3) np.`

```
array([[0.696, 0.124, 0.923, 0.996],
[0.18 , 0. , 0.017, 0.893],
[0.995, 0.888, 0.023, 0.822],
[0.947, 0.026, 0. , 0.122]])
```

`round(np.append(layer_2, y, axis = 1), 3) np.`

```
array([[0.004, 0. ],
[0.995, 1. ],
[0.996, 1. ],
[0.006, 0. ]])
```

`round(synapse_0, 3) np.`

```
array([[ 4.407, 4.028, -6.225, -4.092],
[-2.342, -5.706, -6.53 , -3.506],
[ 0.826, -1.959, 2.488, 5.623]])
```

`round(synapse_1, 3) np.`

```
array([[ -6.611],
[ 6.857],
[-10.069],
[ 7.501]])
```

### R

In general, different weights can ultimately produce similar predictions, and combinations for a particular node are arbitrary. For example, node 1 in one network might become more like node 3 on a separate run with different starting points. So things are somewhat more difficult to compare across runs, let alone across different tools. However, what matters more is the final prediction, and you should see essentially the same result for R and Python. Again we create a function for repeated use.

```
= matrix(
X c(0, 0, 1,
0, 1, 1,
1, 0, 1,
1, 1, 1),
nrow = 4,
ncol = 3,
byrow = TRUE
)
= matrix(as.integer(xor(X[,1], X[,2])), ncol = 1) # make the relationship explicit
y
set.seed(1)
# or do randomly in same fashion
= matrix(runif(12, -1, 1), 3, 4)
synapse_0 = matrix(runif(12, -1, 1), 4, 1)
synapse_1
# synapse_0
# synapse_1
<- function(
nn_2
X,
y,
synapse_0_start,
synapse_1_start,maxiter = 30000,
verbose = TRUE
) {
= synapse_0_start
synapse_0 = synapse_1_start
synapse_1
for (j in 1:maxiter) {
= plogis(X %*% synapse_0) # 4 x 4
layer_1 = plogis(layer_1 %*% synapse_1) # 4 x 1
layer_2
# how much did we miss the target value?
= y - layer_2
layer_2_error
if (verbose && (j %% 10000) == 0) {
message(glue::glue("Error: {mean(abs(layer_2_error))}"))
}
# in what direction is the target value?
# were we really sure? if so, don't change too much.
= (y - layer_2) * (layer_2 * (1 - layer_2))
layer_2_delta
# how much did each l1 value contribute to the l2 error (according to the weights)?
= layer_2_delta %*% t(synapse_1)
layer_1_error
# in what direction is the target l1?
# were we really sure? if so, don't change too much.
= tcrossprod(layer_2_delta, synapse_1) * (layer_1 * (1 - layer_1))
layer_1_delta
# update
= synapse_1 + crossprod(layer_1, layer_2_delta)
synapse_1 = synapse_0 + crossprod(X, layer_1_delta)
synapse_0
}
list(
layer_1_error = layer_1_error,
layer_2_error = layer_2_error,
synapse_0 = synapse_0,
synapse_1 = synapse_1,
layer_1 = layer_1,
layer_2 = layer_2
) }
```

With function in place, we’re ready to try it out.

```
= nn_2(
fit_nn
X,
y,synapse_0_start = synapse_0,
synapse_1_start = synapse_1,
maxiter = 30000
)
```

`Error: 0.0105538166393651`

`Error: 0.00729252475321203`

`Error: 0.0058973637409426`

`::glue('Final error: {round(mean(abs(fit_nn$layer_2_error)), 5)}') glue`

`Final error: 0.0059`

`round(fit_nn$layer_1, 3)`

```
[,1] [,2] [,3] [,4]
[1,] 0.259 0.889 0.364 0.445
[2,] 0.000 0.037 0.978 0.030
[3,] 0.946 1.000 0.984 0.020
[4,] 0.016 0.984 1.000 0.001
```

`round(cbind(fit_nn$layer_2, y), 3)`

```
[,1] [,2]
[1,] 0.002 0
[2,] 0.993 1
[3,] 0.994 1
[4,] 0.008 0
```

`round(fit_nn$synapse_0, 3)`

```
[,1] [,2] [,3] [,4]
[1,] 3.915 7.364 4.705 -3.669
[2,] -6.970 -5.351 4.337 -3.242
[3,] -1.050 2.079 -0.559 -0.221
```

`round(fit_nn$synapse_1, 3)`

```
[,1]
[1,] 10.988
[2,] -10.733
[3,] 5.576
[4,] -2.987
```

### Comparison

Let’s compare our results to the nnet package that comes with the base R installation. Note that it will include intercepts (a.k.a. *biases*) for each node, so it’s estimating more parameters in total.

`= nnet::nnet(X, y, size = 4) fit_nnet `

```
# weights: 21
initial value 1.152323
iter 10 value 0.998588
iter 20 value 0.558917
final value 0.000000
converged
```

`data.frame(coef(fit_nnet))`

```
coef.fit_nnet.
b->h1 5320.330
i1->h1 -49025.055
i2->h1 -15547.859
i3->h1 5320.491
b->h2 12329.141
i1->h2 -13914.491
i2->h2 -13782.790
i3->h2 12328.763
b->h3 22292.322
i1->h3 -18601.227
i2->h3 34270.828
i3->h3 22292.567
b->h4 22269.956
i1->h4 -3504.603
i2->h4 30916.959
i3->h4 22269.420
b->o 34687.508
h1->o -37130.643
h2->o 32075.031
h3->o -36018.926
h4->o -18446.081
```

`fitted(fit_nnet)`

```
[,1]
[1,] 0
[2,] 1
[3,] 1
[4,] 0
```

## Example 3

The next example follows the code from the second post listed in the introduction. A couple of changes are seen here. We have a notably larger hidden layer (size 32). We also split the previous nonlin function into two parts. And finally, we add an `alpha`

parameter, which is akin to the *learning rate* in gradient descent (stepsize in our previous implementation). It basically puts a control on the gradient so that we hopefully don’t make too large a jump in our estimated weights on each update. However, we can show what will happen as a result of setting the parameter to small and large values.

### Python

As before, I make names more explicit, and other very minor updates to the original code.

```
import numpy as np
= [0.001, 0.01, 0.1, 1, 10, 100, 1000]
alphas = 32
hidden_size
# compute sigmoid nonlinearity
def sigmoid(x):
= 1/(1 + np.exp(-x))
output return output
# convert output of sigmoid function to its derivative
def sigmoid_output_to_derivative(output):
return output*(1 - output)
= np.array([
X 0, 0, 1],
[0, 1, 1],
[1, 0, 1],
[1, 1, 1]
[
])
= np.array([
y 0],
[1],
[1],
[0]
[
])
for alpha in alphas:
print("\nTraining With Alpha:" + str(alpha))
1)
np.random.seed(
# randomly initialize our weights with mean 0
= 2*np.random.random((3, hidden_size)) - 1
synapse_0 = 2*np.random.random((hidden_size, 1)) - 1
synapse_1
for j in np.arange(30000):
# Feed forward through layers input, 1, and 2
= sigmoid(np.dot(X, synapse_0))
layer_1 = sigmoid(np.dot(layer_1, synapse_1))
layer_2
# how much did we miss the target value?
= layer_2 - y
layer_2_error
if (j% 10000) == 0:
print("Error after "+ str(j) +" iterations:" +
str(np.mean(np.abs(layer_2_error))))
# in what direction is the target value?
# were we really sure? if so, don't change too much.
= layer_2_error*sigmoid_output_to_derivative(layer_2)
layer_2_delta
# how much did each l1 value contribute to the l2 error (according to the weights)?
= layer_2_delta.dot(synapse_1.T)
layer_1_error
# in what direction is the target l1?
# were we really sure? if so, don't change too much.
= layer_1_error * sigmoid_output_to_derivative(layer_1)
layer_1_delta
-= alpha * (layer_1.T.dot(layer_2_delta))
synapse_1 -= alpha * (X.T.dot(layer_1_delta)) synapse_0
```

```
Training With Alpha:0.001
Error after 0 iterations:0.49643992250078794
Error after 10000 iterations:0.49104946812904954
Error after 20000 iterations:0.4849763070274596
Training With Alpha:0.01
Error after 0 iterations:0.49643992250078794
Error after 10000 iterations:0.35637906164802124
Error after 20000 iterations:0.14693984546475994
Training With Alpha:0.1
Error after 0 iterations:0.49643992250078794
Error after 10000 iterations:0.030540490838554993
Error after 20000 iterations:0.019063872533418427
Training With Alpha:1
Error after 0 iterations:0.49643992250078794
Error after 10000 iterations:0.007360522342493724
Error after 20000 iterations:0.004972517050388164
Training With Alpha:10
Error after 0 iterations:0.49643992250078794
Error after 10000 iterations:0.0022489756615412743
Error after 20000 iterations:0.0015326357139237015
Training With Alpha:100
Error after 0 iterations:0.49643992250078794
Error after 10000 iterations:0.5
Error after 20000 iterations:0.5
Training With Alpha:1000
Error after 0 iterations:0.49643992250078794
Error after 10000 iterations:0.5
Error after 20000 iterations:0.5
```

Since `alpha = 10`

seems reasonable, let’s inspect the results just for that. Note that this value being the ‘best’ likely won’t hold on another random run, and in general we would assess this *hyperparameter* using cross-validation or other means.

```
import numpy as np
= 10
alpha
# randomly initialize our weights with mean 0
1)
np.random.seed(
= 2*np.random.random((3, hidden_size)) - 1
synapse_0 = 2*np.random.random((hidden_size, 1)) - 1
synapse_1
for j in np.arange(30000):
# Feed forward through layers input, 1, and 2
= sigmoid(np.dot(X, synapse_0))
layer_1 = sigmoid(np.dot(layer_1, synapse_1))
layer_2
# how much did we miss the target value?
= layer_2 - y
layer_2_error
# in what direction is the target value?
# were we really sure? if so, don't change too much.
= layer_2_error*sigmoid_output_to_derivative(layer_2)
layer_2_delta
# how much did each l1 value contribute to the l2 error (according to the weights)?
= layer_2_delta.dot(synapse_1.T)
layer_1_error
# in what direction is the target l1?
# were we really sure? if so, don't change too much.
= layer_1_error * sigmoid_output_to_derivative(layer_1)
layer_1_delta
-= alpha * (layer_1.T.dot(layer_2_delta))
synapse_1 -= alpha * (X.T.dot(layer_1_delta))
synapse_0
round(layer_2, 4), y, axis = 1) np.append(np.
```

```
array([[0.0011, 0. ],
[0.9989, 1. ],
[0.9989, 1. ],
[0.0015, 0. ]])
```

### R

The following has little to no modification, but again creates a function for easier manipulation of inputs.

```
# input dataset
= matrix(
X c(0, 0, 1,
0, 1, 1,
1, 0, 1,
1, 1, 1),
nrow = 4,
ncol = 3,
byrow = TRUE
)
# output dataset
= matrix(c(0, 1, 1, 0), ncol = 1)
y
= c(0.001, 0.01, 0.1, 1, 10, 100, 1000)
alphas = 32
hidden_size
# compute sigmoid nonlinearity
= plogis # already part of base R, no function needed
sigmoid
# convert output of sigmoid function to its derivative
<- function(output) {
sigmoid_output_to_derivative * (1 - output)
output
}
<- function(
nn_3
X,
y,
hidden_size,
alpha,maxiter = 30000,
show_messages = FALSE
) {
for (val in alpha) {
if(show_messages)
message(glue::glue("Training With Alpha: {val}"))
set.seed(1)
# randomly initialize our weights with mean 0
= matrix(runif(3 * hidden_size, -1, 1), 3, hidden_size)
synapse_0 = matrix(runif(hidden_size), hidden_size, 1)
synapse_1
for (j in 1:maxiter) {
# Feed forward through layers input, 1, and 2
= sigmoid(X %*% synapse_0)
layer_1 = sigmoid(layer_1 %*% synapse_1)
layer_2
# how much did we miss the target value?
= layer_2 - y
layer_2_error
if ((j %% 10000) == 0 & show_messages) {
message(glue::glue("Error after {j} iterations: {mean(abs(layer_2_error))}"))
}
# in what direction is the target value?
# were we really sure? if so, don't change too much.
= layer_2_error * sigmoid_output_to_derivative(layer_2)
layer_2_delta
# how much did each l1 value contribute to the l2 error (according to the weights)?
= layer_2_delta %*% t(synapse_1)
layer_1_error
# in what direction is the target l1?
# were we really sure? if so, don't change too much.
= layer_1_error * sigmoid_output_to_derivative(layer_1)
layer_1_delta
= synapse_1 - val * crossprod(layer_1, layer_2_delta)
synapse_1 = synapse_0 - val * crossprod(X, layer_1_delta)
synapse_0
}
}
list(
layer_1_error = layer_1_error,
layer_2_error = layer_2_error,
synapse_0 = synapse_0,
synapse_1 = synapse_1,
layer_1 = layer_1,
layer_2 = layer_2
) }
```

With function in place, we are now ready to try it out. You can change whether you show the messages or not to compare with the Python. I don’t in this case for the sake of the document, but it won’t be overwhelming if you do interactively, and is recommended. We will also see that `alpha = 10`

is the better option under the default settings as it was with the Python code. Note that this simple demo code will probably not work well beyond a hidden size of 50 (add more iterations even then).

```
set.seed(1)
= nn_3(
fit_nn
X,
y,hidden_size = 32,
maxiter = 30000,
alpha = alphas,
show_messages = FALSE
)
```

Let’s rerun at `alpha = 10`

.

```
set.seed(1)
= nn_3(
fit_nn
X,
y,hidden_size = 32,
alpha = 10,
show_messages = TRUE
)
```

`Training With Alpha: 10`

`Error after 10000 iterations: 0.00483502508464005`

`Error after 20000 iterations: 0.00207985328313795`

`Error after 30000 iterations: 0.00152092933542719`

`cbind(round(fit_nn$layer_2, 4), y)`

```
[,1] [,2]
[1,] 0.0013 0
[2,] 0.9985 1
[3,] 0.9986 1
[4,] 0.0018 0
```

### Comparison

Again we can compare to nnet.

`= nnet::nnet(X, y, size = 32) fit_nnet `

```
# weights: 161
initial value 1.851187
iter 10 value 0.936564
iter 20 value 0.008594
final value 0.000065
converged
```

`fitted(fit_nnet)`

```
[,1]
[1,] 0.0001274475
[2,] 0.9933000550
[3,] 0.9959851931
[4,] 0.0019116939
```

## Source

This was never on the original repo but I may put it there eventually.