BHHH Optimization From Individual Scores

Why Use Individual Scores?

Maximum-likelihood optimization requires information about both the direction in which the log likelihood increases and its curvature. Newton–Raphson uses the observed curvature; Fisher scoring uses expected information. BHHH, named for Berndt, Hall, Hall, and Hausman, offers another route: it approximates curvature using the outer products of the individual score contributions. This makes the update possible using scores alone, but the approximation can be singular or unstable when those scores do not provide enough independent directions.

Suppose the log likelihood is written as a sum of contributions, \(\ell(\theta)=\sum_{i=1}^n \ell_i(\theta)\), where \(\theta\) is a \(k\)-dimensional parameter vector. The contribution \(\ell_i\) may be the log likelihood for observation \(i\), or a conditional log-likelihood term in a model whose likelihood has been factored into such terms. All contributions and derivatives in an update are evaluated at the same current parameter value.

Definition: The individual score is the column vector \(s_i(\theta)=\partial \ell_i(\theta)/\partial\theta\). The total score is \(S(\theta)=\sum_{i=1}^n s_i(\theta)\). The BHHH matrix, also called the outer-product-of-gradients matrix, is \(B(\theta)=\sum_{i=1}^n s_i(\theta)s_i(\theta)'\). It is a sample approximation to information, not the observed information itself.

For any vector \(a\), \(a'Ba=\sum_i(a's_i)^2\), which cannot be negative. Thus \(B\) is positive semidefinite. It is positive definite, and invertible, only when the individual score vectors span all \(k\) parameter directions. This geometric condition is central to deciding whether a BHHH update is available and numerically useful.

From Score Outer Products to an Update

