Homework #5 – Getting Started Guide

1 Logistic Regression and the MNIST Dataset

1.1 Binary Classification with Logistic Regression

Logistic regression extends linear models to classification by applying a nonlinear transformation to the output. For logistic regression specifically, this transformation maps the linear predictor to conditional probabilities \(P(Y=1|X)\) between 0 and 1, providing a probabilistic framework for making binary decisions.

1.1.1 The Sigmoid Function as a Nonlinear Detector

The logistic (sigmoid) function transforms any real-valued input into the range (0,1):

\[\sigma(z) = \frac{1}{1 + e^{-z}}\]

This transformation functions as a soft threshold detector. The function’s gradient provides critical information about sensitivity to input changes:

\[\frac{d\sigma}{dz} = \sigma(z)(1-\sigma(z))\]

The gradient reaches its maximum value of 0.25 at z=0 (where \(\sigma(z)=0.5\)), making the function most sensitive precisely at the decision boundary. This mathematical property is fundamental to the learning process in logistic regression.

Code

def sigmoid(z):
    """Numerically stable implementation of sigmoid function"""
    return np.where(z >= 0, 
                  1 / (1 + np.exp(-z)),
                  np.exp(z) / (1 + np.exp(z)))

def sigmoid_derivative(z):
    """Derivative of the sigmoid function"""
    s = sigmoid(z)
    return s * (1 - s)

# Visualize sigmoid and its derivative
z = np.linspace(-6, 6, 100)
sig = sigmoid(z)
dsig = sigmoid_derivative(z)

# Plot sigmoid and derivative
plt.figure(figsize=(10, 5))
plt.plot(z, sig, 'b-', linewidth=2, label='Sigmoid function')
plt.plot(z, dsig, 'r-', linewidth=2, label='Derivative')

# Illustrate tangent line at a specific point
point_z = 2
point_s = sigmoid(point_z)
point_d = sigmoid_derivative(point_z)

# Mark the point
plt.plot(point_z, point_s, 'ko', markersize=8)

# Draw tangent line
x_tangent = np.array([point_z - 2, point_z + 2])
y_tangent = point_s + point_d * (x_tangent - point_z)
plt.plot(x_tangent, y_tangent, 'g--', linewidth=2, 
         label=f'Tangent at z={point_z}, slope={point_d:.3f}')

plt.grid(True)
plt.legend()
plt.title('Sigmoid Function and Its Derivative')
plt.xlabel('z')
plt.ylabel('σ(z)')
plt.axhline(y=0, color='k', linestyle='-', alpha=0.3)
plt.axhline(y=1, color='k', linestyle='-', alpha=0.3)
plt.axvline(x=0, color='k', linestyle='-', alpha=0.3)

Sigmoid function with derivative and tangent line

The forward propagation in logistic regression combines the linear model with the sigmoid function:

Code

def predict_proba(X, w, b):
    """
    Forward propagation for logistic regression
    
    Parameters:
    X: input features matrix (n_samples, n_features)
    w: weight vector (n_features,)
    b: bias term (scalar)
    
    Returns:
    Predicted probabilities for class 1
    """
    z = X @ w + b
    return sigmoid(z)

1.1.2 Binary Cross-Entropy Loss

Binary cross-entropy loss has foundations in information theory and maximum likelihood estimation:

\[\mathcal{L}(w) = -\frac{1}{N}\sum_{i=1}^{N}[y_i\log(p(x_i)) + (1-y_i)\log(1-p(x_i))]\]

This loss function comes directly from maximum likelihood estimation. When modeling p(x) as the probability that y=1 given model parameters w, the likelihood of our dataset is:

\[L(w) = \prod_{i=1}^{N} p(x_i)^{y_i} (1-p(x_i))^{1-y_i}\]

Taking the negative logarithm yields the binary cross-entropy loss. This direct derivation from probability theory makes binary cross-entropy the mathematically optimal loss function for logistic regression.

Code

def binary_cross_entropy(y_true, y_pred):
    """
    Compute binary cross-entropy loss with numerical stability
    
    Parameters:
    y_true: True binary labels (0 or 1)
    y_pred: Predicted probabilities for class 1
    
    Returns:
    Average binary cross-entropy loss
    """
    epsilon = 1e-15  # Prevent log(0)
    y_pred = np.clip(y_pred, epsilon, 1 - epsilon)
    return -np.mean(y_true * np.log(y_pred) + (1 - y_true) * np.log(1 - y_pred))

1.2 Regularization

Regularization imposes constraints on the system response (weights) to prevent high-gain solutions that would amplify noise in the input data.

1.2.1 L2 Regularization (Ridge)

\[\mathcal{L}_{L2}(w) = \mathcal{L}(w) + \lambda\|w\|_2^2\]

