Newton-Raphson Coefficient Updates in GLMs

From the Score Equation to a Coefficient Update

A fitted generalized linear model is often found by solving the score equation $U(\beta)=0$, where $\beta$ is the vector of regression coefficients and $U(\beta)$ is the gradient of the log-likelihood. In many models this equation cannot be solved directly. Newton–Raphson replaces it, near the current coefficient vector, by a linear approximation that can be solved for a step.

The sign in this calculation matters. Let $H(\beta)$ denote the Hessian matrix of the log-likelihood, the matrix of second derivatives. When $H(\beta)$ is nonsingular, the Newton–Raphson update for solving the score equation is $\beta-H(\beta)^{-1}U(\beta)$. Equivalently, if $J(\beta)=-H(\beta)$ is the observed information, the update is $\beta+J(\beta)^{-1}U(\beta)$. These forms are the same update written with different matrices.

Definition: At a current coefficient vector $\beta$, define the Newton–Raphson step $\delta$ by solving $H(\beta)\delta=-U(\beta)$. If $H(\beta)$ is nonsingular, the next iterate is $$ \beta_{\mathrm{new}}=\beta+\delta=\beta-H(\beta)^{-1}U(\beta). $$ Writing $J(\beta)=-H(\beta)$, the same update is $\beta_{\mathrm{new}}=\beta+J(\beta)^{-1}U(\beta)$. A unique ordinary Newton step requires the relevant matrix to be nonsingular.

For a GLM with independent observations and fixed design, write the linear predictor for observation $i$ as $\eta_i=x_i^\mathsf{T}\beta+o_i$, where $x_i$ is its design vector and $o_i$ is a fixed offset. Let $\ell_i(\eta_i)$ be its log-likelihood contribution, with any other model quantities held fixed. Define $u_i(\eta_i)=\partial\ell_i/\partial\eta_i$ and $h_i(\eta_i)=\partial^2\ell_i/\partial\eta_i^2$. The chain rule gives the following useful expressions.

Result: If the log-likelihood is $\ell(\beta)=\sum_{i=1}^n\ell_i(\eta_i)$, and each contribution is twice differentiable in $\eta_i$, then $$ U(\beta)=\sum_{i=1}^n x_i u_i(\eta_i), \qquad H(\beta)=\sum_{i=1}^n x_i x_i^\mathsf{T}h_i(\eta_i). $$ Equivalently, with $X$ the design matrix, $u$ the vector of $u_i$'s and $D$ the diagonal matrix with entries $h_i$, $U=X^\mathsf{T}u$ and $H=X^\mathsf{T}DX$.

To verify these expressions, differentiate each contribution with respect to $\beta$. Since $\partial\eta_i/\partial\beta=x_i$, the gradient of $\ell_i(\eta_i)$ is $x_i u_i(\eta_i)$. Differentiating that gradient again gives $x_i x_i^\mathsf{T}h_i(\eta_i)$, because $x_i$ is fixed. Summing the contributions yields the stated score and Hessian. This derivation uses the observed second derivatives at the current coefficient values; it does not replace them with an expectation.

Why the Newton Formula Has This Sign

The Newton step comes from applying Newton’s method to the equation $U(\beta)=0$. At the current value $\beta$, approximate the score at $\beta+\delta$ by its first-order expansion, $U(\beta+\delta)\approx U(\beta)+H(\beta)\delta$. Setting this linear approximation equal to zero gives $H(\beta)\delta=-U(\beta)$. If $H(\beta)$ is nonsingular, solving this equation gives the update in the definition.

The same step can be understood through the local quadratic approximation to the log-likelihood. Define $$ q(\delta)=\ell(\beta)+U(\beta)^\mathsf{T}\delta+\frac{1}{2}\delta^\mathsf{T}H(\beta)\delta. $$ Its gradient with respect to $\delta$ is $U(\beta)+H(\beta)\delta$. Thus its stationary point solves the Newton step equation. If, in addition, $H(\beta)$ is negative definite, this quadratic has a unique global maximum at that step: for any nonzero vector $v$, its quadratic curvature in direction $v$ is $v^\mathsf{T}H(\beta)v<0$. This statement concerns the quadratic approximation, not a guarantee that the actual log-likelihood increases after the update.

