Solving Information Systems Without Matrix Inversion

Why Solve a System Instead of Forming an Inverse?

A Newton or Fisher-scoring update for a generalized linear model is often written using the inverse of an information matrix. That notation is compact, but it does not mean that a computer should calculate and store the inverse. The computational task is to find the vector that solves a linear system. Solving that system directly typically takes less work and avoids introducing the extra rounding error that can come from forming an inverse and then multiplying by a vector.

Let \(\beta\) be the current coefficient vector, \(U(\beta)=\nabla\ell(\beta)\) the score, and \(J(\beta)=-\nabla^2\ell(\beta)\) the observed information. A Newton step for maximizing the log-likelihood solves a system with \(J(\beta)\). Fisher scoring instead uses the expected information \(I(\beta)\). Provided the chosen matrix is nonsingular, the step can be expressed as follows.

J(\beta)\,\delta=U(\beta) \qquad\text{or}\qquad I(\beta)\,\delta=U(\beta), \qquad \beta_{\mathrm{new}}=\beta+\delta.

The familiar expression \(\delta=I(\beta)^{-1}U(\beta)\) describes the same mathematical solution when the inverse exists. In an implementation, however, the preferred operation is to solve \(I(\beta)\delta=U(\beta)\) directly. This tutorial develops the method for a symmetric positive-definite information matrix, where Cholesky factorization provides an especially efficient route.

Cholesky Factorization and the Two Substitutions

Definition: A real symmetric matrix \(A\) is positive definite if \(x^{T}Ax>0\) for every nonzero real vector \(x\). The Cholesky factorization of a positive-definite matrix is a representation \(A=LL^{T}\), where \(L\) is lower triangular and has strictly positive diagonal entries.

If the information matrix \(A\) is positive definite, factor it as \(A=LL^{T}\). To solve \(A\delta=b\), first solve the lower-triangular system \(Lz=b\), then solve the upper-triangular system \(L^{T}\delta=z\). Each triangular solve uses substitution rather than a matrix inverse.

Factor.
Compute \(A=LL^{T}\), with \(L\) lower triangular and positive on its diagonal.

Forward substitution.
Compute \(z\) from \(Lz=b\), proceeding from the first row to the last.

Back substitution.
Compute \(\delta\) from \(L^{T}\delta=z\), proceeding from the last row to the first.

Check.
Evaluate the residual \(r=b-A\delta\) to assess whether the computed vector satisfies the original system to the desired numerical accuracy.

For a dense \(p\times p\) matrix, Cholesky factorization requires about \(p^{3}/3\) arithmetic operations, followed by two triangular solves requiring about \(p^{2}\) operations each. If the same information matrix is used with several right-hand sides, its factorization can be reused. By contrast, computing an entire inverse does work that is unnecessary when only one or a few solutions are needed.

Why the Factorization Exists

Theorem: Every real symmetric positive-definite matrix \(A\) has a unique factorization \(A=LL^{T}\) in which \(L\) is lower triangular with strictly positive diagonal entries.

Proof. Use induction on the dimension. In dimension one, positive definiteness means \(A=[a]\) with \(a>0\), and \(A=LL^{T}\) for \(L=[\sqrt{a}]\); the positive diagonal requirement makes this factor unique.

For the induction step, write a symmetric positive-definite \(p\times p\) matrix in block form,

A= \begin{pmatrix} a & r^{T}\\ r & B \end{pmatrix}.

Positive definiteness applied to the first coordinate vector gives \(a>0\). Define \(v=r/\sqrt{a}\) and \(S=B-vv^{T}=B-rr^{T}/a\). To show \(S\) is positive definite, take any nonzero vector \(y\) of dimension \(p-1\), and set \(x=-r^{T}y/a\). Then \((x,y^{T})^{T}\) is nonzero, so positive definiteness of \(A\) gives

0< \begin{pmatrix}x\\y\end{pmatrix}^{T} A \begin{pmatrix}x\\y\end{pmatrix} = y^{T}\left(B-\frac{rr^{T}}{a}\right)y =y^{T}Sy.

By the induction hypothesis, \(S=CC^{T}\) for a unique lower-triangular \(C\) with positive diagonal. Direct block multiplication now gives

