Tracing Successive Logistic IRLS Iterations

Introduction

A single IRLS update can be written as a weighted least-squares calculation, but an iteration is more than solving that one system. Once the coefficients change, the linear predictors change; the fitted probabilities, working responses, and weights must all be recalculated from the new predictors. Tracing these quantities in order makes the algorithm concrete and helps catch arithmetic errors.

This tutorial follows successive updates for logistic regression. The examples use binomial observations, including Bernoulli data as the special case with one trial per observation. We use no offsets. At each step, all working quantities are computed from the current coefficient vector, and the next coefficient vector is obtained by solving the resulting weighted least-squares problem.

Quantities Recomputed at Every Iteration

Definition: For binomial observation $i$, let $m_i$ be the number of trials, $y_i$ the observed proportion of successes, and $x_i^\mathsf{T}$ the $i$th row of the design matrix $X$. At coefficient vector $\beta^{(k)}$, the logistic linear predictor and fitted probability are $$ \eta_i^{(k)}=x_i^\mathsf{T}\beta^{(k)}, \qquad p_i^{(k)}=\frac{1}{1+\exp(-\eta_i^{(k)})}. $$ The working response and working weight are $$ z_i^{(k)}=\eta_i^{(k)}+ \frac{y_i-p_i^{(k)}}{p_i^{(k)}(1-p_i^{(k)})}, \qquad w_i^{(k)}=m_i p_i^{(k)}(1-p_i^{(k)}). $$ The next vector $\beta^{(k+1)}$ minimizes $\sum_i w_i^{(k)}(z_i^{(k)}-x_i^\mathsf{T}b)^2$ over coefficient vectors $b$.

These expressions assume finite predictors, so each fitted probability lies strictly between zero and one, and positive trial counts. If $X$ has full column rank, the working weights are positive and the weighted least-squares solution is unique. The key bookkeeping rule is that the superscript $k$ applies to every working quantity: do not update the coefficients while leaving probabilities, responses, or weights at their previous values.

Result: Let $W^{(k)}$ be the diagonal matrix with entries $w_i^{(k)}$, and let $z^{(k)}$ collect the working responses. If $X$ has full column rank, the next IRLS vector is $$ \beta^{(k+1)} = \left(X^\mathsf{T}W^{(k)}X\right)^{-1} X^\mathsf{T}W^{(k)}z^{(k)}. $$

To verify this expression, differentiate the weighted sum of squares with respect to $b$. Its gradient is $-2X^\mathsf{T}W^{(k)}(z^{(k)}-Xb)$. Setting that gradient equal to zero gives the weighted normal equations $X^\mathsf{T}W^{(k)}Xb=X^\mathsf{T}W^{(k)}z^{(k)}$. Full column rank and positive weights make the matrix on the left invertible, giving the stated update. This explains why each iteration has two stages: first form the working quantities at the current vector, then solve the weighted system for the next vector. $\square$

Worked Trace: Two Symmetric Binomial Groups

Consider two groups, each with four trials. The first group has $x=-1$ and one success, so its observed proportion is $y_1=1/4$. The second has $x=1$ and three successes, so $y_2=3/4$. The model includes an intercept and a slope, with design rows $(1,-1)$ and $(1,1)$. Start at $\beta^{(0)}=(0,0)^\mathsf{T}$.

Worked Example: Three Successive Logistic Updates

Iteration 0: calculate the first working system. Both predictors are zero, so $p_1^{(0)}=p_2^{(0)}=0.5$. The weights are $4(0.5)(0.5)=1$ for each group. The working responses are $$ z_1^{(0)}=0+\frac{0.25-0.5}{0.25}=-1,\qquad z_2^{(0)}=0+\frac{0.75-0.5}{0.25}=1. $$ With equal weights, the fitted working values at $x=-1$ and $x=1$ are $-1$ and $1$. Thus the intercept is zero and the slope is one: $\beta^{(1)}=(0,1)^\mathsf{T}$.