Derivation: If $H(\beta)$ is nonsingular, the linearized score equation has the unique solution $\delta=-H(\beta)^{-1}U(\beta)$. Since $J(\beta)=-H(\beta)$, this is also $\delta=J(\beta)^{-1}U(\beta)$. If $H(\beta)$ is negative definite, the quadratic $q$ has negative-definite Hessian, so its unique stationary point is its unique global maximum. Therefore the Newton update both solves the linearized score equation and maximizes this local quadratic approximation under the stated curvature condition. $\square$

Worked Coefficient Updates

Worked Example: One Newton Step for an Intercept-Only Logistic Model

Suppose four independent Bernoulli observations contain three successes and one failure. Fit an intercept-only logistic model, so $p=\exp(\beta)/(1+\exp(\beta))$. Apart from terms not involving $\beta$, the log-likelihood is $\ell(\beta)=3\beta-4\log(1+\exp(\beta))$. Start at $\beta^{(0)}=0$.

The probability at the starting value is $p=1/2$. Differentiating the log-likelihood gives the score $U(\beta)=3-4p$, and differentiating again gives the Hessian $H(\beta)=-4p(1-p)$. At zero, $U(0)=3-4(1/2)=1$ and $H(0)=-4(1/2)(1/2)=-1$.

The step is $\delta=-H(0)^{-1}U(0)=-(-1)^{-1}(1)=1$. Therefore $\beta^{(1)}=0+1=1$. The updated probability is $p=\exp(1)/(1+\exp(1))\approx0.7311$. This is a single Newton step, not the exact maximum-likelihood estimate.

For a check on the update, the score at the new value is $3-4\exp(1)/(1+\exp(1))\approx0.0758$, not zero. The step has reduced the score’s magnitude from $1$ to about $0.0758$, but the score equation still needs to be solved more closely.

Worked Example: A Second Logistic Step and the Sign Convention

Continue the same model from $\beta^{(1)}=1$. At this value $p=\exp(1)/(1+\exp(1))\approx0.7310586$, so $$ U(1)=3-\frac{4\exp(1)}{1+\exp(1)}\approx0.0757657, \qquad H(1)=-4p(1-p)\approx-0.7864477. $$

Using the Hessian form, the increment is $\delta=-U(1)/H(1)\approx-0.0757657/(-0.7864477)\approx0.09634$. Thus $\beta^{(2)}\approx1.09634$. Using the observed-information form gives the same calculation: $J(1)=0.7864477$, so $\delta=U(1)/J(1)\approx0.09634$.

The positive increment is consistent with the positive score and negative Hessian. The exact logistic maximum-likelihood estimate here is $\log(3)\approx1.09861$, obtained by setting the fitted probability equal to the observed success proportion $3/4$. The second iterate is close to, but not exactly, that value.

Worked Example: One Poisson Log-Link Update

Consider two independent Poisson observations with a common mean $\lambda=\exp(\beta)$, so this is an intercept-only Poisson model with a log link. The observed counts are $1$ and $3$. Ignoring terms constant in $\beta$, the log-likelihood is $\ell(\beta)=4\beta-2\exp(\beta)$. Start at $\beta^{(0)}=0$.

Differentiating gives $U(\beta)=4-2\exp(\beta)$ and $H(\beta)=-2\exp(\beta)$. At zero, $U(0)=4-2=2$ and $H(0)=-2$. Hence $\delta=-2/(-2)=1$, and the first updated coefficient is $\beta^{(1)}=1$, corresponding to fitted mean $\exp(1)\approx2.7183$.

At $\beta=1$, the score is $4-2\exp(1)\approx-1.4366$, and the Hessian is $-2\exp(1)\approx-5.4366$. The next increment is $$ \delta=-\frac{4-2\exp(1)}{-2\exp(1)} =\frac{4-2\exp(1)}{2\exp(1)} \approx-0.26424. $$ Thus the next coefficient is approximately $1-0.26424=0.73576$. The exact maximum-likelihood estimate is $\log(2)\approx0.69315$, since the fitted common mean at the maximum is the sample average $(1+3)/2=2$. The Newton iterates move toward that solution in these calculations.