L= \begin{pmatrix} \sqrt{a} & 0\\ r/\sqrt{a} & C \end{pmatrix}, \qquad LL^{T}= \begin{pmatrix} a & r^{T}\\ r & rr^{T}/a+CC^{T} \end{pmatrix} =A.

This construction also proves uniqueness: the upper-left entry of \(L\) must be the positive square root of \(a\), its first column below that entry must be \(r/\sqrt{a}\), and its remaining block must be the unique positive-diagonal Cholesky factor of \(S\). The induction is complete. \(\square\)

Once the factorization exists, the solution is unique. Indeed, if \(A=LL^{T}\) and \(L\) has positive diagonal, both triangular systems have unique solutions. This gives a precise reason the method applies under the positive-definiteness condition; it is not merely a convenient formula.

Worked GLM Updates

Worked Example: A Poisson Intercept Update

Suppose \(n=5\) independent Poisson responses share an intercept-only mean \(\mu=\exp(\beta)\). At the current value \(\beta=\log 2\), the fitted mean is \(\mu=2\). If the observed counts sum to 14, the score for \(\beta\) is \(\sum_i y_i-n\mu\), and the information is \(n\mu\).

The score is \(U=14-5(2)=4\), and the information is \(I=5(2)=10\). The update is found by solving the one-dimensional information system:

10\delta=4, \qquad \delta=0.4.

Thus the updated intercept is \(\beta_{\mathrm{new}}=\log 2+0.4\). In this scalar case, writing \(I^{-1}U\) gives the same arithmetic, but the system formulation is the one that generalizes to many coefficients: compute the solution of \(I\delta=U\), not an inverse as an intermediate object.

Worked Example: Fisher Scoring for Grouped Logistic Data

Consider a logistic model with two coefficients and two groups. The design rows are \(x_1^{T}=(1,0)\) and \(x_2^{T}=(1,1)\); each group has four trials and respectively three and one successes. At \(\beta=(0,0)^{T}\), both fitted success probabilities are \(1/2\), so each group has variance weight \(4(1/2)(1/2)=1\).

The score is the design matrix transposed times the vector of observed-minus-fitted successes. The fitted success counts are two in each group, so the residuals are \(1\) and \(-1\). Therefore

U= \begin{pmatrix}1&1\\0&1\end{pmatrix} \begin{pmatrix}1\\-1\end{pmatrix} = \begin{pmatrix}0\\-1\end{pmatrix}.

The Fisher information is the sum of the two weighted design-row outer products:

I= \begin{pmatrix}1&0\\0&0\end{pmatrix} + \begin{pmatrix}1&1\\1&1\end{pmatrix} = \begin{pmatrix}2&1\\1&1\end{pmatrix}.

A Cholesky factor is \(L=\begin{pmatrix}\sqrt{2}&0\\1/\sqrt{2}&1/\sqrt{2}\end{pmatrix}\), since multiplying \(LL^{T}\) gives \(I\). First solve \(Lz=U\): the first row gives \(z_1=0\), and the second gives \(z_2=-\sqrt{2}\). Next solve \(L^{T}\delta=z\). Its second row gives \(\delta_2=-2\); its first row gives \(\sqrt{2}\delta_1-\sqrt{2}=0\), so \(\delta_1=1\). The result is \(\delta=(1,-2)^{T}\). Substitution checks both equations of the original system:

I\delta= \begin{pmatrix}2&1\\1&1\end{pmatrix} \begin{pmatrix}1\\-2\end{pmatrix} = \begin{pmatrix}0\\-1\end{pmatrix} =U.

The Fisher-scoring proposal is consequently \(\beta_{\mathrm{new}}=(1,-2)^{T}\). This calculation determines the update; it does not, by itself, guarantee that a full step improves every practical fitting objective.

Worked Example: A Three-Coefficient System and a Residual Check

Suppose a GLM iteration produces the symmetric information system \(A\delta=b\), with

A= \begin{pmatrix} 4&2&2\\ 2&5&1\\ 2&1&3 \end{pmatrix}, \qquad b= \begin{pmatrix}8\\8\\6\end{pmatrix}.

A proposed Cholesky factor is

L= \begin{pmatrix} 2&0&0\\ 1&2&0\\ 1&0&\sqrt{2} \end{pmatrix}.