Iteration 1: recompute at $\beta^{(1)}$. The predictors are $-1$ and $1$. The fitted probabilities are $p_1^{(1)}\approx0.268941421$ and $p_2^{(1)}\approx0.731058579$. Both variance factors $p_i^{(1)}(1-p_i^{(1)})$ are approximately $0.196611933$, so both weights are approximately $4(0.196611933)=0.786447733$. The working responses are $$ z_1^{(1)} =-1+\frac{0.25-0.268941421}{0.196611933} \approx -1.0963391, \qquad z_2^{(1)} =1+\frac{0.75-0.731058579}{0.196611933} \approx 1.0963391. $$ The weights are equal, and the two fitted working values are opposite. The new intercept is zero and the new slope is approximately $1.0963391$, giving $\beta^{(2)}\approx(0,1.0963391)^\mathsf{T}$.

Iteration 2: recompute again, without rounding the preceding calculations internally. The predictors are approximately $-1.0963391$ and $1.0963391$. For the positive predictor, the recomputed probability and variance factor are $$ p_2^{(2)}\approx0.749573535,\qquad p_2^{(2)}(1-p_2^{(2)})\approx0.187713051. $$ The negative-predictor probability is the complement, so its variance factor is the same. The working responses are therefore $$ z_2^{(2)} \approx 1.0963391+ \frac{0.75-0.749573535}{0.187713051} \approx 1.0986110,\qquad z_1^{(2)}\approx-1.0986110. $$ Equal weights and opposite working responses again give a zero intercept. Thus $\beta^{(3)}\approx(0,1.0986110)^\mathsf{T}$.

One further update. The target probability $0.75$ has logit $\log(0.75/0.25)=\log(3)\approx1.0986123$. At a predictor of about $1.0986110$, the fitted probability is about $0.74999976$, so the next working response is about $1.0986123$. The coefficients move to approximately $\beta^{(4)}=(0,1.0986123)^\mathsf{T}$. This numerical approach to the target illustrates the iteration; it is not a general guarantee of convergence for every data set or starting value.

Why the Symmetry Simplifies the Trace

In the preceding example, symmetry is useful not just for shortening arithmetic but also for checking it. If the coefficient vector has zero intercept and slope $t$, the two predictors are $-t$ and $t$. Their fitted probabilities add to one. Since the observed proportions are $1/4$ and $3/4$, their working responses are opposites, and their weights are equal.

Result: In the two-group example above, if an iteration starts at $\beta=(0,t)^\mathsf{T}$ for finite $t$, its next weighted least-squares update also has zero intercept. The new slope is the working response for the group at $x=1$.

Proof. At $(0,t)^\mathsf{T}$, the predictors are $-t$ and $t$. Write $p=1/(1+\exp(-t))$; the two probabilities are $1-p$ and $p$. Their variance factors are both $p(1-p)$, and the group sizes are equal, so their weights are equal. Substitution into the working-response formula, using $y_1=1/4$ and $y_2=3/4$, gives $z_1=-z_2$. The two design rows are $(1,-1)$ and $(1,1)$. With equal positive weights, their weighted least-squares fit matches the two working responses: the intercept is $(z_1+z_2)/2=0$, and the slope is $(z_2-z_1)/2=z_2$. The design matrix has full rank, so this solution is unique. $\square$

Start with the current coefficients.
Use $\beta^{(k)}$ to calculate every $\eta_i^{(k)}$ and $p_i^{(k)}$.

Build the working system.
Use those probabilities to calculate $z_i^{(k)}$ and $w_i^{(k)}$, then form the weighted least-squares problem.

Update the coefficients.
Solve the weighted normal equations to obtain $\beta^{(k+1)}$. Only then begin the next iteration.

Worked Trace: An Intercept-Only Model

A simpler design shows that the same recomputation rule applies even when there is only one coefficient. Suppose there are ten trials with six successes, so $y=0.6$, and the model has an intercept only. Start at $\beta^{(0)}=0$. The design value for every trial is one, and the grouped weight includes all ten trials.

Worked Example: Updating a Binomial Proportion

At $\beta^{(0)}=0$, the predictor is zero and $p^{(0)}=0.5$. The working response and weight are $$ z^{(0)}=0+\frac{0.6-0.5}{0.5(0.5)}=0.4,\qquad w^{(0)}=10(0.5)(0.5)=2.5. $$ With an intercept-only design, the weighted least-squares estimate is the working response itself. Hence $\beta^{(1)}=0.4$.