L2 regularization applies a quadratic penalty on weights, analogous to energy constraints in signal processing. It produces a smoother decision boundary by limiting the energy in the weight vector. From a Bayesian perspective, this represents a Gaussian prior on the weights.

1.2.2 L1 Regularization (Lasso)

\[\mathcal{L}_{L1}(w) = \mathcal{L}(w) + \lambda\|w\|_1\]

L1 regularization promotes sparsity in the weight vector, similar to compressed sensing in signal reconstruction. It enforces a constraint that drives many weights to precisely zero, effectively performing feature selection by identifying the most salient pixels for detecting the target digit.

1.2.3 L1 vs L2 Regularization: Geometric Interpretation

Code

from matplotlib.patches import Circle

# Create geometric visualization of L1 vs L2 regularization
fig, ax = plt.subplots(1, 2, figsize=(12, 5))

# Set up centered coordinates
x = np.linspace(-2, 2, 100)
y = np.linspace(-2, 2, 100)
X, Y = np.meshgrid(x, y)

# True optimal point shifted from origin to show tradeoff
# Modified to be not at 45 degrees from origin for better visualization
optimal_point = np.array([1.2, 0.6])  # Intentionally asymmetric

# Standard loss function (elliptical to avoid 45-degree symmetry)
Z_loss = 2*(X-optimal_point[0])**2 + (Y-optimal_point[1])**2

# Contour levels for loss function
loss_levels = np.array([0.2, 0.5, 1.0, 2.0, 3.0])

# Regularization constraint values
l2_constraint = 0.8  # Radius of circle for L2
l1_constraint = 0.8  # Radius of diamond for L1

# L2 regularization plot
cs0 = ax[0].contour(X, Y, Z_loss, loss_levels, colors='blue', alpha=0.8, linestyles='-')
l2_circle = Circle((0, 0), l2_constraint, fill=False, color='red', 
                  linestyle='-', linewidth=2)
ax[0].add_patch(l2_circle)

# L1 regularization plot
cs1 = ax[1].contour(X, Y, Z_loss, loss_levels, colors='blue', alpha=0.8, linestyles='-')

# L1 diamond shape
l1_points = np.array([
    [l1_constraint, 0], [0, l1_constraint], 
    [-l1_constraint, 0], [0, -l1_constraint], [l1_constraint, 0]
])
ax[1].plot(l1_points[:,0], l1_points[:,1], 'r-', linewidth=2)

# Mark the true optimal point (standard loss only)
for a in ax:
    a.plot(optimal_point[0], optimal_point[1], 'b*', markersize=10, label='Loss optimum')

# Analytically determine the constrained optimal points
# For L2: Point on circle closest to optimal point (projection)
l2_direction = optimal_point / np.linalg.norm(optimal_point)  # Unit vector toward optimal
l2_optimal = l2_constraint * l2_direction  # Scale to lie on circle

# For L1: Point on diamond closest to optimal point
# Due to the asymmetric loss, the optimal point will be on the x-axis
l1_optimal = np.array([l1_constraint, 0])  # Sparse solution with w₂=0

# Plot optimal points
ax[0].plot(l2_optimal[0], l2_optimal[1], 'ko', markersize=8, label='Regularized optimum')
ax[1].plot(l1_optimal[0], l1_optimal[1], 'ko', markersize=8, label='Regularized optimum')

# Calculate loss at optimal points
l2_loss_value = 2*(l2_optimal[0]-optimal_point[0])**2 + (l2_optimal[1]-optimal_point[1])**2
l1_loss_value = 2*(l1_optimal[0]-optimal_point[0])**2 + (l1_optimal[1]-optimal_point[1])**2

# Add custom contours for the exact loss values at optimal points
ax[0].contour(X, Y, Z_loss, [l2_loss_value], colors='blue', linewidths=2)
ax[1].contour(X, Y, Z_loss, [l1_loss_value], colors='blue', linewidths=2)

# Add contour labels
plt.clabel(cs0, inline=1, fontsize=8, fmt='%.1f')
plt.clabel(cs1, inline=1, fontsize=8, fmt='%.1f')

# Labels and formatting
ax[0].set_title('L2 Regularization')
ax[0].text(-1.8, 1.7, 'Blue: Loss contours', color='blue', fontsize=9)
ax[0].text(-1.8, 1.4, 'Red: L2 constraint (‖w‖₂ ≤ c)', color='red', fontsize=9)
ax[0].text(0.9, 0.3, 'Both weights\nnon-zero', color='black', fontsize=9)
ax[0].legend(loc='upper right', fontsize=9)

ax[1].set_title('L1 Regularization')
ax[1].text(-1.8, 1.7, 'Blue: Loss contours', color='blue', fontsize=9)
ax[1].text(-1.8, 1.4, 'Red: L1 constraint (‖w‖₁ ≤ c)', color='red', fontsize=9)
ax[1].text(l1_constraint + 0.1, -0.2, 'w₂=0\n(Sparse solution)', color='black', ha='center', fontsize=9)

