Fisher Scoring as Weighted Least Squares

From a Scoring Step to a Regression Problem

Fisher scoring updates a coefficient vector by moving in a direction determined by the score and the expected information. In a generalized linear model, that step has a useful second interpretation: it is the solution to a weighted least-squares problem. The response and weights in that problem are recalculated from the current fit, so this is an algebraic description of one iteration, not a claim that the original data follow a least-squares model.

We will show the equivalence directly. At the current coefficient vector, the working response adjusts each fitted linear predictor by a scaled residual. The working weight measures how strongly that adjustment contributes to the Fisher scoring equations. Together, they make the weighted least-squares normal equations exactly the Fisher scoring update.

Model and Iteration Quantities

Suppose $Y_1,\ldots,Y_n$ are independent observations from a GLM with mean $\mu_i=E[Y_i]$ and variance $\operatorname{Var}(Y_i)=\phi V(\mu_i)$, where $\phi>0$ is common to the observations. Let $x_i^T$ be row $i$ of the fixed $n\times p$ design matrix $X$, and suppose there is no offset, so $\eta_i=x_i^T\beta$. The link satisfies $\eta_i=g(\mu_i)$.

At a current coefficient vector $\beta$, assume each fitted mean is in the interior of its parameter space, $g$ is differentiable there, and $a_i=d\mu_i/d\eta_i$ exists, is finite, and is nonzero. Write $\mu_i$, $a_i$, and $V_i=V(\mu_i)$ for quantities evaluated at this current fit. Assume also that $0<V_i<\infty$.

Definition: The working response and working weight for observation $i$, evaluated at the current coefficient vector, are $$ z_i=\eta_i+\frac{y_i-\mu_i}{a_i}, \qquad w_i=\frac{a_i^2}{V_i}. $$ Let $z=(z_1,\ldots,z_n)^T$ and $W=\operatorname{diag}(w_1,\ldots,w_n)$. These are iteration-specific quantities: changing $\beta$ generally changes $z$ and $W$.

The weights here omit the common dispersion factor $\phi$. Including $1/\phi$ in every weight would multiply the least-squares objective by the same positive constant and would not change its minimizer. The update below requires $X^TWX$ to be invertible; for example, this holds when $X$ has full column rank and all the working weights are positive.

Why the Fisher Scoring Update Has This Form

For a regular GLM likelihood, the score for the mean parameter of observation $i$ is $(y_i-\mu_i)/\{\phi V(\mu_i)\}$. Applying the chain rule through $\mu_i$ and $\eta_i=x_i^T\beta$ gives the coefficient score $$ U(\beta)=\frac{1}{\phi}X^T\operatorname{diag}\left(\frac{a_i}{V_i}\right)(y-\mu). $$ The expected information for $\beta$ is $$ I(\beta)=\frac{1}{\phi}X^T\operatorname{diag}\left(\frac{a_i^2}{V_i}\right)X =\frac{1}{\phi}X^TWX. $$ Here $y=(y_1,\ldots,y_n)^T$ and $\mu=(\mu_1,\ldots,\mu_n)^T$. Fisher scoring uses the current score and expected information, rather than the observed Hessian, in its update.

Theorem: If $X^TWX$ is invertible, the Fisher scoring update from $\beta$ is the unique solution $\beta_{\mathrm{new}}$ of $$ X^TWX\beta_{\mathrm{new}}=X^TWz. $$ Equivalently, it is the weighted least-squares fit of the working response $z$ on the columns of $X$, using the current working weights $w_i$.

Proof. By definition, the Fisher scoring update is $\beta_{\mathrm{new}}=\beta+I(\beta)^{-1}U(\beta)$. Substituting the score and information expressions, the common factor $1/\phi$ cancels:

\beta_{\mathrm{new}} =\beta+(X^TWX)^{-1}X^T \operatorname{diag}\left(\frac{a_i}{V_i}\right)(y-\mu).

