From One Neuron to a Network

import numpy as np

rng = np.random.default_rng(3)

def relu(x):
    return np.maximum(0.0, x)

weights = np.array([0.6, -1.2, 0.3])
bias = 0.1
inputs = np.array([1.0, 0.5, -2.0])

weighted_sum = float(inputs @ weights + bias)
output = float(relu(weighted_sum))

print(f'inputs        {inputs}')
print(f'weights       {weights}')
print(f'weighted sum  {weighted_sum:.3f}')
print(f'after ReLU    {output:.3f}')

import numpy as np

rng = np.random.default_rng(5)

W1 = rng.normal(0, 0.5, (4, 3))
W2 = rng.normal(0, 0.5, (3, 2))
x = rng.normal(0, 1, 4)

# Two layers with NO activation.
stacked = (x @ W1) @ W2

# One layer with the combined weight matrix.
combined = x @ (W1 @ W2)

print(f'two linear layers: {stacked.round(6)}')
print(f'one layer:         {combined.round(6)}')
print(f'identical: {bool(np.allclose(stacked, combined))}')
print()
print('A hundred linear layers are exactly one linear layer. The')
print('non-linearity is what makes depth mean anything.')

import numpy as np
import pandas as pd

x = np.array([-3.0, -1.0, -0.1, 0.0, 0.1, 1.0, 3.0])

def sigmoid(v):
    return 1 / (1 + np.exp(-v))

table = pd.DataFrame({
    'x': x,
    'relu': np.maximum(0, x).round(3),
    'leaky_relu': np.where(x > 0, x, 0.01 * x).round(3),
    'gelu': (0.5 * x * (1 + np.tanh(np.sqrt(2 / np.pi) * (x + 0.044715 * x ** 3)))).round(3),
    'tanh': np.tanh(x).round(3),
    'sigmoid': sigmoid(x).round(3),
})
print(table.to_string(index=False))
print()
print('relu     - the default for hidden layers: cheap, and its')
print('           gradient is 1 or 0, which does not vanish')
print('gelu     - a smooth relu; the default inside transformers')
print('sigmoid  - the OUTPUT layer for a binary probability')
print('tanh     - centred on zero; mostly historical for hidden layers')

import numpy as np

rng = np.random.default_rng(7)
n = 2000

# The input to a hidden layer is the PREVIOUS ReLU's output, so it is
# never negative. That matters: against a zero-mean input, a negative
# weight mean cancels out and nothing dies. It is the bias that kills,
# and a single too-large update is what drives it down.
inputs = np.maximum(0, rng.normal(0, 1, (n, 50)))
weights = rng.normal(0, 0.3, (50, 100))
bias = rng.normal(-2.0, 1.5, 100)      # after one oversized step

activations = np.maximum(0, inputs @ weights + bias)

is_dead = activations.max(axis=0) == 0
print(f'{is_dead.mean():.0%} of the 100 units never fire on any of {n} examples')
print(f'mean bias, dead units {bias[is_dead].mean():+.2f}')
print(f'mean bias, live units {bias[~is_dead].mean():+.2f}')
print('Their gradient is permanently zero. Leaky ReLU or GELU avoids')
print('this, and so does sensible initialisation.')

import numpy as np
import pandas as pd
from sklearn.neural_network import MLPClassifier
from sklearn.model_selection import train_test_split
from sklearn.metrics import roc_auc_score
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline

rng = np.random.default_rng(11)
n = 6000

# A problem that needs composition: a checkerboard in two dimensions
# plus a third interacting feature.
X = rng.uniform(-3, 3, (n, 3))
y = ((np.sin(X[:, 0] * 2) * np.sin(X[:, 1] * 2) > 0).astype(int)
     ^ (X[:, 2] > 0).astype(int))
X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.3, random_state=0, stratify=y)

rows = []
for shape in [(2,), (16,), (256,), (16, 16), (16, 16, 16), (64, 64)]:
    parameters = sum(a * b for a, b in zip((3,) + shape, shape + (1,)))
    model = Pipeline([('s', StandardScaler()),
                      ('m', MLPClassifier(hidden_layer_sizes=shape, max_iter=600,
                                          random_state=0))]).fit(X_tr, y_tr)
    rows.append({
        'hidden layers': str(shape),
        'approx params': parameters,
        'test AUC': round(roc_auc_score(y_te, model.predict_proba(X_te)[:, 1]), 3),
    })
print(pd.DataFrame(rows).to_string(index=False))
print()
print('(16, 16) has far fewer parameters than (256,) and does better,')
print('because the problem is a composition of simple pieces - which')
print('is what depth is good at.')

import numpy as np

rng = np.random.default_rng(seed=13)

# ---- the data: XOR, which no single linear layer can solve --------
n = 2000
X = rng.uniform(-1, 1, (n, 2))
y = ((X[:, 0] > 0) ^ (X[:, 1] > 0)).astype(float).reshape(-1, 1)

