Machine Learning · Regression · Logistic Regression · Part 5 of 10

Multiclass Logistic and Softmax Regression

July 2026 · 13 min read

Binary logistic regression chooses between two classes. A problem with $K$ classes needs a probability for every class. Two common approaches are available. One-vs-rest builds $K$ separate yes-or-no models. Softmax regression fits all $K$ classes together and makes them share one unit of probability.

One-vs-Rest

For three classes, red, green, and blue, one-vs-rest fits these binary problems:

red versus not red,
green versus not green,
blue versus not blue.

The class with the largest output wins. This approach is simple, reuses any binary classifier, and allows the $K$ models to be fitted independently.

The models do not communicate during training. Their three outputs therefore are not forced to sum to $1$, and each subproblem often has one small class against a much larger rest group. Renormalizing the outputs can produce a usable probability vector, but it is an extra heuristic rather than part of the fitted likelihood. One-vs-rest is often effective when only the winning class matters. A joint softmax model is usually cleaner when consistent multiclass probabilities matter.

Multinomial Softmax Regression

Softmax regression gives class $k$ its own coefficient vector $\beta_k$. For one row, each class receives a linear score $z_k=\beta_k^\top\tilde{x}$. Softmax converts those scores into probabilities:

$$\Pr(Y=k\mid x)=\operatorname{softmax}(z)_k=\frac{e^{z_k}}{\sum_{l=1}^{K}e^{z_l}}.$$

Exponentiation makes every value positive. Dividing by the same sum makes all class probabilities add to $1$. A larger score means a larger share of the probability, but only score differences matter.

With two classes, softmax becomes ordinary logistic regression:

$$p_1=\frac{e^{z_1}}{e^{z_1}+e^{z_2}}=\frac{1}{1+e^{-(z_1-z_2)}}=\sigma(z_1-z_2).$$

The binary logit is simply the difference between the two class scores.

The Softmax and Log-Sum-Exp

Adding the same constant $c$ to every score does not change softmax. The common factor $e^c$ cancels:

$$\operatorname{softmax}(z+c\mathbf{1})_k=\frac{e^{z_k+c}}{\sum_l e^{z_l+c}}=\frac{e^{c}e^{z_k}}{e^{c}\sum_l e^{z_l}}=\operatorname{softmax}(z)_k.$$

This identity prevents numerical overflow. Directly computing $e^{1000}$ produces infinity in ordinary floating-point arithmetic. Stable softmax uses four operations:

The largest score $m=\max_l z_l$ is found.
The same $m$ is subtracted from every score.
The shifted scores are exponentiated. Their largest exponential is now $e^0=1$.
The exponentials are divided by their sum.

The loss uses the related log-sum-exp function, computed with the same shift:

$$\operatorname{logsumexp}(z)=\ln\sum_{l}e^{z_l}=m+\ln\sum_{l}e^{z_l-m}.$$

For $z=(12,10,8)$, the stable probabilities are $(0.8668,0.1173,0.0159)$ and log-sum-exp is $12.1429$. For $(1000,1001,1002)$, direct exponentiation overflows, while the shifted calculation safely returns $(0.0900,0.2447,0.6652)$.

▶ Stable Softmax by the Max-Shift

Each step performs one operation: read the scores, subtract the maximum, exponentiate, normalize, then compare the naive and stable routes on very large scores. Every value is computed at runtime.

Categorical Likelihood, Gradient, and Hessian

A one-hot label has one entry per class. The true class receives $1$, and every other class receives $0$. For example, the second of three classes is $(0,1,0)$. The loss for row $i$ is categorical cross-entropy:

$$\ell_i=-\sum_{k=1}^{K}y_{ik}\ln p_{ik}=-\ln p_{i,c_i}=\operatorname{logsumexp}(z_i)-z_{i,c_i},$$

Here, $c_i$ is the true class. The loss is just the negative log-probability assigned to that class. The final expression is numerically stable because it uses log-sum-exp instead of taking the logarithm of an already rounded probability.

The gradient keeps the same simple pattern as binary logistic regression: predicted probability minus observed label, multiplied by the features.

