ML Factor-Analysis Score Equations

Why the Likelihood Equations Matter

In Gaussian factor analysis, the covariance model is simple to state but generally difficult to fit by direct algebra. The covariance combines products of loadings with variable-specific uniquenesses, and the likelihood depends on the inverse and determinant of that covariance. Differentiating the likelihood gives equations that describe an interior maximum-likelihood solution, even when those equations cannot be solved in closed form.

This tutorial derives those equations and interprets their structure. The derivation uses the common-factor model with standardized, mutually uncorrelated factors; it does not impose a particular rotation or loading pattern.

Definition: For \(p\) observed variables and \(k\) factors, let \(X_i\) be independent observations from \(N_p(\mu,\Sigma)\), where \(\Sigma=\Lambda\Lambda^\mathsf{T}+\Psi\). Here \(\Lambda\) is a \(p\times k\) loading matrix and \(\Psi=\operatorname{diag}(\psi_1,\ldots,\psi_p)\) is the uniqueness matrix, with \(\psi_j>0\). Define \(S=n^{-1}\sum_{i=1}^n(X_i-\bar X)(X_i-\bar X)^\mathsf{T}\). The loadings and uniquenesses are the parameters whose likelihood score equations are derived below.

The mean \(\mu\) is unrestricted and is estimated by \(\bar X\). After substituting that estimate, the log-likelihood, up to an additive constant not depending on \(\Lambda\) or \(\Psi\), is

\ell(\Lambda,\Psi) =-\frac{n}{2}\left\{\log|\Sigma|+\operatorname{tr}(S\Sigma^{-1})\right\}.

The factor \(n\) appears because \(S\) uses divisor \(n\). The equations below are for an interior point of the parameter space: in particular, every \(\psi_j\) is strictly positive. If a solution is constrained to a boundary or to a particular identification surface, its constrained stationarity conditions must also respect those constraints.

Deriving the Score Equations

Write \(B=\Sigma^{-1}(S-\Sigma)\Sigma^{-1}\). This symmetric matrix summarizes the covariance discrepancy in the metric imposed by the model covariance. To differentiate the likelihood, use the matrix differentials \(d\log|\Sigma|=\operatorname{tr}(\Sigma^{-1}d\Sigma)\) and \(d\Sigma^{-1}=-\Sigma^{-1}(d\Sigma)\Sigma^{-1}\). The second identity follows by differentiating \(\Sigma\Sigma^{-1}=I\); the first is the differential of the log determinant for nonsingular \(\Sigma\).

Applying these identities and using the cyclic property of the trace gives

\begin{aligned} d\ell &=-\frac{n}{2}\left\{ \operatorname{tr}(\Sigma^{-1}d\Sigma) -\operatorname{tr}\big(S\Sigma^{-1}(d\Sigma)\Sigma^{-1}\big) \right\}\\ &=\frac{n}{2}\operatorname{tr}\left[ \Sigma^{-1}(S-\Sigma)\Sigma^{-1}d\Sigma \right]\\ &=\frac{n}{2}\operatorname{tr}(B\,d\Sigma). \end{aligned}

Now \(d\Sigma=d\Lambda\,\Lambda^\mathsf{T}+\Lambda\,d\Lambda^\mathsf{T}+d\Psi\). Since \(B\) is symmetric, the two loading terms contribute equally:

\operatorname{tr}\left[B(d\Lambda\,\Lambda^\mathsf{T} +\Lambda\,d\Lambda^\mathsf{T})\right] =2\operatorname{tr}\big[(B\Lambda)^\mathsf{T}d\Lambda\big].

Also, \(d\Psi\) is diagonal, with entries \(d\psi_j\), so \(\operatorname{tr}(B\,d\Psi)=\sum_{j=1}^pB_{jj}\,d\psi_j\). These identities yield the two score matrices.

Theorem: At an interior stationary point of the unconstrained Gaussian factor-analysis likelihood, the loading and uniqueness scores are \(\partial\ell/\partial\Lambda=nB\Lambda\) and \(\partial\ell/\partial\psi_j=(n/2)B_{jj}\). Thus the score equations are \[ \Sigma^{-1}(S-\Sigma)\Sigma^{-1}\Lambda=0, \qquad \operatorname{diag}\!\left\{\Sigma^{-1}(S-\Sigma)\Sigma^{-1}\right\}=0. \]

The equations are necessary conditions for an interior maximum, not a guarantee of a maximum. The loading equation requires \(B\Lambda=0\), whereas the uniqueness equations require only the diagonal of \(B\) to vanish. In particular, they do not generally require \(B=0\), or \(S=\Sigma\).

Worked Examples

Worked Example: Evaluating Scores in a Two-Variable Model

Take \(p=2\), \(k=1\), with loading vector \(\lambda=(0.5,0.4)^\mathsf{T}\), uniquenesses \(\psi_1=0.75\), \(\psi_2=0.84\), and sample covariance \(S=\begin{pmatrix}1.1&0.25\\0.25&0.9\end{pmatrix}\). Compute the scores.

The model covariance and its inverse are

\Sigma=\lambda\lambda^\mathsf{T}+\Psi =\begin{pmatrix}1&0.2\\0.2&1\end{pmatrix}, \qquad \Sigma^{-1}=\frac{1}{0.96} \begin{pmatrix}1&-0.2\\-0.2&1\end{pmatrix}.

Here \(S-\Sigma=\begin{pmatrix}0.1&0.05\\0.05&-0.1\end{pmatrix}\). Multiplying on both sides by \(\Sigma^{-1}\) gives

