Signature
linkedinemailgithubgoogleScholar

Multinomial Logistic Regression

On this page

Multinomial Logistic Regression: From Binary to Multiclass Outcomes

Binary logistic regression assigns a probability to each of two outcomes. For an outcome YY with KK classes, we need KK nonnegative probabilities that sum to one. Softmax supplies this normalization. We derive its likelihood and gradient, then explain why fitting the coefficients requires a reference-class constraint.

Recap: The binary logistic score

In binary logistic regression, observation ii contributes

∇βℓi=xi(yi−πi).\nabla_{\boldsymbol\beta}\ell_i =\mathbf x_i(y_i-\pi_i).

The predictor xi\mathbf x_i and observed binary label yiy_i are fixed during optimization; πi\pi_i is the current model probability. Below, each free class coefficient has the same observed-minus-predicted form, with class membership replacing the binary response.

1. Model Setup

For each observation i=1,…,ni = 1,\dots,n, let xi∈R(p+1)×1\mathbf{x}_i \in \mathbb{R}^{(p+1) \times 1} denote the predictor vector (including the intercept), defined as

xi=(1, xi1, xi2, …, xip)T.\mathbf{x}_i = (1,\, x_{i1},\, x_{i2},\, \dots,\, x_{ip})^T.

The response Yi∈{1,…,K}Y_i\in\{1,\ldots,K\} is a class label. For likelihood calculations, we encode its observed value as a vector yi∈{0,1}K\mathbf{y}_i\in\{0,1\}^K, with yik=1y_{ik}=1 when observation ii belongs to class kk and zero otherwise. Thus exactly one entry equals 1. For example, when K=3K=3 and the observed class is 2, the encoding is yi=(0,1,0)T\mathbf{y}_i=(0,1,0)^T. The label and the vector carry the same information but have different mathematical roles.

We assign a coefficient vector to each class k=1,…,Kk = 1,\dots,K, with βk∈R(p+1)×1\boldsymbol{\beta}_k \in \mathbb{R}^{(p+1) \times 1}, and collect them in B=[β1, β2, …, βK]∈R(p+1)×K.\boldsymbol{B} = \bigl[\boldsymbol{\beta}_1,\, \boldsymbol{\beta}_2,\, \dots,\, \boldsymbol{\beta}_K\bigr] \in \mathbb{R}^{(p+1) \times K}. The softmax function converts their linear scores into class probabilities:

Pr⁡(Yi=k∣xi)=exp⁡ ⁣(xi⊤βk)∑j=1Kexp⁡ ⁣(xi⊤βj).\Pr(Y_i = k \mid \mathbf{x}_i) = \frac{\exp\!\bigl(\mathbf{x}_i^\top \boldsymbol{\beta}_k\bigr)}{\sum_{j=1}^{K} \exp\!\bigl(\mathbf{x}_i^\top \boldsymbol{\beta}_j\bigr)}.

Softmax Probability for Class kk

Pr⁡(Yi=k∣xi)=πik=σk(xi;B)=softmaxk(xi;B)=exp⁡ ⁣(xi⊤βk)∑j=1Kexp⁡ ⁣(xi⊤βj)\Pr(Y_i = k \mid \mathbf{x}_i) = \pi_{ik} = \sigma_k(\mathbf{x}_i; \boldsymbol{B}) = \text{softmax}_k(\mathbf{x}_i; \boldsymbol{B}) = \frac{\exp\!\bigl(\mathbf{x}_i^\top \boldsymbol{\beta}_k\bigr)}{\sum_{j=1}^K \exp\!\bigl(\mathbf{x}_i^\top \boldsymbol{\beta}_j\bigr)}

Adding the same vector to every class coefficient leaves these probabilities unchanged: the extra exponential factor cancels between numerator and denominator. To identify the coefficients, we use class KK as the reference and fix βK=0\boldsymbol{\beta}_K=\mathbf{0}. We retain all KK columns in the notation for convenience, but estimate only the first K−1K-1. This reference-class convention is also described in the Stan User's Guide.

Conditional on the fixed design X\mathbf X and coefficients B\boldsymbol B, the labels are independent with Yi∼Categorical⁡(πi1,…,πiK)Y_i\sim\operatorname{Categorical}(\pi_{i1},\ldots,\pi_{iK}). Their random one-hot encodings satisfy Yi∼Multinomial⁡(1,πi)\mathbf Y_i\sim\operatorname{Multinomial}(1,\boldsymbol\pi_i); lowercase yi\mathbf y_i denotes the observed encoding.

The likelihood treats the predictors and coefficients as fixed. In a Bayesian analysis, we instead assign a prior to the free coefficients and update it using the observed data to obtain their posterior distribution.

2. Likelihood and Estimation

