Vanishing Gradients — Sigmoid Freezes Neural Networks

import numpy as np

def sigmoid(raw_output):
    """Maps any real number to the open interval (0, 1).
    Values far below zero approach 0; far above zero approach 1.
    Useful for expressing confidence as a probability.
    Note: derivative maxes at 0.25 — a critical limitation for deep hidden layers."""
    return 1 / (1 + np.exp(-raw_output))

def relu(raw_output):
    """Rectified Linear Unit — the default choice for hidden layers.
    Returns the input if positive, zero otherwise.
    Derivative is 1 for positive inputs, so gradients do not shrink."""
    return np.maximum(0, raw_output)

def single_neuron_forward(inputs, weights, bias, activation='sigmoid'):
    """One complete forward pass through a single neuron.

    inputs  : numpy array of input values (features)
    weights : numpy array of learned weights, one per input
    bias    : scalar bias term
    """
    # Step 1: linear combination — the weighted vote
    weighted_sum = np.dot(inputs, weights) + bias

    # Step 2: apply the non-linear gate
    if activation == 'sigmoid':
        output = sigmoid(weighted_sum)
    elif activation == 'relu':
        output = relu(weighted_sum)
    else:
        raise ValueError(f'Unknown activation: {activation}')

    return weighted_sum, output

# --- Example: predicting whether a house is 'expensive' ---
# Features are normalised to roughly the same scale before being passed in.
# Skipping normalisation is the #1 cause of erratic training — do not skip it.
house_inputs    = np.array([0.85, 0.60])   # normalised square footage and bedroom count
initial_weights = np.array([0.40, 0.35])   # relative importance learned during training
bias_term       = -0.20                    # shifts the decision threshold

raw, output = single_neuron_forward(house_inputs, initial_weights, bias_term, 'sigmoid')

print(f"Weighted sum (before activation): {raw:.4f}")
print(f"Sigmoid output (prediction):      {output:.4f}")
print(f"Interpretation: {output*100:.1f}% confidence the house is expensive")
print()

# --- Demonstrating why activation choice matters for hidden layers ---
raw_relu, out_relu = single_neuron_forward(house_inputs, initial_weights, bias_term, 'relu')
print(f"Same neuron with ReLU: raw={raw_relu:.4f}, output={out_relu:.4f}")
print("ReLU output is not squished to (0,1) — it preserves scale in hidden layers,")
print("which keeps gradients alive during backpropagation through deep networks.")

import numpy as np

np.random.seed(42)  # Reproducibility is non-negotiable for debugging

# ── Activation functions ──────────────────────────────────────────────────────

def sigmoid(z):
    """Maps any real number to (0, 1). Use ONLY at the output layer for binary tasks."""
    return 1 / (1 + np.exp(-z))

def sigmoid_derivative(z):
    """Gradient of sigmoid with respect to its input z.
    Maximum value is 0.25 — this is the root cause of vanishing gradients in deep networks."""
    s = sigmoid(z)
    return s * (1 - s)

def relu(z):
    """Rectified Linear Unit. Default choice for hidden layers.
    Derivative is 1 for positive inputs — gradients do not shrink through this function."""
    return np.maximum(0, z)

def relu_derivative(z):
    """Gradient of ReLU. Dead neurons (z <= 0) have zero gradient and stop learning.
    He initialisation and small learning rates help keep most neurons alive."""
    return (z > 0).astype(float)

# ── Network initialisation ────────────────────────────────────────────────────

def initialise_network(layer_sizes):
    """Creates weight matrices and bias vectors for a network of arbitrary depth.

    layer_sizes: e.g. [2, 4, 4, 1] means 2 inputs -> 4 neurons -> 4 neurons -> 1 output.

    He initialisation scales weights by sqrt(2/fan_in).
    This keeps activation variance stable through ReLU layers so gradients
    do not vanish or explode before training has a chance to do anything useful.
    Xavier initialisation (sqrt(1/fan_in)) is the alternative for sigmoid/tanh.
    """
    parameters = {}
    for layer_idx in range(1, len(layer_sizes)):
        fan_in  = layer_sizes[layer_idx - 1]
        fan_out = layer_sizes[layer_idx]
        parameters[f'W{layer_idx}'] = np.random.randn(fan_out, fan_in) * np.sqrt(2 / fan_in)
        parameters[f'b{layer_idx}'] = np.zeros((fan_out, 1))
    return parameters

# ── Forward propagation ───────────────────────────────────────────────────────

def forward_pass(input_data, parameters, num_layers):
    """Passes input through every layer in sequence.
    Caches both Z (pre-activation) and A (post-activation) at every layer.
    These cached values are required by backpropagation — do not discard them.
    """
    cache     = {'A0': input_data}
    current_A = input_data

    for idx in range(1, num_layers + 1):
        W = parameters[f'W{idx}']
        b = parameters[f'b{idx}']
        Z = np.dot(W, current_A) + b
        cache[f'Z{idx}'] = Z

        # Hidden layers use ReLU; output layer uses sigmoid for [0,1] probability
        current_A = relu(Z) if idx < num_layers else sigmoid(Z)
        cache[f'A{idx}'] = current_A

    return current_A, cache

# ── Loss function ─────────────────────────────────────────────────────────────

def binary_cross_entropy_loss(predictions, true_labels):
    """Cross-entropy penalises confident wrong answers very heavily,
    which is why it converges faster than MSE for classification tasks.
    The epsilon prevents log(0) from crashing training."""
    epsilon     = 1e-15
    predictions = np.clip(predictions, epsilon, 1 - epsilon)
    return -np.mean(
        true_labels * np.log(predictions) +
        (1 - true_labels) * np.log(1 - predictions)
    )

