Multi-Layer Perceptrons

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188

import numpy as np
from typing import List, Tuple, Callable, Dict
from dataclasses import dataclass
 
@dataclass
class ForwardCache:
    """
    Cache storing all intermediate computations during forward pass.
    Required for backpropagation during training.
    """
    inputs: np.ndarray          # Original input x
    pre_activations: List[np.ndarray]  # z^(l) for each layer
    activations: List[np.ndarray]      # a^(l) for each layer, a^(0) = x
 
 
def sigmoid(z: np.ndarray) -> np.ndarray:
    """Numerically stable sigmoid activation."""
    return np.where(
        z >= 0,
        1 / (1 + np.exp(-z)),
        np.exp(z) / (1 + np.exp(z))
    )
 
 
def relu(z: np.ndarray) -> np.ndarray:
    """ReLU activation: max(0, z)."""
    return np.maximum(0, z)
 
 
def softmax(z: np.ndarray) -> np.ndarray:
    """Numerically stable softmax activation."""
    z_shifted = z - np.max(z, axis=-1, keepdims=True)
    exp_z = np.exp(z_shifted)
    return exp_z / np.sum(exp_z, axis=-1, keepdims=True)
 
 
def forward_propagation(
    x: np.ndarray,
    weights: List[np.ndarray],
    biases: List[np.ndarray],
    activations: List[Callable],
    return_cache: bool = True
) -> Tuple[np.ndarray, ForwardCache]:
    """
    Complete forward propagation through an MLP.
    
    This implementation is designed for clarity and correctness,
    demonstrating the exact computation at each step.
    
    Args:
        x: Input array of shape (n_features,) or (batch_size, n_features)
        weights: List of weight matrices W^(l), each shape (n_out, n_in)
        biases: List of bias vectors b^(l), each shape (n_out,)
        activations: List of activation functions σ^(l)
        return_cache: Whether to return intermediate values for backprop
    
    Returns:
        output: Network output of shape (n_output,) or (batch_size, n_output)
        cache: ForwardCache with all intermediate values
    """
    # Ensure input is 2D: (batch_size, n_features)
    if x.ndim == 1:
        x = x.reshape(1, -1)
    
    L = len(weights)  # Number of layers
    
    # Initialize cache
    pre_activations = []
    layer_activations = [x]  # a^(0) = x
    
    a = x  # Current activation
    
    for l in range(L):
        # Pre-activation: z^(l) = a^(l-1) @ W^(l).T + b^(l)
        # Note: We use (batch, features) convention, so we transpose W
        z = a @ weights[l].T + biases[l]
        pre_activations.append(z)
        
        # Activation: a^(l) = σ^(l)(z^(l))
        a = activations[l](z)
        layer_activations.append(a)
    
    cache = ForwardCache(
        inputs=x,
        pre_activations=pre_activations,
        activations=layer_activations
    ) if return_cache else None
    
    return a, cache
 
 