The definition of $z_i$ implies $z_i-\eta_i=(y_i-\mu_i)/a_i$. Since $\eta=X\beta$ and $w_i=a_i^2/V_i$, the $i$th component of $W(z-X\beta)$ is $$ w_i(z_i-x_i^T\beta) =\frac{a_i^2}{V_i}\frac{y_i-\mu_i}{a_i} =\frac{a_i(y_i-\mu_i)}{V_i}. $$ Consequently, $$ X^TW(z-X\beta)=X^T\operatorname{diag}\left(\frac{a_i}{V_i}\right)(y-\mu). $$ Substitute this identity into the scoring update and collect the terms involving $\beta$: $$ \beta_{\mathrm{new}} =\beta+(X^TWX)^{-1}X^TW(z-X\beta) =(X^TWX)^{-1}X^TWz. $$ Multiplying by $X^TWX$ gives the stated normal equations. Invertibility makes their solution unique. $\square$

There is also a direct optimization interpretation. For the $z$ and $W$ computed at the current fit, define the weighted least-squares criterion for a candidate coefficient vector $b$ by $$ Q(b)=(z-Xb)^TW(z-Xb)=\sum_{i=1}^n w_i(z_i-x_i^Tb)^2. $$

Result: If $W$ has positive diagonal entries and $X^TWX$ is positive definite, the unique minimizer of $Q(b)$ is the Fisher scoring update.

Proof. Expanding $Q(b)$ and differentiating with respect to $b$ gives $$ \nabla_b Q(b)=-2X^TW(z-Xb), \qquad \nabla_b^2 Q(b)=2X^TWX. $$ A minimizer must therefore satisfy $X^TWXb=X^TWz$. Positive definiteness of the Hessian makes $Q$ strictly convex, so this equation has exactly one minimizer. The theorem shows that this solution is $\beta_{\mathrm{new}}$. $\square$

This second result clarifies the role of the weights: they define the quadratic criterion being minimized for this one iteration. They are not generally fixed across iterations, and $Q$ is not the original likelihood.

Worked Updates

Worked Example: Normal Responses With a Log Link

Let one observation satisfy $Y\sim N(\mu,\phi)$, with the log link $\eta=\log\mu$, so $\mu=\exp(\eta)$. Use an intercept-only model, take the observed value $y=6$, and start at $\beta=\log 2$. Here $V(\mu)=1$, $\eta=\beta$, and $a=d\mu/d\eta=\mu$.

At the current fit, $\eta=\log 2$, $\mu=2$, and $a=2$. Thus the working weight and response are $$ w=\frac{2^2}{1}=4, \qquad z=\log 2+\frac{6-2}{2}=\log 2+2. $$ For an intercept-only model with one observation, the weighted least-squares equation is $4\beta_{\mathrm{new}}=4z$. Hence $\beta_{\mathrm{new}}=z=\log 2+2$.

Check the scoring calculation directly. The score is $U=(1/\phi)(2/1)(6-2)=8/\phi$, and the expected information is $I=4/\phi$. The increment is $I^{-1}U=2$, giving $\beta_{\mathrm{new}}=\log 2+2$, as before. The fitted mean after the update is $\exp(\log 2+2)$; this is not the observed value 6, nor is it required to equal it.

Worked Example: A Bernoulli Probit Update

Suppose $Y_i$ are independent Bernoulli observations with $\mu_i=\Phi(\eta_i)$, where $\Phi$ is the standard normal distribution function. Use a probit link, the design rows $x_1^T=(1,-1)$ and $x_2^T=(1,1)$, observations $y_1=1,y_2=0$, and current coefficients $\beta=(0,0)^T$.

