III-16 · Multiclass Classification and Softmax Regression

Inspect a gradient with one column per class

Edit a multiclass weight matrix, normalize shifted logits, inspect cross entropy and gradient shapes, and replay real softmax updates.

Probability and threshold laboratory

Class logits → shifted probabilities → gradient matrix

Hand-authored three-class toy batch; no seed. Labels are class indices 0, 1 and 2.

Class probabilities from the active weight matrix

s1: actual 2, predicted 2; row sum ≈ 1

Class 0≈ 0.09
Class 1≈ 0.2447
Class 2≈ 0.6652

s2: actual 1, predicted 0; row sum ≈ 1

Class 0≈ 0.3333
Class 1≈ 0.3333
Class 2≈ 0.3333

s3: actual 0, predicted 0; row sum ≈ 1

Class 0≈ 0.6652
Class 1≈ 0.2447
Class 2≈ 0.09

Change the active experiment

W row 1 class 0-1000 to 1000 · Enter to applyW row 1 class 1-1000 to 1000 · Enter to applyW row 1 class 2-1000 to 1000 · Enter to applyW row 2 class 0-1000 to 1000 · Enter to applyW row 2 class 1-1000 to 1000 · Enter to applyW row 2 class 2-1000 to 1000 · Enter to apply

Common logit shift-1000 to 1000 · Enter to applySelected multiclass rowMulticlass feature x-3 to 3 · Enter to applyClass labelLearning rate η0.001 to 1 · Enter to apply

X · rows are samplesIntercept column is fixed at 1.: 3 × 2
W, ∇W: 2 × 3
Z, P, one-hot Y: 3 × 3
Class-index labels: 3 × 1
Cross entropy: mean / n

Gradient and prospective weight matrix; update has not yet been applied
Row	∇W	Next W
1	[≈ 0.0295, ≈ −0.0591, ≈ 0.0295]	[≈ −0.0059, ≈ 0.0118, ≈ −0.0059]
2	[≈ −0.1416, ≈ 0, ≈ 0.1416]	[≈ 1.0283, ≈ 0, ≈ −1.0283]

Live state

Mean cross entropy: ≈ 0.6379

Z=XW+c\mathbf1_{3\times3},\quad c\approx0,\quad P_{ic}=\frac{e^{Z_{ic}-m_i}}{\sum_j e^{Z_{ij}-m_i}}

L=\frac1{3}\sum_i\left(\log\sum_c e^{Z_{ic}-m_i}-(Z_{i,y_i}-m_i)\right)\approx 0.6379

\nabla_W L=\frac1{3}X^\mathsf{T}(P-Y),\quad W^+=W-\eta\nabla_W L

Rows of Z, shifted Z, P and Y; labels are class indices, Y is one-hot
ID	Logits	Shifted logits	Probabilities	Y	CE
s1	[≈ −1, ≈ 0, ≈ 1]	[≈ −2, ≈ −1, ≈ 0]	[≈ 0.09, ≈ 0.2447, ≈ 0.6652]	[0, 0, 1]	≈ 0.4076
s2	[≈ 0, ≈ 0, ≈ 0]	[≈ 0, ≈ 0, ≈ 0]	[≈ 0.3333, ≈ 0.3333, ≈ 0.3333]	[0, 1, 0]	≈ 1.0986
s3	[≈ 1, ≈ 0, ≈ −1]	[≈ 0, ≈ −1, ≈ −2]	[≈ 0.6652, ≈ 0.2447, ≈ 0.09]	[1, 0, 0]	≈ 0.4076

Change the controls first, then inspect the active equation, numeric rows and shapes. Exact inputs and approximate computed values are distinguished; Reset restores the hand-authored experiment.

Change one class score and inspect all probabilities

Edit a weight-matrix cell in the opening sheet. Each row of logits is normalized across classes, so increasing one class score changes every class probability in that sample. Press Apply one update and inspect the gradient and prospective next matrix before replaying three updates.

Set Common logit shift to $1000$ . The logits change, but the probabilities, cross entropy and gradient stay the same up to binary64 rounding. The mechanism subtracts a row maximum before exponentiating. Reset restores the original hand-authored three-class batch and weights.

Define class columns and one-hot targets

With $n$ samples, $d$ augmented features and $k$ mutually exclusive classes,

X\in\mathbb R^{n\times d},\quad W\in\mathbb R^{d\times k},\quad Z=XW\in\mathbb R^{n\times k}.