For independent observations, the one-hot encoding selects the probability of each observed class: yik=1y_{ik} = 1 for that class and 0 for the others, with ∑k=1Kyik=1\sum_{k=1}^K y_{ik} = 1. The conditional likelihood is therefore:

L(B)=∏i=1n∏k=1Kπik yikL(\boldsymbol{B}) = \prod_{i=1}^{n} \prod_{k=1}^{K} \pi_{ik}^{\,y_{ik}}

Taking the logarithm, we have the log-likelihood:

ℓ(B)=∑i=1n∑k=1Kyikln⁡ ⁣(exp⁡ ⁣(xi⊤βk)∑j=1Kexp⁡ ⁣(xi⊤βj))\ell(\boldsymbol{B}) = \sum_{i=1}^{n} \sum_{k=1}^{K} y_{ik} \ln\!\left(\frac{\exp\!\bigl(\mathbf{x}_i^\top \boldsymbol{\beta}_k\bigr)}{\sum_{j=1}^K \exp\!\bigl(\mathbf{x}_i^\top \boldsymbol{\beta}_j\bigr)}\right)

Expanding the logarithm yields:

ℓ(B)=∑i=1n∑k=1Kyik[xi⊤βk−ln⁡ ⁣(∑j=1Kexp⁡ ⁣(xi⊤βj))]\ell(\boldsymbol{B}) = \sum_{i=1}^{n} \sum_{k=1}^{K} y_{ik} \Bigl[ \mathbf{x}_i^\top \boldsymbol{\beta}_k - \ln\!\Bigl(\sum_{j=1}^K \exp\!\bigl(\mathbf{x}_i^\top \boldsymbol{\beta}_j\bigr)\Bigr) \Bigr]

Distributing the sum over classes gives:

ℓ(B)=∑i=1n[∑k=1Kyik xi⊤βk  −  ∑k=1Kyik ln⁡ ⁣(∑j=1Kexp⁡ ⁣(xi⊤βj))].\ell(\boldsymbol{B}) = \sum_{i=1}^{n} \left[ \sum_{k=1}^{K} y_{ik}\,\mathbf{x}_i^\top \boldsymbol{\beta}_k \;-\; \sum_{k=1}^{K} y_{ik}\,\ln\!\left(\sum_{j=1}^{K} \exp\!\bigl(\mathbf{x}_i^\top \boldsymbol{\beta}_j\bigr)\right) \right].

The log-normalizing term ln⁡ ⁣(∑j=1Kexp⁡(xi⊤βj))\ln\!\left(\sum_{j=1}^{K} \exp(\mathbf{x}_i^\top \boldsymbol{\beta}_j)\right) does not depend on kk, so we can factor it out of the sum over classes:

∑k=1Kyik ln⁡ ⁣(∑j=1Kexp⁡(xi⊤βj))=ln⁡ ⁣(∑j=1Kexp⁡(xi⊤βj))∑k=1Kyik.\sum_{k=1}^{K} y_{ik}\,\ln\!\left(\sum_{j=1}^{K} \exp(\mathbf{x}_i^\top \boldsymbol{\beta}_j)\right) = \ln\!\left(\sum_{j=1}^{K} \exp(\mathbf{x}_i^\top \boldsymbol{\beta}_j)\right) \sum_{k=1}^{K} y_{ik}.

Because ∑k=1Kyik=1\sum_{k=1}^{K} y_{ik} = 1, the log-likelihood simplifies to:

ℓ(B)=∑i=1n[∑k=1Kyik xi⊤βk−ln⁡ ⁣(∑j=1Kexp⁡(xi⊤βj))].\ell(\boldsymbol{B}) = \sum_{i=1}^{n} \left[ \sum_{k=1}^{K} y_{ik}\,\mathbf{x}_i^\top \boldsymbol{\beta}_k - \ln\!\left(\sum_{j=1}^{K} \exp(\mathbf{x}_i^\top \boldsymbol{\beta}_j)\right) \right].

3. Assumptions

The softmax formula defines probabilities, but it does not by itself guarantee identifiable coefficients or a finite maximum likelihood estimate. Those questions require the following conditions:

  • Linear Log-Probability Ratios: Relative to the reference class, log⁡(πik/πiK)=xi⊤βk\log(\pi_{ik}/\pi_{iK})=\mathbf{x}_i^\top\boldsymbol{\beta}_k for k<Kk<K. These are not the one-versus-all logits log⁡(πik/(1−πik))\log(\pi_{ik}/(1-\pi_{ik})).
  • Independence: Each observation is assumed to be independent.
  • Finite Estimates: Complete or quasi-complete separation can prevent a finite unpenalized maximum likelihood estimate from existing.
  • Identifiability and Design Rank: We fix the reference coefficient vector and require a full-column-rank design for uniquely identified free coefficients. Strong, but not exact, collinearity can still make estimates unstable.
  • Probability Model: The softmax specification must adequately describe the class probabilities for the intended use. It models an unordered categorical response; an ordered outcome may call for a different model.

