On Asymptotic Covariance of Unrotated Factor Solutions

1. Introduction

Let $X_1,X_2,\dots ,X_p$ be random variables, collectively denoted by the random vector Y, Common or exploratory factor analysis postulates the existence of k independent common factors, with $k<p$, that underlie the observed variables Y. Let these factors be jointly represented by the random vector F. The relationship between the observed variables and the common factors is specified by

$$Y=c+\Lambda F+ϵ,$$

(1)

where c is a constant vector, $\Lambda$ is a $p\times k$ factor-loading matrix, and $ϵ$ is a vector of p mutually independent idiosyncratic disturbances.

Let $\Sigma$ denote the variance-covariance (or correlation) matrix of Y, and assume that the common factors have unit variance. Under these assumptions, the covariance structure of Y is given by

$$\Sigma =\Lambda {\Lambda }^{\top }+\mathrm{\Psi },$$

(2)

where the superscript ⊤ denotes matrix transpose, and $\mathrm{\Psi }$ is a $p\times p$ diagonal matrix whose entries correspond to the variances of the idiosyncratic disturbances, often referred to as uniquenesses.

Eq. (2) implies that the off-diagonal covariances of Y are entirely determined by the factor loadings. In contrast, the diagonal variances reflect contributions from both the common and idiosyncratic components: the common factors capture the variation shared across variables Y, while the idiosyncratic disturbances represent variable-specific noise. Let ${\lambda }_{it}$ denote the element in the i-th row and t-th column of $\Lambda$. The quantity

$$h_i=\sum _{t=1}^k{\lambda }_{it}^2,$$

(3)

represents the proportion of the variance of $X_i$ explained by the common factors and is known as the communality of variable $X_i$.

A wide range of factor-extraction methods has been proposed in the literature (e.g., Harman [1976] [7]; Lawley and Maxwell [1971] [15]; SAS [2015] [19]), including principal-factor, maximum-likelihood, and least-squares approaches. These methods aim to recover the factor model in Eq. (2) from a given sample covariance or correlation matrix $\widehat{\Sigma }$. From a practical perspective, it is important not only to obtain point estimates of factor loadings and uniquenesses, but also to assess the precision of these estimates and to identify substantively meaningful—often rotated—factor structures. Such assessments require knowledge of the asymptotic covariance matrices (acov) of the estimators, as well as an understanding of how these covariances are affected by the chosen rotation criterion.

The asymptotic covariance of the estimators depends on both the factor-extraction method and the covariance of $\widehat{\Sigma }$. Let n denote the sample size, and suppose that the estimator of the loading matrix can be expressed as $\Lambda =f\left(\Sigma \right)$ for a given extraction method f. Assume that f is differentiable at $\Sigma$, with differential $df$, which defines a linear mapping from the space of $p\times p$ symmetric matrices to the space of $p\times k$ matrices. By the delta method (see Rao [2009] [17]), the difference between $n^{1/2}d\Lambda$ and $df\left(n^{1/2}d\Sigma \right)$ converges to zero in probability as $n\to \infty$. In this setting, $d\Lambda \approx \widehat{\Lambda }-\Lambda$ and $d\Sigma \approx \widehat{\Sigma }-\Sigma$ denote deviations from the unobservable population parameters $\Sigma$ and $\Lambda$, respectively.

Let $\Sigma =\left[{\sigma }_{ij}\right]$ and define $d{\sigma }_{ij}={\widehat{\sigma }}_{ij}-{\sigma }_{ij}$, where ${\widehat{\sigma }}_{ij}$ denotes the sample covariance (or correlation) between $X_i$ and $X_j$. The asymptotic covariance between elements of the estimated loading matrix $\widehat{\Lambda }$ can then be expressed as

$$\mathrm{acov}({\widehat{\lambda }}_{ir},{\widehat{\lambda }}_{js})=\sum _{m,n,x,y=1}^p{\displaystyle \frac{\partial {\lambda }_{ir}}{\partial {\sigma }_{mx}}}\mathrm{acov}({\widehat{\sigma }}_{mx},{\widehat{\sigma }}_{ny}){\displaystyle \frac{\partial {\lambda }_{js}}{\partial {\sigma }_{ny}}}.$$

Similarly, let $\psi ={({\psi }_1,\dots ,{\psi }_p)}^{\top }$ denote the diagonal elements of the uniqueness matrix $\mathrm{\Psi }$. The asymptotic covariance of the uniqueness estimators can be written as

$$\mathrm{acov}({\widehat{\psi }}_i,{\widehat{\psi }}_j)=\sum _{m,n,x,y=1}^p{\displaystyle \frac{\partial {\psi }_i}{\partial {\sigma }_{mx}}}\mathrm{acov}({\widehat{\sigma }}_{mx},{\widehat{\sigma }}_{ny}){\displaystyle \frac{\partial {\psi }_j}{\partial {\sigma }_{ny}}}.$$

In this paper, we derive explicit expressions for the differentials $d\Lambda$ and $d\psi$ in terms of $d\Sigma$ for several widely used factor extraction methods. The motivation for this research is twofold. First, given the covariance of $\widehat{\Sigma }$, our results enable practitioners to compute standard errors for estimated factor loadings $\widehat{\Lambda }$ and uniquenesses $\widehat{\psi }$. When the data are drawn from a multivariate normal population, for example, the covariance of $\widehat{\Sigma }$ is well-known (see Girshick [1939] [6]). Second, delta-method-based inference for rotated factor loadings (e.g., Archer and Jennrich [1973] [2]; Jennrich [1973] [10]; Yung and Hayashi [2001] [21]) requires these differentials through the chain rule of differentiation. At present, readily available standard error formulas exist primarily for maximum-likelihood factor analysis, as implemented in software such as SAS [2015] [19], reflecting a general property of maximum-likelihood estimators.