B=\frac{1}{0.96^2} \begin{pmatrix}0.076&0.052\\0.052&-0.116\end{pmatrix} \approx \begin{pmatrix} 0.08247&0.05642\\ 0.05642&-0.12587 \end{pmatrix}.

The loading score, divided by \(n\), is

B\lambda\approx \begin{pmatrix} 0.08247(0.5)+0.05642(0.4)\\ 0.05642(0.5)-0.12587(0.4) \end{pmatrix} = \begin{pmatrix}0.06380\\-0.02214\end{pmatrix}.

The uniqueness scores, divided by \(n\), are \(B_{11}/2\approx0.04123\) and \(B_{22}/2\approx-0.06293\). None of these scores is zero, so this parameter choice does not satisfy the interior score equations. The calculations show both how to evaluate the equations and why a discrepancy in one covariance entry can affect several parameter scores.

Worked Example: A One-Variable Model Has Many Score Solutions

For one observed variable and one factor, write \(\Sigma=\lambda^2+\psi=v\), where \(\psi>0\), and let the sample variance be \(s>0\). Find the score equations.

In this case \(B=(s-v)/v^2\). The loading and uniqueness scores are therefore

\frac{\partial\ell}{\partial\lambda} =n\lambda\frac{s-v}{v^2}, \qquad \frac{\partial\ell}{\partial\psi} =\frac{n}{2}\frac{s-v}{v^2}.

For the uniqueness score to vanish, \(v=s\). Every choice \(0<\psi<s\), paired with either \(\lambda=\sqrt{s-\psi}\) or \(\lambda=-\sqrt{s-\psi}\), has \(v=s\) and makes both scores zero. The likelihood identifies the total variance here, not its decomposition into a factor contribution and a uniqueness. Thus, even exact satisfaction of the score equations does not by itself identify all parameters.

Worked Example: Zero Scores Without an Exact Covariance Fit

Let \(p=2\), \(k=1\), set \(\lambda=(0,0)^\mathsf{T}\) and \(\Psi=I\), and suppose \(S=\begin{pmatrix}1&0.2\\0.2&1\end{pmatrix}\). This \(S\) is positive definite because its determinant is \(1-0.2^2=0.96>0\). Evaluate the scores.

The model covariance is \(\Sigma=I\), so

B=\Sigma^{-1}(S-\Sigma)\Sigma^{-1} =\begin{pmatrix}0&0.2\\0.2&0\end{pmatrix}.

Its diagonal is zero, and \(B\lambda=0\). Consequently, every loading and uniqueness score is zero, even though \(S\ne\Sigma\): the sample has a nonzero off-diagonal covariance that the fitted covariance does not reproduce. This is a stationary point of the unconstrained parametrization, not evidence that the model fits exactly or that the point is a global maximum.

Why Closed-Form Solutions Are Generally Unavailable

The score equations are coupled through \(\Sigma^{-1}\), while \(\Sigma\) itself contains the products \(\Lambda\Lambda^\mathsf{T}\) and the diagonal parameters \(\psi_j\). Changing a loading affects multiple covariance entries; changing a uniqueness affects one diagonal entry but also changes the inverse and determinant used in every score. The equations are therefore nonlinear and must generally be solved jointly by numerical methods.

There is an additional structural reason not to expect a unique algebraic answer. With an unrestricted loading matrix, different factor orientations can represent the same covariance: for an orthogonal matrix \(Q\), \((\Lambda Q)(\Lambda Q)^\mathsf{T}=\Lambda\Lambda^\mathsf{T}\). Identification constraints can select a representation, but they do not turn the general coupled likelihood equations into a universal closed-form solution. The stationary point in the third example also illustrates that solving score equations is not, on its own, a global optimization check.

Common Mistake: Requiring the Whole Matrix \(B\) to Vanish

A tempting but incorrect argument is: “At a maximum, the model covariance must equal the sample covariance, so \(B=0\).” That would be a sufficient condition for zero scores, but it is not what the derivative requires. The loading score tests \(B\Lambda\), and the uniqueness scores test only the diagonal of \(B\). The third example has zero scores while \(B\) has nonzero off-diagonal entries.

The correction is to impose the parameter-specific equations, not an unnecessarily strong matrix equation. Then assess whether a solution is an appropriate maximum, whether the parameters are identified, and whether any constraints or boundaries affect stationarity. For a boundary value such as \(\psi_j=0\), the interior equality score condition for that uniqueness need not apply; boundary optimization requires the corresponding constrained conditions.

Key takeaway: For an interior Gaussian factor-analysis solution, the loading score equation is \(B\Lambda=0\) and the uniqueness score equations are \(\operatorname{diag}(B)=0\), where \(B=\Sigma^{-1}(S-\Sigma)\Sigma^{-1}\). These coupled equations are generally nonlinear; their satisfaction is necessary for an interior optimum, but does not alone establish fit, identification, or global optimality.

Check Your Understanding

Use the definitions and score equations above to answer the following.

Why does the uniqueness score involve the diagonal of \(B\), rather than every entry of \(B\)?
If \(B\Lambda=0\) and \(\operatorname{diag}(B)=0\), must \(B\) be the zero matrix? Explain.
In the one-variable example, what feature of the covariance is determined by the likelihood, and what feature is not separately identified?
Why does a zero score at an interior parameter value not prove that the value is a global maximum?
How would the stationarity conditions need to change if a uniqueness parameter were constrained to the boundary \(\psi_j=0\)?