Activation Functions

Stack two linear transformations: y = W2(W1x + b1) + b2. Expand it and you get y = W2W1x + W2b1 + b2, which is just y = Ax + c: a single linear transformation. No matter how many linear layers you stack, the result collapses to one matrix multiply. Your 100-layer network has the exact same representational power as a single layer.

That's not a theoretical curiosity. A deep linear network literally cannot learn XOR, cannot classify a spiral dataset, cannot recognize a face. Without activation functions, depth is an illusion.

Activation functions break the linearity. They warp the output of each layer through a nonlinear function, giving the network room to bend decision boundaries, approximate arbitrary functions, and actually learn. Pick the wrong one, though, and something breaks: gradients vanish to zero (sigmoid in deep networks), explode to infinity (unbounded activations without careful initialization), or neurons die permanently (ReLU with large negative biases). The activation function you choose determines whether your network learns at all.

Why Nonlinearity Is Necessary

Matrix multiplication is composable. Multiplying a vector by matrix A then matrix B is identical to multiplying by AB. Stack ten linear layers and the result is mathematically equivalent to one linear layer with one big matrix. All those parameters, all that depth. Wasted. You need something to break the chain. That's what activation functions do.

Here's the proof. A linear layer computes f(x) = Wx + b. Stack two:

Layer 1: h = W1 * x + b1
Layer 2: y = W2 * h + b2

Substitute:

y = W2 * (W1 * x + b1) + b2
y = (W2 * W1) * x + (W2 * b1 + b2)
y = A * x + c

One layer. Insert a nonlinear activation g() between layers:

h = g(W1 * x + b1)
y = W2 * h + b2

Now the substitution breaks. W2 * g(W1 * x + b1) + b2 cannot be reduced to a single linear transformation. The network can represent nonlinear functions. Each additional layer with an activation adds representational capacity.

Sigmoid

Sigmoid squeezes any real number into a smooth, differentiable S-curve: a probability-like output, and the original activation function for neural networks.

sigmoid(x) = 1 / (1 + e^(-x))

Output range: (0, 1). The derivative:

sigmoid'(x) = sigmoid(x) * (1 - sigmoid(x))

The maximum value of this derivative is 0.25, occurring at x = 0. In backpropagation, gradients multiply through layers. Ten layers of sigmoid means the gradient gets multiplied by at most 0.25 ten times:

0.25^10 = 0.000000953674

Less than one millionth of the original signal. This is the vanishing gradient problem. Gradients in early layers become so small that weights barely update. The network appears to learn (loss decreases in later layers), but the first layers stay frozen. Deep sigmoid networks simply do not train.

One more problem: sigmoid outputs are always positive (0 to 1), which means gradients on weights are always the same sign. This causes zig-zagging during gradient descent.

Tanh

Tanh is sigmoid recentered: the same S-curve, but zero-centered, which eliminates the zig-zag problem.

tanh(x) = (e^x - e^(-x)) / (e^x + e^(-x))

Output range: (-1, 1). The derivative:

tanh'(x) = 1 - tanh(x)^2

Maximum derivative is 1.0 at x = 0: four times better than sigmoid. But the vanishing gradient problem still exists. For large positive or negative inputs, the derivative approaches zero. Ten layers still crush the gradient, just less aggressively.

ReLU: The Breakthrough

ReLU (Rectified Linear Unit) zeroes out every negative input and passes positive inputs straight through, unchanged:

relu(x) = max(0, x)

Nair and Hinton popularized it for deep learning in 2010, though the function itself dates to Fukushima's 1969 work.

Output range: [0, infinity). The derivative is trivially simple:

relu'(x) = 1  if x > 0
            0  if x <= 0

No vanishing gradient for positive inputs: the gradient is exactly 1, passed straight through. This is why deep networks became trainable: ReLU preserves gradient magnitude across layers.

But there is a failure mode: the dead neuron problem. If a neuron's weighted input is always negative (due to a large negative bias or unfortunate weight initialization), its output is always zero, its gradient is always zero, and it never updates. It is permanently dead. In practice, 10-40% of neurons in a ReLU network can die during training.

ReLU's output is also unbounded on the positive side, which cuts the other way from vanishing gradients. With careless weight initialization, activations and their gradients can grow layer over layer instead of shrinking, exploding toward NaN. Vanishing and exploding are the same compounding-multiplication mechanism, just pointed in opposite directions: a per-layer multiplier under 1 shrinks a signal to nothing over enough layers, and a multiplier over 1 blows it up.

Leaky ReLU

Leaky ReLU gives the negative side a small slope instead of zero: the simplest fix for dead neurons.