When the assumption of normality is violated, the results developed here can be combined with existing asymptotic results for sample covariances or correlations (e.g., Hsu [1949] [9]; de Leeuw [1983] [4]) to obtain more general asymptotic covariance formulas for unrotated factor loadings. While previous studies have examined the asymptotic properties of principal component analysis (e.g., Girshick [1939] [6]; Anderson [1963] [1]) and derived standard error formulas for maximum-likelihood factor analysis (Lawley [1967] [14], Jennrich and Thayer [1973] [11]), our focus is on least-squares, iterative principal component, principal, alpha, and image factor solutions, as implemented in SAS [2015] [19].

Despite the availability of maximum-likelihood (ML) factor analysis, non–ML extraction methods remain widely used due to their computational simplicity, robustness, distribution-free restrictions, and broad implementation in standard statistical software (e.g., Lin et al. [2023] [16]; Raykov et al. [2026] [18]). In many empirical applications, factor analysis is applied to binary, categorical, or qualitative data that may be high-frequency, high-dimensional, or heavy-tailed. Although recent research has emphasized robustness, high-dimensional asymptotics, and software-based inference for rotated solutions, explicit closed-form asymptotic covariance formulas for non-ML factor extraction methods remain largely unavailable. This gap substantially limits formal statistical inference and model comparison. The present paper addresses this limitation by deriving analytically tractable asymptotic covariance expressions, thereby enabling delta-method inference under both normal and non-normal settings.

Our approach is based on multivariate implicit differentiation of the estimating equations that define the factor loadings $\widehat{\Lambda }$ and uniquenesses $\widehat{\Sigma }$ under each extraction method, following techniques widely used in factor-analysis literature (e.g., Archer and Jennrich [1973] [2]; Girshick [1939] [6]; Jennrich [1973] [10]; Yung and Hayashi [2001] [21]). We assess the accuracy of the proposed standard errors for the uniqueness estimator $\widehat{\psi }$ through a simulation study and present an empirical application illustrating their use for the least-squares factor solution.

Accordingly, the main contributions of this paper are as follows:

(i)

Derivation of explicit asymptotic covariance formulas for several non-maximum-likelihood factor solutions;

(ii)

A unified implicit-differentiation framework applicable across multiple factor extraction methods;

(iii)

Closed-form standard error formulas for uniqueness estimators;

(iv)

A practical extension to rotated factor loadings via the delta method; and

(v)

Simulation evidence demonstrating finite-sample accuracy.

2. Iterative Principal Component and Principal Factor Analysis

We begin by introducing notation. For a $p\times k$ matrix $\Lambda$, let $\overrightarrow{\Lambda }$ denote its vectorization obtained by stacking the columns of $\Lambda$, that is,

$$\overrightarrow{\Lambda }={({\lambda }_{11},\dots ,{\lambda }_{p1},{\lambda }_{12},\dots ,{\lambda }_{p2},\dots ,{\lambda }_{1k},\dots ,{\lambda }_{pk})}^{\top }.$$

When no confusion arises, we continue to use double subscripts to refer to individual elements of $\overrightarrow{\Lambda }$. In particular, the $(i,r)$-th element of $\Lambda$, denoted by ${\lambda }_{ir}$, corresponds to the $\left(\right(r-1)p+i)$-th element of $\overrightarrow{\Lambda }$. For a square matrix M, let $\mathrm{diag}\left(M\right)$ denote the vector formed by the diagonal elements of M, and let $\mathrm{Diag}\left(M\right)$ denote the diagonal matrix whose diagonal is given by $\mathrm{diag}\left(M\right)$. We also use ${\delta }_{ij}$ to denote the Kronecker delta, which equals 1 if $i=j$ and 0 otherwise.

Both iterative principal component factor analysis and principal factor analysis construct the loading matrix $\Lambda$ from the principal eigenvectors of $\Sigma -\mathrm{\Psi }$:

$$\Sigma -\mathrm{\Psi }=\Lambda {\Lambda }^{\top }=\sum _{r=1}^k{\lambda }_r{\lambda }_r^{\top },$$

(4)

where ${\lambda }_r$ denotes the r-th column of $\Lambda =[{\lambda }_1,{\lambda }_2,\dots ,{\lambda }_k]$. These columns are ordered according to the eigenvalues of $\Sigma -\mathrm{\Psi }$, with ${\lambda }_r$ corresponding to the r-th largest eigenvalue. Without loss of generality, we assume that these eigenvalues are distinct, which is typically the case for a sample variance-covariance matrix. Consequently, the inner product between any two distinct columns of $\Lambda$ is zero. In particular, Eq. (4) implies

$${\sigma }_{ii}-{\psi }_i=\sum _{t=1}^k{\lambda }_{it}^2,\forall i=1,...,p.$$

(5)

For the sample variance-covariance matrix $\widehat{\Sigma }$, Eq. (4) may, in general, include additional items of the form ${\widehat{\lambda }}_r{\widehat{\lambda }}_r^{\top }$ for $r>k$. These items are omitted in accordance with the model specification in Eq. (2).

Using Eq. (4), we obtain

$$\left(\Sigma -\mathrm{\Psi }\right){\lambda }_r=\Lambda {\Lambda }^{\top }{\lambda }_r=\sum _{s=1}^k{\lambda }_s{\lambda }_s^{\top }{\lambda }_r=\sum _{s=1}^k{\lambda }_s\left({\lambda }_s^{\top }{\lambda }_r\right)=\left({\lambda }_r^{\top }{\lambda }_r\right){\lambda }_r.$$

(6)

This result holds for any $1\le r\le k$. Now fix indices $1\le i\le p$ and $1\le r\le k$. Applying Eq. (5) to the i-th components of both sides of Eq. (6) yields

$$\sum _{j\ne i}{\sigma }_{ij}{\lambda }_{jr}+{\lambda }_{ir}\sum _{t=1}^k{\lambda }_{it}^2={\lambda }_{ir}\sum _{j=1}^p{\lambda }_{jr}^2.$$

Taking differentials on both sides, we obtain

$$\sum _{j\ne i}\left({\sigma }_{ij}d{\lambda }_{jr}+{\lambda }_{jr}d{\sigma }_{ij}\right)+\left(\sum _{t=1}^k{\lambda }_{it}^2\right)d{\lambda }_{ir}+2{\lambda }_{ir}\sum _{t=1}^k{\lambda }_{it}d{\lambda }_{it}=\left(\sum _{j=1}^p{\lambda }_{jr}^2\right)d{\lambda }_{ir}+2{\lambda }_{ir}\sum _{j=1}^p{\lambda }_{jr}d{\lambda }_{jr}.$$