# ── Backpropagation ───────────────────────────────────────────────────────────

def backward_pass(predictions, true_labels, cache, parameters, num_layers):
    """Assigns blame to every weight by traversing the network in reverse.
    Uses the cached Z and A values from the forward pass to compute each gradient.
    """
    gradients   = {}
    num_samples = true_labels.shape[1]

    # Combined derivative of binary cross-entropy loss and sigmoid output activation
    dA_current = -(true_labels / predictions) + (1 - true_labels) / (1 - predictions)

    for idx in reversed(range(1, num_layers + 1)):
        Z      = cache[f'Z{idx}']
        A_prev = cache[f'A{idx - 1}']
        W      = parameters[f'W{idx}']

        dZ             = dA_current * (sigmoid_derivative(Z) if idx == num_layers else relu_derivative(Z))
        gradients[f'dW{idx}'] = np.dot(dZ, A_prev.T) / num_samples
        gradients[f'db{idx}'] = np.sum(dZ, axis=1, keepdims=True) / num_samples
        dA_current            = np.dot(W.T, dZ)  # propagate gradient to previous layer

    return gradients

# ── Gradient descent update ───────────────────────────────────────────────────

def update_weights(parameters, gradients, learning_rate, num_layers):
    """Nudges every weight in the direction that reduces the loss.
    The learning rate controls step size — too large overshoots, too small crawls.
    """
    for idx in range(1, num_layers + 1):
        parameters[f'W{idx}'] -= learning_rate * gradients[f'dW{idx}']
        parameters[f'b{idx}'] -= learning_rate * gradients[f'db{idx}']
    return parameters

# ── Full training loop ────────────────────────────────────────────────────────

def train(input_data, true_labels, layer_sizes, learning_rate=0.01, num_epochs=1000):
    num_layers = len(layer_sizes) - 1
    parameters = initialise_network(layer_sizes)

    for epoch in range(num_epochs):
        predictions, cache = forward_pass(input_data, parameters, num_layers)
        loss               = binary_cross_entropy_loss(predictions, true_labels)
        gradients          = backward_pass(predictions, true_labels, cache, parameters, num_layers)
        parameters         = update_weights(parameters, gradients, learning_rate, num_layers)

        if epoch % 200 == 0:
            print(f'Epoch {epoch:>5} | Loss: {loss:.4f}')

    return parameters

# ── XOR: the classic proof that depth enables non-linear learning ──────────────
# XOR (exclusive OR) cannot be separated by a single straight line.
# No single-layer model can solve it. A two-hidden-layer network solves it reliably.
# If your implementation solves XOR cleanly, backpropagation and layer structure work.

XOR_inputs  = np.array([[0, 0, 1, 1], [0, 1, 0, 1]])  # 2 features, 4 samples
XOR_outputs = np.array([[0, 1, 1, 0]])                 # 1 output, 4 samples

print('=== Training on XOR Problem ===')
trained_params = train(
    input_data    = XOR_inputs,
    true_labels   = XOR_outputs,
    layer_sizes   = [2, 4, 4, 1],
    learning_rate = 0.1,
    num_epochs    = 1001
)

final_preds, _ = forward_pass(XOR_inputs, trained_params, len([2, 4, 4, 1]) - 1)
print('\n=== Final Predictions ===')
for i in range(4):
    inp  = XOR_inputs[:, i]
    pred = final_preds[0, i]
    exp  = XOR_outputs[0, i]
    print(f'Input: {inp} | Expected: {exp} | Predicted: {pred:.4f}')

import numpy as np

# ── Gradient descent on a simple 1D loss surface ──────────────────────────────
# In a real network, the loss surface is millions of dimensions.
# The same principle applies: move opposite to the gradient, scaled by the learning rate.
# This 1D version lets you watch exactly what is happening without the complexity.

def simplified_loss(weight_value):
    """A bowl-shaped loss surface. The minimum is at weight_value = 2.0.
    Chosen because the gradient has a simple closed form — easy to verify by hand."""
    return (weight_value - 2.0) ** 2 + 0.5

def loss_gradient(weight_value):
    """The derivative of the loss with respect to the weight.
    For (w-2)^2 + 0.5, this is 2*(w-2).
    Positive gradient => increasing w increases loss => decrease w.
    Negative gradient => increasing w decreases loss => increase w."""
    return 2 * (weight_value - 2.0)

current_weight = -3.0   # Start far from optimal
learning_rate  = 0.15   # How large each step is
num_steps      = 20     # How many updates to apply

print('=== Gradient Descent: Watching a Weight Find Its Optimal Value ===')
print(f'{"Step":>5} | {"Weight":>10} | {"Loss":>10} | {"Gradient":>10}')
print('-' * 47)

for step in range(num_steps):
    current_loss = simplified_loss(current_weight)
    gradient     = loss_gradient(current_weight)

    # Core update rule: w = w - learning_rate * gradient
    # Moving OPPOSITE to the gradient reduces the loss.
    # This is gradient descent — the same operation applied to millions of weights in parallel.
    current_weight -= learning_rate * gradient

    if step % 4 == 0 or step == num_steps - 1:
        print(f'{step:>5} | {current_weight:>10.4f} | {current_loss:>10.4f} | {gradient:>10.4f}')

print(f'\nFinal weight: {current_weight:.4f}  (target: 2.0)')
print(f'Final loss:   {simplified_loss(current_weight):.6f}  (minimum: 0.5)')
print()
print('Notice: the gradient shrinks as the weight approaches the minimum.')
print('This is why learning rate decay is useful — steps should be smaller')
print('as the gradient becomes smaller and we get closer to the optimum.')