split = int(n * 0.8)
X_tr, X_te, y_tr, y_te = X[:split], X[split:], y[:split], y[split:]

# ---- the network: 2 -> 8 -> 1 -------------------------------------
HIDDEN = 8
# He initialisation: scale by sqrt(2 / fan_in), which keeps the
# activations from shrinking or exploding as depth grows.
W1 = rng.normal(0, np.sqrt(2 / 2), (2, HIDDEN))
b1 = np.zeros(HIDDEN)
W2 = rng.normal(0, np.sqrt(2 / HIDDEN), (HIDDEN, 1))
b2 = np.zeros(1)

def sigmoid(z):
    return 1 / (1 + np.exp(-np.clip(z, -60, 60)))

def forward(X, W1, b1, W2, b2):
    z1 = X @ W1 + b1
    a1 = np.maximum(0, z1)              # ReLU
    z2 = a1 @ W2 + b2
    a2 = sigmoid(z2)
    return z1, a1, z2, a2

def loss(predicted, actual):
    p = np.clip(predicted, 1e-9, 1 - 1e-9)
    return float(-np.mean(actual * np.log(p) + (1 - actual) * np.log(1 - p)))

LEARNING_RATE = 0.5
history = []

for epoch in range(1, 2001):
    z1, a1, z2, a2 = forward(X_tr, W1, b1, W2, b2)

    # ---- backpropagation: the chain rule, one layer at a time ----
    # dL/dz2 for binary cross-entropy with a sigmoid output.
    dz2 = (a2 - y_tr) / len(X_tr)
    dW2 = a1.T @ dz2
    db2 = dz2.sum(axis=0)

    da1 = dz2 @ W2.T
    dz1 = da1 * (z1 > 0)                # the ReLU derivative
    dW1 = X_tr.T @ dz1
    db1 = dz1.sum(axis=0)

    # ---- the update ---------------------------------------------
    W2 -= LEARNING_RATE * dW2
    b2 -= LEARNING_RATE * db2
    W1 -= LEARNING_RATE * dW1
    b1 -= LEARNING_RATE * db1

    if epoch % 250 == 0 or epoch == 1:
        train_loss = loss(a2, y_tr)
        test_loss = loss(forward(X_te, W1, b1, W2, b2)[3], y_te)
        accuracy = float(((forward(X_te, W1, b1, W2, b2)[3] > 0.5) == y_te).mean())
        history.append((epoch, train_loss, test_loss, accuracy))

print(f'{"epoch":>7} {"train loss":>11} {"test loss":>10} {"test acc":>9}')
for epoch, train_loss, test_loss, accuracy in history:
    print(f'{epoch:>7} {train_loss:>11.4f} {test_loss:>10.4f} {accuracy:>9.3f}')
print()

# What the hidden layer learned: each unit is a half-plane.
print('hidden layer weights (each column is one unit):')
print(np.round(W1, 2))
print()
print('Each hidden unit learned a line through the input space. The')
print('output layer combines them into the XOR pattern, which no')
print('single line could ever produce.')

import numpy as np
import pandas as pd

rng = np.random.default_rng(17)

def propagate(scale, layers=12, width=200, seed=0):
    """Push a signal through many layers and watch its size."""
    generator = np.random.default_rng(seed)
    activation = generator.normal(0, 1, (500, width))
    sizes = []
    for _ in range(layers):
        W = generator.normal(0, scale, (width, width))
        activation = np.maximum(0, activation @ W)
        sizes.append(float(activation.std()))
    return sizes

rows = []
for name, scale in [('too small (0.01)', 0.01),
                    ('too large (0.20)', 0.20),
                    ('He (sqrt(2/n))', np.sqrt(2 / 200))]:
    sizes = propagate(scale)
    rows.append({'initialisation': name,
                 'layer 1 std': round(sizes[0], 4),
                 'layer 6 std': round(sizes[5], 4),
                 'layer 12 std': round(sizes[11], 6)})
print(pd.DataFrame(rows).to_string(index=False))
print()
print('Too small and the signal vanishes by layer 12; too large and it')
print('explodes. He initialisation keeps it roughly constant, which is')
print('what makes a deep network trainable at all.')

import numpy as np
import pandas as pd

def parameters(shape):
    return sum(a * b + b for a, b in zip(shape[:-1], shape[1:]))

rows = []
for shape in [(100, 64, 1), (100, 256, 1), (100, 64, 64, 1),
              (100, 256, 256, 1), (100, 512, 512, 512, 1)]:
    count = parameters(shape)
    rows.append({'architecture': ' -> '.join(map(str, shape)),
                 'parameters': f'{count:,}',
                 'memory (float32)': f'{count * 4 / 1024:,.0f} KiB'})
print(pd.DataFrame(rows).to_string(index=False))
print()
print('Parameters grow with the PRODUCT of adjacent layer widths, so')
print('doubling the width roughly quadruples the layer. Depth is the')
print('cheaper way to add capacity.')