Limitations of Linear Models in Complex Classification

Traditional logistic regression performs well for simple classification tasks with limited features. However, when dealing with problems involving numerous input features, the number of polynomial terms required to capture non-linear relationships grows exponentially. Consider a classification problem with just two features: constructing a decision boundary may require adding quadratic, cubic, and higher-order interaction terms to the hypothesis function.

For datasets with hundreds of features, the combinatorial explosion of polynomial terms makes logistic regression computationally impractical. Neural networks provide an elegant solution by learning hierarchical feature representations automatically.

Neural Network Architecture Fundamentals

A biological neuron receives signals through dendrites, processes information in the cell body, and transmits output through axons. Artificial neurons model this behavior using mathematical activation functions. A single artificial neuron with a sigmoid activation function accepts multiple inputs, applies weighted summation with bias, and produces an output between 0 and 1.

Neural networks organize multiple neurons into layers. The architecture typically consists of three types of layers:

• Input Layer: Receives raw feature values from the dataset
• Hidden Layer(s): Intermediate processing units that learn feature transformations
• Output Layer: Produces final prediction values

Consider a network with three layers where the input layer has n units, the hidden layer has m units, and the output layer has k units. The notation a_i^(j) represents the activation of unit i in layer j, while Theta^(j) denotes the weight matrix controlling the transformation from layer j to layer j+1.

Forward Propagation: Vectorized Implementation

Forward propagation computes activations layer by layer from input to output. For a network with one hidden layer, the process involves:

First, compute the weighted sum for the hidden layer:

z_hidden = Theta_input @ x_with_bias
a_hidden = sigmoid(z_hidden)

Add a bias unit to the hidden layer activations:

a_hidden_with_bias = insert_bias_unit(a_hidden)

Compute the output layer activations:

z_output = Theta_hidden @ a_hidden_with_bias
hypothesis = sigmoid(z_output)

The dimensionality of weight matrices follows a predictable pattern. If layer j contains s_j units and layer j+1 contains s_{j+1} units, then Theta^(j) has dimensions s_{j+1} x (s_j + 1), where the +1 accounts for the bias term.

Feature Learning in Neural Networks

Neural networks differ fundamentally from logistic regression in how they handle features. In logistic regression, features are fixed inputs from the dataset. In neural networks, hidden layers learn to construct new features from the original inputs through the learned weight matrices. Each hidden unit learns to detect specific patterns or combinations of input features, creating a hierarchy of increasingly abstract representations.

Cost Function for Multi-Class Classification

For a multi-class classification problem with K output classes, the network produces a K-dimensional output vector. The cost function extends the logistic regression cost to handle multiple outputs:

J(Theta) = (-1/m) * sum(sum(y_k * log(h_k) + (1-y_k) * log(1-h_k)))
           + (lambda/2m) * sum(sum(sum(Theta_ij^2)))

The regularization term sums the squares of all weights, excluding bias terms, across all layers. This penalizes large weight values and helps prevent overfitting.

Backpropagation Algorithm

Backpropagation efficiently computes partial derivatives of the cost function with respect to each weight. The key insight is propagating error signals backward through the network.

Define delta^(l) as the error term for layer l. For the output layer:

delta_output = hypothesis - y_true

For hidden layers, propagate the error backward:

delta_hidden = (Theta_hidden.T @ delta_output) * sigmoid_gradient(z_hidden)

The sigmoid gradient computes the derivative of the activation function:

sigmoid_gradient(z) = sigmoid(z) * (1 - sigmoid(z))

Accumulate the gradient contributions from all training examples, then compute the partial derivatives.

Gradient Verification

Backpropagation implementation is error-prone. Numerical gradient checking verifies correctness by comparing analytical gradients with finite difference approximations:

numerical_gradient = (J(theta + epsilon) - J(theta - epsilon)) / (2 * epsilon)

For parameter vectors, compute numerical gradients for each element individually. Verify that the analytical gradients from backpropagation closely match the numerical approximation. After verification, disable gradient checking during actual training to avoid computational overhead.

Training Procedure Summary

Training a neural network involves the following steps:

1. Initialize all weight matrices with small random values to break symmetry
2. Implement forward propagation to compute predictions
3. Compute the cost function value
4. Implement backpropagation to compute gradients
5. Verify gradients using numerical checking
6. Use optimization algorithms (gradient descent, L-BFGS, etc.) to minimize the cost function

