Softmax Regression for Multi-Class Classification

1. The Missing Piece: The Gradient

To train our model using gradient descent, we need to know how to adjust the weights $W$ and biases $b$ to minimize the loss. This requires calculating the gradient of the loss function, $J$, with respect to these parameters: $\frac{\partial J}{\partial W}$ and $\frac{\partial J}{\partial b}$.

The derivation involves the chain rule, as the loss is a function of the probabilities ($\hat{y}$), which are a function of the logits ($z$), which are a function of the weights and biases ($W,b$).

The full derivation can get quite involved with matrix calculus. However, it leads to a surprisingly simple and elegant result. The gradient of the cross-entropy loss with respect to the pre-activation logits $z$ is:

$\frac{\partial J}{\partial z}=\hat{y}-y_{true}$

Here, $\hat{y}$ is the vector of predicted probabilities from the softmax function, and $y_{true}$ is the one-hot encoded vector of the true label. This simple difference, predictions - actuals, is the error signal that drives learning.

To understand how we arrive at this clean result, the following video provides a detailed, step-by-step walkthrough. The presenter's approach of "pulling back the curtain" on how these derivations are often done in practice—by simplifying to scalars and then re-arranging dimensions—will likely resonate with your engineering background.

2. From Logit Gradient to Parameter Gradients

Now that we have the gradient with respect to the logits, $\frac{\partial J}{\partial z}$, we can use the chain rule to find the gradients for our parameters $W$ and $b$. Let's denote the gradient $\frac{\partial J}{\partial z}$ as d_scores.

Recall that our logits are calculated by a linear transformation: $z=XW+b$.

Gradient with respect to Weights $W$:
Using the chain rule, $\frac{\partial J}{\partial W}=\frac{\partial z}{\partial W}\cdot \frac{\partial J}{\partial z}$.
The derivative of $XW$ with respect to $W$ is $X^T$. Therefore:
$\frac{\partial J}{\partial W}=X^T\cdot (\hat{y}-y_{true})=X^T\cdot \text{d\_scores}$
The dimensions work out: If $X$ is $(N,D)$ and d_scores is $(N,C)$, then $X^T$ is $(D,N)$, and the resulting gradient $\frac{\partial J}{\partial W}$ is $(D,C)$, which matches the shape of our weight matrix $W$.

Gradient with respect to Biases $b$:
Similarly, $\frac{\partial J}{\partial b}=\frac{\partial z}{\partial b}\cdot \frac{\partial J}{\partial z}$.
The derivative of $XW+b$ with respect to $b$ is simply 1. We just need to sum the gradients from d_scores across all $N$ samples in the batch.
$\frac{\partial J}{\partial b}={\sum }_{i=1}^N({\hat{y}}^{(i)}-y_{true}^{(i)})=\text{np.sum(d\_scores, axis=0)}$
The result is a vector of shape $(1,C)$, matching our bias vector $b$.

The following guide provides a complete from-scratch implementation in NumPy that puts these gradient calculations into practice within a training loop.

3. Implementing Softmax Regression from Scratch

We're ready to combine everything. We'll create a SoftmaxRegression class that encapsulates the model's logic. This class will use the softmax and categorical_cross_entropy functions you built in the previous lesson.

Let's walk through the complete code.

import numpy as np
from sklearn.datasets import make_blobs
from sklearn.model_selection import train_test_split
import matplotlib.pyplot as plt

# Reuse the functions from our previous lesson
def softmax(z):
    if z.ndim == 1:
        z = z.reshape(1, -1)
    z_stable = z - np.max(z, axis=1, keepdims=True)
    exp_z = np.exp(z_stable)
    return exp_z / np.sum(exp_z, axis=1, keepdims=True)

def categorical_cross_entropy(y_pred_proba, y_true_indices):
    n_samples = y_pred_proba.shape[0]
    correct_class_probabilities = y_pred_proba[range(n_samples), y_true_indices]
    epsilon = 1e-15
    correct_class_probabilities = np.clip(correct_class_probabilities, epsilon, 1 - epsilon)
    log_likelihoods = -np.log(correct_class_probabilities)
    return np.sum(log_likelihoods) / n_samples