Checking the factorization gives \(LL^{T}=A\): its diagonal entries are \(4,5,3\), and its off-diagonal entries are \(2,2,1\), as required. Forward substitution in \(Lz=b\) gives \(2z_1=8\), hence \(z_1=4\); \(z_1+2z_2=8\), hence \(z_2=2\); and \(z_1+\sqrt{2}z_3=6\), hence \(z_3=\sqrt{2}\).

Back substitution in \(L^{T}\delta=z\) gives \(\sqrt{2}\delta_3=\sqrt{2}\), so \(\delta_3=1\); \(2\delta_2=2\), so \(\delta_2=1\); and \(2\delta_1+\delta_2+\delta_3=4\), so \(\delta_1=1\). Finally,

A\delta= \begin{pmatrix} 4+2+2\\ 2+5+1\\ 2+1+3 \end{pmatrix} = \begin{pmatrix}8\\8\\6\end{pmatrix} =b.

The residual \(b-A\delta\) is exactly zero in this arithmetic. In floating-point computation it will generally not be exactly zero, so its size should be assessed relative to the scale of the system and the accuracy required.

What the Update Means, and When Cholesky Applies

Result: If a symmetric positive-definite matrix \(A\) solves \(A\delta=U\), then \(U^{T}\delta>0\) whenever \(U\ne 0\). Thus \(\delta\) is an ascent direction for a differentiable objective with gradient \(U\).

Proof. Since \(A\delta=U\), we have \(U^{T}\delta=(A\delta)^{T}\delta=\delta^{T}A\delta\), using symmetry of \(A\). If \(U\ne 0\), then \(\delta\ne 0\), because \(A\delta=0\) would imply \(\delta^{T}A\delta=0\), which positive definiteness permits only for \(\delta=0\), and then \(U=0\). Positive definiteness therefore gives \(\delta^{T}A\delta>0\), hence \(U^{T}\delta>0\). By the definition of differentiability, the directional derivative of the objective in direction \(\delta\) is \(U^{T}\delta\), so it is positive. \(\square\)

This result explains why a positive-definite information system produces a locally sensible direction for maximizing the log-likelihood. It does not assert that any chosen step length will improve the objective, nor that the method will converge from every starting value. Those are separate questions from solving the linear system accurately.

In practice, a Cholesky routine should report whether factorization succeeded, and the computed solution should be checked against the original system. A small residual is useful evidence that the equations were solved accurately, though it is not by itself a guarantee that the solution is insensitive to small perturbations in the input. Cholesky is appropriate when the matrix supplied is symmetric positive definite. If it is not, the assumptions behind this method do not hold; do not treat a failed factorization as permission to silently alter the system.

Common Mistakes

Computing an inverse just to obtain an update. Forming \(A^{-1}\) and multiplying by \(U\) adds unnecessary work. Solve \(A\delta=U\) directly.
Using the wrong sign for the update. With \(U=\nabla\ell(\beta)\) and positive information \(J=-\nabla^{2}\ell(\beta)\), the Newton equation is \(J\delta=U\), followed by \(\beta_{\mathrm{new}}=\beta+\delta\).
Assuming every information matrix can be Cholesky-factorized. The method described here requires symmetry and positive definiteness. These are mathematical conditions, not optional numerical preferences.
Checking only the factorization, not the solution. After substitution, verify that the computed \(\delta\) makes \(A\delta\) close to \(b\); this tests the original system directly.

Key takeaway: For a suitable GLM information matrix, compute an update by factoring \(A=LL^{T}\), solving two triangular systems, and checking \(A\delta=b\). The inverse notation expresses a solution; it is usually not the computation to perform.

Check Your Understanding

Try these questions without looking back.

For a Newton update with score \(U(\beta)\) and observed information \(J(\beta)=-\nabla^{2}\ell(\beta)\), write the linear system and the coefficient update using the sign convention in this tutorial.
Given a Cholesky factorization \(A=LL^{T}\), what are the two triangular systems that yield the solution to \(A\delta=b\)?
Why does positive definiteness guarantee a unique Cholesky factor with positive diagonal entries?
In the grouped logistic example, multiply the computed update by the information matrix to verify the score equation.
What does a small residual tell you about a computed solution, and what does it not establish on its own?