Python Implementation for MNIST Digit Recognition

The following implementation demonstrates handwritten digit classification using a three-layer neural network.

Import required libraries:

import numpy as np
from scipy.io import loadmat
from scipy.optimize import minimize
from sklearn.preprocessing import OneHotEncoder

Load and preprocess the dataset. Each digit image contains 20x20 pixels, flattened into a 400-dimensional feature vector:

def load_dataset(filepath):
    raw_data = loadmat(filepath)
    features = raw_data['X']
    labels = raw_data['y']
    return features, labels

features, labels = load_dataset('ex4data1.mat')

Encode labels using one-hot encoding for multi-class classification:

encoder = OneHotEncoder(sparse_output=False)
labels_encoded = encoder.fit_transform(labels)

Define activation functions and their derivatives:

def sigmoid(z):
    return 1.0 / (1.0 + np.exp(-z))

def sigmoid_derivative(z):
    s = sigmoid(z)
    return s * (1 - s)

Implement forward propagation through the network:

def forward_pass(weights, input_data, layer_sizes):
    theta1, theta2 = reshape_weights(weights, layer_sizes)
    
    m = input_data.shape[0]
    a1 = np.c_[np.ones(m), input_data]
    
    z2 = a1 @ theta1.T
    a2 = np.c_[np.ones(m), sigmoid(z2)]
    
    z3 = a2 @ theta2.T
    predictions = sigmoid(z3)
    
    return a1, z2, a2, z3, predictions

Compute the regularized cost function:

def compute_cost(weights, input_data, targets, layer_sizes, reg_factor=1.0):
    m = input_data.shape[0]
    _, _, _, _, predictions = forward_pass(weights, input_data, layer_sizes)
    
    cross_entropy = -np.sum(targets * np.log(predictions) + (1 - targets) * np.log(1 - predictions)) / m
    
    theta1, theta2 = reshape_weights(weights, layer_sizes)
    regularization = (reg_factor / (2 * m)) * (np.sum(theta1[:, 1:]**2) + np.sum(theta2[:, 1:]**2))
    
    return cross_entropy + regularization

Implement backpropagation to compute gradients:

def compute_gradients(weights, input_data, targets, layer_sizes, reg_factor=1.0):
    m = input_data.shape[0]
    theta1, theta2 = reshape_weights(weights, layer_sizes)
    
    a1, z2, a2, z3, predictions = forward_pass(weights, input_data, layer_sizes)
    
    delta3 = predictions - targets
    delta2 = (delta3 @ theta2)[:, 1:] * sigmoid_derivative(z2)
    
    theta1_grad = (delta2.T @ a1) / m
    theta2_grad = (delta3.T @ a2) / m
    
    theta1_grad[:, 1:] += (reg_factor / m) * theta1[:, 1:]
    theta2_grad[:, 1:] += (reg_factor / m) * theta2[:, 1:]
    
    return np.concatenate([theta1_grad.ravel(), theta2_grad.ravel()])

Initialize weights with small random values:

def initialize_weights(input_size, hidden_size, output_size):
    total_params = hidden_size * (input_size + 1) + output_size * (hidden_size + 1)
    return np.random.uniform(-0.12, 0.12, total_params)

Train the network using scipy's optimization routine:

input_dim = 400
hidden_dim = 25
output_dim = 10
layer_config = (input_dim, hidden_dim, output_dim)

initial_weights = initialize_weights(input_dim, hidden_dim, output_dim)

result = minimize(
    fun=compute_cost,
    x0=initial_weights,
    args=(features, labels_encoded, layer_config, 1.0),
    method='TNC',
    jac=compute_gradients,
    options={'maxiter': 400}
)

trained_weights = result.x

Evaluate the trained model on the training set:

_, _, _, _, final_predictions = forward_pass(trained_weights, features, layer_config)
predicted_classes = np.argmax(final_predictions, axis=1) + 1

accuracy = np.mean(predicted_classes.reshape(-1, 1) == labels)
print(f'Training Accuracy: {accuracy * 100:.2f}%')

This implementation achieves approximately 97-99% accuracy on the training set for handwritten digit recognition. The network architecture consists of 400 input units corresponding to pixel intensities, 25 hidden units that learn intermediate feature representations, and 10 output units representing digit classes 0-9.