4. Derivation of the Log-Likelihood Derivative

To derive the fitting direction, return to the simplified log-likelihood:

ℓ(B)=∑i=1n[∑k=1Kyik xi⊤βk−ln⁡ ⁣(∑j=1Kexp⁡ ⁣(xi⊤βj))],\ell(\boldsymbol{B}) = \sum_{i=1}^{n} \left[ \sum_{k=1}^{K} y_{ik}\,\mathbf{x}_i^\top \boldsymbol{\beta}_k - \ln\!\left(\sum_{j=1}^{K} \exp\!\bigl(\mathbf{x}_i^\top \boldsymbol{\beta}_j\bigr)\right) \right],

We differentiate this expression with respect to a free coefficient vector βr\boldsymbol{\beta}_r, where r<Kr<K. Only the class-rr contribution in the linear term depends on this vector. For the logarithm of the normalizing sum, we apply the chain rule to the class-rr exponential:

∂ℓ(B)∂βr=∑i=1n[∂∂βr(yir xi⊤βr)−∂∂βrln⁡ ⁣(∑j=1Kexp⁡(xi⊤βj))].\frac{\partial \ell(\boldsymbol{B})}{\partial \boldsymbol{\beta}_r} =\sum_{i=1}^{n}\left[\frac{\partial}{\partial \boldsymbol{\beta}_r}\Bigl(y_{ir}\,\mathbf{x}_i^\top \boldsymbol{\beta}_r\Bigr) -\frac{\partial}{\partial \boldsymbol{\beta}_r}\ln\!\left(\sum_{j=1}^{K} \exp(\mathbf{x}_i^\top \boldsymbol{\beta}_j)\right)\right].

For the first term, we differentiate

yir xi⊤βr,y_{ir}\,\mathbf{x}_i^\top \boldsymbol{\beta}_r,

with respect to βr\boldsymbol{\beta}_r. Applying the basic rule of matrix calculus—that is, for a constant vector a\mathbf{a},

∂∂β(a⊤β)=a,\frac{\partial}{\partial \boldsymbol{\beta}} (\mathbf{a}^\top \boldsymbol{\beta}) = \mathbf{a},

we obtain:

∂∂βr(yir xi⊤βr)=yir xi.\frac{\partial}{\partial \boldsymbol{\beta}_r}\Bigl(y_{ir}\,\mathbf{x}_i^\top \boldsymbol{\beta}_r\Bigr) = y_{ir}\,\mathbf{x}_i.

The normalization term also depends on the coefficients. Applying the chain rule gives:

∂∂βrln⁡ ⁣(∑j=1Kexp⁡(xi⊤βj))=1∑j=1Kexp⁡(xi⊤βj)⋅∂∂βrexp⁡(xi⊤βr)=exp⁡(xi⊤βr)∑j=1Kexp⁡(xi⊤βj) xi.\frac{\partial}{\partial \boldsymbol{\beta}_r}\ln\!\left(\sum_{j=1}^{K} \exp(\mathbf{x}_i^\top \boldsymbol{\beta}_j)\right) =\frac{1}{\sum_{j=1}^{K} \exp(\mathbf{x}_i^\top \boldsymbol{\beta}_j)} \cdot \frac{\partial}{\partial \boldsymbol{\beta}_r}\exp(\mathbf{x}_i^\top \boldsymbol{\beta}_r) =\frac{\exp(\mathbf{x}_i^\top \boldsymbol{\beta}_r)}{\sum_{j=1}^{K} \exp(\mathbf{x}_i^\top \boldsymbol{\beta}_j)}\,\mathbf{x}_i.

Here, only the j=rj=r term depends on βr\boldsymbol{\beta}_r; defining πir=exp⁡(xi⊤βr)∑j=1Kexp⁡(xi⊤βj)\pi_{ir} = \frac{\exp(\mathbf{x}_i^\top \boldsymbol{\beta}_r)}{\sum_{j=1}^{K} \exp(\mathbf{x}_i^\top \boldsymbol{\beta}_j)} gives the result πir xi\pi_{ir}\,\mathbf{x}_i.

∂∂βrln⁡ ⁣(∑j=1Kexp⁡ ⁣(xi⊤βj))=πir xi.\frac{\partial}{\partial \boldsymbol{\beta}_r}\ln\!\left(\sum_{j=1}^{K} \exp\!\bigl(\mathbf{x}_i^\top \boldsymbol{\beta}_j\bigr)\right) = \pi_{ir}\,\mathbf{x}_i.

Subtracting the log-normalizing derivative from the linear-term derivative gives:

