Fisher Scoring and Iteratively Reweighted Least Squares

From a Likelihood Score to a Regression Update

In a generalized linear model, the coefficient vector determines each observation's linear predictor, while the link function translates that predictor into a mean. The likelihood score therefore changes with the current coefficient estimate. Fisher scoring uses the expected information at that estimate to turn the score into a coefficient update. For generalized linear models, the update has a useful form: it is a weighted least-squares fit to a constructed response, with both the response and the weights recalculated at each iteration.

The phrase iteratively reweighted least squares (IRLS) describes this computational form. It does not mean that the original observations are repeatedly given new least-squares errors. Instead, each step constructs a local weighted least-squares problem whose solution is exactly one Fisher-scoring update.

Setup: Suppose \(Y_1,\ldots,Y_n\) are independent observations from a regular one-parameter exponential-dispersion model. Let \(E[Y_i]=\mu_i\), \(\operatorname{Var}(Y_i)=\phi V(\mu_i)\), and \(\eta_i=g(\mu_i)=\mathbf{x}_i^T\boldsymbol\beta\). Treat the dispersion parameter \(\phi>0\) as fixed. Assume the means are in the interior of their parameter space, \(g\) is differentiable with nonzero derivative there, \(V(\mu_i)>0\), and the likelihood is differentiable. Write \(\mathbf{X}\) for the matrix whose \(i\)th row is \(\mathbf{x}_i^T\).

For the usual exponential-dispersion likelihood, the derivative of the \(i\)th log-likelihood contribution with respect to \(\mu_i\) is \((y_i-\mu_i)/\{\phi V(\mu_i)\}\). Define \(q_i=d\mu_i/d\eta_i\). Applying the chain rule, the coefficient score is

\mathbf{U}(\boldsymbol\beta) = \sum_{i=1}^n \mathbf{x}_i \frac{(y_i-\mu_i)q_i}{\phi V(\mu_i)}.

All quantities on the right, including \(\mu_i\) and \(q_i\), depend on the current coefficient vector. The expected information uses the model variance and is

\mathcal{I}(\boldsymbol\beta) = \mathbf{X}^T\mathbf{W}\mathbf{X}, \qquad \mathbf{W}=\operatorname{diag}(w_1,\ldots,w_n), \qquad w_i=\frac{q_i^2}{\phi V(\mu_i)}.

When this information matrix is nonsingular, a Fisher-scoring step from \(\boldsymbol\beta^{(t)}\) is

\boldsymbol\beta^{(t+1)} = \boldsymbol\beta^{(t)} + \mathcal{I}(\boldsymbol\beta^{(t)})^{-1} \mathbf{U}(\boldsymbol\beta^{(t)}).

Why Expected Information Has This Form

The weights are not arbitrary numerical devices. They come from the variability of the score. The following derivation establishes the information formula used in the update.

Result: Under the setup above, the covariance matrix of the score equals the expected information, and that matrix is \(\mathbf{X}^T\mathbf{W}\mathbf{X}\), with \(w_i=q_i^2/\{\phi V(\mu_i)\}\).

Proof. The score contribution from observation \(i\) is \(\mathbf{u}_i=\mathbf{x}_i (Y_i-\mu_i)q_i/\{\phi V(\mu_i)\}\). Since \(E[Y_i-\mu_i]=0\), \(E[\mathbf{u}_i]=\mathbf{0}\). Its covariance is

\begin{aligned} \operatorname{Cov}(\mathbf{u}_i) &= \mathbf{x}_i\mathbf{x}_i^T \frac{q_i^2}{\phi^2 V(\mu_i)^2} \operatorname{Var}(Y_i)\\ &= \mathbf{x}_i\mathbf{x}_i^T \frac{q_i^2}{\phi^2 V(\mu_i)^2} \phi V(\mu_i)\\ &= \mathbf{x}_i\mathbf{x}_i^T \frac{q_i^2}{\phi V(\mu_i)}. \end{aligned}

Independence makes the covariance of the total score the sum of these contribution covariances. Summing gives \(\sum_i w_i\mathbf{x}_i\mathbf{x}_i^T=\mathbf{X}^T\mathbf{W}\mathbf{X}\), as claimed. \(\square\)

Deriving the Working Response

At iteration \(t\), calculate the current linear predictors \(\eta_i^{(t)}=\mathbf{x}_i^T\boldsymbol\beta^{(t)}\), means \(\mu_i^{(t)}\), and derivatives \(q_i^{(t)}\). Define a working response \(z_i^{(t)}\) by

z_i^{(t)} = \eta_i^{(t)} + \frac{y_i-\mu_i^{(t)}}{q_i^{(t)}}, \qquad w_i^{(t)} = \frac{(q_i^{(t)})^2}{\phi V(\mu_i^{(t)})}.