Write \(H(\theta)=\partial^2\ell(\theta)/\partial\theta\partial\theta'\) for the Hessian of the summed log likelihood. A second-order Taylor approximation around the current value \(\theta\), for a proposed change \(d\), is

\ell(\theta+d)\approx \ell(\theta)+S(\theta)'d+\frac{1}{2}d'H(\theta)d.

Newton–Raphson uses this local quadratic approximation with the actual Hessian. BHHH replaces the Hessian by \(-B(\theta)\), treating the outer-product matrix as an approximation to negative curvature. The resulting approximate objective is maximized by setting its derivative with respect to \(d\) to zero.

q(d)=\ell(\theta)+S'd-\frac{1}{2}d'Bd, \qquad \frac{\partial q(d)}{\partial d}=S-Bd.

If \(B\) is invertible, the stationary point satisfies \(Bd=S\), so the BHHH change is \(d=B^{-1}S\). Since \(B\) is positive definite in this case, \(q\) is concave and this stationary point maximizes the approximate quadratic. The proposed parameter value is the current value plus this change.

\theta_{\mathrm{new}}=\theta+B(\theta)^{-1}S(\theta).

The sign is for maximizing a log likelihood. The total score points toward local increase, and the positive information approximation scales that direction to account for curvature. If instead the objective were a negative log likelihood to minimize, the corresponding signs would be written for that minimization problem.

Justification: The BHHH matrix is motivated by information equality. Under a correctly specified, regular likelihood with suitable differentiation and integrability conditions, the expected score is zero and the expected outer product of a score equals the expected negative Hessian for the corresponding contribution. Summing across independent observations gives \(E[B(\theta)]=E[-H(\theta)]\). In a finite sample, however, \(B\) need not equal the observed information or expected information; it is an approximation based on the scores at the current parameter value.

This justification requires care in dependent-data settings. The relevant individual contributions must be the terms in the likelihood factorization being optimized, and the information equality requires the model and regularity conditions to support it. The fact that a matrix can be computed as a sum of score outer products does not, by itself, establish that it is a good curvature approximation for every model or sample.

Worked Example: A Two-Parameter BHHH Step

Consider three hypothetical likelihood contributions evaluated at a current two-parameter value. Suppose their score vectors are \(s_1=(1,0)'\), \(s_2=(0,1)'\), and \(s_3=(1,1)'\). Calculate the BHHH update.

Total score. Add the three vectors:

$$ S=s_1+s_2+s_3=(2,2)'. $$

Outer-product matrix. Each score contributes its outer product:

B=s_1s_1'+s_2s_2'+s_3s_3' =\begin{pmatrix}1&0\\0&0\end{pmatrix} +\begin{pmatrix}0&0\\0&1\end{pmatrix} +\begin{pmatrix}1&1\\1&1\end{pmatrix} =\begin{pmatrix}2&1\\1&2\end{pmatrix}.

Its determinant is \(2(2)-1(1)=3\), so it is invertible, with \(B^{-1}=\frac{1}{3}\begin{pmatrix}2&-1\\-1&2\end{pmatrix}\). Therefore,

d=B^{-1}S =\frac{1}{3}\begin{pmatrix}2&-1\\-1&2\end{pmatrix} \begin{pmatrix}2\\2\end{pmatrix} =\begin{pmatrix}2/3\\2/3\end{pmatrix}.

The BHHH proposal adds \(2/3\) to each parameter. Under the local quadratic approximation, its predicted log-likelihood change is \(S'd-\frac{1}{2}d'Bd=\frac{8}{3}-\frac{4}{3}=\frac{4}{3}\). This is a prediction from the approximation, not a guarantee that the actual log likelihood increases by that amount.

When the BHHH Matrix Is Singular

Because \(B\) is a sum of score outer products, its rank cannot exceed the number of linearly independent score vectors. If all scores lie in a proper subspace of the \(k\)-dimensional parameter space, then there is a nonzero direction \(a\) orthogonal to every score. In that direction, \(Ba=0\), so \(B\) has no inverse. A BHHH update requiring \(B^{-1}\) is then undefined.

Worked Example: Collinear Scores Produce a Singular Matrix

Suppose three hypothetical observations have two-parameter score vectors \(s_1=(1,2)'\), \(s_2=(2,4)'\), and \(s_3=(-1,-2)'\). These scores are all multiples of \((1,2)'\). Construct \(B\) and determine whether the update exists.

The total score is \(S=(1+2-1,\;2+4-2)'=(2,4)'\). Each outer product is a multiple of \(\begin{pmatrix}1&2\\2&4\end{pmatrix}\): the multiplier is \(1^2\), \(2^2\), or \((-1)^2\). Thus

B=(1+4+1)\begin{pmatrix}1&2\\2&4\end{pmatrix} =\begin{pmatrix}6&12\\12&24\end{pmatrix}.

The determinant is \(6(24)-12(12)=0\), so the inverse does not exist. In fact, \(B(2,-1)'=(0,0)'\), confirming that the direction \((2,-1)'\) has zero curvature according to this approximation. The equation \(Bd=S\) reduces to \(d_1+2d_2=1/3\), which has infinitely many solutions rather than a unique update. At this iteration, the observed score vectors provide information in only one direction. This fact alone does not prove that the model is globally unidentified; it shows why this BHHH matrix cannot produce a unique two-parameter step.

Poor Conditioning and Sensitivity

An invertible matrix can still be nearly singular. When the score vectors are close to lying in a lower-dimensional subspace, the smallest eigenvalue of \(B\) is positive but very small. The inverse then magnifies components of the score in the associated weakly represented direction. As a result, small numerical errors in the scores or their computation can produce large changes in the proposed update.

Worked Example: Nearly Collinear Scores and an Unstable Update

Suppose two hypothetical observations have score vectors \(s_1=(1,1)'\) and \(s_2=(1,1.01)'\). Their scores are not exactly collinear, so construct the matrix and assess its conditioning.

The total score and outer-product matrix are

S=\begin{pmatrix}2\\2.01\end{pmatrix}, \qquad B=\begin{pmatrix}1&1\\1&1\end{pmatrix} +\begin{pmatrix}1&1.01\\1.01&1.0201\end{pmatrix} =\begin{pmatrix}2&2.01\\2.01&2.0201\end{pmatrix}.

The determinant is \(2(2.0201)-(2.01)^2=0.0001\), positive but very small. The eigenvalues are approximately \(4.0201\) and \(0.00002488\), giving a condition number of about \(161{,}600\). The matrix is invertible, but the large condition number signals substantial sensitivity in the direction associated with the small eigenvalue.

For this exact example, the inverse is \(10{,}000\begin{pmatrix}2.0201&-2.01\\-2.01&2\end{pmatrix}\). Multiplying by the total score gives

B^{-1}S =10{,}000\begin{pmatrix} 2.0201(2)-2.01(2.01)\\ -2.01(2)+2(2.01) \end{pmatrix} =10{,}000\begin{pmatrix}0.0001\\0\end{pmatrix} =\begin{pmatrix}1\\0\end{pmatrix}.

The resulting step is \((1,0)'\), but the near-collinearity means small perturbations to the scores could change it substantially. The apparent simplicity of this particular answer does not make the matrix well-conditioned; the eigenvalues and condition number reveal the numerical risk.

Interpreting and Checking a BHHH Proposal

A BHHH update is a local optimization proposal. It uses the current total score for direction and the current score outer products as an information approximation for scale. It is distinct from Fisher scoring, which uses expected information, and from Newton–Raphson, which uses the observed Hessian. The methods can coincide under particular conditions, but their curvature inputs are not interchangeable in general.

The outer-product matrix is positive semidefinite by construction, which is useful for maximizing the local quadratic when it is positive definite. Yet positive semidefiniteness does not ensure invertibility, good conditioning, a reliable approximation to the actual likelihood curvature, or improvement in the actual likelihood. In practice, one should inspect the matrix's rank or eigenvalues and the size of the proposed change, and then assess the likelihood at the proposed value. Detailed safeguards for choosing or modifying a step are separate from the BHHH derivation.

If the likelihood is represented by an average rather than a sum, both the total score and the outer-product matrix are divided by \(n\). Their ratio in the update is unchanged, provided the same convention is used consistently. Mixing a summed score with an averaged matrix, or vice versa, changes the step by a factor of \(n\).

Common Mistake

A frequent error is to treat \(B\) as the observed information and conclude that its inverse is automatically a dependable curvature correction. The observed information is \(-H(\theta)\), formed from second derivatives; \(B\) is formed from score outer products. Information equality motivates the latter as an approximation under suitable conditions, but it does not make the two matrices equal in a finite sample. A second error is to assume that a positive determinant always means a safe update: a positive but tiny determinant, as in the preceding example, can accompany severe ill-conditioning. Check rank and conditioning, and interpret the resulting step as a proposal rather than a guaranteed improvement.

Key takeaway: BHHH replaces the negative Hessian with the sum of individual score outer products and proposes \(\theta_{\mathrm{new}}=\theta+B^{-1}S\). The update requires enough independent score directions for \(B\) to be invertible; near-collinearity can make an otherwise defined update highly sensitive.

Check Your Understanding

Use the definitions and calculations above to answer these questions.

For a summed log likelihood, define the individual score, total score, and BHHH matrix.
Starting from the quadratic approximation with Hessian replaced by \(-B\), derive the equation that gives the BHHH change.
Why does a nonzero vector orthogonal to every individual score imply that \(B\) is singular?
In the nearly collinear example, why does a positive determinant not establish that the BHHH update is numerically reliable?
How does the update change if both the total score and BHHH matrix are consistently computed using average rather than summed contributions?