# Draw tangent lines at optimal points to emphasize tangency
def get_tangent_points(center, radius, point):
    # Get two points on tangent line through optimal point
    dx, dy = point[0] - center[0], point[1] - center[1]
    norm = np.sqrt(dx**2 + dy**2)
    nx, ny = -dy/norm, dx/norm  # Normal vector
    return [[point[0] - nx*0.5, point[1] - ny*0.5], 
            [point[0] + nx*0.5, point[1] + ny*0.5]]

# Add tangent line for L2
tangent_l2 = get_tangent_points([0, 0], l2_constraint, l2_optimal)
ax[0].plot([tangent_l2[0][0], tangent_l2[1][0]], 
           [tangent_l2[0][1], tangent_l2[1][1]], 
           'k--', linewidth=1, alpha=0.7)

# Set equal axes limits for better comparison
for a in ax:
    a.set_xlim(-2, 2)
    a.set_ylim(-2, 2)
    a.grid(True, linestyle='--', alpha=0.5)
    a.set_xlabel('w₁')
    a.set_ylabel('w₂')
    a.axhline(y=0, color='k', linestyle='-', alpha=0.3)
    a.axvline(x=0, color='k', linestyle='-', alpha=0.3)
    a.set_aspect('equal')  # Equal aspect ratio

plt.tight_layout()
plt.suptitle('L1 vs L2 Regularization: Geometric Interpretation', y=1.05)

Geometric comparison of L1 and L2 regularization constraints

The figure illustrates the geometric difference between L1 and L2 regularization constraints:

Left (L2): The circular constraint tangentially intersects the loss contour at a point where both weights are non-zero. The smooth boundary of the L2 constraint rarely coincides with coordinate axes.
Right (L1): The diamond-shaped constraint touches the loss contour precisely at a vertex on the w₁ axis (w₂=0). The corners of the L1 constraint align with coordinate axes, promoting sparse solutions.

For MNIST classification, L1 regularization yields sparse solutions where many pixel weights equal zero, effectively performing feature selection. L2 regularization shrinks all weights but rarely nullifies them completely.

1.2.4 Regularization as Bayesian Priors

From a Bayesian perspective, regularization terms correspond directly to prior distributions on model weights. This connection provides additional insight into why different regularization methods produce different behaviors.

In Bayesian inference, we seek to estimate the posterior distribution of weights given the data:

\[p(w|y) = \frac{p(y|w)p(w)}{p(y)}\]

Where \(p(y|w)\) is the likelihood of observing the data given the weights, and \(p(w)\) is the prior probability of the weights. Taking the negative logarithm of the posterior (and dropping the constant denominator), we get:

\[-\log p(w|y) = -\log p(y|w) - \log p(w)\]

For logistic regression with a binary cross-entropy loss function, the negative log-likelihood term is:

\[-\log p(y|w) = \sum_{i=1}^{N}[-y_i\log(p(x_i)) - (1-y_i)\log(1-p(x_i))] = \mathcal{L}(w)\]

The regularization term arises from our prior assumptions about the weights. Consider two common prior distributions:

Gaussian (Normal) Prior: When we assume weights are normally distributed around zero:

\[p(w) = \prod_{j=1}^{d} \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left(-\frac{w_j^2}{2\sigma^2}\right)\]

Taking the negative logarithm:

\[-\log p(w) = \frac{1}{2\sigma^2}\sum_{j=1}^{d} w_j^2 + \text{const} = \frac{\lambda}{2}||w||_2^2 + \text{const}\]

Where \(\lambda = 1/\sigma^2\). This gives us the L2 regularization term.
Laplace Prior: When we assume weights follow a Laplace distribution:

\[p(w) = \prod_{j=1}^{d} \frac{1}{2b} \exp\left(-\frac{|w_j|}{b}\right)\]

Taking the negative logarithm:

\[-\log p(w) = \frac{1}{b}\sum_{j=1}^{d} |w_j| + \text{const} = \lambda||w||_1 + \text{const}\]

Where \(\lambda = 1/b\). This gives us the L1 regularization term.

The critical difference lies in the shape of these distributions near zero. The Laplace distribution has a sharp peak at zero, placing higher probability density on exactly zero values than the Gaussian distribution does. This mathematical property explains why L1 regularization promotes sparsity while L2 does not.

In MNIST digit classification, using L1 regularization effectively encodes a prior belief that most pixel weights should be exactly zero, with only a few important pixels contributing to the classification decision.

1.2.5 Computing the Gradient with Regularization

Code