Both current linear predictors are zero, so $\mu_1=\mu_2=\Phi(0)=1/2$. For a Bernoulli response, $V(\mu)=\mu(1-\mu)$, which is $1/4$ here. Also $a_i=d\Phi(\eta_i)/d\eta_i$ is the standard normal density at zero: $$ a_1=a_2=\frac{1}{\sqrt{2\pi}}. $$ Therefore each working weight is $$ w_i=\frac{(1/\sqrt{2\pi})^2}{1/4}=\frac{2}{\pi}. $$ The working responses are $$ z_1=0+\frac{1-1/2}{1/\sqrt{2\pi}}=\sqrt{\frac{\pi}{2}}, \qquad z_2=0+\frac{0-1/2}{1/\sqrt{2\pi}}=-\sqrt{\frac{\pi}{2}}. $$

The weights are equal, so the weighted least-squares fit is the ordinary least-squares fit to these two working responses. The two equations $b_0-b_1=z_1$ and $b_0+b_1=z_2$ give $$ \beta_{\mathrm{new},0}=\frac{z_1+z_2}{2}=0, \qquad \beta_{\mathrm{new},1}=\frac{z_2-z_1}{2}=-\sqrt{\frac{\pi}{2}}. $$ Thus this scoring step leaves the intercept at zero and gives a negative slope. In the original model, the current fitted probabilities were both $1/2$; the working responses are local quantities used to compute the update, not Bernoulli outcomes.

Worked Example: A Poisson Model With a Square-Root Link

Let $Y\sim\operatorname{Poisson}(\mu)$, so $V(\mu)=\mu$, and use the link $\eta=\sqrt{\mu}$ on the domain $\eta>0$. In an intercept-only model, take $y=8$ and current coefficient $\beta=2$. Then $\eta=2$, $\mu=4$, and $a=d\mu/d\eta=2\eta=4$.

The current working weight and response are $$ w=\frac{4^2}{4}=4, \qquad z=2+\frac{8-4}{4}=3. $$ With a single intercept and observation, the weighted least-squares fit is $\beta_{\mathrm{new}}=z=3$. Directly, the score is $U=(4/4)(8-4)=4$, and the expected information is $I=4^2/4=4$. Thus the scoring increment is $I^{-1}U=1$, which also gives $\beta_{\mathrm{new}}=3$.

What the Equivalence Does—and Does Not—Say

The equivalence is useful because a scoring update can be computed with familiar weighted least-squares machinery. At each iteration, first calculate the fitted means and linear predictors from the current coefficients; then calculate $z$ and $W$; finally solve the weighted normal equations. The fit is a regression of the working response on the design matrix, which includes an intercept only if an intercept column is present in $X$.

Evaluate the current fit.
Compute $\eta_i=x_i^T\beta$, $\mu_i$, $a_i=d\mu_i/d\eta_i$, and $V_i=V(\mu_i)$.

Form the working quantities.
Set $z_i=\eta_i+(y_i-\mu_i)/a_i$ and $w_i=a_i^2/V_i$.

Update the coefficients.
Solve $X^TWX\beta_{\mathrm{new}}=X^TWz$. The resulting coefficients are the Fisher scoring step from the current $\beta$.

A common mistake is to weight by $1/V_i$ alone. That omits the squared derivative $a_i^2=(d\mu_i/d\eta_i)^2$, which is part of the expected information and is essential for the normal equations to match the score update. Another mistake is to treat $z_i$ as a newly observed response. It is instead constructed from the current residual and local link derivative; once the coefficients change, both the working response and weights are generally recomputed.

The weighted least-squares interpretation also depends on having a solvable system. If the design is rank deficient, or if a weight is zero, $X^TWX$ may fail to be invertible; in that case the update need not be unique, and the inverse-based formula does not apply as written. Under the stated positive-weight and full-rank conditions, the equivalence is exact, not an approximation to the Fisher scoring step.

Check Your Understanding

Use the definitions and equations above to answer the following questions.

Why does the common dispersion factor $\phi$ cancel from the Fisher scoring coefficient update?
For a given observation, explain how $w_i(z_i-\eta_i)$ reduces to $a_i(y_i-\mu_i)/V_i$.
What condition on $X^TWX$ ensures that the weighted least-squares update is unique?
Why is the working response not generally an observed value of the original response variable?
What factor is missing if someone uses $1/V_i$ as the working weight?