Solving the IRLS Weighted Least-Squares System

From an IRLS Step to a Linear System

An IRLS iteration turns the current GLM fit into a weighted least-squares problem. The working response and working weights are recalculated from the current coefficients; the computational task in this tutorial is to solve the resulting weighted regression efficiently. The central principle is simple: factor the system and solve it; do not form an inverse of the information matrix.

Let $X$ be the $n\times p$ design matrix, with one row per observation and one column per coefficient. Let $z$ be the current IRLS working-response vector, $W=\operatorname{diag}(w_1,\ldots,w_n)$ the diagonal matrix of working weights, and $o$ a known offset vector. If there is no offset, take $o=0$. The weighted least-squares step seeks coefficients whose linear predictors $X\beta+o$ are close to $z$, with observations weighted according to $W$.

Definition: The IRLS weighted least-squares criterion is $$ Q(\beta)=(z-o-X\beta)^\mathsf{T}W(z-o-X\beta) =\sum_{i=1}^n w_i(z_i-o_i-x_i^\mathsf{T}\beta)^2. $$ Here $x_i^\mathsf{T}$ is row $i$ of $X$. The weighted normal equations are $$ X^\mathsf{T}WX\beta=X^\mathsf{T}W(z-o). $$

For the factorization methods below, assume first that every $w_i>0$ and that $X$ has full column rank $p$. These conditions make $X^\mathsf{T}WX$ positive definite and the coefficient solution unique. If a weight is zero, its observation contributes nothing and can be omitted from the weighted calculation. Rank deficiency requires separate care, discussed below.

Why the Weighted Normal Equations Hold

The normal equations express a basic least-squares condition: at the minimizing coefficients, the weighted residuals are orthogonal to every column of the design matrix. The result follows directly by differentiating the criterion.

Theorem: If $W$ is diagonal with positive entries and $X$ has full column rank, $Q(\beta)$ has a unique minimizer. It is the solution of $X^\mathsf{T}WX\beta=X^\mathsf{T}W(z-o)$.

Proof. Write $r(\beta)=z-o-X\beta$. For any coefficient vector $h$, $$ Q(\beta+h)-Q(\beta) =-2h^\mathsf{T}X^\mathsf{T}Wr(\beta)+h^\mathsf{T}X^\mathsf{T}WXh. $$ At a minimizer, the linear term must vanish for every $h$, so $X^\mathsf{T}Wr(\beta)=0$. Substituting $r(\beta)=z-o-X\beta$ gives the weighted normal equations. Conversely, if these equations hold, then for any $h$, $$ Q(\beta+h)-Q(\beta)=h^\mathsf{T}X^\mathsf{T}WXh. $$ For nonzero $h$, full column rank implies $Xh\ne0$, and positive weights imply $$ h^\mathsf{T}X^\mathsf{T}WXh=(Xh)^\mathsf{T}W(Xh)>0. $$ Thus a solution of the equations is a strict minimizer, and it is unique. $\square$

Set $A=X^\mathsf{T}WX$ and $b=X^\mathsf{T}W(z-o)$. The equations are $A\beta=b$, but the intended computation is not $\beta=A^{-1}b$. Explicitly forming $A^{-1}$ requires extra work, can amplify numerical error, and is unnecessary: a factorization reduces the task to triangular systems.

Cholesky Factorization and Triangular Solves

Under the full-rank and positive-weight assumptions, $A$ is symmetric positive definite. It therefore has a Cholesky factorization $A=R^\mathsf{T}R$, where $R$ is upper triangular with positive diagonal entries. Solving the weighted normal equations then takes two triangular solves.

Form the system.
Compute $A=X^\mathsf{T}WX$ and $b=X^\mathsf{T}W(z-o)$.

Factor the matrix.
Compute the Cholesky factor $R$ satisfying $A=R^\mathsf{T}R$.

Solve forward.
Find $u$ from $R^\mathsf{T}u=b$.

Solve backward.
Find the coefficient update $\beta$ from $R\beta=u$.

This procedure is correct because the two triangular equations imply $$ R^\mathsf{T}R\beta=R^\mathsf{T}u=b, $$ which is exactly $A\beta=b$. Each triangular solve can be carried out by substitution. In an implementation, it is often useful to compute the weighted design and response directly: let $W^{1/2}$ be diagonal with entries $\sqrt{w_i}$, set $B=W^{1/2}X$, and set $c=W^{1/2}(z-o)$. Then $A=B^\mathsf{T}B$ and $b=B^\mathsf{T}c$.