Each column of $W$ is one class’s coefficient column; the intercept is the first feature row. The integer label column has shape $n\times1$ . Its one-hot encoding $Y$ has shape $n\times k$ : row $i$ has a one at class index $y_i$ and zeros elsewhere.

Softmax normalizes each logit row:

P_{ic}=\frac{e^{Z_{ic}}}{\sum_j e^{Z_{ij}}},\qquad \sum_c P_{ic}=1.

The class prediction is the largest probability, equivalently the largest logit. Here exact maximum ties choose the lowest class index, an explicit deterministic decision convention. These are mutually exclusive classes; independent binary sigmoids would answer a different multilabel problem.

Derive stable cross entropy

For a row maximum $m_i=\max_c Z_{ic}$ , factor out the common exponential:

P_{ic}=\frac{e^{Z_{ic}-m_i}}{\sum_j e^{Z_{ij}-m_i}}.

All shifted logits are nonpositive and at least one is zero, so the exponential sum lies between one and $k$ . This avoids overflow. The class-index cross entropy can be calculated from shifted logits without taking a logarithm of a rounded probability:

\ell_i=\log\sum_c e^{Z_{ic}-m_i}-(Z_{i,y_i}-m_i),\qquad L=\frac1n\sum_i\ell_i.

Adding a common scalar to a logit row changes neither probabilities nor this loss. It does not make the score matrix uniquely identifiable; class probabilities depend on relative scores.

Derive the gradient matrix

Differentiate the row log-sum-exp and subtract the selected class’s logit. The result is

\frac{\partial L}{\partial Z}=\frac1n(P-Y)\in\mathbb R^{n\times k}, \qquad \nabla_WL=\frac1nX^\mathsf T(P-Y)\in\mathbb R^{d\times k}.

There is no MSE factor of two. The $1/n$ belongs exactly once in this mean-loss gradient. Summing each row of $P-Y$ gives zero; each feature row of the weight gradient therefore sums to zero up to rounding. This checks class normalization and the one-hot target representation.

Gradient descent updates the entire old-state matrix simultaneously: $W^+=W-\eta\nabla_WL$ . The UI displays the gradient shape, each numeric row and the prospective next row. Replay frames contain actual detached weight matrices, not pictures of a changing bar chart.

Work the three-row batch

The hand-authored batch has augmented rows $[1,-1]$ , $[1,0]$ , $[1,1]$ , class labels $[2,1,0]^\mathsf T$ and

W=\begin{bmatrix}0&0&0\\1&0&-1\end{bmatrix}.

There is no random seed, downloaded dataset or measured-population claim. Logits are $[-1,0,1]$ , $[0,0,0]$ and $[1,0,-1]$ . Their probability rows are approximately

[0.0900306,0.2447285,0.6652410],\quad [1/3,1/3,1/3],\quad [0.6652410,0.2447285,0.0900306].

The per-row losses are approximately $0.407606$ , $1.098612$ and $0.407606$ , giving mean $0.6379414059$ . The middle row predicts class zero by the tie convention although its actual label is one. A nonzero correct-class probability still yields a finite loss.

The mean gradient is approximately

\nabla_WL=\begin{bmatrix}0.029535&-0.059070&0.029535\\-0.141597&0&0.141597\end{bmatrix}.

At rate $0.2$ , the next first row is approximately $[-0.00590699,0.01181398,-0.00590699]$ and the next second row $[1.02831931,0,-1.02831931]$ . Shapes are $X:3\times2$ , $W:2\times3$ , $Z,P,Y:3\times3$ , and gradient $2\times3$ .

Extreme logits and common mistakes

For logits $[1000,0,-1000]$ and actual class two, probability for class two underflows to zero, yet stable cross entropy is approximately $2000$ and the logit gradient is approximately $[1,0,-1]$ . Taking $-\log$ of a displayed zero probability would lose this finite logit-based calculation.

Do not normalize down the sample axis. Each sample needs its own probability distribution across classes. Do not multiply $X^\mathsf T$ by an integer label column; the gradient subtracts the one-hot matrix. Do not average twice. Do not mistake a common row shift for increased confidence; confidence changes with score differences.

Numerical limits and connections

The pure kernel supports 1–64 rows, 1–7 features and 2–6 classes with finite inputs and valid integer labels. The sheet fixes three classes and three samples, lets you edit the feature, label and six weights, and bounds weights/shift to ±1000. Nonfinite logit differences, losses or updates are rejected with a reason. Underflow follows binary64; a displayed approximate zero does not rewrite the stored value.

