10  Logistic Regression — Additional Material

From Newton’s method to iteratively reweighted least squares

NoteOptional depth after the common experiment

The main chapter derives Bernoulli likelihood, the gradient, and Newton’s method. This page connects the new method to the least-squares solution from Week 1. It assumes you can read a matrix product and solve a small linear system. Work through the scalar interpretation even if the matrix derivation is new.

Try → understand → improve: try one ordinary least-squares fit to binary labels; identify what it does not model; derive a sequence of weighted least-squares problems that does optimize logistic likelihood.

10.1 Try: apply the Week 1 solution directly

Use exactly the chapter’s two features and training split. The scalar targets are zeros and ones. Least squares can fit these numbers, but an affine output is not constrained to be a probability.

import numpy as np
import matplotlib.pyplot as plt
from dlbook.plotting import show_plot

from dlbook.logistic import (
    design, gradient, hessian, load_wdbc, loss, newton, sigmoid, stratified_split,
)

X_all, y_all = load_wdbc()
train, valid, test = stratified_split(y_all, seed=8)
X, y = X_all[train], y_all[train]
mu, scale = X.mean(axis=0), X.std(axis=0)
A = design((X - mu) / scale)
theta_ls = np.linalg.lstsq(A, y, rcond=None)[0]
linear_outputs = A @ theta_ls
print("Least-squares output range:", linear_outputs.min(), linear_outputs.max())
print("Training outputs outside [0, 1]:", np.sum((linear_outputs < 0) | (linear_outputs > 1)))
Least-squares output range: -0.1423821544587424 1.7423487490396152
Training outputs outside [0, 1]: 52

Predict. Would clipping the outputs to \([0,1]\) turn this estimate into the logistic MLE? Would applying a sigmoid after fitting least squares solve the logistic stationarity equation? Neither transformation changes the objective that was optimized. Test the latter prediction:

print("Logistic gradient at the least-squares coefficients:", gradient(A, y, theta_ls))
theta_mle, mle_history, _ = newton(A, y)
print("Logistic loss at LS coefficients:", loss(A, y, theta_ls))
print("Logistic loss at MLE:", loss(A, y, theta_mle))
Logistic gradient at the least-squares coefficients: [ 0.21596935 -0.27471109 -0.15596082]
Logistic loss at LS coefficients: 0.6403967294105761
Logistic loss at MLE: 0.262413276051727

Ordinary least squares on binary responses is sometimes called a linear probability model. It can be studied on its own terms; it is not logistic maximum likelihood.

10.2 Understand: derive an iteratively reweighted problem

At the current iterate \(\theta_t\), define

\[ \boldsymbol z_t=A\theta_t,\quad \boldsymbol p_t=\sigma(\boldsymbol z_t),\quad r_i=p_{it}(1-p_{it}),\quad R_t=\operatorname{diag}(r_i). \]

For finite scores in exact arithmetic, every \(r_i>0\). The unpenalized mean negative log likelihood has

\[ \boldsymbol g_t=\frac1nA^\top(\boldsymbol p_t-\boldsymbol y),\qquad H_t=\frac1nA^\top R_t A. \]

A full Newton step solves \(H_t(\theta_{t+1}-\theta_t)=-\boldsymbol g_t\). Multiply by \(n\) and rearrange:

\[ A^\top R_t A\theta_{t+1} =A^\top R_t A\theta_t+A^\top(\boldsymbol y-\boldsymbol p_t). \]

Now define a working response

\[ \boldsymbol u_t=A\theta_t+R_t^{-1}(\boldsymbol y-\boldsymbol p_t), \qquad u_{it}=z_{it}+\frac{y_i-p_{it}}{p_{it}(1-p_{it})}. \]

Then the same equation becomes

\[ A^\top R_t A\theta_{t+1}=A^\top R_t\boldsymbol u_t. \]

These are the normal equations for the weighted least-squares problem

\[ \theta_{t+1} =\underset{\theta}{\operatorname{argmin}}\; \frac{1}{2n}\sum_i r_i\left(u_{it}-\boldsymbol a_i^\top\theta\right)^2. \]

We have not assumed Gaussian binary observations. This weighted quadratic is a local approximation used to optimize the Bernoulli objective. Its responses and weights depend on the current parameters, so they must be recomputed after each step. This is iteratively reweighted least squares (IRLS). For logistic regression, a full IRLS update is algebraically the full Newton update.

10.2.1 What is being approximated?

For one example, expand its negative log likelihood locally around the current score \(z_t\):

\[ \ell(z,y)\approx \ell(z_t,y)+(p_t-y)(z-z_t)+\frac12 r(z-z_t)^2. \]

Complete the square in \(z\). Apart from a term independent of the candidate \(z\), the expression is

\[ \frac12r\left[z-\left(z_t+\frac{y-p_t}{r}\right)\right]^2 =\frac12r(z-u_t)^2. \]

The working response \(u_t\) lies on the score/log-odds scale; it is not a new label or a probability. It may be negative or greater than one. Its weight \(r\) comes from likelihood curvature, not from a claim that some patients matter less.

Math bridge. At initialization \(z=0\), all \(p=1/2\), all \(r=1/4\), and \(u=4y-2\). The first IRLS step is therefore ordinary least squares on responses \(-2\) and \(2\), with identical weights. After that step, probabilities differ and the weights change. This explains why one least-squares solve is not the logistic MLE.