After algebraic simplification, this expression reduces to

$$\sum _{j\ne i}\left({\sigma }_{ij}-2{\lambda }_{ir}{\lambda }_{jr}\right)d{\lambda }_{jr}+\left(\sum _{t=1}^k{\lambda }_{it}^2-\sum _{j=1}^p{\lambda }_{jr}^2\right)d{\lambda }_{ir}+2{\lambda }_{ir}\sum _{t\ne r}{\lambda }_{it}d{\lambda }_{it}=-\sum _{j\ne i}{\lambda }_{jr}d{\sigma }_{ij}.$$

Collecting all linear equations derived above, the system can be written compactly in matrix form as

$$Ad\overrightarrow{\Lambda }=Bd\overrightarrow{\Sigma },$$

where the coefficient of $d{\lambda }_{jt}$ in the $(i,r)$-th equation is given by

$$A_{ir,jt}=\left\{\begin{array}{cc}\sum _{s=1}^k{\lambda }_{is}^2-\sum _{s=1}^p{\lambda }_{sr}^2,\hfill & \mathrm{if}j=i\mathrm{and}t=r,\hfill \\ 2{\lambda }_{ir}{\lambda }_{it},\hfill & \mathrm{if}j=i\mathrm{but}t\ne r,\hfill \\ {\sigma }_{ij}-2{\lambda }_{ir}{\lambda }_{jr},\hfill & \mathrm{if}j\ne i\mathrm{but}t=r,\hfill \\ 0,\hfill & \mathrm{otherwise},\hfill \end{array}\right.$$

for all $1\le i,j\le p$ and $1\le r,t\le k$. Similarly, the coefficient of $d{\sigma }_{xy}$ in the same equation is

$$B_{ir,xy}=-{\delta }_{ix}(1-{\delta }_{xy}){\lambda }_{yr},\forall 1\le i,x,y\le p,1\le r\le k.$$

Finally, we obtain

$$d\overrightarrow{\Lambda }=A^{-1}Bd\overrightarrow{\Sigma }$$

and, by the delta method, the asymptotic covariance matrix of $\overrightarrow{\widehat{\Lambda }}$ is given by

$$\mathrm{acov}\left(\overrightarrow{\widehat{\Lambda }}\right)=\left(A^{-1}B\right)\mathrm{acov}\left(\overrightarrow{\widehat{\Sigma }}\right){\left(A^{-1}B\right)}^{\top }.$$

3. Least-Squares Factor Analysis

In least-squares factor analysis (see Harman [1976] [7]), the objective is to minimize the sum of squared off-diagonal elements of the residual covariance matrix $\widehat{\Sigma }-\widehat{\Lambda }{\widehat{\Lambda }}^{\top }$. Define the objective function

$$\mathrm{g}(\Lambda )=\sum _{1\le y<x\le p}{\left({\sigma }_{xy}-\sum _{t=1}^k{\lambda }_{xt}{\lambda }_{yt}\right)}^2,$$

that is, the sum of squares of the strictly lower-triangular elements of $\Sigma -\Lambda {\Lambda }^{\top }$.

The first-order optimality conditions yield $p\times k$ constraints on the loading matrix$\Lambda$. Specifically, we require

$${\displaystyle \frac{\partial \mathrm{g}(\Lambda )}{\partial {\lambda }_{ir}}}=0$$

for all $1\le i\le p$ and $1\le r\le k$. Straightforward differentiation gives

$$\begin{array}{ccc}\hfill -{\displaystyle \frac{1}{2}}{\displaystyle \frac{\partial \mathrm{g}}{\partial {\lambda }_{ir}}}(\Lambda )& =& {\displaystyle \sum _{1\le y<x\le p}}\left({\sigma }_{xy}-{\displaystyle \sum _{t=1}^k}{\lambda }_{xt}{\lambda }_{yt}\right){\displaystyle \sum _{t=1}^k}({\delta }_{ix}{\delta }_{rt}{\lambda }_{yt}+{\delta }_{iy}{\delta }_{rt}{\lambda }_{xt})\hfill \\ & =& {\displaystyle \sum _{1\le y<x\le p}}({\delta }_{ix}{\lambda }_{yr}+{\delta }_{iy}{\lambda }_{xr})\left({\sigma }_{xy}-{\displaystyle \sum _{t=1}^k}{\lambda }_{xt}{\lambda }_{yt}\right)\hfill \\ & =& {\displaystyle \sum _{1\le y<i}}{\lambda }_{yr}\left({\sigma }_{iy}-{\displaystyle \sum _{t=1}^k}{\lambda }_{it}{\lambda }_{yt}\right)+{\displaystyle \sum _{i<x\le p}}{\lambda }_{xr}\left({\sigma }_{xi}-{\displaystyle \sum _{t=1}^k}{\lambda }_{xt}{\lambda }_{it}\right)\hfill \\ & =& {\displaystyle \sum _{j\ne i}}{\lambda }_{jr}\left({\sigma }_{ij}-{\displaystyle \sum _{t=1}^k}{\lambda }_{it}{\lambda }_{jt}\right).\hfill \end{array}$$

For any fixed $1\le i\le p$ and $1\le r\le k$, the first-order condition reduces to

$$\sum _{j\ne i}{\lambda }_{jr}{\sigma }_{ij}=\sum _{j\ne i}{\lambda }_{jr}\sum _{t=1}^k{\lambda }_{it}{\lambda }_{jt}.$$

(7)

The differential form of Eq. (7) is given by

$$\sum _{j\ne i}\left({\lambda }_{jr}d{\sigma }_{ij}+{\sigma }_{ij}d{\lambda }_{jr}\right)=\sum _{j\ne i}\left(\sum _{t=1}^k{\lambda }_{it}{\lambda }_{jt}\right)d{\lambda }_{jr}+\sum _{j\ne i}{\lambda }_{jr}\sum _{t=1}^k\left({\lambda }_{jt}d{\lambda }_{it}+{\lambda }_{it}d{\lambda }_{jt}\right).$$

Grouping terms and simplifying yield

$$\sum _{j\ne i}{\lambda }_{jr}d{\sigma }_{ij}=\sum _{j\ne i}\left(\sum _{t=1}^k{\lambda }_{it}{\lambda }_{jt}-{\sigma }_{ij}+{\lambda }_{jr}{\lambda }_{ir}\right)d{\lambda }_{jr}+\sum _{t=1}^k\left(\sum _{j\ne i}{\lambda }_{jr}{\lambda }_{jt}\right)d{\lambda }_{it}+\sum _{j\ne i}\sum _{t\ne r}{\lambda }_{jr}{\lambda }_{it}d{\lambda }_{jt},$$

which completes the differential characterization of the first-order condition.

We rewrite the preceding system of equations in matrix form as

$$-Bd\overrightarrow{\Sigma }=Dd\overrightarrow{\Lambda },$$

where the coefficient of $d{\lambda }_{jt}$ in the $(i,r)$-th equation is defined as

$$D_{ir,jt}=\left\{\begin{array}{cc}\sum _{z\ne i}{\lambda }_{zr}{\lambda }_{zt},\hfill & \mathrm{if}j=i,\hfill \\ \sum _{z=1}^k{\lambda }_{iz}{\lambda }_{jz}-{\sigma }_{ij}+{\lambda }_{ir}{\lambda }_{jr},\hfill & \mathrm{if}j\ne i\mathrm{and}t=r,\hfill \\ {\lambda }_{jr}{\lambda }_{it},\hfill & \mathrm{if}j\ne i\mathrm{and}t\ne r.\hfill \end{array}\right.$$

These definitions apply for all $1\le i,j\le p$ and $1\le r,t\le k$. Finally, solving for $d\overrightarrow{\Lambda }$ yields

$$d\overrightarrow{\Lambda }=-D^{-1}Bd\overrightarrow{\Sigma }.$$

4. Alpha Factor Analysis

Let H be the diagonal matrix whose i-th diagonal element is $\sqrt{h_i}$. Clearly, $h_i>0$; otherwise, Eq. (3) would imply that the i-th row of the loading matrix $\Lambda$ is a zero vector, rendering the corresponding factor redundant. Moreover, $H^{-1}(\Sigma -\mathrm{\Psi })H^{-1}$ is the correlation matrix for $\Lambda F$ and Eq (3) implies that

$$d{h}_i=2\sum _{t=1}^k{\lambda }_{it}d{\lambda }_{it}.$$

(8)

To obtain the spectral decomposition of this correlation matrix, let ${\theta }_r$ denote its r-th largest eigenvalue, and let ${\beta }_r$ be the corresponding eigenvector with length ${\theta }_r^{1/2}$. Thus,

$${\theta }_r={\beta }_r^{\top }{\beta }_r,$$

and

$$H^{-1}(\Sigma -\mathrm{\Psi })H^{-1}{\beta }_r={\theta }_r{\beta }_r,\forall 1\le r\le k.$$

In alpha factor analysis (Kaiser and Caffrey [1965] [13]), the loading matrix $\Lambda =[{\lambda }_1,{\lambda }_2,\dots ,{\lambda }_k]$ is defined as

$$\Lambda =H[{\beta }_1,{\beta }_2,\dots ,{\beta }_k].$$

Additionally, let $\Gamma =[{\gamma }_1,{\gamma }_2,\dots ,{\gamma }_k]$, where ${\gamma }_r=H^{-1}{\beta }_r$. Then, for each $r=1,\dots ,k$, the following relationships—corresponding to Eqs. (9)-(11)—hold:

$${\lambda }_r=H{\beta }_r=H^2{\gamma }_r,$$

(9)

$${\theta }_r={\beta }_r^{\top }{\beta }_r={\gamma }_r^{\top }{H}^2{\gamma }_r={\gamma }_r^{\top }{\lambda }_r,$$

(10)

$$(\Sigma -\mathrm{\Psi }){\gamma }_r=(\Sigma -\mathrm{\Psi })H^{-1}{\beta }_r=H{H}^{-1}(\Sigma -\mathrm{\Psi })H^{-1}{\beta }_r={\theta }_{r}H{\beta }_r={\gamma }_r^{\top }{\lambda }_r{\lambda }_r.$$

(11)

In the remainder of this section, we first express $d\overrightarrow{\Gamma }$ in terms of $d\overrightarrow{\Lambda }$ using Eq. (9). We then derive a matrix equation relating $d\overrightarrow{\Sigma }$, $d\overrightarrow{\Gamma }$, and $d\overrightarrow{\Lambda }$ from Eq. (11).

From Eq. (9), for any $1\le i\le p$ and $1\le r\le k$, we have

$${\lambda }_{ir}=h_i{\gamma }_{ir},$$

whose differential is given by

$$d{\lambda }_{ir}=h_{i}d{\gamma }_{ir}+2{\gamma }_{ir}\sum _{t=1}^k{\lambda }_{it}d{\lambda }_{it}.$$

Collecting these expressions yields the following matrix representation:

$$d\overrightarrow{\Gamma }=Ed\overrightarrow{\Lambda },$$

(12)

where the elements of the matrix E are defined by

$$E_{ir,jt}={\displaystyle \frac{\partial {\gamma }_{ir}}{\partial {\lambda }_{jt}}}={\delta }_{ij}({\delta }_{rt}-2{\gamma }_{ir}{\lambda }_{it})/h_i,$$

for all $1\le i,j\le p$ and $1\le r,t\le k$.

For any fixed $1\le i\le p$ and $1\le r\le k$, the i-th component of Eq. (11) can be written as

$$\sum _{j\ne i}{\sigma }_{ij}{\gamma }_{jr}+h_i{\gamma }_{ir}={\lambda }_{ir}\sum _{j=1}^p{\lambda }_{jr}{\gamma }_{jr}.$$

(13)

Taking differentials of Eq. (13) yields

$$\sum _{j\ne i}\left({\sigma }_{ij}d{\gamma }_{jr}+{\gamma }_{jr}d{\sigma }_{ij}\right)+h_{i}d{\gamma }_{ir}+{\gamma }_{ir}d{h}_i={\lambda }_{ir}\sum _{j=1}^p\left({\lambda }_{jr}d{\gamma }_{jr}+{\gamma }_{jr}d{\lambda }_{jr}\right)+\left(\sum _{j=1}^p{\lambda }_{jr}{\gamma }_{jr}\right)d{\lambda }_{ir}.$$

Substituting Eq. (8) into the above expression and rearranging terms, we obtain

$$\sum _{j\ne i}{\gamma }_{jr}d{\sigma }_{ij}=\sum _{j\ne i}({\lambda }_{ir}{\lambda }_{jr}-{\sigma }_{ij})d{\gamma }_{jr}+({\lambda }_{ir}^2-h_i)d{\gamma }_{ir}+\left(\sum _{j\ne i}{\lambda }_{jr}{\gamma }_{jr}\right)d{\lambda }_{ir}-2{\gamma }_{ir}\sum _{t\ne r}{\lambda }_{it}d{\lambda }_{it}+{\lambda }_{ir}\sum _{j\ne i}{\gamma }_{jr}d{\lambda }_{jr}.$$

We express all the differential forms above in matrix notation as

$$Fd\overrightarrow{\Sigma }=Gd\overrightarrow{\Gamma }+Jd\overrightarrow{\Lambda }.$$

(14)

The coefficient matrices F, G, and J are defined as follows.

• Matrix F: The coefficient of $d{\sigma }_{xy}$ in the $(i,r)$-th equation of Eq. (14) is

  $$F_{ir,xy}={\delta }_{ix}(1-{\delta }_{iy}){\gamma }_{yr}.$$
• Matrix G: The coefficient of $d{\gamma }_{jt}$ is given by

  $$G_{ir,jt}=\left\{\begin{array}{cc}{\lambda }_{ir}^2-h_i,\hfill & \mathrm{if}j=i\mathrm{and}t=r,\hfill \\ {\lambda }_{ir}{\lambda }_{jr}-{\sigma }_{ij},\hfill & \mathrm{if}t=r\mathrm{but}j\ne i,\hfill \\ 0,\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$
• Matrix J: The coefficient of $d{\lambda }_{jt}$ is

  $$J_{ir,jt}=\left\{\begin{array}{cc}{\displaystyle \sum _{s\ne i}}{\lambda }_{sr}{\gamma }_{sr},\hfill & \mathrm{if}j=i\mathrm{and}t=r,\hfill \\ -2{\gamma }_{ir}{\lambda }_{it},\hfill & \mathrm{if}j=i\mathrm{but}t\ne r,\hfill \\ {\lambda }_{ir}{\gamma }_{jr},\hfill & \mathrm{if}t=r\mathrm{but}j\ne i,\hfill \\ 0,\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

Combining Eqs. (12) and (14), we obtain

$$Fd\overrightarrow{\Sigma }=Gd\overrightarrow{\Gamma }+Jd\overrightarrow{\Lambda }=G\left(Ed\overrightarrow{\Lambda }\right)+Jd\overrightarrow{\Lambda }=\left(GE+J\right)\overrightarrow{\Lambda }.$$

Consequently,

$$d\overrightarrow{\Lambda }={(GE+J)}^{-1}Fd\overrightarrow{\Sigma }.$$

5. Image Factor Analysis

Let $\Delta ={\left(\mathrm{Diag}\left({\Sigma }^{-1}\right)\right)}^{-1}$. In image factor analysis (see Jöreskog [1969] [12]), the covariance matrix $\Sigma$ is modeled as

$$\Sigma =\Lambda {\Lambda }^{\top }+\tau \Delta ,$$

(15)

where $\tau >0$ is an unknown scalar parameter. Within this framework, one may write $\mathrm{\Psi }=\tau \Delta$. Jöreskog [1969] [12] proposed a maximum-likelihood procedure for estimating this model, under which the covariance matrix of the estimators follows directly from standard maximum-likelihood theory.

For the principal-factor procedure (see SAS [2015] [19]), the loading matrix $\Lambda$ is obtained from the k largest eigenvalues and their corresponding eigenvectors of $\Sigma -\tau \Delta$. The parameter $\tau$ is simultaneously estimated as the coefficient in the regression of $\mathrm{diag}(\Sigma -\Lambda {\Lambda }^{\top })$ on $\mathrm{diag}(\Delta )$, namely,

$$\tau ={\displaystyle \frac{\langle \mathrm{diag}(\Delta ),\mathrm{diag}(\Sigma -\Lambda {\Lambda }^{\top })\rangle }{\langle \mathrm{diag}(\Delta ),\mathrm{diag}(\Delta )\rangle }}$$

(16)

where $\langle ·,·\rangle$ denotes the standard inner product.

Let ${\sigma }^{xy}$ denote the element at the x-th row and y-th column of ${\Sigma }^{-1}$. Since $\Sigma {\Sigma }^{-1}$ is the $p\times p$ identity matrix, differentiating this identity using the product rule yields

$$d{\Sigma }^{-1}=-{\Sigma }^{-1}\left(d\Sigma \right){\Sigma }^{-1}.$$

In particular, for each $i=1,...,p$,

$$d{\sigma }^{ii}=-\sum _{x,y=1}^p{\sigma }^{ix}\left(d{\sigma }_{xy}\right){\sigma }^{yi}=-\sum _{x,y=1}^p{\sigma }^{ix}{\sigma }^{iy}d{\sigma }_{xy}.$$

Consequently,

$$d\left({\displaystyle \frac{1}{{\sigma }^{ii}}}\right)=-{\displaystyle \frac{1}{{\left({\sigma }^{ii}\right)}^2}}d{\sigma }^{ii}=\sum _{x,y=1}^p{\displaystyle \frac{{\sigma }^{ix}{\sigma }^{iy}}{{\left({\sigma }^{ii}\right)}^2}}d{\sigma }_{xy}.$$

(17)

Eq. (16) implies that

$$\tau \langle \mathrm{diag}(\Delta ),\mathrm{diag}(\Delta )\rangle =\langle \mathrm{diag}(\Delta ),\mathrm{diag}(\Sigma -\Lambda {\Lambda }^{\top })\rangle ,$$

and therefore

$$\tau \sum _{i=1}^p{\left({\displaystyle \frac{1}{{\sigma }^{ii}}}\right)}^2=\sum _{i=1}^p{\displaystyle \frac{1}{{\sigma }^{ii}}}({\sigma }_{ii}-h_i).$$

Differentiating both sides yields

$$\left[\sum _{i=1}^p{\left({\displaystyle \frac{1}{{\sigma }^{ii}}}\right)}^2\right]d\tau +\sum _{i=1}^p{\displaystyle \frac{2\tau }{{\sigma }^{ii}}}d\left({\displaystyle \frac{1}{{\sigma }^{ii}}}\right)=\sum _{i=1}^p\left[{\displaystyle \frac{1}{{\sigma }^{ii}}}(d{\sigma }_{ii}-d{h}_i)+({\sigma }_{ii}-h_i)d\left({\displaystyle \frac{1}{{\sigma }^{ii}}}\right)\right].$$

Substituting Eq. (17) into the above expression, we obtain

$$\left[\sum _{i=1}^p{\left({\displaystyle \frac{1}{{\sigma }^{ii}}}\right)}^2\right]d\tau =\sum _{x,y=1}^p\left[\sum _{i=1}^p\left({\sigma }_{ii}-h_i-{\displaystyle \frac{2\tau }{{\sigma }^{ii}}}\right){\displaystyle \frac{{\sigma }^{ix}{\sigma }^{iy}}{{\left({\sigma }^{ii}\right)}^2}}\right]d{\sigma }_{xy}+\sum _{i=1}^p{\displaystyle \frac{1}{{\sigma }^{ii}}}\left(d{\sigma }_{ii}-2\sum _{t=1}^k{\lambda }_{it}d{\lambda }_{it}\right).$$

These expressions can be written compactly as

$$d\tau ={\mu }^{\top }d\overrightarrow{\Sigma }+{\eta }^{\top }d\overrightarrow{\Lambda },$$

(18)

where $\mu$ and $\eta$ are vectors of dimensions $p^2$ and $pk$, respectively. The coefficient of $d{\sigma }_{xy}$ is given by

$${\mu }_{xy}={\displaystyle \frac{\frac{{\delta }_{xy}}{{\sigma }^{xx}}+{\displaystyle \sum _{i=1}^p}\left({\sigma }_{ii}-h_i-\frac{2\tau }{{\sigma }^{ii}}\right)\frac{{\sigma }^{ix}{\sigma }^{iy}}{{\left({\sigma }^{ii}\right)}^2}}{{\displaystyle \sum _{i=1}^p}{\left(\frac{1}{{\sigma }^{ii}}\right)}^2}}$$

and the coefficient of $d{\lambda }_{ir}$ is

$${\eta }_{ir}=-{\displaystyle \frac{2{\lambda }_{ir}}{{\sigma }^{ii}{\displaystyle \sum _{z=1}^p}{\left(\frac{1}{{\sigma }^{zz}}\right)}^2}}.$$

Next, let ${\theta }_r$ denote the r-th largest eigenvalue of $\Sigma -\tau \Delta$, and let ${\lambda }_r$ be the associated eigenvector with length ${\theta }_r^{1/2}$. Then,

$${\lambda }_r^{\top }{\lambda }_r={\theta }_r,$$

and

$$\left(\Sigma -\tau \Delta \right){\lambda }_r={\theta }_r{\lambda }_r={\lambda }_r^{\top }{\lambda }_r{\lambda }_r.$$

For any fixed $1\le i\le p$, the i-th component of this equation is given by

$$\sum _{j=1}^p{\sigma }_{ij}{\lambda }_{jr}-{\displaystyle \frac{\tau {\lambda }_{ir}}{{\sigma }^{ii}}}={\lambda }_{ir}\sum _{j=1}^p{\lambda }_{jr}^2.$$

Differentiating both sides with respect to the underlying parameters yields

$$\sum _{j=1}^p\left({\sigma }_{ij}d{\lambda }_{jr}+{\lambda }_{jr}d{\sigma }_{ij}\right)-{\displaystyle \frac{\tau }{{\sigma }^{ii}}}d{\lambda }_{ir}-{\displaystyle \frac{{\lambda }_{ir}}{{\sigma }^{ii}}}d\tau -\tau {\lambda }_{ir}d\left({\displaystyle \frac{1}{{\sigma }^{ii}}}\right)=2{\lambda }_{ir}\sum _{j=1}^p{\lambda }_{jr}d{\lambda }_{jr}+\left(\sum _{j=1}^p{\lambda }_{jr}^2\right)d{\lambda }_{ir}.$$

Applying Eq. (17) and collecting terms involving $d\Sigma ,d\Gamma$, and $d\tau$, we obtain

$$\sum _{j=1}^p{\lambda }_{jr}d{\sigma }_{ij}-\tau {\lambda }_{ir}\sum _{x,y=1}^p{\displaystyle \frac{{\sigma }^{ix}{\sigma }^{iy}}{{\left({\sigma }^{ii}\right)}^2}}d{\sigma }_{xy}=\sum _{j\ne i}(2{\lambda }_{ir}{\lambda }_{jr}-{\sigma }_{ij})d{\lambda }_{jr}+\left(2{\lambda }_{ir}^2-{\sigma }_{ii}+{\displaystyle \frac{\tau }{{\sigma }^{ii}}}+\sum _{j=1}^p{\lambda }_{jr}^2\right)d{\lambda }_{ir}+{\displaystyle \frac{{\lambda }_{ir}}{{\sigma }^{ii}}}d\tau .$$

The matrix representation of the above system of equations can be written as

$$Ld\overrightarrow{\Sigma }=Pd\overrightarrow{\Lambda }+\pi d\tau$$

(19)

where L, P, and $\pi$ are appropriately defined coefficient matrices. In particular, the coefficient of $d{\sigma }_{xy}$ in the $(i,r)$-th equation is given by

$$L_{ir,xy}={\delta }_{ix}{\lambda }_{yr}-{\displaystyle \frac{\tau {\lambda }_{ir}{\sigma }^{ix}{\sigma }^{iy}}{{\left({\sigma }^{ii}\right)}^2}}.$$

The coefficient of $d{\lambda }_{jt}$ is

$$P_{ir,jt}=\left\{\begin{array}{cc}2{\lambda }_{ir}^2-{\sigma }_{ii}+{\displaystyle \frac{\tau }{{\sigma }^{ii}}}+\sum _{z=1}^p{\lambda }_{zr}^2,\hfill & \mathrm{if}j=i\mathrm{and}t=r,\hfill \\ 2{\lambda }_{ir}{\lambda }_{jr}-{\sigma }_{ij},\hfill & \mathrm{if}t=r\mathrm{but}j\ne i,\hfill \\ 0,\hfill & \mathrm{otherwise},\hfill \end{array}\right.$$

and the coefficient of $d\tau$ is

$${\pi }_{ir}={\displaystyle \frac{{\lambda }_{ir}}{{\sigma }^{ii}}},$$

for any $1\le i,j,x,y\le p$ and $1\le r,t\le k$. Note that $\pi$ is a vector of dimension $pk$.

Combining Eqs. (18) and (19) yields

$$d\overrightarrow{\Lambda }={\left(P+\pi {\eta }^{\top }\right)}^{-1}\left(L-\pi {\mu }^{\top }\right)d\overrightarrow{\Sigma }.$$

Substituting this expression back into Eq. (18) gives the differential of $\tau$ in terms of $d\overrightarrow{\Sigma }$ :

$$d\tau =\left[{\mu }^{\top }+{\eta }^{\top }{\left(P+\pi {\eta }^{\top }\right)}^{-1}\left(L-\pi {\mu }^{\top }\right)\right]d\overrightarrow{\Sigma }.$$

(20)

6. Standard Error of Uniqueness Estimator $\widehat{\psi }$

Assume that $d\overrightarrow{\Lambda }=Md\overrightarrow{\Sigma }$ for some $pk\times p^2$ matrix M. To derive the standard errors of the uniqueness estimator $\widehat{\psi }$, we begin with the generic relationship from Eq. (5): ${\psi }_i={\sigma }_{ii}-\sum _{t=1}^k{\lambda }_{it}^2.$ Differentiating both sides yields

$$d{\psi }_i=d{\sigma }_{ii}-2\sum _{t=1}^k{\lambda }_{it}d{\lambda }_{it},i=1,2,\dots ,p.$$

(21)

Eq. (21) can be written more compactly in matrix form as

$$d\psi =Td\overrightarrow{\Sigma }+Zd\overrightarrow{\Lambda },$$

(22)

where the coefficient of $d{\sigma }_{xy}$ in the i-th equation of Eq. (22) is given by

$$T_{i,xy}={\delta }_{ix}{\delta }_{iy}$$

and the coefficient of $d{\lambda }_{jt}$ in the same equation is

$$Z_{i,jt}=-2{\delta }_{ij}{\lambda }_{jt},$$

for all $1\le i,j,x,y\le p$ and $1\le t\le k$. Consequently, we obtain

$$d\psi =(T+ZM)d\overrightarrow{\Sigma }.$$

When $\Sigma$ is a correlation matrix, $d{\sigma }_{ii}=0$, and the expression simplifies further to $d\psi =ZMd\overrightarrow{\Sigma }$.

In image factor analysis, the relationship ${\widehat{\psi }}_i={\widehat{\sigma }}_{ii}-\sum _{t=1}^k{\widehat{\lambda }}_{it}^2$ may no longer hold. Instead, we rewrite Eq. (20) as $d\tau ={\varphi }^{\top }d\overrightarrow{\Sigma }$ for some vector $\varphi$. Using the identity ${\psi }_i={\displaystyle \frac{\tau }{{\sigma }^{ii}}}$ and applying Eq. (17), we obtain

$$d{\psi }_i={\displaystyle \frac{d\tau }{{\sigma }^{ii}}}+\tau d\left({\displaystyle \frac{1}{{\sigma }^{ii}}}\right)={\displaystyle \frac{1}{{\sigma }^{ii}}}\sum _{x,y=1}^p{\varphi }_{xy}d{\sigma }_{xy}+\tau \sum _{x,y=1}^p{\displaystyle \frac{{\sigma }^{ix}{\sigma }^{iy}}{{\left({\sigma }^{ii}\right)}^2}}d{\sigma }_{xy},$$

(23)

for all $1\le i\le p$. In matrix form, this expression can be written compactly as

$$d\psi =Qd\overrightarrow{\Sigma },$$

where the coefficient of $d{\sigma }_{xy}$ in the i-th equation is given by

$$Q_{i,xy}={\displaystyle \frac{\partial {\psi }_i}{\partial {\sigma }_{xy}}}={\displaystyle \frac{{\sigma }^{ii}{\varphi }_{xy}+\tau {\sigma }^{ix}{\sigma }^{iy}}{{\left({\sigma }^{ii}\right)}^2}},\forall 1\le i,x,y\le p.$$

7. An Empirical and Simulation Study

In this section, we apply the least-squares factor analysis method to a sample correlation matrix consisting of nine variables measured on 211 subjects. The data were originally reported by Lawley and Maxwell [1971, p.43] [15], and the sample is assumed to be drawn from an approximately multivariate normal population. For simplicity, we set the number of common factors to $k=2$. This choice is supported by both Horn [1965] [8] ’s parallel analysis and Velicer [1976] [20] ’s minimum average partial (MAP) test. Estimates of the factor loading matrix $\widehat{\Lambda }$ and the uniqueness vector $\widehat{\psi }$ are obtained using the iterative procedure described by Harman [1976] [7].

To compute the asymptotic covariances of these estimators, we first apply Girshick [1939] [6]’s formula (3.23) to derive the covariances of the sample correlation coefficients. We then use the results developed in Section 3 and Section 6 to compute the asymptotic covariance matrix of $\widehat{\Lambda }$, as well as the standard errors of $\widehat{\psi }$. To illustrate the effect of factor rotation, we orthogonally rotate $\widehat{\Lambda }$ using the varimax criterion and apply the method of Archer and Jennrich [1973] [11] to obtain the standard errors of the rotated factor loadings. Table 1 reports the parameter estimates for the two-factor least-squares model; standard errors are reported in parentheses.

To verify the theoretical derivation of the asymptotic covariance, $\mathrm{acov}\left(\widehat{\psi }\right)$, presented in Section 6, we conduct a simulation study. The sample correlation matrix reported above is treated as the population correlation matrix $\Sigma$. Using this matrix, we generate a large number m of random correlation matrices from a Wishart distribution with parameter $(\Sigma ,211)$. For each simulated correlation matrix, the uniqueness parameters are estimated using the least-squares method with $k=2$.

Based on the resulting collection of simulated estimates of $\psi$, we compute empirical standard errors for the uniqueness estimates. For $m=50,100,500,1000,2000$, the corresponding empirical standard errors are reported in Table 2. The final column of the table presents the theoretical asymptotic standard errors computed using the formulas developed in Section 6. The results demonstrate clear convergence of the empirical standard errors toward their theoretical counterparts, although the convergence occurs at a relatively slow rate.

8. Discussions

Recent years have seen renewed interest in formal statistical inference for factor-analytic models. This resurgence is driven by a range of practical considerations, including applications involving non-normal data, concerns about robustness, and the growing complexity of covariance structures encountered in empirical work. Broadly speaking, the contemporary literature can be grouped into three main strands. The first strand emphasizes likelihood-based approaches to factor analysis and covariance structure modeling, building on classical frameworks and their subsequent extensions (e.g., Lawley [1967] [14]; Jennrich and Thayer [1973] [11]; Raykov et al. [2026] [18]). A second line of research focuses on high-dimensional and dynamic factor models, in which both the cross-sectional and temporal dimensions grow, posing new challenges for estimation and inference (e.g., Bai [2003] [3]; Fan et al. [2021] [5]). The third strand concentrates on software-oriented implementations of standard error estimation for rotated factor solutions, typically relying on numerical differentiation or sandwich-type approximations (e.g., Archer and Jennrich [1973] [2]; Yung and Hayashi [2001] [21]).

In contrast, explicit asymptotic covariance formulas for widely used non-maximum-likelihood extraction method—such as least-squares, principal, alpha, and image factor analysis—remain relatively scarce, despite their continued prevalence in empirical research and standard statistical software (e.g., Harman [1976] [7]; SAS [2015] [19]). The results developed in this paper help bridge this gap by providing closed-form, analytically tractable asymptotic covariance expressions for several classical no-ML factor solutions. By formulating these results within a unified implicit-differentiation framework, the proposed approach accommodates non-normal data and can be readily combined with existing delta-method results for factor rotation (e.g., Archer and Jennrich [1973] [2]; Jennrich [1973] [10]; Yung and Hayashi [2001] [21]). From a practical perspective, this framework enables transparent computation of standard errors without reliance on numerical differentiation or resampling. From a methodological standpoint, it complements recent advances in robust and high-dimensional factor modeling by providing a rigorous inferential foundation for classical extraction methods that remain central to applied factor analysis.

In the preceding sections, the derived asymptotic covariances of $\widehat{\Lambda }$ and $\widehat{\psi }$ are conditional on the assumed number of factors. Consequently, abnormal standard errors may reflect model misspecification or systematic errors—such as an incorrect choice of factor dimensionality—rather than genuine sampling variability. In particular, such distortions may arise from either over-factoring (when k is chosen too large) or under-factoring (when k is chosen too small).

For expositional simplicity, we have not explicitly imposed the symmetry condition ${\sigma }_{ij}={\sigma }_{ji}$. In practical implementations, however, exploiting this property is important not only for reducing memory requirements but also for improving computational efficiency, especially when working with large-scale datasets.

The matrix $\Sigma$ need not represent the variance-covariance matrix of a single group of observed data series; it may also correspond to the variance-covariance matrix of residuals in a system of equations. As an illustration, consider a structural vector autoregressive (SVAR) model with p time series:

$$Y_t=c+\sum _{s=1}^{\zeta }{A}_s{Y}_{t-s}+\Lambda F_t+{ϵ}_t,$$

where $Y_t$ is a p-dimensional vector of observed time series at time t, $\zeta$ denotes the lag length, $F_t$ represents k common structural shocks, and ${ϵ}_t$ captures p idiosyncratic shocks. The composite disturbance $\Lambda F_t+{ϵ}_t$ therefore exhibits a two-component error structure, while the common factors F in Eq. (1) are latent and may possess a nonzero mean.

In much of the SVAR literature, idiosyncratic shocks are frequently ignored, effectively imposing the restriction $k=p$. Allowing for idiosyncratic shocks, however, substantially reduces both the number of structural shocks and the number of identification conditions required. From a computational perspective, incorporating an error-component structure—such as a common-factor model—into SVAR systems is straightforward. The formulas developed in this paper offer practical guidance for implementing such extensions and for assessing the precision of estimated coefficients associated with structural shocks.

This form of dimension reduction is particularly valuable in environments characterized by highly interconnected data. In macroeconomic modeling, structural shocks are typically assumed to be independent. Although an economy may generate a large number of interrelated observable time series, its dynamics are generally driven by only a small number of underlying shocks, such as demand, supply, and technological disturbances.

Finally, an important advantage of the derived formulas is their applicability under non-normality without requiring modification. For example, when the data are multivariate binary, the asymptotic covariance $\mathrm{acov}\left(\widehat{\Lambda }\right)$ remains well defined. More generally, for categorical data—which are common in practical factor-analysis applications—the computation of $\mathrm{acov}\left(\widehat{\Lambda }\right)$ is unlikely to pose significant difficulties.