def compute_gradient(X, y_true, y_pred, w, lambda_reg=0, reg_type=None):
    """
    Compute gradient for binary logistic regression with regularization
    
    Parameters:
    X: Input features matrix (n_samples, n_features)
    y_true: True binary labels (0 or 1)
    y_pred: Predicted probabilities
    w: Current weights
    lambda_reg: Regularization strength
    reg_type: Type of regularization ('l1', 'l2', or None)
    
    Returns:
    Gradient of the loss with respect to weights
    """
    n_samples = X.shape[0]
    grad_loss = X.T @ (y_pred - y_true) / n_samples
    
    if reg_type == 'l2' and lambda_reg > 0:
        grad_reg = lambda_reg * w
    elif reg_type == 'l1' and lambda_reg > 0:
        grad_reg = lambda_reg * np.sign(w)
    else:
        grad_reg = 0
        
    return grad_loss + grad_reg

The gradient computation incorporates the regularization term’s derivative. For L2 regularization, this is simply \(\lambda w\), while for L1 regularization, it’s \(\lambda \cdot \text{sign}(w)\). This difference in gradients explains why L1 regularization can drive weights exactly to zero, while L2 regularization only approaches zero asymptotically.

1.3 Learning Rate Selection Strategy

The learning rate determines the step size in gradient descent. The weight update follows:

\[w^{(t+1)} = w^{(t)} - \alpha \nabla \mathcal{L}(w^{(t)})\]

For high-dimensional inputs like MNIST (784 features), smaller learning rates typically provide better stability. Common approaches include:

Fixed step size starting with a conservative value (e.g., 0.01)
Line search to adaptively find the optimal step size at each iteration
Scheduled decay that reduces the learning rate as training progresses

Code

def gradient_descent(X, y, learning_rate=0.01, max_iter=1000, 
                     lambda_reg=0, reg_type=None, tol=1e-6):
    """
    Simple gradient descent implementation for logistic regression
    
    Parameters:
    X: Input features matrix (n_samples, n_features)
    y: True binary labels (0 or 1)
    learning_rate: Step size for gradient updates
    max_iter: Maximum number of iterations
    lambda_reg: Regularization strength
    reg_type: Type of regularization ('l1', 'l2', or None)
    tol: Convergence tolerance (minimum change in loss)
    
    Returns:
    w: Learned weights
    b: Learned bias
    losses: History of loss values during training
    """
    n_features = X.shape[1]
    # Initialize weights and bias
    w = np.zeros(n_features)
    b = 0
    
    losses = []
    
    for iteration in range(max_iter):
        # Forward pass
        y_pred = predict_proba(X, w, b)
        
        # Compute loss
        loss = binary_cross_entropy(y, y_pred)
        if reg_type == 'l2':
            loss += lambda_reg * 0.5 * np.sum(w**2)
        elif reg_type == 'l1':
            loss += lambda_reg * np.sum(np.abs(w))
            
        losses.append(loss)
        
        # Check convergence
        if iteration > 0 and abs(losses[-1] - losses[-2]) < tol:
            break
            
        # Compute gradients
        grad_w = compute_gradient(X, y, y_pred, w, lambda_reg, reg_type)
        grad_b = np.mean(y_pred - y)
        
        # Update parameters
        w -= learning_rate * grad_w
        b -= learning_rate * grad_b
        
    return w, b, losses

1.4 Convergence Criteria

Establishing proper convergence criteria prevents both premature termination and unnecessary computation. Consider monitoring:

Absolute loss threshold: Stop when loss falls below a small value
Relative improvement: Stop when improvements between iterations become minimal
Validation performance: Stop when validation performance plateaus or degrades

1.5 Weight Interpretation as Filter Response

The learned weight vector w in logistic regression represents a spatial filter matched to the target pattern (digit “2”). Areas with positive weights increase the activation when pixel values are high, while negative weights suppress the response.

Visualizing this weight vector as a 28×28 image typically reveals a pattern resembling the target digit, providing insight into the features the model has learned to identify.

Code

def visualize_weights(w):
    """
    Visualize the learned weights as an image
    
    Parameters:
    w: Weight vector of length 784 (for 28x28 MNIST images)
    """
    plt.figure(figsize=(6, 6))
    plt.imshow(w.reshape(28, 28), cmap='RdBu_r')
    plt.colorbar()
    plt.title('Weight Vector Visualization')
    plt.show()

1.6 Evaluation Beyond Accuracy

The sensitivity-specificity tradeoff provides an important perspective for binary classifiers:

Accuracy: Overall correctness (potentially misleading with imbalanced classes)
Precision: Proportion of true positives among positive predictions
Recall: Proportion of true positives that are correctly identified
F1 score: Harmonic mean of precision and recall

These metrics connect to detection theory concepts of probability of detection and false alarm rates in signal processing.

Code