Worked Example: A Two-Coefficient Grouped Logistic Step

Two groups have covariate values $x=0$ and $x=1$. Their numbers of trials are $m_1=2$ and $m_2=3$, and their success counts are $y_1=1$ and $y_2=2$. Fit a logistic model with $\eta_i=\beta_0+\beta_1x_i$, and start at $\beta=(0,0)^\mathsf{T}$. Both fitted probabilities are $1/2$.

For grouped binomial data, the score contribution from group $i$ is $x_i(y_i-m_ip_i)$, with the intercept included in $x_i$. The residuals $y_i-m_ip_i$ are $1-2(1/2)=0$ and $2-3(1/2)=1/2$. Therefore $$ U=\begin{pmatrix}0+1/2\\0(0)+1(1/2)\end{pmatrix} =\begin{pmatrix}1/2\\1/2\end{pmatrix}. $$

The observed information for a logistic model at these probabilities is $X^\mathsf{T}WX$, where the group weights are $m_ip_i(1-p_i)$. They are $2(1/2)(1/2)=1/2$ and $3(1/2)(1/2)=3/4$. With design rows $(1,0)$ and $(1,1)$, this gives $$ J=\begin{pmatrix}1/2+3/4&3/4\\3/4&3/4\end{pmatrix} =\begin{pmatrix}5/4&3/4\\3/4&3/4\end{pmatrix}. $$

Solve $J\delta=U$. The equations are $(5/4)\delta_0+(3/4)\delta_1=1/2$ and $(3/4)\delta_0+(3/4)\delta_1=1/2$. Subtracting the second equation from the first gives $(1/2)\delta_0=0$, so $\delta_0=0$. Substitution gives $(3/4)\delta_1=1/2$, hence $\delta_1=2/3$. The updated coefficient vector is $(0,2/3)^\mathsf{T}$.

What the Update Does—and Common Errors

A Newton update uses the current score and the local second derivatives to choose a coefficient step. It is not simply a step in the score’s direction: the Hessian scales and mixes the coefficient components. In several dimensions, computing the step means solving a linear system such as $J\delta=U$; explicitly forming a matrix inverse is usually unnecessary.

Reversing the sign. With $H$ defined as the log-likelihood Hessian, the correct update is $\beta-H^{-1}U$. The expression $\beta+H^{-1}U$ has the wrong sign. With $J=-H$, the correct equivalent is $\beta+J^{-1}U$.
Mixing up the Hessian and information. The Hessian $H$ and observed information $J$ have opposite signs. Substituting one for the other without changing the update sign produces an incorrect step.
Assuming every update solves the score equation. Newton’s method solves a linear approximation to the score equation at the current value. Except in special cases, the new score is not exactly zero.
Assuming a step must increase the likelihood. The local quadratic calculation alone does not guarantee that the actual log-likelihood rises after a full step. Whether and how to modify a step for acceptance is a separate issue.
Ignoring a singular Hessian. If $H$ is singular, the ordinary Newton formula does not define a unique update. A numerical solver cannot make the underlying linear system uniquely solvable merely by changing notation.

Key takeaway: Evaluate the score and Hessian at the current coefficients, solve $H\delta=-U$ (or $J\delta=U$), and add the resulting increment. Keep the sign tied to the matrix definition, and interpret the result as a local Newton step—not automatically as a completed fit or a guaranteed likelihood increase.

Check Your Understanding

Use the score and curvature at the current coefficients to reason through each question.

If $U(\beta)$ is positive in a one-parameter model and $H(\beta)$ is negative, what sign does the Newton increment have?
Write the Newton update once in terms of $H$ and once in terms of $J=-H$.
For an independent-observation GLM with fixed design, explain why the Hessian has the form $X^\mathsf{T}DX$ when each contribution depends on coefficients through its linear predictor.
In the grouped logistic example, why is the intercept increment zero even though the score vector is nonzero?
Does solving the quadratic approximation’s maximization problem prove that the actual log-likelihood increases after the full Newton step? Explain.