def trace_forward_pass(
    x: np.ndarray,
    weights: List[np.ndarray],
    biases: List[np.ndarray],
    activations: List[Callable],
    activation_names: List[str]
) -> Dict:
    """
    Trace forward pass with detailed output at each step.
    Useful for debugging and understanding information flow.
    """
    print("=" * 60)
    print("FORWARD PROPAGATION TRACE")
    print("=" * 60)
    
    if x.ndim == 1:
        x = x.reshape(1, -1)
    
    print(f"
Input x: shape {x.shape}")
    print(f"  values: {x[0, :5]}..." if x.shape[1] > 5 else f"  values: {x[0]}")
    
    a = x
    trace = {"input": x, "layers": []}
    
    for l in range(len(weights)):
        print(f"
--- Layer {l + 1} ({activation_names[l]}) ---")
        print(f"W^({l+1}) shape: {weights[l].shape}")
        print(f"b^({l+1}) shape: {biases[l].shape}")
        
        # Compute pre-activation
        z = a @ weights[l].T + biases[l]
        print(f"z^({l+1}) = a^({l}) @ W^({l+1}).T + b^({l+1})")
        print(f"  z shape: {z.shape}")
        print(f"  z range: [{z.min():.4f}, {z.max():.4f}]")
        print(f"  z mean: {z.mean():.4f}, std: {z.std():.4f}")
        
        # Compute activation
        a = activations[l](z)
        print(f"a^({l+1}) = {activation_names[l]}(z^({l+1}))")
        print(f"  a shape: {a.shape}")
        print(f"  a range: [{a.min():.4f}, {a.max():.4f}]")
        
        if activation_names[l] == "ReLU":
            sparsity = np.mean(a == 0)
            print(f"  sparsity (fraction of zeros): {sparsity:.2%}")
        
        trace["layers"].append({
            "z": z.copy(),
            "a": a.copy(),
            "W_shape": weights[l].shape,
            "activation": activation_names[l]
        })
    
    print(f"
{'=' * 60}")
    print(f"OUTPUT: shape {a.shape}")
    print(f"  values: {a[0]}")
    print("=" * 60)
    
    trace["output"] = a
    return trace
 
 
# Demonstration
if __name__ == "__main__":
    np.random.seed(42)
    
    # Create a simple network: 4 → 8 → 6 → 3
    architecture = [4, 8, 6, 3]
    
    weights = []
    biases = []
    
    for i in range(len(architecture) - 1):
        n_in = architecture[i]
        n_out = architecture[i + 1]
        # He initialization
        W = np.random.randn(n_out, n_in) * np.sqrt(2.0 / n_in)
        b = np.zeros(n_out)
        weights.append(W)
        biases.append(b)
    
    # Activation functions: ReLU for hidden, softmax for output
    activations = [relu, relu, softmax]
    activation_names = ["ReLU", "ReLU", "Softmax"]
    
    # Single sample
    x = np.random.randn(4)
    
    # Trace the forward pass
    trace = trace_forward_pass(x, weights, biases, activations, activation_names)
    
    # Verify output sums to 1 (softmax property)
    print(f"
Output sum (should be 1.0): {trace['output'].sum():.6f}")

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104

import numpy as np
 
def sigmoid_naive(z):
    """Naive sigmoid - can overflow for negative z."""
    return 1 / (1 + np.exp(-z))
 
def sigmoid_stable(z):
    """Numerically stable sigmoid."""
    return np.where(
        z >= 0,
        1 / (1 + np.exp(-z)),
        np.exp(z) / (1 + np.exp(z))
    )
 
def softmax_naive(z):
    """Naive softmax - overflows for large z."""
    exp_z = np.exp(z)
    return exp_z / np.sum(exp_z, axis=-1, keepdims=True)
 
def softmax_stable(z):
    """Numerically stable softmax."""
    z_shifted = z - np.max(z, axis=-1, keepdims=True)
    exp_z = np.exp(z_shifted)
    return exp_z / np.sum(exp_z, axis=-1, keepdims=True)
 
def log_softmax_stable(z):
    """
    Numerically stable log-softmax.
    More accurate than log(softmax(z)) for cross-entropy.
    """
    max_z = np.max(z, axis=-1, keepdims=True)
    log_sum_exp = max_z + np.log(np.sum(np.exp(z - max_z), axis=-1, keepdims=True))
    return z - log_sum_exp
 
def cross_entropy_stable(z, targets):
    """
    Stable cross-entropy from logits.
    
    Args:
        z: Logits (pre-softmax), shape (batch, classes)
        targets: Class indices, shape (batch,)
    
    Returns:
        Cross-entropy loss per sample
    """
    log_probs = log_softmax_stable(z)
    # Select log probability of correct class
    return -log_probs[np.arange(len(targets)), targets]
 
 
# Demonstrate the difference
def demonstrate_numerical_stability():
    print("Numerical Stability Demonstration")
    print("=" * 50)
    
    # Test sigmoid with extreme values
    z_extreme = np.array([-1000, -100, 0, 100, 1000])
    
    print("
Sigmoid with extreme values:")
    print(f"Input z: {z_extreme}")
    
    # Naive (will have issues)
    with np.errstate(over='ignore', invalid='ignore'):
        naive_result = sigmoid_naive(z_extreme)
    print(f"Naive sigmoid: {naive_result}")
    
    # Stable
    stable_result = sigmoid_stable(z_extreme)
    print(f"Stable sigmoid: {stable_result}")
    
    # Test softmax with large values
    print("
Softmax with large values:")
    z_large = np.array([1000, 1001, 1002])
    print(f"Input z: {z_large}")
    
    with np.errstate(over='ignore', invalid='ignore'):
        naive_softmax = softmax_naive(z_large)
    print(f"Naive softmax: {naive_softmax}")
    
    stable_softmax = softmax_stable(z_large)
    print(f"Stable softmax: {stable_softmax}")
    
    # Log-softmax precision
    print("
Log-softmax precision:")
    z = np.array([[0.1, 0.2, 100.0]])  # Very confident prediction
    
    # Naive: log(softmax(z))
    with np.errstate(divide='ignore'):
        naive_log = np.log(softmax_stable(z))
    print(f"Naive log(softmax): {naive_log[0]}")
    
    # Stable log-softmax
    stable_log = log_softmax_stable(z)
    print(f"Stable log_softmax: {stable_log[0]}")
    
    print("
(Stable version preserves precision in the low-probability values)")
 
 
if __name__ == "__main__":
    demonstrate_numerical_stability()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123

import numpy as np
import time
from typing import List, Tuple
 
def batch_forward(
    X: np.ndarray,  # Shape: (batch_size, n_input)
    weights: List[np.ndarray],
    biases: List[np.ndarray],
    activations: List[callable]
) -> Tuple[np.ndarray, List[np.ndarray]]:
    """
    Efficient batch forward propagation.
    
    Uses matrix-matrix multiplication for all samples simultaneously.
    
    Args:
        X: Input batch, shape (batch_size, n_features)
        weights: List of weight matrices
        biases: List of bias vectors  
        activations: List of activation functions
    
    Returns:
        output: Predictions, shape (batch_size, n_output)
        all_activations: List of activation matrices for each layer
    """
    A = X  # Current activation matrix
    all_activations = [X]
    
    for W, b, sigma in zip(weights, biases, activations):
        # Z = A @ W.T + b
        # With broadcasting, b of shape (n_out,) is added to each row
        Z = A @ W.T + b
        A = sigma(Z)
        all_activations.append(A)
    
    return A, all_activations
 
 
def benchmark_batch_sizes():
    """
    Compare computation time for different batch sizes.
    Demonstrates the efficiency gains from batching.
    """
    np.random.seed(42)
    
    # Network: 784 → 512 → 256 → 128 → 64 → 10
    sizes = [784, 512, 256, 128, 64, 10]
    weights = [np.random.randn(sizes[i+1], sizes[i]) * 0.01 
               for i in range(len(sizes)-1)]
    biases = [np.zeros(sizes[i+1]) for i in range(len(sizes)-1)]
    activations = [lambda z: np.maximum(0, z)] * 4 + [lambda z: z]  # ReLU + Linear
    
    print("Batch Size Benchmark")
    print("=" * 60)
    print(f"Network: {' → '.join(map(str, sizes))}")
    print("-" * 60)
    
    # Test different batch sizes
    n_samples = 10000
    batch_sizes = [1, 8, 32, 128, 512, 2048, n_samples]
    
    for batch_size in batch_sizes:
        n_batches = n_samples // batch_size
        
        # Generate all data
        X_all = np.random.randn(n_samples, 784).astype(np.float32)
        
        # Time the processing
        start = time.time()
        
        for i in range(n_batches):
            X_batch = X_all[i*batch_size:(i+1)*batch_size]
            output, _ = batch_forward(X_batch, weights, biases, activations)
        
        elapsed = time.time() - start
        samples_per_sec = n_samples / elapsed
        
        print(f"Batch size {batch_size:5d}: {elapsed:.3f}s | "
              f"{samples_per_sec:,.0f} samples/sec | "
              f"{n_batches} batches")
    
    print("-" * 60)
    print("Larger batches = better hardware utilization = higher throughput")
 
 
def verify_batch_equivalence():
    """
    Verify that batch processing gives same results as processing
    samples individually.
    """
    np.random.seed(42)
    
    # Simple network
    sizes = [10, 20, 5]
    weights = [np.random.randn(sizes[i+1], sizes[i]) * 0.1 
               for i in range(len(sizes)-1)]
    biases = [np.zeros(sizes[i+1]) for i in range(len(sizes)-1)]
    activations = [lambda z: np.maximum(0, z), lambda z: z]
    
    # Batch of 5 samples
    X = np.random.randn(5, 10)
    
    # Process as batch
    batch_output, _ = batch_forward(X, weights, biases, activations)
    
    # Process individually
    individual_outputs = []
    for i in range(5):
        x = X[i:i+1]  # Keep 2D
        out, _ = batch_forward(x, weights, biases, activations)
        individual_outputs.append(out[0])
    individual_output = np.array(individual_outputs)
    
    # Compare
    max_diff = np.max(np.abs(batch_output - individual_output))
    print(f"Max difference: {max_diff}")
    print(f"Batch and individual processing are {'equivalent' if max_diff < 1e-10 else 'NOT equivalent'}!")
 
 
if __name__ == "__main__":
    verify_batch_equivalence()
    print()
    benchmark_batch_sizes()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139

import torch
import torch.nn as nn
import torch.nn.functional as F
 
class MLPWithHooks(nn.Module):
    """
    MLP with hooks for observing forward pass internals.
    
    Demonstrates:
    - Standard forward() implementation
    - Forward hooks for debugging/analysis
    - Intermediate activation extraction
    """
    
    def __init__(self, layer_sizes: list, activation: str = "relu"):
        super().__init__()
        
        self.layer_sizes = layer_sizes
        self.activation_name = activation
        
        # Create layers
        self.layers = nn.ModuleList()
        for i in range(len(layer_sizes) - 1):
            self.layers.append(nn.Linear(layer_sizes[i], layer_sizes[i+1]))
        
        # For storing activations during forward pass
        self.activations = {}
        self._register_hooks()
    
    def _get_activation_fn(self):
        """Return activation function."""
        if self.activation_name == "relu":
            return F.relu
        elif self.activation_name == "gelu":
            return F.gelu
        elif self.activation_name == "tanh":
            return torch.tanh
        else:
            raise ValueError(f"Unknown activation: {self.activation_name}")
    
    def _register_hooks(self):
        """Register forward hooks to capture activations."""
        def get_hook(name):
            def hook(module, input, output):
                self.activations[name] = output.detach()
            return hook
        
        for i, layer in enumerate(self.layers):
            layer.register_forward_hook(get_hook(f"layer_{i}"))
    
    def forward(self, x: torch.Tensor) -> torch.Tensor:
        """
        Forward pass through the network.
        
        Args:
            x: Input tensor of shape (batch_size, input_dim)
        
        Returns:
            Output tensor of shape (batch_size, output_dim)
        """
        activation_fn = self._get_activation_fn()
        
        # Forward through all layers except last
        for layer in self.layers[:-1]:
            x = activation_fn(layer(x))
        
        # Last layer without activation (raw logits)
        x = self.layers[-1](x)
        
        return x
    
    def forward_with_intermediates(self, x: torch.Tensor) -> dict:
        """
        Forward pass returning all intermediate values.
        Useful for analysis and debugging.
        """
        intermediates = {"input": x.detach().clone()}
        activation_fn = self._get_activation_fn()
        
        for i, layer in enumerate(self.layers[:-1]):
            z = layer(x)  # Pre-activation
            x = activation_fn(z)  # Post-activation
            intermediates[f"layer_{i}_pre"] = z.detach()
            intermediates[f"layer_{i}_post"] = x.detach()
        
        # Final layer
        z = self.layers[-1](x)
        intermediates["output"] = z.detach()
        
        return z, intermediates
 
 
def demonstrate_forward_pass():
    """Demonstrate forward pass with analysis."""
    
    # Create model
    model = MLPWithHooks([784, 256, 128, 64, 10], activation="relu")
    
    # Random batch
    batch_size = 32
    x = torch.randn(batch_size, 784)
    
    # Forward pass
    output = model(x)
    
    print("Forward Pass Analysis")
    print("=" * 50)
    print(f"Input shape: {x.shape}")
    print(f"Output shape: {output.shape}")
    
    print("
Layer-wise activation statistics:")
    for name, activation in model.activations.items():
        mean = activation.mean().item()
        std = activation.std().item()
        sparsity = (activation == 0).float().mean().item()
        print(f"  {name}: mean={mean:.4f}, std={std:.4f}, sparsity={sparsity:.2%}")
    
    # Forward with intermediates
    output, intermediates = model.forward_with_intermediates(x)
    
    print("
Intermediate tensor shapes:")
    for name, tensor in intermediates.items():
        print(f"  {name}: {tensor.shape}")
    
    # Verify gradient flow
    print("
Gradient flow test:")
    loss = output.sum()
    loss.backward()
    
    for i, layer in enumerate(model.layers):
        grad_norm = layer.weight.grad.norm().item()
        print(f"  Layer {i} weight gradient norm: {grad_norm:.4f}")
 
 
if __name__ == "__main__":
    demonstrate_forward_pass()