Logistic Regression and Softmax

1. \(0.5+1.2x_1-0.8x_2=0\), i.e. \(x_2 = \frac{0.5+1.2x_1}{0.8} = 0.625+1.5x_1\). A straight line with the normal vector \(\vw=(1.2,-0.8)\) pointing towards the positive class.

2. Holding \(x_2\) fixed, a one-unit increase in \(x_1\) increases the log-odds by 1.2. It does not increase the probability by any fixed amount — the effect on \(\hat p\) depends on where you start, being largest near \(\hat p=0.5\) (where \(\sigma'\) peaks at 0.25) and negligible in the tails. This is why "a one-unit change in \(x_1\) raises the probability by 30%" is almost always a misstatement.

3. Log-odds increase by 1.2, so the odds multiply by \(e^{1.2}=3.32\). Odds ratios are the natural currency of logistic regression precisely because they are constant across the range, unlike probability differences. If the odds were 1:4 (\(\hat p = 0.2\)) they become 3.32:4, i.e. \(\hat p=0.454\).

4. \(\hat p=\sigma(0.5)=0.6225\). Doubling everything gives \(\sigma(1.0)=0.7311\).

The boundary does not move — \(2(0.5+1.2x_1-0.8x_2)=0\) has exactly the same solution set — but every prediction becomes more extreme. Accuracy is unchanged; cross-entropy and calibration change completely. That is the \(\|\vw\|\) effect from the second figure, and it is why regularisation on \(\|\vw\|\) is a statement about confidence, not about the decision rule.

Proof. At the optimum \(\grad\mathcal{L}=\mX^\top(\hat{\vp}-\vy)=\vzero\), i.e. the residual vector \(\hat{\vp}-\vy\) is orthogonal to every column of \(\mX\). Applying this to the intercept column \(\mathbf{1}\):

Exactly the structure of "residuals sum to zero" for OLS (Lecture 02, exercise 2) — same normal equations, same consequence.

For a general feature \(j\): \(\sum_i x_{ij}\hat p_i = \sum_i x_{ij}y_i\). The probability-weighted total of feature \(j\) matches its label-weighted total. For a binary feature this says within that subgroup, average predicted probability equals the observed rate — calibration holds within every group defined by a feature in the model.

Why "self-calibrating". Calibration in the aggregate is not something you have to check for — it is forced by the first-order conditions, on the training set, for every direction the model can see. This is a real advantage of logistic regression over, say, an SVM or a boosted tree, whose raw scores have no such guarantee and genuinely need a post-hoc Platt scaling or isotonic step.

Two caveats. It holds on training data, not automatically on test data — though for a well-specified model with enough data the two agree. And it holds only in directions spanned by the features: a subgroup not encoded as a feature can be badly miscalibrated while every aggregate looks perfect. That is precisely the mechanism behind many fairness failures.

Why L2 breaks it. The regularised optimality condition is \(\mX^\top(\hat{\vp}-\vy) + \lambda\begin{bmatrix}0\\ \vw_{1:}\end{bmatrix} = \vzero\). The intercept row is unpenalised, so \(\sum_i(\hat p_i - y_i)=0\) survives — aggregate calibration is preserved. But the feature rows now read \(\sum_i x_{ij}(\hat p_i-y_i) = -\lambda w_j \neq 0\), so per-feature calibration is deliberately sacrificed. The shrinkage pulls predictions towards \(\hat p = \sigma(w_0)\), making the model systematically under-confident. Usually a good trade; worth knowing you made it.

This also explains why you should never regularise the intercept. Doing so would break aggregate calibration too, biasing every prediction towards \(0.5\) for no benefit.

1. Write \(S=\sum_j e^{z_j}\), so \(\hat\pi_k = e^{z_k}/S\) and \(\partial S/\partial z_m = e^{z_m}\). Quotient rule:

(\(\delta_{km}=1\) if \(k=m\), else 0.) Note this is the softmax Jacobian \(\operatorname{diag}(\hat{\boldsymbol\pi})-\hat{\boldsymbol\pi}\hat{\boldsymbol\pi}^\top\), which is symmetric positive semi-definite with the all-ones vector in its null space — the redundancy of part 3, in matrix form.

2. Chain rule through every \(k\):

using \(\sum_k y_k=1\) for a one-hot target. ∎ The \(\hat\pi_k\) cancels exactly, as in the binary case, and for the same exponential-family reason.

3. Replacing \(z_k\) with \(z_k+c\):

Softmax is shift-invariant, hence not injective — the observation from Phase A Lecture 1, now proved.

Consequences. With \(K\) classes there are only \(K-1\) free directions: the loss has a flat direction along \(\vw_1=\cdots=\vw_K\), so the Hessian is singular and the MLE is not unique. You can fix the redundancy by setting \(\vw_K=\vzero\) (which recovers binary logistic regression exactly when \(K=2\)), or simply leave it — L2 regularisation selects the minimum-norm representative, and gradient descent from zero initialisation never leaves the identifiable subspace anyway.

The practical payoff is the numerical one: shift-invariance lets you subtract \(\max_k z_k\) before exponentiating, which is exactly the log-sum-exp trick and is what makes softmax safe for large logits (Phase A Lecture 2).

import numpy as np
rng = np.random.default_rng(0)
sig = lambda z: 1/(1+np.exp(-np.clip(z,-500,500)))

def fit(X, y, lam=1e-6, T=30000, lr=0.05):
    w = np.zeros(X.shape[1])
    for _ in range(T):
        g = X.T @ (sig(X@w) - y) / len(y)
        g[1:] += lam*w[1:]
        w -= lr*g
    return w

def ece(p, y, B=10):
    e, idx = 0.0, np.clip(np.digitize(p, np.linspace(0,1,B+1))-1, 0, B-1)
    for b in range(B):
        m = idx==b
        if m.sum(): e += m.sum()/len(y)*abs(p[m].mean()-y[m].mean())
    return e

def make(kind, n):
    Z = rng.normal(size=(n,2))
    X = np.column_stack([np.ones(n), Z])
    z = 1.5*Z[:,0] - 2*Z[:,1] + 0.3 if kind=="linear" else 2.0*(Z[:,0]**2 + Z[:,1]**2 - 1.6)
    return X, (rng.random(n) < sig(z)).astype(float)

for kind in ("linear", "quadratic"):
    Xtr, ytr = make(kind, 4000); Xte, yte = make(kind, 8000)
    for lam, tag in ((1e-6, "unreg"), (100.0, "lam=100")):
        w = fit(Xtr, ytr, lam)
        p = sig(Xte@w)
        acc = ((p>=0.5)==yte).mean()
        print(f"{kind:10s} {tag:8s} acc {acc:.4f}  ECE {ece(p,yte):.4f}  "
              f"mean p {p.mean():.4f} vs base rate {yte.mean():.4f}  "
              f"p range [{p.min():.3f},{p.max():.3f}]")

# Platt scaling on a calibration split
Xtr, ytr = make("quadratic", 4000); Xca, yca = make("quadratic", 4000)
Xte, yte = make("quadratic", 8000)
w = fit(Xtr, ytr)
lg = lambda p: np.log(np.clip(p,1e-9,1-1e-9)/np.clip(1-p,1e-9,1))
Zc = np.column_stack([np.ones(len(yca)), lg(sig(Xca@w))])
ab = fit(Zc, yca, lam=1e-8)
p_raw = sig(Xte@w)
p_cal = sig(ab[0] + ab[1]*lg(p_raw))
print(f"\\nquadratic, raw    ECE {ece(p_raw,yte):.4f}  acc {((p_raw>=.5)==yte).mean():.4f}")
print(f"quadratic, Platt  ECE {ece(p_cal,yte):.4f}  acc {((p_cal>=.5)==yte).mean():.4f}")
print("Platt parameters a, b =", np.round(ab[::-1], 4))

1. The well-specified linear case gives ECE around \(0.005\) — essentially perfect, for the reason proved in exercise 2: the first-order conditions force it.

2. The quadratic truth with a linear model is a disaster, and much worse than accuracy suggests. The true boundary is a circle; the best linear separator is nearly useless, so accuracy lands near the base rate. But the calibration failure is qualitatively different: the model has no direction in which it can express the real structure, so its predicted probabilities collapse into a narrow band around the base rate. It is simultaneously uninformative and, in the aggregate, still calibrated — because the intercept condition still holds. Aggregate calibration and usefulness are independent properties: a model that always outputs the base rate is perfectly calibrated and worthless. This is why ECE alone is never sufficient; report a proper scoring rule (log loss or Brier score), which penalises both miscalibration and lack of resolution.

3. With \(\lambda=100\) the coefficients are crushed towards zero, so \(\hat p\to\sigma(w_0)\approx\) the base rate, and the predicted range narrows drastically (something like \([0.42,0.58]\) instead of \([0.01,0.99]\)). The model is under-confident — it never commits. Aggregate ECE may still look small because the mean is right; the failure is in the spread.

4. Platt scaling fits a one-dimensional logistic recalibration on the model's own logit, using a held-out split (using the training split would fit \(a=1,b=0\) and do nothing). It repairs case 3 essentially completely — the fitted \(a\) comes out well above 1, undoing the shrinkage — and improves case 2's calibration.

What it cannot repair: Platt scaling is a monotone transformation of the score, so it cannot change the ranking of examples, and therefore cannot change accuracy, AUC, or which points fall on which side of any given quantile. If the model's ordering is wrong — as it is for the quadratic case, where the linear score is nearly uninformative — no post-hoc calibration can help. Calibration fixes the numbers, never the model. The fix for case 2 is to add quadratic features, i.e. change the inductive bias.