10.3 Improve: verify the same step by two routes

Use square roots of weights to solve the least-squares system. Forming normal equations can square the condition number; lstsq on the weighted design is preferable to explicitly inverting \(A^\top R A\).

theta_current = np.array([-0.2, 0.5, 0.1])
z = A @ theta_current
p = sigmoid(z)
r = p * (1 - p)
assert r.min() > 1e-6  # This demonstration is deliberately away from saturation.
working_response = z + (y - p) / r
sqrt_r = np.sqrt(r)
theta_irls = np.linalg.lstsq(sqrt_r[:, None] * A, sqrt_r * working_response, rcond=None)[0]
theta_newton = theta_current - np.linalg.solve(hessian(A, theta_current), gradient(A, y, theta_current))
np.testing.assert_allclose(theta_irls, theta_newton, atol=1e-12)
print("Maximum coordinate difference:", np.max(np.abs(theta_irls - theta_newton)))
print("Loss before / after one full step:", loss(A, y, theta_current), loss(A, y, theta_irls))
Maximum coordinate difference: 1.7763568394002505e-15
Loss before / after one full step: 0.5108228911207965 0.324071585561428

This verifies one full step. If backtracking selects \(\alpha<1\), use \(\theta_t+\alpha(\theta_{\mathrm{IRLS}}-\theta_t)\) and check the actual logistic loss. Replacing the whole fit by one weighted solve, or keeping the first iteration’s weights forever, changes the algorithm.

When \(p\) is numerically close to zero or one, the working response divides by tiny numbers. The main chapter’s gradient/Hessian implementation avoids explicitly forming that response, though its linear system can still become ill-conditioned. Scaling, stable logits, rank checks, and damping remain useful.

10.4 Why Newton’s trajectory survives a change of coordinates

The chapter found identical loss trajectories for raw and standardized Newton fits. This is more than an accident of one random seed.

Represent an invertible affine change of features, including the intercept, by a nonsingular matrix \(T\) such that \(A_{\mathrm{new}}=AT\). The same score is obtained with \(\theta=T\beta\). By the chain rule,

\[ \boldsymbol g_\beta=T^\top\boldsymbol g_\theta,\qquad H_\beta=T^\top H_\theta T. \]

If \(H_\beta\boldsymbol s_\beta=\boldsymbol g_\beta\), multiply through to see that \(T\boldsymbol s_\beta=\boldsymbol s_\theta\). A full Newton update therefore produces the same scores in either coordinate system. The decrement and Armijo condition are also unchanged if initial scores and damping rules match. Finite-precision arithmetic and nearly singular transformations can spoil that agreement.

Plain GD generally does not have this invariance: mapping its step back multiplies the gradient by \(TT^\top\). This is why changing feature units can change a GD trajectory even when the model family stays the same.

10.5 Why a finite solution is not guaranteed

Here is a separate synthetic, perfectly separable example. The intercept is zero and the one feature’s sign exactly determines the label. Increase a positive weight and inspect the loss.

x_separable = np.array([-2.0, -1.0, 1.0, 2.0])
y_separable = np.array([0.0, 0.0, 1.0, 1.0])
weights = np.linspace(0, 20, 120)
# For every example here, its correctly signed margin is w * abs(x).
# This form avoids cancellation for large positive, correct scores.
separable_losses = np.array([np.logaddexp(0.0, -w * np.abs(x_separable)).mean() for w in weights])
fig, ax = plt.subplots(figsize=(6, 3.5))
ax.semilogy(weights, separable_losses)
ax.set(xlabel="weight", ylabel="mean logistic loss", title="No finite unpenalized MLE")
fig.tight_layout()
show_plot(fig)

Loss versus positive weight on a logarithmic vertical scale; the curve keeps decreasing as the weight grows for a perfectly separable four-point dataset.

Loss versus positive weight on a logarithmic vertical scale; the curve keeps decreasing as the weight grows for a perfectly separable four-point dataset.

Figure 10.1: A separate, perfectly separable toy dataset: increasing the weight drives logistic loss toward zero without attaining it at any finite weight. This is not the overlapping two-feature WDBC training problem.

For any finite positive weight, every assigned probability is still short of perfect certainty, so the loss remains positive. Larger weights always improve it. A report that Newton “converged” merely because an improvement fell below a tolerance does not prove a finite MLE exists.

An L2 penalty changes the problem and can prevent this divergence. For the small two-feature, overlapping training example used in the chapter, the finite fitted parameters have a tiny gradient and a positive-definite Hessian. Do not transfer that observation to all 30 UCI features or to every classification dataset.

10.6 Check your understanding

  1. At zero initialization, derive the working responses for \(y=0\) and \(y=1\).
  2. Which objects must be recomputed after every IRLS step, and why?
  3. Why can we solve each weighted subproblem in closed form while logistic MLE has no general closed-form solution?
  4. Does a curvature weight near zero mean the example can safely be discarded? Consider a confidently wrong prediction and its working response.
  5. What is preserved under an invertible feature transformation: the coefficients, the representable logits, or both?
  6. Give one reason a Hessian might be numerically singular even when the data matrix has full column rank.

Return to the lab and use one of these results to explain your optimizer comparison.