$$\frac{\partial \ell_i}{\partial z_{ik}}=p_{ik}-y_{ik},\qquad \nabla_{\beta_k}J=\frac{1}{n}\sum_{i=1}^{n}(p_{ik}-y_{ik})\,\tilde{x}_i,$$

There is one gradient block for each class. The curvature couples pairs of classes because increasing one class probability must reduce probability elsewhere. With $\delta_{kl}=1$ when $k=l$ and $0$ otherwise, the Hessian block is

$$\nabla^2_{\beta_k,\beta_l}J=\frac{1}{n}\sum_{i=1}^{n}p_{ik}(\delta_{kl}-p_{il})\,\tilde{x}_i\tilde{x}_i^\top.$$

Each row contributes the positive-semidefinite class matrix $\operatorname{diag}(p_i)-p_ip_i^\top$. Therefore the full cross-entropy is convex in the coefficient vectors. Newton and IRLS still apply, but their coupled linear system grows with the number of classes.

Reference-Class Identifiability

Only differences between class scores affect the probabilities. Adding the same vector $v$ to every class coefficient adds the same value $v^\top\tilde{x}$ to every score:

$$\beta_k\mapsto\beta_k+v\quad\Longrightarrow\quad p_k\text{ remains unchanged for every class}.$$

The likelihood therefore cannot distinguish coefficient sets that differ only by this common vector. Without a constraint, the Hessian contains a redundant direction and is singular.

A common solution sets one reference class coefficient to zero, such as $\beta_K=0$. Every remaining vector then describes log-odds against class $K$:

$$\ln\frac{p_k}{p_K}=\beta_k^\top\tilde{x}.$$

Another solution imposes a sum-to-zero constraint across class vectors. An L2 penalty covering all redundant components also selects the centered, minimum-norm representative. These choices produce the same probabilities; they only choose how the same model is written. With $K=2$ and one reference vector fixed to zero, the remaining vector is exactly the binary logistic coefficient.

Decision Regions

The class with the largest score also has the largest probability:

$$\arg\max_k p_k=\arg\max_k z_k.$$

Softmax changes scores into probabilities but does not change their order. Classes $k$ and $l$ tie where $\beta_k^\top\tilde{x}=\beta_l^\top\tilde{x}$. This equation defines a line in two feature dimensions and a hyperplane in higher dimensions. The model therefore divides feature space into convex regions with straight boundaries. The animation shows those regions for three classes.

▶ Piecewise-Linear Decision Regions for Three Classes

The plane is shaded by the argmax class of a three-class softmax with class 3 as reference; each step then classifies a sample point by evaluating the softmax there. Regions and probabilities are computed at runtime.

Calibration and Cost-Sensitive Decisions

Softmax probabilities sum to $1$, but that does not guarantee calibration. If class $k$ receives probability near $0.7$ on many rows, good calibration means roughly $70\%$ of those rows truly belong to class $k$. This property must be checked for each class or with a multiclass reliability analysis.

Probability estimation and the final decision are separate. Choosing the most probable class is appropriate when every kind of mistake has the same cost. With unequal costs, the selected action minimizes expected cost:

$$\hat{c}(x)=\arg\min_{a}\sum_{k=1}^{K}C(a,k)\,\Pr(Y=k\mid x),$$

Here, $C(a,k)$ is the cost of taking action $a$ when the true class is $k$. A rare but costly class can justify an action even when it is not the most probable class. Calibrated probabilities matter because expected cost depends on their numerical values, not only their ranking.

Edge Cases and the Implementation

Redundant coefficients. All $K$ unconstrained vectors contain one shared direction that probabilities cannot identify. A reference class, centering constraint, or suitable penalty removes it.
Overflow. Direct exponentiation can overflow on large scores. Max-shifted softmax and log-sum-exp keep the calculation finite.
Separation. A class perfectly separated from the others can still send coefficients toward infinity. Regularization provides the same finite remedy as in the binary case.
One-vs-rest probabilities. Independent binary outputs need not sum to $1$. Renormalization is a useful heuristic, while joint softmax models the probability competition directly.

The C++ multiclass section implements stable softmax, log-sum-exp loss, and the categorical gradient by reusing the binary building blocks. A dedicated multiclass softmax implementation chapter is also planned.