Worked Example: A Poisson Log-Link Step by Cholesky

Consider three independent Poisson observations with log means $\eta_i=\beta_0+\beta_1x_i$, covariate values $x=(0,1,2)$, and counts $y=(1,2,4)$. Start the iteration at $\beta=(0,0)^\mathsf{T}$. For a Poisson log-link model, the current means are $\mu_i=1$, the working weights are $w_i=\mu_i=1$, and the working responses are $z_i=\eta_i+(y_i-\mu_i)/\mu_i=y_i-1$. Thus $z=(0,1,3)$.

The design matrix has rows $(1,0),(1,1),(1,2)$. With unit weights and no offset, $$ A=X^\mathsf{T}X= \begin{pmatrix}3&3\\3&5\end{pmatrix}, \qquad b=X^\mathsf{T}z= \begin{pmatrix}4\\7\end{pmatrix}. $$ A Cholesky factor is $$ R=\begin{pmatrix}\sqrt{3}&\sqrt{3}\\0&\sqrt{2}\end{pmatrix}, \qquad R^\mathsf{T}R= \begin{pmatrix}3&3\\3&5\end{pmatrix}. $$

First solve $R^\mathsf{T}u=b$. The first equation gives $\sqrt{3}u_1=4$, hence $u_1=4/\sqrt{3}$. The second gives $\sqrt{3}u_1+\sqrt{2}u_2=7$, so $4+\sqrt{2}u_2=7$ and $u_2=3/\sqrt{2}$. Next solve $R\beta=u$. Its second equation gives $\sqrt{2}\beta_1=3/\sqrt{2}$, so $\beta_1=3/2$. Its first gives $\sqrt{3}(\beta_0+\beta_1)=4/\sqrt{3}$, so $\beta_0+\beta_1=4/3$ and $\beta_0=-1/6$. The weighted least-squares update is therefore $(-1/6,3/2)^\mathsf{T}$.

As a check, the fitted working responses at $x=0,1,2$ are $-1/6,4/3,17/6$. The residuals are $1/6,-1/3,1/6$; they sum to zero, and their products with the centered covariate values $-1,0,1$ also sum to zero, as the least-squares condition requires.

QR Factorization: Solving Without Forming the Normal Matrix

Cholesky is convenient when the weighted normal matrix is already available and well-conditioned. A QR factorization instead works directly with the weighted design $B=W^{1/2}X$. When numerical stability matters, this can be preferable: forming $B^\mathsf{T}B$ squares the condition number of $B$, whereas QR avoids that step.

Theorem: Suppose $B=W^{1/2}X$ has full column rank. Let $B=QR$ be a reduced QR factorization, where $Q^\mathsf{T}Q=I_p$ and $R$ is an invertible upper-triangular matrix. Then the unique minimizer of $\lVert c-B\beta\rVert^2$, with $c=W^{1/2}(z-o)$, is obtained by solving $$ R\beta=Q^\mathsf{T}c. $$

Proof. Since the columns of $Q$ are orthonormal, decompose $c-B\beta$ into its projection onto the column space of $Q$ and its orthogonal component: $$ \lVert c-QR\beta\rVert^2 =\lVert Q^\mathsf{T}c-R\beta\rVert^2 +\lVert (I-QQ^\mathsf{T})c\rVert^2. $$ The second term does not depend on $\beta$. The first is nonnegative and equals zero when $R\beta=Q^\mathsf{T}c$. Because $R$ is invertible, that equation has exactly one solution, which is the unique minimizer. $\square$

Worked Example: A Grouped Binomial Step by QR

Suppose two groups have binomial counts $s=(1,2)$ out of $m=(2,2)$, with a logit model $\operatorname{logit}(p_i)=\beta_0+\beta_1x_i$ and $x=(0,1)$. At the current value $\beta=(0,0)^\mathsf{T}$, both fitted probabilities are $p_i=1/2$. For grouped proportions $r_i=s_i/m_i$, the working response is $z_i=\eta_i+(r_i-p_i)/[p_i(1-p_i)]$, and the working weight is $w_i=m_ip_i(1-p_i)$. Hence $z=(0,2)$ and $w_1=w_2=1/2$.