Native controls support numeric, keyboard and touch edits; immutable replay follows the displayed weight matrix, and editing a historical frame ends replay. Independent 80-digit Decimal unshifted probabilities verify the stable results, and centered finite differences separately verify every gradient entry. The canonical Python multiclass uses explicit NumPy matrices and one-hot targets from scratch.

Softmax generalizes logistic probability modeling and leads toward the neuron/matrix composition chapters. You should be able to state all shapes, derive $P-Y$ and $X^\mathsf T(P-Y)/n$ , explain shift invariance and retain finite loss when a wrong class probability underflows.

Execute the actual Python and NumPy

The imported source below is the runnable file used in canonical comparisons. Read the named mechanism from this chapter; the other functions coordinate the same Part experiment. From the book root:

.venv/bin/python -B code/python/p02/classification_reference.py

code/python/p02/classification_reference.py

"""III-13–16: stable logit losses, exact thresholds/ROC and softmax gradients.

Run: .venv/bin/python -B code/python/p02/classification_reference.py
Samples occupy rows; binary W is d×1, multiclass W is d×k.
"""
import json
import math
import numpy as np
from fit_reference import column, design, finite

SCORE_FIXTURE = {"scores": [.9, .8, .8, .4, .1], "labels": [1, 0, 1, 0, 1]}


def class_labels(y, n, k=2):
    column(y, n, "Labels")
    if any(int(v) != v or not 0 <= v < k for v in y):
        raise ValueError(f"Labels must be class indices 0–{k-1}")


def sigmoid(z):
    if isinstance(z, bool) or not isinstance(z, (int, float)) or not math.isfinite(z):
        raise ValueError("Logit must be finite")
    return 1 / (1 + math.exp(-z)) if z >= 0 else math.exp(z) / (1 + math.exp(z))


def bce(z, y):
    sigmoid(z)
    if isinstance(y, bool) or y not in (0, 1):
        raise ValueError("Binary label must be 0 or 1")
    return float(finite(max(z, 0) - z * y + math.log1p(math.exp(-abs(z))), "BCE"))


def learning_rate(rate):
    if isinstance(rate, bool) or not isinstance(rate, (int, float)) or not math.isfinite(rate) or not 0 < rate <= 1:
        raise ValueError("Learning rate must be in (0,1]")


def logistic(X, y, weights, rate):
    A, _ = design(X, y)
    W = column(weights, A.shape[1], "Weights")
    class_labels(y, len(A))
    learning_rate(rate)
    logits = finite(A @ W, "Logits")[:, 0].tolist()
    scores = [sigmoid(z) for z in logits]
    losses = [bce(z, label) for z, label in zip(logits, y)]
    loss = float(finite(sum(losses) / len(A), "Mean BCE"))
    gradient = finite(A.T @ (np.array(scores) - np.array(y)).reshape(-1, 1) / len(A), "Gradient")
    next_weights = finite(W - rate * gradient, "Updated weights")
    n, d = A.shape
    return {"logits": logits, "scores": scores, "losses": losses, "loss": loss,
            "gradient": gradient[:, 0].tolist(), "next": next_weights[:, 0].tolist(),
            "shapes": {"X": [n, d], "W": [d, 1], "y": [n, 1], "gradient": [d, 1]}}


def ratio(top, bottom, reason):
    return {"value": top / bottom, "reason": None} if bottom else {"value": None, "reason": reason}


def confusion(scores, y, threshold):
    column(scores, len(scores), "Scores")
    if not 1 <= len(scores) <= 64 or any(not 0 <= p <= 1 for p in scores):
        raise ValueError("Scores require 1–64 probabilities in [0,1]")
    class_labels(y, len(scores))
    if threshold is not None and (isinstance(threshold, bool) or not isinstance(threshold, (int, float)) or not math.isfinite(threshold) or not 0 <= threshold <= 1):
        raise ValueError("Threshold must be in [0,1] or the explicit above-all None sentinel")
    rows = []
    cells = dict(TP=0, FP=0, TN=0, FN=0)
    for p, label in zip(scores, y):
        prediction = int(threshold is not None and p >= threshold)
        cell = ("TP" if label else "FP") if prediction else ("FN" if label else "TN")
        cells[cell] += 1
        rows.append(dict(score=p, label=label, prediction=prediction, cell=cell))
    tp, fp, tn, fn = (cells[k] for k in ("TP", "FP", "TN", "FN"))
    return {"rows": rows, "tp": tp, "fp": fp, "tn": tn, "fn": fn,
            "precision": ratio(tp, tp + fp, "No predicted positives (TP+FP=0)."),
            "recall": ratio(tp, tp + fn, "No actual positives (TP+FN=0)."),
            "f1": ratio(2 * tp, 2 * tp + fp + fn, "No actual or predicted positives (2TP+FP+FN=0)."),
            "accuracy": (tp + tn) / len(scores)}