The working response is the current linear predictor plus a first-order correction based on the observed residual. It is generally not an observed response on the original measurement scale. The key identity is

w_i^{(t)}\bigl(z_i^{(t)}-\eta_i^{(t)}\bigr) = \frac{q_i^{(t)}(y_i-\mu_i^{(t)})}{\phi V(\mu_i^{(t)})}.

The right side is the scalar multiplying \(\mathbf{x}_i\) in the score. Consequently, the weighted least-squares normal equations for fitting \(z_i^{(t)}\) to \(\mathbf{x}_i^T\boldsymbol\beta\) have exactly the score-correction term as their right-hand side.

Theorem: If \(\mathbf{X}^T\mathbf{W}^{(t)}\mathbf{X}\) is nonsingular, the Fisher-scoring update is the unique minimizer over \(\boldsymbol\beta\) of \(\sum_{i=1}^n w_i^{(t)}(z_i^{(t)}-\mathbf{x}_i^T\boldsymbol\beta)^2\).

Proof. Differentiating this weighted sum of squares with respect to \(\boldsymbol\beta\), and setting the gradient to zero, gives \(\mathbf{X}^T\mathbf{W}^{(t)}\mathbf{X}\boldsymbol\beta=\mathbf{X}^T\mathbf{W}^{(t)}\mathbf{z}^{(t)}\). Because the matrix is nonsingular, the unique solution is \[ \boldsymbol\beta_{\mathrm{WLS}} = (\mathbf{X}^T\mathbf{W}^{(t)}\mathbf{X})^{-1} \mathbf{X}^T\mathbf{W}^{(t)}\mathbf{z}^{(t)}. \] Now use \(\eta_i^{(t)}=\mathbf{x}_i^T\boldsymbol\beta^{(t)}\) and the identity above:

\begin{aligned} \mathbf{X}^T\mathbf{W}^{(t)}\mathbf{z}^{(t)} &= \mathbf{X}^T\mathbf{W}^{(t)}\mathbf{X}\boldsymbol\beta^{(t)} + \mathbf{U}(\boldsymbol\beta^{(t)}). \end{aligned}

Substitution into the weighted least-squares solution yields \(\boldsymbol\beta_{\mathrm{WLS}}=\boldsymbol\beta^{(t)}+ (\mathbf{X}^T\mathbf{W}^{(t)}\mathbf{X})^{-1}\mathbf{U}(\boldsymbol\beta^{(t)})\), which is exactly the Fisher-scoring update. Positive weights and nonsingularity make the weighted sum strictly convex, so this solution is its unique minimizer. \(\square\)

Worked Examples

Worked Example: A Poisson Intercept-Only Update

Take independent Poisson observations \(y_1=1\) and \(y_2=3\), with an intercept-only model and log link. Thus \(\eta_i=\beta\), \(\mu_i=e^\beta\), \(V(\mu_i)=\mu_i\), \(\phi=1\), and \(q_i=d\mu_i/d\eta_i=\mu_i\). Start at \(\beta^{(0)}=0\).

The starting means are \(\mu_1^{(0)}=\mu_2^{(0)}=1\), and both weights are \(w_i^{(0)}=1^2/1=1\). The score is \(U(\beta^{(0)})=(1-1)+(3-1)=2\), and the information is \(1+1=2\). Therefore the Fisher-scoring update is \(\beta^{(1)}=0+2/2=1\).

The working responses are \(z_1^{(0)}=0+(1-1)/1=0\) and \(z_2^{(0)}=0+(3-1)/1=2\). Their weighted least-squares fit with an intercept is their mean, \((0+2)/2=1\), matching the coefficient update. This is one scoring step, not the exact maximum likelihood estimate: the latter is \(\log\{(1+3)/2\}=\log 2\), which differs from \(1\).

Worked Example: A Bernoulli Logistic Update

Suppose three independent Bernoulli outcomes are \(0,1,1\), and use an intercept-only logistic model. At \(\beta^{(0)}=0\), each fitted probability is \(p_i^{(0)}=1/2\). For this model, \(\mu_i=p_i\), \(V(\mu_i)=p_i(1-p_i)\), \(\phi=1\), and \(q_i=p_i(1-p_i)\).

Each weight is \(w_i^{(0)}=(1/4)^2/(1/4)=1/4\). The score is \(U(\beta^{(0)})=\sum_i(y_i-p_i^{(0)})=-1/2+1/2+1/2=1/2\); the information is \(3/4\). Hence \(\beta^{(1)}=0+(1/2)/(3/4)=2/3\).

