Softmax / Multinomial Regression
Logistic regression generalized to K classes at once.
On this page
Beginner: section 2.4.1's logistic regression only ever answers a yes/no question — one score, squashed through a sigmoid, into one probability. Softmax regression answers a "which one of K?" question instead. Give every class its own linear score, then run all K scores through the softmax function, which turns them into K probabilities that are all non-negative and add up to exactly 1 — a proper probability distribution over the classes, ready to compare against a one-hot encoded true label.
Intermediate: softmax is not a new idea bolted onto logistic regression — it's the same idea, generalized. Where 2.4.1 fit one weight vector θ by minimizing binary cross-entropy, softmax regression fits K weight vectors θ₁, …, θ_K at once by minimizing categorical cross-entropy, and (as this lesson derives concretely below) setting K = 2 in the softmax formula reproduces 2.4.1's sigmoid exactly, term for term.
Advanced: two facts drive most of softmax regression's practical behavior. First, softmax is invariant to shifting every logit by the same constant, which means the K weight vectors are only identified up to a shared offset — a genuine redundancy, not a bug, that's usually resolved by fixing one class as a reference. Second, the gradient of categorical cross-entropy with respect to the logits collapses to the strikingly simple p̂ − y (predicted probabilities minus the one-hot truth) — the same "prediction error times features" pattern seen for OLS and for 2.4.1's sigmoid gradient, and (per the Going Deeper note below) an instance of a much more general exponential-family phenomenon from section 2.3.9.
One linear score z_k per class, generalizing 2.4.1's single z = θᵀx; softmax turns the K scores into a probability vector p̂; categorical cross-entropy is the negative log-probability the model assigned to each example's true, one-hot encoded class y_i — the direct multinomial generalization of 2.4.1's binary cross-entropy.
In 2.4.1, a single linear score z = θᵀx was fed through the sigmoid, mapping ℝ → (0,1), read as P(y=1∣x). Generalizing to K classes means giving every class its own weight vector and its own linear score, or logit:
Stack these into a vector z ∈ ℝ^K. What's needed now is a function turning K arbitrary real numbers into a valid probability distribution — every output non-negative, all of them summing to 1. The softmax function does exactly that, by exponentiating (which guarantees positivity, since e^x > 0 for every real x) and then dividing by the total (which guarantees the sum is exactly 1):
So for any real-valued vector z whatsoever, softmax(z) lands somewhere inside the probability simplex Δ^{K-1} = \{p ∈ ℝ^K : p_k ≥ 0, Σp_k = 1\} — a valid categorical distribution over the K classes, by construction, for every possible input.
Translation invariance. Add the same constant c to every logit at once, z + c𝟙, and factor e^{z_k+c} as e^c·e^{z_k}:
The e^c factor is identical across every term, top and bottom, and cancels — so softmax literally cannot tell z apart from any shifted version z + c𝟙. That has a real consequence for the parameters, not just the logits: if θ_k' = θ_k + v for the same vector v added to every class's weights, then z_k' = z_k + vᵀx for every k — a per-example constant shift, identical across classes — and softmax's output is completely unchanged for every input x. The model has K·D raw parameters, but only (K−1)·D of them actually affect any prediction; the rest is a genuine, provable redundancy.
The standard fix is to pick one class as a fixed reference and remove its redundant freedom entirely — most commonly, fix θ_K = 0 so only K−1 weight vectors remain free. (Derivation 4 below shows this is exactly how K = 2 softmax turns back into 2.4.1's single-θ logistic regression.) In practice, many implementations — including the diagrams and code in this lesson — simply keep the full, redundant K·D parameterization anyway, since L2 regularization or the optimizer's own initialization quietly breaks the symmetry without changing any prediction.
Computed literally, softmax needs e^{z_k} for every class — and any moderately large logit (a raw, unbounded linear score, easily in the hundreds or thousands for an untrained or poorly-scaled model) makes e^{z_k} overflow to infinity in floating point long before the ratio that would have been a perfectly ordinary probability is ever formed. Section 1.21 (Numerical Stability & the Log-Sum-Exp Trick) derives the general fix in full; here is the one-line specialization to softmax specifically. Subtracting the max logit m = max_j z_j before exponentiating is an exact identity, not an approximation — the same e^{-m} cancellation shown in Derivation 1 above, just with c = −m instead of an arbitrary shift:
Choosing m to be the max specifically (rather than any other shift) is what makes this useful: every exponent computed is now ≤ 0, so the largest term is e⁰ = 1 and nothing can overflow. See section 1.21 for the full derivation and its C++/Python implementations — it is not repeated here.
One more stability point specific to training softmax regression: computing softmax and then separately taking its log (as categorical cross-entropy needs) reintroduces a risk, since a probability that rounds to exactly 0 makes log(0) = −∞. The fix is to compute log-softmax directly, without ever forming the intermediate probability:
Where this is used: every framework's built-in log_softmax/cross_entropy function combines both steps into one numerically stable operation for exactly this reason.
Encode each label as a one-hot vector y_i ∈ {0,1}^K (exactly one entry is 1, the rest are 0). A single draw from a categorical distribution with class probabilities p̂ = softmax(z_i) has likelihood P(y_i∣x_i;θ) = ∏_k p̂_{ik}^{y_{ik}} — every factor is 1 except the true class's, since raising anything to the power 0 gives 1. Following section 2.1.6's MLE recipe exactly (log-likelihood, summed over an i.i.d. sample):
(the last step holds because the one-hot y_{ik} zeroes out every term except the true class). Minimizing a loss rather than maximizing a likelihood, flip the sign — this is exactly the categorical cross-entropy from the Formula box above:
The gradient, worked example by example. Fix one example, drop the i subscript, and write c for its true class. Using \log \hat p_c = \log \text{softmax}(z)_c = z_c - \log\sum_j e^{z_j} (the log-sum-exp identity from Derivation 2), the per-example loss is:
Differentiate with respect to a single logit z_k, in two cases. If k = c (the true class):
If k ≠ c (any other class), the −z_c term doesn't involve z_k at all, so only the log-sum-exp term contributes:
Both cases are the same formula once written with the one-hot indicator y_k (1 for the true class, 0 otherwise): ∂L/∂z_k = p̂_k − y_k in every case. Stacked over all K classes, this is the clean closed form the plain-English section promised:
— "predicted probabilities minus the one-hot truth," nothing more. This is exactly why softmax and cross-entropy are so often paired: no matter how complicated softmax and log look individually, their combined gradient with respect to the logits is this trivial vector subtraction. One more chain-rule step reaches the actual parameters, using ∂z_k/∂θ_k = x:
Where this is used: this is the exact gradient the from-scratch gradient-descent code below implements, one (p̂ − y) update per class, and it's what the training diagram animates.
Set K = 2, classes {0, 1}, and apply the reference-class fix from Derivation 1: pin the redundant parameter θ₀ = 0, keeping only θ₁ = θ free. Then z₀ = 0 and z₁ = θᵀx. Plugging directly into the softmax definition for class 1:
— the sigmoid, exactly, matching 2.4.1's P(y=1∣x) = σ(θᵀx) term for term. Since softmax's two outputs must sum to 1, softmax(z)₀ = 1 − σ(θᵀx) automatically, matching P(y=0∣x) from 2.4.1 with no extra work. The gradient agrees too: Derivation 3's general result ∂L/∂z_1 = p̂_1 − y becomes σ(θᵀx) − y — precisely 2.4.1's binary cross-entropy gradient. Softmax regression isn't merely analogous to logistic regression; for K = 2 it degenerates into the identical model, loss, and gradient.
The intro sweeps a shared shift c across all three logits -- the bars don't move, proving softmax(z+c)=softmax(z). Afterward, drag z1..z3 individually and watch how fast probability mass concentrates as one logit pulls ahead of the others.
Three Gaussian blobs, softmax regression fit with the exact grad = X^T(softmax(Z)-onehot(y)) update derived above. Shaded regions are the model's current argmax prediction across the plane; dots are the true labels. Watch the three linear decision regions carve themselves out as the loss readout falls.
#include <algorithm>
#include <cmath>
#include <cstdio>
#include <vector>
using Vec = std::vector<double>;
using Mat = std::vector<Vec>;
Mat softmaxRows(const Mat& Z) {
Mat P(Z.size(), Vec(Z[0].size()));
for (size_t i = 0; i < Z.size(); ++i) {
double m = *std::max_element(Z[i].begin(), Z[i].end());
double sum = 0.0;
for (size_t k = 0; k < Z[i].size(); ++k) {
P[i][k] = std::exp(Z[i][k] - m);
sum += P[i][k];
}
for (size_t k = 0; k < Z[i].size(); ++k) P[i][k] /= sum;
}
return P;
}
int main() {
// Small 3-class toy dataset; feature 0 is the bias column (always 1).
Mat X = {{1, -1.4, 1.1}, {1, -1.6, 0.8}, {1, 1.5, 1.0},
{1, 1.7, 0.9}, {1, 0.1, -1.5}, {1, -0.2, -1.7}};
std::vector<int> y = {0, 0, 1, 1, 2, 2};
int n = (int)X.size(), d = (int)X[0].size(), k = 3;
Mat theta(k, Vec(d, 0.0));
double lr = 0.5;
for (int iter = 0; iter < 500; ++iter) {
Mat Z(n, Vec(k, 0.0));
for (int i = 0; i < n; ++i)
for (int c = 0; c < k; ++c)
for (int j = 0; j < d; ++j) Z[i][c] += X[i][j] * theta[c][j];
Mat P = softmaxRows(Z);
Mat grad(k, Vec(d, 0.0));
for (int i = 0; i < n; ++i)
for (int c = 0; c < k; ++c) {
double diff = P[i][c] - (y[i] == c ? 1.0 : 0.0); // softmax(z) - onehot(y)
for (int j = 0; j < d; ++j) grad[c][j] += diff * X[i][j] / n;
}
for (int c = 0; c < k; ++c)
for (int j = 0; j < d; ++j) theta[c][j] -= lr * grad[c][j];
}
for (int c = 0; c < k; ++c) {
std::printf("theta[%d] = ", c);
for (double v : theta[c]) std::printf("%.3f ", v);
std::printf("\n");
}
return 0;
}- Handwritten digit recognition (MNIST-style) — a softmax over 10 classes is the textbook final layer for "which digit is this," the simplest possible case of the training diagram above scaled up from 3 classes to 10.
- Document / topic categorization — routing a news article, support ticket, or email into one of several mutually exclusive topics or departments.
- Image classification baselines — almost every deep convolutional or transformer classifier ends in a linear layer followed by softmax and categorical cross-entropy, making this exact lesson's math the final stage of a much larger model.
- Product categorization — assigning a single catalog category to a product listing from a fixed, exhaustive taxonomy.
- Language identification — one probability per candidate language for a piece of text, normalized to sum to 1 across a closed set of languages.
- Computing
softmaxnaively (rawexpdivided by a raw sum) instead of shifting by the max first — the single most common source ofNaNlosses when logits get even moderately large, as derived in Derivation 2 above and covered in full in section 1.21. - Forgetting the parameter redundancy from Derivation 1 — comparing raw weight magnitudes across two separately-trained softmax models (or assuming a specific
θ_kis individually "the" weight for classk) is meaningless without first pinning the same reference class in both, since any shared shiftvadded to everyθ_kleaves every prediction unchanged. - Reaching for softmax on a multi-label problem, where an example can genuinely belong to more than one class at once (e.g. a movie can be both "comedy" and "romance"). Softmax's
Σp_k = 1normalization actively assumes the classes are mutually exclusive — forcing probability mass to compete across labels that aren't actually competing. Section 2.4.5 (Multi-label Classification) covers the right alternative: independent per-label sigmoids, not one shared softmax.
Going deeper
The clean p̂ − y gradient derived above isn't a coincidence specific to softmax — it's the multi-class instance of a general pattern. Section 2.3.9 showed that any exponential-family GLM fit under its canonical link has log-likelihood gradient X^\top(y-\mu), where μ is the distribution's mean under the model. The categorical/multinomial distribution is itself an exponential-family member, and softmax is precisely its canonical link — so ∇L = X^\top(y-\hat p) here (up to the loss-vs-likelihood sign flip) is the exact same phenomenon as 2.3.9's Bernoulli-with-sigmoid case, just with a categorical distribution and K-way mean vector standing in for a Bernoulli and its scalar mean.
A second, more practical aside: dividing every logit by a temperature T before applying softmax, softmax(z/T), changes nothing about the ordering of the probabilities but a great deal about their sharpness. As T → 0, softmax sharpens toward a one-hot arg-max (a near-certain prediction); as T → ∞, it flattens toward the uniform distribution (maximal uncertainty). This single knob underlies both temperature-scaled model calibration and knowledge distillation, where a "student" model is trained against a softened, higher-temperature version of a "teacher" model's output distribution.
If you add the same constant c to every one of the K logits before applying softmax, how do the output probabilities change — and what does that imply about how many independent parameters a K-class softmax model with D features per class actually has?
The probabilities don't change at all -- softmax(z+c1) = softmax(z) exactly, because the shared e^c factor cancels between every term of the numerator and denominator. That means only the DIFFERENCES between logits matter, never their absolute scale, so out of the model's K·D raw weights, only (K-1)·D actually influence any prediction. The remaining D degrees of freedom are a genuine redundancy -- which is why it's standard to fix one class's weight vector (commonly θ_K = 0) as a reference and treat only the other K-1 as free parameters.
Softmax regression is 2.4.1's logistic regression generalized to K classes at once: same linear scores, same cross-entropy-via-MLE derivation, same clean p̂ − y gradient, fit the same way — by smooth, convex, gradient-based optimization. The next lesson, 2.4.3 (Perceptron), is a deliberate change of subject: an older, mistake-driven linear classifier with no probability, no loss surface, and no gradient at all, whose update rule fires only when a prediction is wrong. Seeing that contrast directly is the fastest way to appreciate why the field largely moved from mistake-driven rules toward the probabilistic, gradient-based fitting this and the previous lesson both use.