leaky_relu(x) = x        if x > 0
                alpha * x if x <= 0

Alpha is a small constant, typically 0.01. That small slope means dead neurons still get a gradient signal and can recover.

GELU: The Modern Default

GELU (Gaussian Error Linear Unit) weights each input by how likely it is to be positive under a Gaussian distribution: a probabilistic gate instead of a hard cutoff. Hendrycks and Gimpel introduced it in 2016, and it's now the default activation in BERT, GPT, and most modern transformers.

gelu(x) = x * Phi(x)

Where Phi(x) is the cumulative distribution function of the standard normal distribution. The approximation used in practice:

gelu(x) ~= 0.5 * x * (1 + tanh(sqrt(2/pi) * (x + 0.044715 * x^3)))

GELU is smooth everywhere, allows small negative values (unlike ReLU, which hard-clips to zero), and this smooth gating outperforms ReLU in transformer architectures because it provides better gradient flow and avoids the dead neuron problem entirely.

Swish / SiLU

Swish gates each input by its own sigmoid: a self-gated activation that Ramachandran et al. found in 2017 through automated search.

swish(x) = x * sigmoid(x)

Google ran that search over activation function space: a neural network designing parts of neural networks.

Like GELU, it is smooth, non-monotonic, and allows small negative values. The difference is subtle: Swish uses sigmoid for gating while GELU uses the Gaussian CDF. In practice, performance is nearly identical. Swish is used in EfficientNet and some vision models. GELU dominates in language models.

Softmax: The Output Activation

Softmax converts a vector of raw scores (logits) into a probability distribution: every output lands between 0 and 1, and they all sum to 1. It's never used in hidden layers, only as the standard final activation for multi-class classification.

softmax(x_i) = e^(x_i) / sum(e^(x_j) for all j)

The largest logit gets the highest probability, but unlike argmax, softmax is differentiable and preserves information about relative confidence.

Comparison of Shapes

Gradient Flow Comparison

Which Activation When

Step 1: Implement All Activation Functions with Derivatives

Each function takes a single float and returns a float. Each derivative function takes the same input and returns the gradient.

import math
 
def sigmoid(x):
    x = max(-500, min(500, x))
    return 1.0 / (1.0 + math.exp(-x))
 
def sigmoid_derivative(x):
    s = sigmoid(x)
    return s * (1 - s)
 
def tanh_act(x):
    return math.tanh(x)
 
def tanh_derivative(x):
    t = math.tanh(x)
    return 1 - t * t
 
def relu(x):
    return max(0.0, x)
 
def relu_derivative(x):
    return 1.0 if x > 0 else 0.0
 
def leaky_relu(x, alpha=0.01):
    return x if x > 0 else alpha * x
 
def leaky_relu_derivative(x, alpha=0.01):
    return 1.0 if x > 0 else alpha
 
def gelu(x):
    return 0.5 * x * (1 + math.tanh(math.sqrt(2 / math.pi) * (x + 0.044715 * x ** 3)))
 
def gelu_derivative(x):
    phi = 0.5 * (1 + math.erf(x / math.sqrt(2)))
    pdf = math.exp(-0.5 * x * x) / math.sqrt(2 * math.pi)
    return phi + x * pdf
 
def swish(x):
    return x * sigmoid(x)
 
def swish_derivative(x):
    s = sigmoid(x)
    return s + x * s * (1 - s)
 
def softmax(xs):
    max_x = max(xs)
    exps = [math.exp(x - max_x) for x in xs]
    total = sum(exps)
    return [e / total for e in exps]

Step 2: Visualize Where Gradients Die

Compute the gradient at 100 evenly-spaced points from -5 to 5. Print a text histogram showing where each activation's gradient is near-zero.

def gradient_scan(name, derivative_fn, start=-5, end=5, n=100):
    step = (end - start) / n
    near_zero = 0
    healthy = 0
    for i in range(n):
        x = start + i * step
        g = derivative_fn(x)
        if abs(g) < 0.01:
            near_zero += 1
        else:
            healthy += 1
    pct_dead = near_zero / n * 100
    print(f"{name:15s}: {healthy:3d} healthy, {near_zero:3d} near-zero ({pct_dead:.0f}% dead zone)")
 
gradient_scan("Sigmoid", sigmoid_derivative)
gradient_scan("Tanh", tanh_derivative)
gradient_scan("ReLU", relu_derivative)
gradient_scan("Leaky ReLU", leaky_relu_derivative)
gradient_scan("GELU", gelu_derivative)
gradient_scan("Swish", swish_derivative)