def evaluate_binary_classifier(y_true, y_pred_prob, threshold=0.5):
    """
    Evaluate binary classifier performance metrics
    
    Parameters:
    y_true: True binary labels (0 or 1)
    y_pred_prob: Predicted probabilities for class 1
    threshold: Classification threshold (default: 0.5)
    
    Returns:
    Dictionary of evaluation metrics
    """
    y_pred = (y_pred_prob >= threshold).astype(int)
    
    # Compute metrics
    accuracy = np.mean(y_pred == y_true)
    
    true_pos = np.sum((y_pred == 1) & (y_true == 1))
    false_pos = np.sum((y_pred == 1) & (y_true == 0))
    false_neg = np.sum((y_pred == 0) & (y_true == 1))
    true_neg = np.sum((y_pred == 0) & (y_true == 0))
    
    precision = true_pos / (true_pos + false_pos) if (true_pos + false_pos) > 0 else 0
    recall = true_pos / (true_pos + false_neg) if (true_pos + false_neg) > 0 else 0
    f1 = 2 * precision * recall / (precision + recall) if (precision + recall) > 0 else 0
    
    return {
        'accuracy': accuracy,
        'precision': precision,
        'recall': recall,
        'f1_score': f1
    }

1.7 Learning Curves Analysis

Learning curves plot the model’s performance on training and validation sets against the number of iterations. These curves help diagnose:

Underfitting: Both training and validation errors remain high
Overfitting: Training error continues to decrease while validation error increases
Optimal stopping point: Where validation error reaches its minimum

Code

def plot_learning_curves(train_losses, val_losses):
    """
    Plot learning curves to analyze training progress
    
    Parameters:
    train_losses: List of loss values on training set
    val_losses: List of loss values on validation set
    """
    plt.figure(figsize=(10, 6))
    plt.plot(train_losses, label='Training Loss')
    plt.plot(val_losses, label='Validation Loss')
    plt.xlabel('Iteration')
    plt.ylabel('Loss')
    plt.title('Learning Curves')
    plt.grid(True)
    plt.legend()
    plt.show()

1.8 Saving Model Parameters

Save your trained weights and bias as specified:

Code

def save_model(weights, bias, outFile):
    """
    Save trained model parameters to HDF5 file
    
    Parameters:
    weights: Trained weight vector (length 784)
    bias: Trained bias term (scalar)
    outFile: Output file path
    """
    with h5py.File(outFile, 'w') as hf:
        hf.create_dataset('w', data=np.asarray(weights))
        hf.create_dataset('b', data=np.asarray(bias))

The weights should be a 784-length vector and the bias a scalar value.

2 Feed-Forward Neural Networks for MNIST

2.1 MNIST Dataset Structure

The MNIST dataset contains 28×28 pixel handwritten digit images. Understanding the data organization is essential:

Code

# Create sample MNIST digit for visualization
def create_sample_digit():
    # Generate a simplified "3" digit
    digit = np.zeros((28, 28))
    
    # Top curve
    for j in range(10, 18):
        digit[5, j] = 1.0
    
    # Right side of top curve
    for i in range(5, 14):
        digit[i, 18] = 1.0
    
    # Middle line
    for j in range(10, 18):
        digit[14, j] = 1.0
    
    # Right side of bottom curve
    for i in range(14, 22):
        digit[i, 18] = 1.0
    
    # Bottom curve
    for j in range(10, 18):
        digit[22, j] = 1.0
    
    return digit

# Visualize the digit
digit = create_sample_digit()
plt.figure(figsize=(5, 5))
plt.imshow(digit, cmap='gray')
plt.axis('off')
plt.title('Sample Digit (3)')
plt.show()

# Show vector representation
flattened = digit.flatten()
print(f"Image shape: {digit.shape}, Flattened length: {len(flattened)}")

Image shape: (28, 28), Flattened length: 784

2.2 Neural Network Architecture

A Multi-Layer Perceptron (MLP) is organized as a sequence of layers, each performing a linear transformation followed by a non-linear activation:

Code

