How does the softmax function generalize logistic regression to multiclass classification, mathematically?

Binary logistic regression uses the sigmoid function to convert a single linear score into a probability between 0 and 1. For k classes, the softmax function generalises this: given k linear scores (logits) z₁,...,zₖ, softmax computes p_i = e^{z_i} / Σⱼ e^{z_j} for each class i. This produces a valid probability distribution — all values are positive and sum to exactly 1 — by exponentiating each score (making them positive) and normalising by the sum of all exponentials.

Softmax reduces exactly to the sigmoid function in the binary case: with two classes and scores z₁, z₂, p_1 = e^{z_1}/(e^{z_1}+e^{z_2}) = 1/(1+e^{-(z_1-z_2)}) — the sigmoid of the score difference. The training objective generalises from binary cross-entropy to categorical cross-entropy: L = -Σᵢ yᵢ log(pᵢ), summed over all classes for each example, which reduces to the negative log of the predicted probability for the true class.

import numpy as np
from sklearn.linear_model import LogisticRegression

def softmax(z):
    z_shifted = z - np.max(z)         # numerical stability trick
    exp_z = np.exp(z_shifted)
    return exp_z / np.sum(exp_z)

logits = np.array([2.0, 1.0, 0.1])
probs = softmax(logits)
print(probs)        # [0.659, 0.242, 0.099] — sums to 1.0
print(probs.sum())  # 1.0

# scikit-learn handles multinomial logistic regression natively
model = LogisticRegression(multi_class='multinomial', solver='lbfgs')
model.fit(X_train, y_train)
probs_sklearn = model.predict_proba(X_test)
print(probs_sklearn.sum(axis=1))  # each row sums to 1.0

Why subtract the max before exponentiating: e^z can overflow to infinity for large z; subtracting the maximum logit before exponentiating keeps all values ≤ 1 inside the exponential while leaving the final softmax output mathematically unchanged (since the normalisation cancels the shift).