Step 3: Vanishing Gradient Experiment

Forward-pass a signal through N layers using sigmoid, ReLU, and GELU. Measure how the activation magnitude changes.

import random
 
def vanishing_gradient_experiment(activation_fn, name, n_layers=10, n_inputs=5):
    random.seed(42)
    values = [random.gauss(0, 1) for _ in range(n_inputs)]
 
    print(f"\n{name} through {n_layers} layers:")
    for layer in range(n_layers):
        weights = [random.gauss(0, 1) for _ in range(n_inputs)]
        z = sum(w * v for w, v in zip(weights, values))
        activated = activation_fn(z)
        magnitude = abs(activated)
        bar = "#" * int(magnitude * 20)
        print(f"  Layer {layer+1:2d}: magnitude = {magnitude:.6f} {bar}")
        values = [activated] * n_inputs
 
vanishing_gradient_experiment(sigmoid, "Sigmoid")
vanishing_gradient_experiment(relu, "ReLU")
vanishing_gradient_experiment(gelu, "GELU")

Step 4: Dead Neuron Detector

Create a ReLU network, pass random inputs through it, count how many neurons never fire.

def dead_neuron_detector(n_inputs=5, hidden_size=20, n_samples=1000):
    random.seed(0)
    weights = [[random.gauss(0, 1) for _ in range(n_inputs)] for _ in range(hidden_size)]
    biases = [random.gauss(0, 1) for _ in range(hidden_size)]
 
    fire_counts = [0] * hidden_size
 
    for _ in range(n_samples):
        inputs = [random.gauss(0, 1) for _ in range(n_inputs)]
        for neuron_idx in range(hidden_size):
            z = sum(w * x for w, x in zip(weights[neuron_idx], inputs)) + biases[neuron_idx]
            if relu(z) > 0:
                fire_counts[neuron_idx] += 1
 
    dead = sum(1 for c in fire_counts if c == 0)
    rarely_fire = sum(1 for c in fire_counts if 0 < c < n_samples * 0.05)
    healthy = hidden_size - dead - rarely_fire
 
    print(f"\nDead Neuron Report ({hidden_size} neurons, {n_samples} samples):")
    print(f"  Dead (never fired):     {dead}")
    print(f"  Barely alive (<5%):     {rarely_fire}")
    print(f"  Healthy:                {healthy}")
    print(f"  Dead neuron rate:       {dead/hidden_size*100:.1f}%")
 
    for i, c in enumerate(fire_counts):
        status = "DEAD" if c == 0 else "WEAK" if c < n_samples * 0.05 else "OK"
        bar = "#" * (c * 40 // n_samples)
        print(f"  Neuron {i:2d}: {c:4d}/{n_samples} fires [{status:4s}] {bar}")
 
dead_neuron_detector()

Step 5: Comparing Sigmoid, ReLU, and GELU

Train the same two-layer network on the circle dataset (points inside a circle = class 1, outside = class 0) with three different activations. Compare convergence speed.

def make_circle_data(n=200, seed=42):
    random.seed(seed)
    data = []
    for _ in range(n):
        x = random.uniform(-2, 2)
        y = random.uniform(-2, 2)
        label = 1.0 if x * x + y * y < 1.5 else 0.0
        data.append(([x, y], label))
    return data
 
 
class ActivationNetwork:
    def __init__(self, activation_fn, activation_deriv, hidden_size=8, lr=0.1):
        random.seed(0)
        self.act = activation_fn
        self.act_d = activation_deriv
        self.lr = lr
        self.hidden_size = hidden_size
 
        self.w1 = [[random.gauss(0, 0.5) for _ in range(2)] for _ in range(hidden_size)]
        self.b1 = [0.0] * hidden_size
        self.w2 = [random.gauss(0, 0.5) for _ in range(hidden_size)]
        self.b2 = 0.0
 
    def forward(self, x):
        self.x = x
        self.z1 = []
        self.h = []
        for i in range(self.hidden_size):
            z = self.w1[i][0] * x[0] + self.w1[i][1] * x[1] + self.b1[i]
            self.z1.append(z)
            self.h.append(self.act(z))
 
        self.z2 = sum(self.w2[i] * self.h[i] for i in range(self.hidden_size)) + self.b2
        self.out = sigmoid(self.z2)
        return self.out
 
    def backward(self, target):
        error = self.out - target
        d_out = error * self.out * (1 - self.out)
 
        for i in range(self.hidden_size):
            d_h = d_out * self.w2[i] * self.act_d(self.z1[i])
            self.w2[i] -= self.lr * d_out * self.h[i]
            for j in range(2):
                self.w1[i][j] -= self.lr * d_h * self.x[j]
            self.b1[i] -= self.lr * d_h
        self.b2 -= self.lr * d_out
 
    def train(self, data, epochs=200):
        losses = []
        for epoch in range(epochs):
            total_loss = 0
            correct = 0
            for x, y in data:
                pred = self.forward(x)
                self.backward(y)
                total_loss += (pred - y) ** 2
                if (pred >= 0.5) == (y >= 0.5):
                    correct += 1
            avg_loss = total_loss / len(data)
            accuracy = correct / len(data) * 100
            losses.append(avg_loss)
            if epoch % 50 == 0 or epoch == epochs - 1:
                print(f"    Epoch {epoch:3d}: loss={avg_loss:.4f}, accuracy={accuracy:.1f}%")
        return losses
 
 
data = make_circle_data()
 
configs = [
    ("Sigmoid", sigmoid, sigmoid_derivative),
    ("ReLU", relu, relu_derivative),
    ("GELU", gelu, gelu_derivative),
]
 
results = {}
for name, act_fn, act_d_fn in configs:
    print(f"\n=== Training with {name} ===")
    net = ActivationNetwork(act_fn, act_d_fn, hidden_size=8, lr=0.1)
    losses = net.train(data, epochs=200)
    results[name] = losses
 
print("\n=== Final Loss Comparison ===")
for name, losses in results.items():
    print(f"  {name:10s}: start={losses[0]:.4f} -> end={losses[-1]:.4f} (improvement: {(1 - losses[-1]/losses[0])*100:.1f}%)")

PyTorch provides all of these as both functional and module forms:

import torch
import torch.nn as nn
import torch.nn.functional as F
 
x = torch.randn(4, 10)
 
relu_out = F.relu(x)
gelu_out = F.gelu(x)
sigmoid_out = torch.sigmoid(x)
swish_out = F.silu(x)
 
logits = torch.randn(4, 5)
probs = F.softmax(logits, dim=1)
 
model = nn.Sequential(
    nn.Linear(10, 64),
    nn.GELU(),
    nn.Linear(64, 32),
    nn.GELU(),
    nn.Linear(32, 5),
)

Hidden layers in a transformer: GELU. Hidden layers in a CNN: ReLU. Output layer for classification: softmax. Output layer for regression: none (linear). Output layer for probabilities: sigmoid. That's it. Start with these defaults. Change them only when you have evidence.

RNNs and LSTMs use tanh for hidden state and sigmoid for gates, but if you're building from scratch today, you're probably not using RNNs. If neurons are dying in your ReLU network, switch to GELU. Don't reach for Leaky ReLU unless you have a specific reason: GELU solves the dead neuron problem and gives better gradient flow.

This lesson produces:

• outputs/prompt-activation-selector.md, a reusable prompt that helps you pick the right activation function for any architecture

1. Implement Parametric ReLU (PReLU) where the negative slope alpha is a learnable parameter. Train it on the circle dataset and compare to fixed Leaky ReLU.

2. Run the vanishing gradient experiment with 50 layers instead of 10. Plot the magnitude at each layer for sigmoid, tanh, ReLU, and GELU. At which layer does each activation's signal effectively reach zero?

3. Implement the ELU (Exponential Linear Unit): elu(x) = x if x > 0, alpha * (e^x - 1) if x <= 0. Compare its dead neuron rate to ReLU on the same network.

4. Build a "gradient health monitor" that runs during training: at each epoch, compute the average gradient magnitude at each layer. Print a warning when any layer's gradient drops below 0.001 or exceeds 100.

5. Modify the training comparison to use the XOR dataset from Lesson 01 instead of circles. Which activation converges fastest on XOR? Why does this differ from the circle results?

• Nair & Hinton, "Rectified Linear Units Improve Restricted Boltzmann Machines" (2010) — the paper that introduced ReLU and enabled training of deep networks
• Hendrycks & Gimpel, "Gaussian Error Linear Units (GELUs)" (2016) — introduced the activation function that became the default for transformers
• Ramachandran et al., "Searching for Activation Functions" (2017) — used automated search to discover Swish, showing that activation design can be automated
• Glorot & Bengio, "Understanding the difficulty of training deep feedforward neural networks" (2010) — the paper that diagnosed vanishing/exploding gradients and proposed Xavier initialization
• Goodfellow, Bengio, Courville, "Deep Learning" Chapter 6.3 (https://www.deeplearningbook.org/) — rigorous treatment of hidden units and activation functions