def visualize_network_architecture():
    # Create figure
    fig, ax = plt.subplots(figsize=(8, 4))
    
    # Layer sizes (simplified)
    layers = [
        {"name": "Input", "size": 784, "color": "lightskyblue", "width": 0.8, "height": 2.8},
        {"name": "Hidden 1", "size": 200, "color": "lightgreen", "width": 0.8, "height": 2.4},
        {"name": "Hidden 2", "size": 100, "color": "lightgreen", "width": 0.8, "height": 2.0},
        {"name": "Output", "size": 10, "color": "salmon", "width": 0.8, "height": 1.6}
    ]
    
    # Position layers horizontally
    spacing = 2.0
    for i, layer in enumerate(layers):
        # Draw layer box
        x = i * spacing
        y = (3 - layer["height"]) / 2  # Center vertically
        rect = plt.Rectangle((x, y), layer["width"], layer["height"], 
                            facecolor=layer["color"], edgecolor="gray", alpha=0.7)
        ax.add_patch(rect)
        
        # Add layer name and size
        ax.text(x + layer["width"]/2, y - 0.3, 
               f"{layer['name']}\n({layer['size']} neurons)", 
               ha='center', va='center', fontsize=9)
        
        # Add nodes (just a few representative ones)
        if layer["size"] <= 10:
            # Show all nodes if few
            for j in range(layer["size"]):
                node_y = y + 0.2 + j * (layer["height"] - 0.4) / max(1, layer["size"]-1)
                ax.plot(x + layer["width"]/2, node_y, 'o', markersize=6, color='white', 
                       markeredgecolor='gray')
        else:
            # Show just 3 nodes with ellipsis
            for j in range(3):
                pos = [0.1, 0.5, 0.9]  # Relative positions
                node_y = y + 0.2 + pos[j] * (layer["height"] - 0.4)
                ax.plot(x + layer["width"]/2, node_y, 'o', markersize=6, color='white', 
                       markeredgecolor='gray')
            ax.text(x + layer["width"]/2, y + layer["height"]/2, "...", fontsize=14)
        
        # Draw connections to next layer
        if i < len(layers) - 1:
            # Connection label (activation function)
            label = "ReLU" if i < len(layers) - 2 else "Softmax"
            ax.text(x + spacing/2 + layer["width"]/2, 3.3, label, ha='center', fontsize=9)
            
            # Connection arrows
            next_x = (i+1) * spacing
            for src_y in [y + 0.3, y + layer["height"]/2, y + layer["height"] - 0.3]:
                for dst_y in [layers[i+1]["height"]/2 + (3 - layers[i+1]["height"])/2]:
                    ax.arrow(x + layer["width"], src_y, 
                            next_x - (x + layer["width"]), dst_y - src_y,
                            head_width=0.05, head_length=0.1, fc='gray', ec='gray', alpha=0.3)
    
    ax.set_xlim(-0.5, (len(layers)-1) * spacing + 1.5)
    ax.set_ylim(-0.8, 3.8)
    ax.set_axis_off()
    ax.set_title("Multi-Layer Perceptron for MNIST", fontsize=12)
    
    plt.tight_layout()
    plt.show()

visualize_network_architecture()

MLP architecture for MNIST classification

2.3 ReLU Activation Function

The Rectified Linear Unit (ReLU) introduces non-linearity by setting negative values to zero:

Code

def plot_relu():
    x = np.linspace(-3, 3, 100)
    y = np.maximum(0, x)
    
    plt.figure(figsize=(6, 4))
    plt.plot(x, y, 'b-', linewidth=2)
    plt.plot(x, x, 'r--', alpha=0.5)
    plt.grid(True, alpha=0.3)
    plt.axhline(y=0, color='k', linestyle='-', alpha=0.3)
    plt.axvline(x=0, color='k', linestyle='-', alpha=0.3)
    plt.xlabel('Input (x)')
    plt.ylabel('Output (ReLU(x))')
    plt.title('ReLU Activation Function')
    plt.legend(['ReLU(x) = max(0, x)', 'Linear'])
    plt.show()
    
    # Simple implementation
    print("ReLU Implementation:")
    print("def relu(x):")
    print("    return np.maximum(0, x)")

plot_relu()

ReLU Implementation:
def relu(x):
    return np.maximum(0, x)

2.4 Softmax Function for Classification

Softmax converts raw network outputs into class probabilities that sum to 1:

Code

def plot_softmax():
    # Example raw outputs from network
    logits = np.array([2.0, 5.0, 1.0, 0.5, 3.0])
    
    # Naive (potentially unstable) implementation
    def naive_softmax(x):
        return np.exp(x) / np.sum(np.exp(x))
    
    # Numerically stable implementation
    def stable_softmax(x):
        shifted_x = x - np.max(x)
        exp_x = np.exp(shifted_x)
        return exp_x / np.sum(exp_x)
    
    # Apply softmax
    probs = stable_softmax(logits)
    
    # Create visualization
    plt.figure(figsize=(8, 4))
    
    # Plot raw values
    plt.subplot(1, 2, 1)
    plt.bar(np.arange(len(logits)), logits, color='blue', alpha=0.7)
    plt.title('Raw Network Outputs')
    plt.xlabel('Class')
    plt.ylabel('Value')
    plt.xticks(np.arange(len(logits)))
    
    # Plot probabilities
    plt.subplot(1, 2, 2)
    plt.bar(np.arange(len(probs)), probs, color='red', alpha=0.7)
    plt.title('Softmax Probabilities')
    plt.xlabel('Class')
    plt.ylabel('Probability')
    plt.xticks(np.arange(len(probs)))
    plt.ylim(0, 1)
    
    plt.tight_layout()
    plt.show()
    
    # Show implementation
    print("Softmax Implementation (Numerically Stable):")
    print("def softmax(x):")
    print("    shifted_x = x - np.max(x)")
    print("    exp_x = np.exp(shifted_x)")
    print("    return exp_x / np.sum(exp_x)")
    
    print("\nInput values:", logits)
    print("Softmax output:", probs)
    print("Sum of probabilities:", np.sum(probs))