∂ℓ(B)∂βr=∑i=1n[yir xi−πir xi]=∑i=1n(yir−πir)xi.\frac{\partial \ell(\boldsymbol{B})}{\partial \boldsymbol{\beta}_r} =\sum_{i=1}^{n}\left[y_{ir}\,\mathbf{x}_i - \pi_{ir}\,\mathbf{x}_i\right] =\sum_{i=1}^{n}\left(y_{ir}-\pi_{ir}\right)\mathbf{x}_i.

The gradient for βr\boldsymbol{\beta}_r is a sum of predictor vectors weighted by the difference between observed class membership and predicted class probability. This is the multiclass counterpart of the binary logistic score.

Let X∈Rn×(p+1)X \in \mathbb{R}^{n \times (p+1)} have rows xiT\mathbf{x}_i^T, let Y∈Rn×KY \in \mathbb{R}^{n \times K} collect the observed one-hot encodings, and let Π∈Rn×K\Pi \in \mathbb{R}^{n \times K} have entries πik\pi_{ik}. The formal gradient with respect to all columns is:

∇Bℓ(B)=X⊤ (Y−Π).\nabla_{\boldsymbol{B}} \ell(\boldsymbol{B}) = X^\top \,(Y - \Pi).

Under the reference-class constraint, we use only its first K−1K-1 columns; the last coefficient vector remains fixed at zero.

5. Gradient Ascent Update

To maximize ℓ(B)\ell(\boldsymbol{B}), we move the free coefficients in the gradient direction. Write B−K\boldsymbol{B}_{-K} for the first K−1K-1 columns. At iteration tt, the gradient ascent update is:

B−K(t+1)=B−K(t)+η ∇B−Kℓ(B(t)),\boldsymbol{B}_{-K}^{(t+1)} = \boldsymbol{B}_{-K}^{(t)} + \eta\, \nabla_{\boldsymbol{B}_{-K}} \ell(\boldsymbol{B}^{(t)}),

where η>0\eta > 0 is the learning rate. Using the vectorized gradient expression

∇Bℓ(B)=X⊤(Y−Π),\nabla_{\boldsymbol{B}} \ell(\boldsymbol{B}) = X^\top (Y - \Pi),

the update becomes:

B−K(t+1)=B−K(t)+η X⊤(Y−K−Π−K(t)),βK(t+1)=0,\boldsymbol{B}_{-K}^{(t+1)} = \boldsymbol{B}_{-K}^{(t)} + \eta\, X^\top \left(Y_{-K} - \Pi_{-K}^{(t)}\right), \qquad \boldsymbol{\beta}_K^{(t+1)}=\mathbf{0},

with Π(t)\Pi^{(t)} denoting the softmax probability matrix at iteration tt, where the (i,k)(i,k)-th entry is given by

πik(t)=exp⁡ ⁣(xi⊤βk(t))∑j=1Kexp⁡ ⁣(xi⊤βj(t)).\pi_{ik}^{(t)} = \frac{\exp\!\bigl(\mathbf{x}_i^\top \boldsymbol{\beta}_k^{(t)}\bigr)}{\sum_{j=1}^{K} \exp\!\bigl(\mathbf{x}_i^\top \boldsymbol{\beta}_j^{(t)}\bigr)}.

Here, Y−KY_{-K} and Π−K(t)\Pi_{-K}^{(t)} retain the first K−1K-1 columns. We choose a step size that increases the log-likelihood and check convergence, for example through the norm of the gradient for the free coefficients. A small change in ℓ(B)\ell(\boldsymbol{B}) alone does not establish a finite optimum, especially when the classes can be separated. Minimizing the negative log-likelihood gives the equivalent gradient-descent procedure.

6. Applications in Real-World Problems

The same model can be used whenever observations belong to a finite set of known, mutually exclusive categories. Examples include:

  • Epidemiology: Predicting a disease subtype from measured risk factors.
  • Marketing: Predicting membership in predefined customer groups from purchasing behavior.
  • Natural Language Processing: Classifying text documents into multiple categories.
  • Customer Choice: Modeling which of several products a customer selects.
  • Medical Diagnosis: Differentiating among various diagnoses using multiple clinical indicators.

In these settings, the output is a probability for each class rather than only a predicted label. Those probabilities still require assessment on data relevant to the intended use. A fitted coefficient describes a change in the log probability ratio of a class to the reference class, holding other predictors fixed—not a causal effect or a guarantee of predictive accuracy.

Conclusion

The softmax function extends binary logistic regression to multiple classes while preserving a simple gradient: observed class indicators minus predicted probabilities, multiplied by the predictors. The reference-class constraint removes a redundancy in the coefficient representation. Together, the probability model, gradient, and constraint specify what we fit; model checking then determines how useful that fit is for the application.