At this new coefficient, $p^{(1)}=1/(1+\exp(-0.4))\approx0.598687660$, and $p^{(1)}(1-p^{(1)})\approx0.240260746$. The new weight is approximately $10(0.240260746)=2.402607457$. Recalculating the working response gives $$ z^{(1)} =0.4+\frac{0.6-0.598687660}{0.240260746} \approx0.40546215. $$ Therefore $\beta^{(2)}\approx0.40546215$. Recomputing once more gives $\beta^{(3)}\approx0.40546511$, close to $\log(0.6/0.4)\approx0.40546511$. The new weight does not alter the fitted value in this one-coefficient example, but it is still part of the working system and would matter when several design rows have different predictors.

Worked Trace: Changing the Design Coordinates

The coefficient values depend on how predictors are encoded, even when the fitted linear predictors are the same. Return to the two groups with proportions $1/4$ and $3/4$, but use predictor values $x=0$ and $x=1$, with an intercept. This example makes the effect of the design rows visible in the coefficient update.

Worked Example: From Working Values to Intercept and Slope

Take four trials in each group and start at $\beta^{(0)}=(0,0)^\mathsf{T}$. Both predictors are zero, both probabilities are $0.5$, and both weights are $4(0.25)=1$. The working responses are $-1$ for the group at $x=0$ and $1$ for the group at $x=1$, by the same calculations $(0.25-0.5)/0.25=-1$ and $(0.75-0.5)/0.25=1$. The two working values are fitted exactly by intercept $-1$ and slope $2$, so $\beta^{(1)}=(-1,2)^\mathsf{T}$.

At this vector, the predictors are $-1$ and $1$. The probabilities are approximately $0.268941421$ and $0.731058579$; the variance factors are both $0.196611933$, and hence the weights are equal at approximately $0.786447733$. The working responses are approximately $-1.0963391$ and $1.0963391$. With design rows $(1,0)$ and $(1,1)$, the fitted working values equal those responses when the intercept is $-1.0963391$ and the slope is their difference, $2.1926782$. Thus $\beta^{(2)}\approx(-1.0963391,2.1926782)^\mathsf{T}$.

The predictor values at this updated vector are again approximately $-1.0963391$ and $1.0963391$, just as in the centered coding. The coefficients differ because the intercept now describes the fitted log-odds at $x=0$, rather than at the midpoint of the two groups.

Common Pitfalls and Practical Checks

Reusing stale working quantities. After replacing $\beta^{(k)}$ by $\beta^{(k+1)}$, recompute predictors, probabilities, working responses, and weights. The next weighted system is not the previous system with a new answer inserted.
Rounding too early. The second and later working responses depend on differences such as $y_i-p_i^{(k)}$. In the symmetric example, the difference at the second update is only about $0.000426465$. Rounding the preceding probability or variance factor too aggressively changes the update. Keep full precision internally and round displayed values only.
Confusing coefficients with predictors. An intercept and slope are not generally the fitted linear predictors themselves. Calculate $X\beta^{(k)}$ before calculating probabilities, and account for the design coding.
Assuming every update is a final estimate. An IRLS step is a new iterate. Its numerical proximity to a target in one example does not establish that an arbitrary sequence will converge.

Key takeaway: One iteration produces a complete new working system, not just a new coefficient vector. Trace the sequence in order: current coefficients, linear predictors, probabilities, working responses and weights, then the weighted least-squares solution.

Check Your Understanding

Use the current coefficient vector at each step; do not carry working quantities forward unchanged.

In the symmetric two-group example, what are the fitted probabilities at $\beta^{(1)}=(0,1)^\mathsf{T}$, and why are the two working weights equal?
At the second update of that example, why can rounding $p_2^{(2)}$ too early have a noticeable effect on the working response?
For an intercept-only binomial model, why does the weighted least-squares update equal the single working response, regardless of the positive weight?
With design rows $(1,0)$ and $(1,1)$, explain how two fitted working values determine the intercept and slope.
State the rank condition that makes the weighted least-squares update unique, and explain why positive logistic working weights matter.