plot_softmax()

Softmax Implementation (Numerically Stable):
def softmax(x):
    shifted_x = x - np.max(x)
    exp_x = np.exp(shifted_x)
    return exp_x / np.sum(exp_x)

Input values: [2.  5.  1.  0.5 3. ]
Softmax output: [0.04099229 0.82335225 0.01508022 0.00914662 0.11142861]
Sum of probabilities: 1.0

2.5 Forward Propagation Through Layers

The feed-forward process passes data through successive transformations:

Code

# Create simplified network components
np.random.seed(42)

# Small example (3 inputs → 2 hidden → 2 outputs)
x = np.array([[0.5], [0.1], [0.9]])         # Input
W1 = np.array([[0.1, 0.2, -0.1], 
               [-0.1, 0.1, 0.3]])           # First layer weights
b1 = np.array([[0.1], [0.2]])               # First layer bias

# First layer computation
z1 = W1 @ x + b1                            # Linear transformation
a1 = np.maximum(0, z1)                      # ReLU activation

W2 = np.array([[0.2, 0.3], [0.1, -0.2]])    # Second layer weights
b2 = np.array([[0.1], [0.2]])               # Second layer bias

# Output layer computation
z2 = W2 @ a1 + b2                           # Linear transformation
e_z2 = np.exp(z2 - np.max(z2))              # Stabilized exponential
a2 = e_z2 / np.sum(e_z2)                    # Softmax normalization

# Visualization
plt.figure(figsize=(8, 4))

# Layer positions
layer_x = [1, 4, 7]
node_y = [[1, 2, 3], [1.5, 2.5], [1.5, 2.5]]
labels = ['Input', 'Hidden\n(ReLU)', 'Output\n(Softmax)']
colors = ['skyblue', 'lightgreen', 'salmon']
values = [x.flatten(), a1.flatten(), a2.flatten()]

# Draw each layer
for i, (x_pos, y_pos, label, color, vals) in enumerate(
    zip(layer_x, node_y, labels, colors, values)):
    
    # Layer label
    plt.text(x_pos, 0.2, label, ha='center')
    
    # Draw nodes with values
    for j, (y, val) in enumerate(zip(y_pos, vals)):
        circle = plt.Circle((x_pos, y), 0.4, color=color, alpha=0.7)
        plt.gca().add_patch(circle)
        plt.text(x_pos, y, f"{val:.2f}", ha='center', va='center')

# Connections between layers
for i in range(len(layer_x)-1):
    for y1 in node_y[i]:
        for y2 in node_y[i+1]:
            plt.plot([layer_x[i], layer_x[i+1]], [y1, y2], 'k-', alpha=0.1)

# Transformation labels
plt.text((layer_x[0] + layer_x[1])/2, 3.5, "W1, b1", ha='center')
plt.text((layer_x[1] + layer_x[2])/2, 3.5, "W2, b2", ha='center')

plt.xlim(0, 8)
plt.ylim(0, 4)
plt.axis('off')
plt.tight_layout()
plt.show()

2.6 Batch Processing and Output Formatting

When processing multiple images at once, efficient matrix operations can be used:

Code

def demonstrate_batch_processing():
    # Create simplified batch of 5 images with 4 pixels each
    batch_size = 5
    image_size = 4
    
    # Random batch of images
    np.random.seed(0)
    images = np.random.rand(batch_size, image_size)
    
    # Simple model weights (4 inputs, 3 outputs)
    W = np.random.randn(3, 4) * 0.1
    b = np.zeros((3, 1))
    
    print("Batch of images shape:", images.shape)
    print("Weights shape:", W.shape)
    print("Bias shape:", b.shape)
    
    # Process one by one
    results_individual = []
    for i in range(batch_size):
        # Get single image and reshape to column vector
        img = images[i].reshape(-1, 1)
        
        # Forward pass
        z = W @ img + b
        # Simplified output (no activation for demo)
        results_individual.append(z.flatten())
    
    # Process as batch
    # Transpose images to have shape (4, 5)
    images_t = images.T
    # Compute W @ images_t to get shape (3, 5)
    z_batch = W @ images_t + b
    # Each column is the result for one image
    results_batch = z_batch.T
    
    print("\nResults match:", np.allclose(results_individual, results_batch))
    
    # Demonstrate the HDF5 output format for submission
    # activations: one row of 10 softmax outputs per image
    # yhat: predicted class (index of the maximum activation) per image
    activations = np.random.rand(batch_size, 10)
    activations /= activations.sum(axis=1, keepdims=True)
    yhat = np.argmax(activations, axis=1)
    
    with h5py.File('mlp_output.hdf5', 'w') as hf:
        hf.create_dataset('activations', data=activations)
        hf.create_dataset('yhat', data=yhat)
    
    with h5py.File('mlp_output.hdf5', 'r') as hf:
        print("\nSaved datasets:", list(hf.keys()))
        print("activations shape:", hf['activations'].shape, hf['activations'].dtype)
        print("yhat shape:", hf['yhat'].shape, hf['yhat'].dtype)