The design and weighted quantities are $$ X=\begin{pmatrix}1&0\\1&1\end{pmatrix}, \qquad B=W^{1/2}X= \begin{pmatrix}1/\sqrt{2}&0\\1/\sqrt{2}&1/\sqrt{2}\end{pmatrix}, \qquad c=W^{1/2}z= \begin{pmatrix}0\\\sqrt{2}\end{pmatrix}. $$ The columns of $B$ are independent. Applying Gram–Schmidt gives $$ Q=\begin{pmatrix}1/\sqrt{2}&-1/\sqrt{2}\\1/\sqrt{2}&1/\sqrt{2}\end{pmatrix}, \qquad R=\begin{pmatrix}1&1/2\\0&1/2\end{pmatrix}. $$ Indeed, $Q^\mathsf{T}Q=I$ and $QR=B$. Also, $$ Q^\mathsf{T}c=\begin{pmatrix}1\\1\end{pmatrix}. $$ Thus $R\beta=Q^\mathsf{T}c$ gives $\beta_1/2=1$, so $\beta_1=2$, and then $\beta_0+\beta_1/2=1$, so $\beta_0=0$. The result fits the two working responses exactly. It can also be checked from $A=B^\mathsf{T}B=\begin{pmatrix}1&1/2\\1/2&1/2\end{pmatrix}$ and $b=B^\mathsf{T}c=(1,1)^\mathsf{T}$, which yield the same equations.

Rank Deficiency and the Limits of a Unique Solve

A factorization does not create information that the design lacks. If columns of $X$ are linearly dependent, $X^\mathsf{T}WX$ is singular: Cholesky cannot produce a nonsingular factor, and a full-rank QR solve is not available. A rank-revealing, often column-pivoted, QR factorization can identify dependent directions. The model may still determine fitted linear predictors, even when it does not determine one unique coefficient vector.

Worked Example: Duplicate Predictors and a Nonunique Update

Take a Poisson log-link model with two observations, both having the same predictor value $x_i=1$, and use an intercept and that predictor as the two columns of $X$. Let the counts be $y=(1,3)$ and start at $\beta=(0,0)^\mathsf{T}$. The current means and weights are both $1$, and the working responses are $z=y-1=(0,2)$. The design has two identical rows: $$ X=\begin{pmatrix}1&1\\1&1\end{pmatrix},\qquad A=X^\mathsf{T}X=\begin{pmatrix}2&2\\2&2\end{pmatrix},\qquad b=X^\mathsf{T}z=\begin{pmatrix}2\\2\end{pmatrix}. $$

The equations reduce to $2\beta_0+2\beta_1=2$, or $\beta_0+\beta_1=1$. There are infinitely many solutions; for example, $(1,0)^\mathsf{T}$ and $(0,1)^\mathsf{T}$ both give fitted working responses $(1,1)$. The columns of $X$ are identical, so $A$ is singular and has no Cholesky factor with positive diagonal. A rank-revealing factorization would expose this dependence. If a constraint such as $\beta_1=0$ is imposed, the constrained solution is $(1,0)^\mathsf{T}$, but that constraint is additional information, not a consequence of the weighted least-squares criterion.

Common Mistakes and Practical Choices

Explicitly inverting the information matrix. Writing $\beta=A^{-1}b$ is an algebraic identity when $A$ is invertible, not a recommended computation. Factor $A$ and use triangular solves, or solve the weighted regression with QR.
Ignoring the offset. With a known offset, the response on the right-hand side is $z-o$, not $z$. Omitting $o$ changes the least-squares problem.
Assuming every system has a unique solution. A successful-looking set of equations does not guarantee full column rank. Check the design; use a rank-revealing factorization when dependence is possible.
Treating Cholesky and QR as equally stable in every setting. Cholesky is efficient for a well-conditioned positive-definite normal matrix. QR works on the weighted design directly and is generally safer when the columns are nearly dependent.

The factorization choice matters because each IRLS step may be repeated many times. A reliable solve helps ensure that the computed update is the solution to the intended weighted regression, while a rank check warns when the coefficients are not separately determined. The working responses and weights define the current subproblem; the factorization is the numerical method for solving it.

Check Your Understanding

Use the weighted least-squares formulation and factorization ideas to answer these questions.

Write the weighted normal equations when the IRLS working regression includes an offset vector $o$.
Why does $X^\mathsf{T}WX$ become positive definite when all weights are positive and $X$ has full column rank?
In a Cholesky solve with $A=R^\mathsf{T}R$, which triangular system is solved first, and which is solved second?
What equation does a reduced QR factorization require you to solve for $\beta$?
If two columns of the design are identical, what does that imply about uniqueness of the coefficient update, and how can a factorization help diagnose the problem?