Since \(q_i^{(0)}=1/4\), the working responses are \(z_i^{(0)}=0+4(y_i-1/2)\), giving \(-2,2,2\). Their weighted mean, with equal weights \(1/4\), is \((-2+2+2)/3=2/3\). Thus the IRLS fit again reproduces the Fisher-scoring step. The updated probability is \(1/(1+e^{-2/3})\), approximately \(0.661\).

Worked Example: Normal Regression as a Least-Squares Step

Let \(Y_i\sim N(\mu_i,1)\) independently, with identity link and design rows \((1,0),(1,1),(1,2)\). The observed responses are \(1,2,2\). Here \(q_i=1\), \(V(\mu_i)=1\), and \(\phi=1\), so every weight is \(1\). The working response is \(z_i=\eta_i+(y_i-\mu_i)=y_i\), regardless of the starting coefficient vector.

The weighted least-squares normal equations are therefore \(\mathbf{X}^T\mathbf{X}\boldsymbol\beta=\mathbf{X}^T\mathbf{y}\), with \(\mathbf{X}^T\mathbf{X}=\begin{pmatrix}3&3\\3&5\end{pmatrix}\) and \(\mathbf{X}^T\mathbf{y}=(5,6)^T\).

The determinant is \(3(5)-3(3)=6\), so the solution is \[ \boldsymbol\beta = \frac{1}{6} \begin{pmatrix}5&-3\\-3&3\end{pmatrix} \begin{pmatrix}5\\6\end{pmatrix} = \begin{pmatrix}7/6\\1/2\end{pmatrix}. \] Thus an IRLS step returns intercept \(7/6\) and slope \(1/2\). In this normal identity-link case, this is the least-squares estimate and does not depend on the starting point.

Carrying Out IRLS and Avoiding Misinterpretation

A practical iteration consists of calculating the current means and derivatives, constructing the working response and weights, solving the weighted least-squares problem, and repeating with the updated coefficients. Each iteration solves a different problem because the weights and working response generally change with \(\boldsymbol\beta^{(t)}\).

Evaluate the current fit.
Compute \(\eta_i^{(t)}=\mathbf{x}_i^T\boldsymbol\beta^{(t)}\), \(\mu_i^{(t)}\), and \(q_i^{(t)}=d\mu_i/d\eta_i\).

Construct the working quantities.
Set \(z_i^{(t)}=\eta_i^{(t)}+(y_i-\mu_i^{(t)})/q_i^{(t)}\) and \(w_i^{(t)}=(q_i^{(t)})^2/\{\phi V(\mu_i^{(t)})\}\).

Fit by weighted least squares.
Solve \((\mathbf{X}^T\mathbf{W}^{(t)}\mathbf{X})\boldsymbol\beta^{(t+1)}=\mathbf{X}^T\mathbf{W}^{(t)}\mathbf{z}^{(t)}\).

Recalculate and assess.
Use the new coefficients to rebuild the means, weights, and working responses before the next iteration; stop according to a specified numerical convergence criterion.

Common mistake: treating the working response as fixed data. A tempting but incorrect argument is that IRLS creates a new response \(z_i\), after which one performs ordinary least squares once and regards the result as the final fit. The correction is that \(z_i^{(t)}\) is defined using the current fitted mean and derivative, and \(w_i^{(t)}\) also depends on that fit. They must be updated at each iteration unless the model makes them constant, as in normal regression with identity link.

Another important distinction is that Fisher scoring uses expected information, not the observed negative Hessian. The weighted least-squares identity follows from that expected-information choice. If \(\mathbf{X}^T\mathbf{W}^{(t)}\mathbf{X}\) is singular, the unique-update formula does not apply; coefficients may be unidentified or the local weighted design may lack full rank. If the matrix is nonsingular, the algebra establishes the form of the step, but it does not by itself guarantee that successive steps converge to a likelihood maximum.

IRLS matters because it connects likelihood-based estimation to a familiar linear-algebra problem while preserving the model's mean-variance relationship through changing weights. The working response is a local device for computing a score update; interpreting it as a new observed outcome, or assuming every iteration is ordinary least squares with fixed weights, loses that distinction.

Check Your Understanding

Use the definitions and derivations above to answer these questions.

For a model with \(q_i=d\mu_i/d\eta_i\), write the Fisher-scoring weight and working response.
Explain why the covariance of the score contribution from observation \(i\) is \(w_i\mathbf{x}_i\mathbf{x}_i^T\).
Show how the weighted least-squares normal equations yield the Fisher-scoring coefficient update.
In a Bernoulli logistic model, what is \(q_i\) in terms of \(p_i\), and why do the IRLS weights change between iterations?
What condition permits the unique weighted least-squares update written in this tutorial, and what does failure of that condition imply for that formula?