def roc(scores, y):
    confusion(scores, y, .5)
    positives, negatives = y.count(1), y.count(0)
    points = []
    for threshold in [None] + sorted(set(scores), reverse=True):
        c = confusion(scores, y, threshold)
        points.append({"threshold": threshold, **{k: c[k] for k in ("tp", "fp", "tn", "fn")},
                       "tpr": ratio(c["tp"], positives, "ROC TPR needs at least one positive."),
                       "fpr": ratio(c["fp"], negatives, "ROC FPR needs at least one negative.")})
    reason = None if positives and negatives else "ROC/AUC requires both actual classes."
    area = None if reason else sum((b["fpr"]["value"] - a["fpr"]["value"]) * (a["tpr"]["value"] + b["tpr"]["value"]) / 2 for a, b in zip(points, points[1:]))
    return {"points": points, "auc": {"value": area, "reason": reason}, "positives": positives, "negatives": negatives}


def softmax(logits):
    if not isinstance(logits, (list, tuple)) or not 2 <= len(logits) <= 6:
        raise ValueError("Softmax requires 2–6 logits")
    column(logits, len(logits), "Logits")
    maximum = max(logits)
    shifted = finite(np.array(logits) - maximum, "Shifted logits")
    exp = np.exp(shifted)
    total = float(np.sum(exp))
    return {"maximum": maximum, "shifted": shifted.tolist(), "probabilities": (exp / total).tolist(), "logSumExpShifted": math.log(total)}


def multiclass(X, y, W, rate, shift=0):
    A, _ = design(X, y)
    weights, _ = design(W, [0.] * len(W))
    n, d = A.shape
    k = weights.shape[1]
    if weights.shape[0] != d or not 2 <= k <= 6:
        raise ValueError("W must be d×k with 2–6 classes")
    class_labels(y, n, k)
    learning_rate(rate)
    if isinstance(shift, bool) or not isinstance(shift, (int, float)) or not math.isfinite(shift):
        raise ValueError("Common shift must be finite")
    logits = finite(A @ weights + shift, "Logits")
    normalized = [softmax(row) for row in logits.tolist()]
    P = np.array([r["probabilities"] for r in normalized])
    Y = np.eye(k)[np.array(y, dtype=int)]
    losses = [r["logSumExpShifted"] - r["shifted"][int(y[i])] for i, r in enumerate(normalized)]
    loss = float(finite(sum(losses) / n, "Mean cross entropy"))
    gradient = finite(A.T @ (P - Y) / n, "Gradient")
    next_weights = finite(weights - rate * gradient, "Updated weights")
    return {"logits": logits.tolist(), "normalized": normalized, "probabilities": P.tolist(), "losses": losses, "loss": loss,
            "gradient": gradient.tolist(), "next": next_weights.tolist(), "predicted": np.argmax(P, axis=1).tolist(), "oneHot": Y.tolist(),
            "shapes": {"X": [n, d], "W": [d, k], "logits": [n, k], "Y": [n, k], "gradient": [d, k], "labels": [n, 1]}}


def demonstration():
    X, y, w = [[1, -2], [1, -1], [1, 0], [1, 1], [1, 2]], [0, 0, 1, 0, 1], [0, 1]
    trace = []
    for i in range(4):
        result = logistic(X, y, w, .2)
        trace.append({"index": i, "weights": w, "derived": result})
        w = result["next"]
    return {"logisticTrace": trace, "threshold": confusion(SCORE_FIXTURE["scores"], SCORE_FIXTURE["labels"], .5),
            "roc": roc(SCORE_FIXTURE["scores"], SCORE_FIXTURE["labels"]),
            "softmax": multiclass([[1, -1], [1, 0], [1, 1]], [2, 1, 0], [[0, 0, 0], [1, 0, -1]], .2),
            "extremeBCE": [[z, label, bce(z, label)] for z in (-1000, 1000) for label in (0, 1)]}


if __name__ == "__main__":
    print(json.dumps(demonstration(), indent=2, allow_nan=False))