class SoftmaxRegression:
    """
    Softmax Regression classifier implemented from scratch.
    
    Parameters:
    -----------
    learning_rate : float
        The step size for parameter updates in gradient descent.
    n_epochs : int
        The number of passes over the training dataset.
    """
    def __init__(self, learning_rate=0.1, n_epochs=100):
        self.learning_rate = learning_rate
        self.n_epochs = n_epochs
        self.weights = None
        self.bias = None

    def fit(self, X, y):
        """
        Train the classifier.
        
        Parameters:
        -----------
        X : np.array, shape (n_samples, n_features)
            Training data.
        y : np.array, shape (n_samples,)
            Target labels (as integer indices).
        """
        n_samples, n_features = X.shape
        self.n_classes = len(np.unique(y))
        
        # 1. Initialize parameters
        self.weights = np.zeros((n_features, self.n_classes))
        self.bias = np.zeros((1, self.n_classes))
        
        # Gradient descent loop
        for epoch in range(self.n_epochs):
            # 2. Compute logits (linear score)
            logits = np.dot(X, self.weights) + self.bias
            
            # 3. Compute probabilities using softmax
            y_pred_proba = softmax(logits)
            
            # Optional: Compute and print loss for monitoring
            loss = categorical_cross_entropy(y_pred_proba, y)
            if epoch % 10 == 0:
                print(f"Epoch: {epoch}, Loss: {loss:.4f}")

            # 4. Compute the gradient of loss w.r.t logits (d_scores)
            # This is the y_hat - y_true part.
            d_scores = y_pred_proba.copy()
            d_scores[range(n_samples), y] -= 1
            d_scores /= n_samples # Average gradient over the batch
            
            # 5. Compute gradients for W and b
            dW = np.dot(X.T, d_scores)
            db = np.sum(d_scores, axis=0, keepdims=True)
            
            # 6. Update parameters
            self.weights -= self.learning_rate * dW
            self.bias -= self.learning_rate * db
            
    def predict_proba(self, X):
        """Returns class probabilities for each sample."""
        logits = np.dot(X, self.weights) + self.bias
        return softmax(logits)

    def predict(self, X):
        """Returns the predicted class label for each sample."""
        probabilities = self.predict_proba(X)
        return np.argmax(probabilities, axis=1)

# --- Let's test our model ---
# 1. Generate synthetic data
X, y = make_blobs(n_samples=500, centers=3, n_features=2, random_state=42)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=123)

# 2. Instantiate and train the model
model = SoftmaxRegression(learning_rate=0.1, n_epochs=100)
model.fit(X_train, y_train)

# 3. Evaluate the model
y_pred = model.predict(X_test)
accuracy = np.mean(y_pred == y_test)
print(f"\nFinal Test Accuracy: {accuracy * 100:.2f}%")

When you run this code, you should see the loss decreasing over the epochs and a final test accuracy of 100%. This is because make_blobs creates linearly separable data, which is the ideal scenario for Softmax Regression.

Conclusion

Congratulations! You have successfully implemented a complete Softmax Regression classifier from the ground up. You have now mastered a cornerstone algorithm for multi-class classification.

Key Takeaways:

• The gradient of the cross-entropy loss with respect to the logits is a simple and powerful error signal: predicted probabilities minus one-hot encoded true labels.
• Using the chain rule, we can easily find the gradients for the weights ($\frac{\partial J}{\partial W}=X^T\cdot \text{d\_scores}$) and biases ($\frac{\partial J}{\partial b}=\text{sum(d\_scores)}$).
• A full training loop consists of a forward pass to get predictions, calculation of gradients (backpropagation), and updating parameters using gradient descent.
• By structuring your implementation in a class, you've created a reusable and organized piece of code, which is excellent engineering practice.

Preview of the Next Lesson:
Our current model works perfectly on clean, linearly separable data. But what happens when the data is noisy or the model is too complex? It can lead to overfitting. In the next lesson, we will explore Regularization. You will learn how to modify the cost function and its gradient to include L1 (Lasso) and L2 (Ridge) penalties, which help prevent overfitting and improve the generalization of your models.