demonstrate_batch_processing()

Batch of images shape: (5, 4)
Weights shape: (3, 4)
Bias shape: (3, 1)

Results match: True

Saved datasets: ['activations', 'yhat']
activations shape: (5, 10) float64
yhat shape: (5,) int64

2.7 Visualizing MNIST Digits

Examination of correctly and incorrectly classified digits gives insights into model performance:

Code

def visualize_classification_examples():
    # Create a few synthetic MNIST-like digits
    def create_digit(digit_type):
        img = np.zeros((28, 28))
        
        if digit_type == 0:    # Zero
            for i in range(8, 20):
                for j in range(8, 20):
                    if ((i-14)**2 + (j-14)**2 <= 36) and ((i-14)**2 + (j-14)**2 >= 16):
                        img[i, j] = 1.0
        
        elif digit_type == 1:  # One
            for i in range(7, 21):
                img[i, 14] = 1.0
        
        elif digit_type == 3:  # Three
            # Top horizontal line
            for j in range(10, 18):
                img[7, j] = 1.0
            # Middle horizontal line
            for j in range(10, 18):
                img[14, j] = 1.0
            # Bottom horizontal line
            for j in range(10, 18):
                img[21, j] = 1.0
            # Right vertical lines
            for i in range(7, 14):
                img[i, 18] = 1.0
            for i in range(14, 21):
                img[i, 18] = 1.0
                
        return img
    
    # Create example digits
    digits = [create_digit(0), create_digit(1), create_digit(3)]
    
    # Random network outputs (predicted probabilities)
    predictions = [
        [0.9, 0.02, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01],  # Correct: 0
        [0.05, 0.8, 0.02, 0.01, 0.01, 0.05, 0.02, 0.02, 0.01, 0.01],  # Correct: 1
        [0.01, 0.01, 0.01, 0.7, 0.01, 0.01, 0.01, 0.20, 0.03, 0.01]   # Correct: 3
    ]
    
    # Create visualization with fixed layout
    fig = plt.figure(figsize=(12, 6))
    
    for i, (digit, pred) in enumerate(zip(digits, predictions)):
        # Main digit display
        ax1 = fig.add_subplot(2, 3, i+1)
        ax1.imshow(digit, cmap='gray')
        ax1.set_title(f"Digit example: {np.argmax(pred)}")
        ax1.axis('off')
        
        # Prediction bars in separate row below
        ax2 = fig.add_subplot(2, 3, i+4)
        ax2.bar(range(10), pred)
        ax2.set_xticks(range(10))
        ax2.set_ylim(0, 1)
        ax2.set_title("Predictions")
    
    plt.tight_layout()
    plt.show()

visualize_classification_examples()

2.8 Classification Output Processing

After forward propagation, converting outputs to a final classification involves:

Code

# Example of processing network outputs
raw_output = np.array([0.3, 1.2, -0.5, 2.1, 0.8, -0.4, 0.2, 0.5, 0.1, -0.2])

# Apply softmax
def softmax(x):
    exp_x = np.exp(x - np.max(x))
    return exp_x / np.sum(exp_x)

probabilities = softmax(raw_output)
predicted_class = np.argmax(probabilities)

# Visualization
plt.figure(figsize=(7, 4))
plt.bar(range(10), probabilities)
plt.axvline(x=predicted_class, color='red', linestyle='--')
plt.text(predicted_class, 0.01, f"Prediction: {predicted_class}", color='red',
         ha='center', bbox=dict(facecolor='white', alpha=0.8))
plt.xlabel('Digit Class')
plt.ylabel('Probability')
plt.title('Network Output Probabilities')
plt.xticks(range(10))
plt.ylim(0, 1)
plt.grid(True, alpha=0.3)
plt.show()

# Output statistics
print(f"Predicted class: {predicted_class}")
print(f"Confidence: {probabilities[predicted_class]:.4f}")

Predicted class: 3
Confidence: 0.3864