from _api_doc_utils import *MultinomialLogit
Multiclass logistic regression
1 Where it fits
Group: Regression
MultinomialLogit generalizes binary logit to \(K\) classes with softmax probabilities:
\[ \Pr(Y_i=k\mid X_i=x_i)=\frac{\exp(\alpha_k+x_i'\beta_k)}{\sum_\ell \exp(\alpha_\ell+x_i'\beta_\ell)}. \]
The summary identifies the last sorted class as reference_class and reports identifiable class-versus-reference coefficient contrasts. Fisher-information standard errors are available only for alpha=0.
2 Likelihood, parameterization, and solver
For one-hot outcomes \(y_{ic}\) and logits \(\eta_{ic}=\alpha_c+x_i'\beta_c\), the native solver minimizes
\[ Q(\alpha,\beta) = -\sum_{i=1}^n\sum_{c=1}^C y_{ic}\log p_{ic} +\frac{\lambda}{2}\sum_{c=1}^C\|\beta_c\|_2^2, \qquad p_{ic}=\frac{\exp(\eta_{ic})}{\sum_\ell\exp(\eta_{i\ell})}. \]
Intercepts are unpenalized. The implementation uses a row-wise log-sum-exp shift for the objective and softmax gradient, then runs ten-vector-memory L-BFGS with a More-Thuente line search from an all-zero parameter vector. L-BFGS operates on all \(C\) coefficient columns. This full softmax parameterization is redundant because adding a common coefficient vector to every class leaves probabilities unchanged. With \(\lambda>0\), the common slope direction is pinned down by the penalty, while the common intercept direction remains flat; fitted probabilities and reported contrasts are invariant to it. The public summary removes the redundancy by reporting
\[ \delta_c=(\alpha_c-\alpha_r,\ \beta_c-\beta_r'), \]
against the last sorted class \(r\). Reaching max_iterations is a failed fit: fit() raises ValueError and clears any previous fitted state. Successful summaries include converged, iterations, termination_reason, and the final penalized objective.
3 Inference
Inference is available only at \(\lambda=0\) and is computed directly in the \((C-1)\) reference-class contrast parameterization. For non-reference classes \(a,b\) and \(\tilde x_i=(1,x_i')'\), the information blocks are
\[ H_{ab} = \sum_i \hat p_{ia} \{\mathbf1(a=b)-\hat p_{ib}\} \tilde x_i\tilde x_i'. \]
The returned covariance is \((H+10^{-8}I)^{-1}\). It is model-based only; there is no robust sandwich or Wald-test method on this class. Penalized fits return no covariance. The pairs bootstrap reports class-versus-reference contrasts, but it raises an error if any resample omits an outcome class.
4 Performance and numerical behavior
An L-BFGS objective or gradient evaluation costs \(O(npC)\). The optimizer stores \((p+1)C\) parameters and limited-memory history; the native cost and gradient use \(O(C)\) row scratch rather than materializing an \(n\times C\) probability matrix. Public prediction does return a dense \(n\times C\) array. Fisher inference has dimension \((p+1)(C-1)\); constructing its blocks is costly for many classes and dense inversion is cubic in that total dimension. The full fitted parameterization is unidentified in its common intercept direction, and is also unidentified in common slope directions when \(\lambda=0\), even though reported contrasts are identified. Rare classes make both the Fisher matrix and pairs bootstrap fragile.
5 Python API
Constructor: cm.MultinomialLogit
Use integer class labels in fit(x, y_int32); at least two distinct classes are required and class order is sorted. predict(x) returns an \(n\times C\) probability matrix, predict_lin(x) returns logits, and predict_label(x) returns original class labels. summary() returns fit diagnostics and contrast rows aligned with class_labels; penalized fits mark inference unavailable and omit se/vcov.
print(inspect.signature(cm.MultinomialLogit))(alpha=0.0, max_iterations=100, gradient_tolerance=0.0001)cls = cm.MultinomialLogit
display(HTML(html_table(["Public method"], public_methods(cls))))| Public method |
|---|
bootstrap(self, /, n_bootstrap, seed=None) |
fit(self, /, x, y) |
predict(self, /, x) |
predict_label(self, /, x) |
predict_lin(self, /, x) |
summary(self, /) |
6 Minimal example
rng = np.random.default_rng(6)
x = rng.normal(size=(240, 2))
logits = x @ np.array([[0.6, -0.3], [-0.4, 0.5], [0.2, 0.2]]).T + np.array([0.1, -0.2, 0.0])
p = np.exp(logits - logits.max(axis=1, keepdims=True))
p = p / p.sum(axis=1, keepdims=True)
y = np.array([rng.choice(3, p=row) for row in p], dtype=np.int32)
model = cm.MultinomialLogit(max_iterations=200)
model.fit(x, y)
fit = model.summary()
print({key: fit[key] for key in ['converged', 'iterations', 'termination_reason', 'objective']})
print(fit['coef'])
print(model.predict(x[:5])){'converged': True, 'iterations': 9, 'termination_reason': 'Solver converged', 'objective': 237.44870578984387}
[[ 0.25831439 0.53237632 -0.55020927]
[-0.15395916 -0.36374742 0.31279128]]
[[0.29713335 0.35470924 0.34815742]
[0.10437552 0.60469433 0.29093015]
[0.35652528 0.30570259 0.33777213]
[0.27898795 0.37429089 0.34672116]
[0.3679027 0.30230937 0.32978794]]7 summary() contract
The table below is generated by fitting the live class in this repository and then inspecting summary(). Shapes are shown because most values are plain NumPy arrays or scalars.
rng = np.random.default_rng(106)
x = rng.normal(size=(100, 2))
logits = x @ np.array([[0.6, -0.3], [-0.4, 0.5], [0.2, 0.2]]).T
p = np.exp(logits - logits.max(1, keepdims=True))
p = p / p.sum(1, keepdims=True)
y = np.array([rng.choice(3, p=row) for row in p], dtype=np.int32)
model = cm.MultinomialLogit(max_iterations=200)
model.fit(x, y)
summary = model.summary()
display(HTML(html_table(["summary() key", "shape"], summary_shape_rows(summary))))| summary() key | shape |
|---|---|
coef |
(2, 3) |
class_labels |
(2,) |
reference_class |
() |
penalty |
() |
inference_available |
() |
converged |
() |
iterations |
() |
termination_reason |
() |
objective |
() |
se |
(2, 3) |
vcov |
(6, 6) |