2.1 Maximum likelihood estimation

Let ${\boldsymbol {X}}=(X_1,\ldots ,X_p)^\top $ be a p-dimensional observed random vector. The factor analysis model is

$$ \begin{align*} {\boldsymbol{X}} =\boldsymbol{\mu} + \Lambda{\boldsymbol{f}} + {\boldsymbol{\varepsilon}} , \end{align*} $$

where $\boldsymbol {\mu }$ is a mean vector, $\Lambda =(\lambda _{ij})$ is a $p \times k$ matrix of factor loadings, and ${\boldsymbol {f}} = (f_1,\ldots ,f_k)^\top $ and ${\boldsymbol {\varepsilon }} = (\varepsilon _1,\ldots , \varepsilon _p)^\top $ are unobservable random vectors. Here, $A^\top $ denotes the transpose of a matrix A. The elements of $\boldsymbol {f}$ and $\boldsymbol {\varepsilon }$ are referred to as common factors and unique factors, respectively. It is assumed that ${\boldsymbol {f}} \sim N_k(\boldsymbol {0}, I_k)$ and ${\boldsymbol {\varepsilon }} \sim N_p(\boldsymbol {0}, \Psi )$ and ${\boldsymbol {f}}$ and ${\boldsymbol {\varepsilon }}$ are independent. Here, $I_k$ is an identity matrix of order k and $\Psi $ is a $p \times p$ diagonal matrix with the i-th diagonal element being $\psi _i$ , referred to as a unique variance. Under these assumptions, the observed vector ${\boldsymbol {X}}$ follows multivariate-normal distribution; ${\boldsymbol {X}} \sim N_p(\boldsymbol {\mu }, \Sigma )$ with $\Sigma =\Lambda \Lambda ^\top +\Psi $ . When the observed variables are standardized, the i-th diagonal element of $ \Lambda \Lambda ^\top $ is called communality, which stands for the proportion of variance in $X_i$ explained by common factors.

It is well-known that factor loadings have a rotational indeterminacy since both $\Lambda $ and $\Lambda T$ generate the same covariance matrix $\Sigma $ , where T is an arbitrary orthonormal matrix. To obtain an identifiable model, we often use a restriction that $\Lambda ^\top \Psi ^{-1}\Lambda $ is a diagonal matrix or the upper triangular elements of $\Lambda $ are fixed by 0.

Suppose that we have a random sample of N observations ${\boldsymbol {x}}_1,\ldots ,{\boldsymbol {x}}_N$ from the p-dimensional normal population $N_p(\boldsymbol {\mu } ,\Sigma )$ . The MLEs of factor loadings and unique variances are obtained by maximizing the log-likelihood function

(1) $$ \begin{align} \log {f}(\boldsymbol{x}_1,\dots,\boldsymbol{x}_N | \Lambda,\Psi ) = -\frac{N}{2} \bigg\{ p\log(2\pi)+\log |\Sigma| + \mathrm{tr}(\Sigma^{-1} S ) \bigg\}, \end{align} $$

or equivalently, minimizing the discrepancy function:

(2) $$ \begin{align} {q}( \Lambda,\Psi ) = \log|\Sigma|+ \mathrm{tr}(\Sigma^{-1}S) -\log |S|-p, \end{align} $$

where $S=(s_{ij})$ is the sample variance–covariance matrix,

$$ \begin{align*} S=\frac{1}{N} \sum_{n=1}^N({\boldsymbol{x}}_n-\bar{{\boldsymbol{x}}})({\boldsymbol{x}}_n-\bar{{\boldsymbol{x}}})^\top, \quad \bar{\boldsymbol{x}} = \frac{1}{N}\sum_{n=1}^N\boldsymbol{x}_n. \end{align*} $$

We remark that it is common to use the unbiased sample covariance matrix with $n = N-1$ , and the choice between N and $N-1$ slightly affects the statistical properties of the estimator. In this study, however, the focus is on computing estimates rather than analyzing estimator properties. Since factor analysis typically uses the sample correlation matrix, the choice between N and $N-1$ does not affect the discrepancy function.

We assume that S is non-singular to ensure the existence of $S^{-1}$ , which plays an important role in constructing an algorithm to find MLEs. This assumption can be violated when the number of variables exceeds the sample size. In such a case, we may employ the ridge regularization on S; that is, we use $S_{\lambda }:=S+\lambda I_p$ with $\lambda>0$ (Yuan & Chan, Reference Yuan and Chan2008).

2.2 Numerical algorithms

We briefly review the algorithms for computing the MLEs (Jöreskog, Reference Jöreskog1967; Lawley & Maxwell, Reference Lawley and Maxwell1971). The candidates of the MLEs, say $(\hat {\Lambda },\hat {\Psi })$ , are given as the solutions of

(3) $$ \begin{align} \frac{\partial {q}(\Lambda,\Psi ) }{ \partial \Lambda} = 0,~~ \frac{\partial {q}(\Lambda,\Psi) }{ \partial \psi_i} = 0. \quad (i=1,\ldots, p) \end{align} $$

The solution to the above equation is given by

(4) $$ \begin{align} \begin{aligned} \Sigma^{-1}(\Sigma-S)\Sigma^{-1}\Lambda &=0,\\ \mathrm{diag}(\Sigma^{-1}(\Sigma-S)\Sigma^{-1}) &=0, \end{aligned} \end{align} $$

where $\mathrm {{diag}}(A)$ denotes the diagonal matrix consisting of the diagonal elements of a matrix A.

Since the solutions are not expressed in a closed form, some iterative algorithm is required. Typically, we use an algorithm based on the gradient of the log-likelihood function, such as the Newton–Raphson method or Quasi-Newton method. When $\Psi $ is not singular, Eq. (3) results in the following expressions (Lawley & Maxwell, Reference Lawley and Maxwell1971, page 28, Eq. (4.10)):

(5) $$ \begin{align} (\Psi^{-\frac{1}{2}} S \Psi^{-\frac{1}{2}}) \Psi^{-\frac{1}{2}} \Lambda&= \Psi^{-\frac{1}{2}} \Lambda(I_k +\Lambda^\top\Psi^{-1}\Lambda), \end{align} $$

(6) $$ \begin{align}\kern38pt \Psi&=\mathrm{diag}(S-\Lambda\Lambda^\top). \end{align} $$

Let $\Lambda ^\top \Psi ^{-1}\Lambda =\Delta $ . The rotational indeterminacy is resolved when $\Delta $ is a diagonal matrix (e.g., Lawley & Maxwell, Reference Lawley and Maxwell1971); in this case, (5) is an eigenvalue/eigenvector problem. Let $(\theta _1,\boldsymbol {e}_1),\ldots ,(\theta _p,\boldsymbol {e}_p)$ be the eigenvalues and eigenvectors of $\Psi ^{-\frac {1}{2}} S \Psi ^{-\frac {1}{2}}$ , respectively, where the eigenvalues are rearranged in a weakly decreasing order; that is, $\theta _1 \geq \theta _2 \geq \cdots \geq \theta _p$ . Let $P =(\boldsymbol {e}_1,\ldots ,\boldsymbol {e}_k)$ and $\Theta =\mathrm {diag}(\theta _1,\theta _2, \ldots , \theta _k)$ . Eq. (5) implies

(7) $$ \begin{align} P (\Theta -I_k )^{\frac{1}{2}}=\Psi^{-\frac{1}{2}} \Lambda \end{align} $$

or equivalently,

(8) $$ \begin{align} \Lambda&= \Psi^{\frac{1}{2}} P (\Theta -I_k )^{\frac{1}{2}}. \end{align} $$

Thus, for given $\Psi $ , $\Lambda $ can be calculated using (8). For given $\Lambda $ , one can update $\Psi $ with (6). However, the update with (6) can be slow when some of the diagonal elements of $\hat {\Psi }$ are small. Thus, the unique variances are often updated by using the Newton–Raphson method or Quasi-Newton method (Jöreskog, Reference Jöreskog1967; Lawley & Maxwell, Reference Lawley and Maxwell1971). The factanal function in R uses the optim function to update the unique variances.

To summarize, the estimation algorithm in Jöreskog (Reference Jöreskog1967) is given as follows:

1. 1. Set an initial value of $\Psi $ .

2. 2. For given $\Psi $ , Eq. (8) implies that $\Lambda $ is a function of $\Psi $ ; that is, $\Lambda $ can be written as $\Lambda =\Lambda (\Psi )$ . Substituting $\Lambda (\Psi )$ into the discrepancy function (2), it is a function of $\Psi $ . The elements of $\Psi $ are then optimized by Newton–Raphson method or Quasi-Newton method.

2.3 A system of algebraic equations for the MLEs

In this study, we will employ a computational algebraic approach to get exact algebraic solutions to the system (4) (i.e., all candidates of the MLEs). However, Eq. (4) may not be directly used as it includes an inverse matrix $\Sigma ^{-1}$ that consists of significantly complicated rational functions. Consequently, a simpler system of algebraic equations is required. One may use (5) and (6) as they do not include $\Sigma ^{-1}$ , but Eq. (5) includes $\Psi ^{-1}$ , which is unavailable when $\Psi $ is singular. As will be demonstrated in the numerical example in Section 4, there exist many estimates whose unique variances are exactly zeros. Thus, it is essential to derive a system of algebraic equations that does not include complicated rational functions and also satisfies (4), even when $\Psi $ is singular.

In the following theorem, we provide a necessary condition for (4) that does not have rational functions and holds even when $\Psi $ is singular. Therefore, all solutions of (4) must satisfy that necessary condition.

Theorem 2.1. Whether $\Psi $ is singular or not, if (4) is satisfied, then we have

(9) $$ \begin{align} \Lambda&=(\Lambda\Lambda^\top + \Psi) S^{-1}\Lambda, \end{align} $$

(6) $$ \begin{align} \kern-6pt \Psi&=\mathrm{diag}(S-\Lambda\Lambda^\top). \end{align} $$

Proof. It is shown that Eq. (9) is equivalent to the upper part of Eq. (4), since we have

$$ \begin{align*} \Sigma^{-1}(\Sigma-S)\Sigma^{-1}\Lambda =0 &\Longleftrightarrow (\Sigma S^{-1} \Sigma) \Sigma^{-1}(\Sigma-S)\Sigma^{-1}\Lambda =0 \\ &\Longleftrightarrow \Sigma S^{-1}\Lambda=\Lambda \Longleftrightarrow (\Lambda\Lambda^\top + \Psi) S^{-1}\Lambda=\Lambda. \end{align*} $$

Therefore, what we need to prove is that Eq. (4) implies Eq. (6). First, we have

$$ \begin{align*} \Sigma - S = (\Lambda\Lambda^\top + \Psi) \Sigma^{-1} (\Sigma - S) \Sigma^{-1} (\Lambda\Lambda^\top + \Psi) \end{align*} $$

by $\Sigma = \Lambda \Lambda ^\top + \Psi $ . Hence, the upper part of (4) implies that

$$ \begin{align*} \Sigma - S = \Psi \Sigma^{-1} (\Sigma - S) \Sigma^{-1} \Psi. \end{align*} $$

To handle the case where $\Psi $ is singular, suppose that $\psi _1, \ldots , \psi _q$ are not equal to zero and $\psi _{q+1}, \ldots , \psi _p$ are equal to zero, where $q \leq p$ . Let $J = \{ 1, \ldots , q\}, K = \{q+1,\ldots ,p\}$ , and let

$$ \begin{align*} T = \mathrm{diag}\left(\frac{1}{\psi_1}, \ldots, \frac{1}{\psi_q}, 1, \ldots, 1\right), ~~ I_{q \leq p} = \mathrm{diag}\left(\overbrace{1, \ldots, 1}^{q}, \overbrace{0, \ldots, 0}^{p-q}\right). \end{align*} $$

Let $M_{JK}$ denote a submatrix of a matrix M from which the rows corresponding to the indices in J and the columns corresponding to the indices in K are selected, respectively. Since T is not singular,

$$ \begin{align*} \Sigma - S = \Psi \Sigma^{-1} (\Sigma - S) \Sigma^{-1} \Psi \Longleftrightarrow T(\Sigma - S) T &= T \Psi \Sigma^{-1} (\Sigma - S) \Sigma^{-1} \Psi T \\ &= I_{q \leq p} \Sigma^{-1} (\Sigma - S) \Sigma^{-1} I_{q \leq p} \\ &= \begin{pmatrix} (\Sigma^{-1} (\Sigma - S) \Sigma^{-1})_{JJ} & 0 \\ 0 & 0 \end{pmatrix}. \end{align*} $$

As $(T (\Sigma - S) T)_{JJ} = (\Sigma ^{-1} (\Sigma - S) \Sigma ^{-1})_{JJ}$ holds, the lower part of (4) implies that $\mathrm {diag}((T (\Sigma - S) T)_{JJ}) = 0$ . Hence, $\Psi _{JJ} = \mathrm {diag}((S - \Lambda \Lambda ^\top )_{JJ})$ holds because T is diagonal. Since $(\Sigma - S)_{KK} = I_{p-q} (\Sigma - S)_{KK} I_{p-q} = (T (\Sigma - S) T)_{KK} = 0$ holds, we also have $s_{ii} - \sum _{j=1}^k \lambda _{ij}^2 = 0 = \psi _i$ for each $i \in K$ . Consequently, we obtain Eq. (6) without an assumption that $\Psi $ is not singular.

Remark 1. As an alternative algorithm to obtain the MLEs given by Lawley & Maxwell (Reference Lawley and Maxwell1971), Jennrich & Robinson (Reference Jennrich and Robinson1969) used the eigenvalue/eigenvector problem of $S^{-\frac {1}{2}} \Psi S^{-\frac {1}{2}}$ instead of $\Psi ^{-\frac {1}{2}} S \Psi ^{-\frac {1}{2}}$ by using the following equation:

(10) $$ \begin{align} (S^{-\frac{1}{2}} \Psi S^{-\frac{1}{2}}) S^{-\frac{1}{2}} \Lambda&= S^{-\frac{1}{2}} \Lambda(I_k -\Lambda^\top S^{-1}\Lambda), \end{align} $$

which is equivalent to Eq. (9). Hence, Jennrich and Robinson’s (Reference Jennrich and Robinson1969) algorithm can produce zero or negative unique variance estimates.

3.1 Gröbner bases and systems of algebraic equations

3.1.1 Affine varieties and systems of algebraic equations

Let $\mathbb {R}[\boldsymbol {z}]$ denote the set of all polynomials in m variables $\boldsymbol {z} = (z_1, \ldots , z_m)$ with coefficients in the real number field $\mathbb {R}$ . Recall that $\mathbb {R}[\boldsymbol {z}]$ is a commutative ring, which is called the polynomial ring in variables $\boldsymbol {z}$ with coefficients in $\mathbb {R}$ . For details of fields and commutative rings, please refer to Appendix A. Let $f_1, \ldots , f_r$ be polynomials in the polynomial ring $\mathbb {R}[\boldsymbol {z}]$ . We consider a system of algebraic equations which has the following general form:

(11) $$ \begin{align} f_1(\boldsymbol{z}) = \cdots = f_r(\boldsymbol{z}) = 0. \end{align} $$

Define

$$ \begin{align*} \mathbb{V}_{W}(f_1, \ldots, f_r) = \{ \boldsymbol{z} \in W^m : f_1(\boldsymbol{z}) = \cdots = f_r(\boldsymbol{z}) = 0\} \subset W^m, \end{align*} $$

where $W = \mathbb {R}$ or $\mathbb {C}$ . The set $ \mathbb {V}_{W}(f_1, \ldots , f_r)$ is referred to as an affine variety defined by $f_1, \ldots , f_r$ in $W^{m}$ (Cox et al., Reference Cox, Little and O’Shea2015, Section 1.2, Definition 1). Note that the solutions to (11) in $W^m$ form the affine variety $\mathbb {V}_{W}(f_1,\ldots , f_r)$ , that is, the solution space to (11) in $W^m$ coincides with the affine variety $\mathbb {V}_{W}(f_1,\ldots , f_r)$ .

Here, let us consider an illustrative example related to the system of algebraic equations (6) and (9).

Example 1. Let

$$ \begin{align*} S = \begin{pmatrix} 1 & \frac{1}{2} & \frac{1}{3} \\ \frac{1}{2} & 1 & \frac{2}{3} \\ \frac{1}{3} & \frac{2}{3} & 1 \end{pmatrix}, ~~ \Lambda = \begin{pmatrix} \lambda_{11} \\ \lambda_{21} \\ \lambda_{31} \end{pmatrix}, ~~ \Psi = \begin{pmatrix} \psi_1 & 0 & 0 \\ 0 & \psi_2 & 0 \\ 0 & 0 & \psi_3 \end{pmatrix}. \end{align*} $$

We consider the system of algebraic equations (6) and (9) for the above $S, \Lambda , \Psi $ , that is,

(12) $$ \begin{align} \left\{\begin{array}{l} 0 = (\Psi - \mathrm{diag}(S - \Lambda\Lambda^\top))_{11} = (\Psi - \mathrm{diag}(S - \Lambda\Lambda^\top))_{22} = (\Psi - \mathrm{diag}(S - \Lambda\Lambda^\top))_{33}, \\ 0 = (\Lambda - (\Lambda\Lambda^\top + \Psi)S^{-1}\Lambda)_{11} = (\Lambda - (\Lambda\Lambda^\top + \Psi)S^{-1}\Lambda)_{21} = (\Lambda - (\Lambda\Lambda^\top + \Psi)S^{-1}\Lambda)_{31}. \end{array}\right. \end{align} $$

Because we have

$$ \begin{align*} S^{-1} = \begin{pmatrix} \frac{4}{3} & -\frac{2}{3} & 0 \\ -\frac{2}{3} & \frac{32}{15} & -\frac{6}{5} \\ 0 & -\frac{6}{5} & \frac{9}{5} \end{pmatrix}, \end{align*} $$

we consider polynomials $f_1, \ldots , f_6 \in \mathbb {R}[\boldsymbol {\psi }, \boldsymbol {\lambda }] = \mathbb {R}[\psi _1, \psi _2, \psi _3, \lambda _{11},\lambda _{21},\lambda _{31}]$ defined by

$$ \begin{align*} f_1 &= (\Psi - \mathrm{diag}(S - \Lambda\Lambda^\top))_{11} = \psi_1 - (1 - \lambda_{11}^2), \\ f_2 &= (\Psi - \mathrm{diag}(S - \Lambda\Lambda^\top))_{22} = \psi_2 - (1 - \lambda_{21}^2), \\ f_3 &= (\Psi - \mathrm{diag}(S - \Lambda\Lambda^\top))_{33} = \psi_3 - (1 - \lambda_{31}^2), \\ f_4 &= (\Lambda - (\Lambda\Lambda^\top + \Psi)S^{-1}\Lambda)_{11} \\ &= \lambda_{11} - \left(\begin{array}{l} \lambda_{11} \left(- \frac{2}{3} \lambda_{11} \lambda_{21} + \frac{4}{3} (\lambda_{11}^2 + \psi_1)\right) \\ + \; \lambda_{21} \left(\frac{32}{15} \lambda_{11} \lambda_{21} - \frac{6}{5} \lambda_{11} \lambda_{31} - \frac{2}{3} (\lambda_{11}^2 + \psi_1)\right) + \; \lambda_{31} \left(-\frac{6}{5} \lambda_{11} \lambda_{21} + \frac{9}{5}\lambda_{11} \lambda_{31}\right) \end{array}\right), \\ f_5 &= (\Lambda - (\Lambda\Lambda^\top + \Psi)S^{-1}\Lambda)_{21} \\ &= \lambda_{21} - \left(\begin{array}{l} \lambda_{11} \left(\frac{4}{3} \lambda_{11} \lambda_{21} - \frac{2}{3} (\lambda_{21}^2 + \psi_2)\right) \\ + \; \lambda_{21} \left(-\frac{2}{3} \lambda_{11} \lambda_{21} - \frac{6}{5} \lambda_{21} \lambda_{31} + \frac{32}{15} (\lambda_{21}^2 + \psi_2)\right) + \; \lambda_{31} \left(\frac{9}{5} \lambda_{21} \lambda_{31} - \frac{6}{5} (\lambda_{21}^2 + \psi_2)\right) \end{array}\right), \\ f_6 &= (\Lambda - (\Lambda\Lambda^\top + \Psi)S^{-1}\Lambda)_{31} \\ &= \lambda_{31} - \left(\begin{array}{l} \lambda_{11} \left(\frac{4}{3} \lambda_{11} \lambda_{31} - \frac{2}{3} \lambda_{21} \lambda_{31}\right) \\ + \; \lambda_{21} \left(-\frac{2}{3} \lambda_{11} \lambda_{31} + \frac{32}{15} \lambda_{21} \lambda_{31} - \frac{6}{5} (\lambda_{31}^2 + \psi_3)\right) + \; \lambda_{31} \left(-\frac{6}{5} \lambda_{21} \lambda_{31} + \frac{9}{5} (\lambda_{31}^2 + \psi_3)\right) \end{array}\right). \end{align*} $$

All solutions of the above equation form the affine variety defined by $f_1,\ldots ,f_6$ . Indeed, there are seven real solutions to (12), and these solutions, $(\boldsymbol {\psi }, \boldsymbol {\lambda }) = (\psi _1, \psi _2, \psi _3, \lambda _{11}, \lambda _{21}, \lambda _{31}) \in \mathbb {R}^6$ , form the following affine variety:

$$ \begin{align*} &\quad\mathbb{V}_{\mathbb{R}}(f_1,\ldots,f_6) \\ &= \left\{ \left(0, \frac{3}{4}, \frac{8}{9}, \pm 1, \pm \frac{1}{2}, \pm \frac{1}{3}\right), \left(\frac{3}{4}, 0, \frac{5}{9}, \pm \frac{1}{2}, \pm 1, \pm \frac{2}{3} \right), \left(\frac{8}{9}, \frac{5}{9}, 0, \pm \frac{1}{3}, \pm \frac{2}{3}, \pm 1\right), \left(1, 1, 1, 0, 0, 0\right) \right\}. \end{align*} $$

3.1.2 Sums and saturations

In Example 1, we obtain only seven solutions to the system of algebraic equations (12). For a larger model, however, the number of solutions of (6) and (9) is considerably large and the solution space (i.e., the affine variety) is too complicated. In some cases, there are an infinite number of solutions. Therefore, it is essential to enhance an efficient computation for obtaining all solutions. A key idea is a decomposition of (6) and (9) into several simple sub-problems. This decomposition is related to the theory of polynomial ideals and Gröbner bases.

We define the following set of polynomials:

$$ \begin{align*} \langle f_1, \ldots, f_r \rangle = \left\{\sum_{i=1}^r q_i f_i : q_i \in \mathbb{R}[\boldsymbol{z}]\right\}, \end{align*} $$

where $f_1, \ldots , f_r \in \mathbb {R}[\boldsymbol {z}]$ . The above set is an ideal in the polynomial ring $\mathbb {R}[\boldsymbol {z}]$ . So, it is called the ideal generated by $f_1,. .. ,f_r$ (Cox et al., Reference Cox, Little and O’Shea2015, Section 1.4, Definition 2). Let

$$ \begin{align*} \mathcal{J} = \langle f_1, \ldots, f_r \rangle \subset \mathbb{R}[\boldsymbol{z}]. \end{align*} $$

For details of the ideal in the commutative ring, please refer to Appendix A. The affine variety defined by $\mathcal {J}$ in $W^m$ is expressed as $\mathbb {V}_{W}(\mathcal {J}) = \{ \boldsymbol {z} \in W^m : f(\boldsymbol {z}) = 0, \forall f \in \mathcal {J}\}.$ Since $\mathbb {V}_{W}(f_1, \ldots , f_r) = \mathbb {V}_{W}(\mathcal {J})$ holds, the solution space of (11) coincides with $\mathbb {V}_{W}(\mathcal {J})$ .

Since ideals are algebraic objects, there exist natural algebraic operations. Essentially, such algebraic operations give us some simple sub-problems of the system of algebraic equations (6) and (9). Let $\mathcal {K}$ be an ideal in $\mathbb {R}[\boldsymbol {z}]$ . We note that any ideal in $\mathbb {R}[\boldsymbol {z}]$ is finitely generated by the Hilbert Basis Theorem (Cox et al., Reference Cox, Little and O’Shea2015, Theorem 4, Section 2.5), that is, there exist $h_1, \ldots , h_s \in \mathcal {K}$ such that $\mathcal {K} = \langle h_1, \ldots , h_s \rangle $ .

Now, we define two operations. As the first operation for ideals, we define the sum of $\mathcal {J}$ and $\mathcal {K}$ , say $\mathcal {J} + \mathcal {K}$ , by

$$ \begin{align*} \mathcal{J} + \mathcal{K} = \{ f + h : f \in \mathcal{J}, h \in \mathcal{K} \}. \end{align*} $$

For details of the definition of sums, see (Cox et al., Reference Cox, Little and O’Shea2015, Section 4.3, Definition 1). Note that any sum of ideals is an ideal in $\mathbb {R}[\boldsymbol {z}]$ (Cox et al., Reference Cox, Little and O’Shea2015, Proposition 2, Section 4.3). Moreover, the affine variety $\mathbb {V}_W(\mathcal {J} + \mathcal {K})$ coincides with $\mathbb {V}_W(\mathcal {J}) \cap \mathbb {V}_W(\mathcal {K})$ (Cox et al., Reference Cox, Little and O’Shea2015, Corollary 3, Section 4.3), that is, $\mathbb {V}_W(\mathcal {J} + \mathcal {K})$ is none other than the solution space of the equations

$$ \begin{align*} (11) \text{ and } h_1(\boldsymbol{z}) = \cdots = h_s(\boldsymbol{z}) = 0. \end{align*} $$

As the second operation, we define the saturation of $\mathcal {J}$ with respect to $\mathcal {K}$ , say $\mathcal {J}:\mathcal {K}^{\infty }$ , by

$$ \begin{align*} \mathcal{J}:\mathcal{K}^{\infty} = \{ f \in \mathbb{R}[\boldsymbol{z}] : \text{for any } h \in \mathcal{K} \text{ there exists } t \in \mathbb{Z}_{\geq 0} \text{ such that } fh^t \in \mathcal{J} \}, \end{align*} $$

where $\mathbb {Z}_{\geq 0} = \{ t \in \mathbb {Z} : t \geq 0\}$ , and $fh^t$ means the product of the polynomial f and the t-th power of h denoted by $h^t$ . For details of the definition of saturations, see (Cox et al., Reference Cox, Little and O’Shea2015, Section 4.4, Definition 8). Note that any saturation is an ideal in $\mathbb {R}[\boldsymbol {z}]$ (Cox et al., Reference Cox, Little and O’Shea2015, Proposition 9, Section 4.4). If $W = \mathbb {C}$ , the affine variety $\mathbb {V}_{\mathbb {C}}(\mathcal {J}:\mathcal {K}^{\infty })$ coincides with the Zariski closure of the subset $\mathbb {V}_{\mathbb {C}}(\mathcal {J}) \setminus \mathbb {V}_{\mathbb {C}}(\mathcal {K}) \subseteq \mathbb {C}^m$ (Cox et al., Reference Cox, Little and O’Shea2015, Theorem 10, Section 4.4), that is, the $\mathbb {V}_{\mathbb {C}}(\mathcal {J}:\mathcal {K}^{\infty })$ is none other than the Zariski closure of the solution space over $\mathbb {C}$ of

$$ \begin{align*} \left((11) \text{ and } h_1(\boldsymbol{z}) \neq 0 \right) \text{ or } \cdots \text{ or } \left((11) \text{ and } h_s(\boldsymbol{z}) \neq 0 \right), \end{align*} $$

where $\mathbb {V}_{\mathbb {C}}(\mathcal {J}) \setminus \mathbb {V}_{\mathbb {C}}(\mathcal {K}) \subseteq \mathbb {C}^m$ is defined by

$$ \begin{align*} \mathbb{V}_{\mathbb{C}}(\mathcal{J}) \setminus \mathbb{V}_{\mathbb{C}}(\mathcal{K}) = \{ \boldsymbol{z} \in \mathbb{V}_{\mathbb{C}}(\mathcal{J}) : \boldsymbol{z} \not\in \mathbb{V}_{\mathbb{C}}(\mathcal{K}) \}. \end{align*} $$

Recall that the Zariski closure of a subset in $W^m = \mathbb {R}^m$ or $\mathbb {C}^m$ is the smallest affine variety containing the subset (see Appendix A for details of Zariski closures).

Now, we provide a proposition that gives theoretical correctness to divide Eqs. (6) and (9) into some sub-problems (Cox et al., Reference Cox, Little and O’Shea2015, Theorem 10, Section 4.4).

Proposition 3.1. We have $\mathbb {V}_W(\mathcal {J}) = \mathbb {V}_W(\mathcal {J} + \mathcal {K}) \cup \mathbb {V}_{W}(\mathcal {J}:\mathcal {K}^{\infty })$ .

Proposition 3.1 tells that if we obtain generator sets of the sum $\mathcal {J}+\mathcal {K}$ and the saturation $\mathcal {J}:\mathcal {K}^{\infty }$ , and both of their affine varieties are not empty, we get meaningful sub-problems of (11).

Now, let us turn our attention to the system of algebraic equations for the maximum likelihood factor analysis. The result of Example 1 shows that the affine variety includes many solutions whose $\Psi $ is singular; 6 out of 7 solutions are singular. In fact, as we will see in the numerical example presented in Section 4, the solution space of (6) and (9) includes a surprisingly large number of solutions whose $\Psi $ is singular. Therefore, we propose to decompose the original problem into the following sub-problems: $\psi _j=0$ and $\psi _j\neq 0$ for $j=1,\dots , p$ . Now, we provide a simple example based on Example 1 to illustrate how the decomposition can be achieved. A rigorous procedure will be provided in Section 3.2.

Example 2. We illustrate the usefulness of Proposition 3.1 using the same example as in Example 1. Let us consider the ideal $\mathcal {J} = \langle f_1, \ldots , f_6 \rangle $ in the polynomial ring $\mathbb {R}[\boldsymbol {\psi }, \boldsymbol {\lambda }],$ where $f_1,\ldots ,f_6$ are polynomials given in Example 1. We construct some ideals sequentially as follows: first, by using the ideal $\mathcal {J}$ , we construct the following sum and saturation:

$$ \begin{align*} \mathcal{I}_0 = \mathcal{J} + \langle \psi_1 \rangle, ~~ \mathcal{I}_1 = \mathcal{J} : \langle \psi_1 \rangle^{\infty}. \end{align*} $$

Next, by using the ideals $\mathcal {I}_{0}$ and $\mathcal {I}_{1}$ , we construct the following sums and saturations:

$$ \begin{align*} \mathcal{I}_{00} = \mathcal{I}_0 + \langle \psi_2 \rangle, ~~ \mathcal{I}_{01} = \mathcal{I}_0 : \langle \psi_2 \rangle^{\infty}, ~~ \mathcal{I}_{10} = \mathcal{I}_{1} + \langle \psi_2 \rangle, ~~ \mathcal{I}_{11} = \mathcal{I}_{1} : \langle \psi_2 \rangle^{\infty}. \end{align*} $$

Last, by using the ideals $\mathcal {I}_{00}$ , $\mathcal {I}_{01}$ , $\mathcal {I}_{10},$ and $\mathcal {I}_{11}$ , we construct the following sums and saturations:

$$ \begin{align*} \mathcal{I}_{000} = \mathcal{I}_{00} + \langle \psi_3 \rangle, ~~ \mathcal{I}_{001} = \mathcal{I}_{00} : \langle \psi_3 \rangle^{\infty}, ~~ \mathcal{I}_{010} = \mathcal{I}_{01} + \langle \psi_3 \rangle, ~~ \mathcal{I}_{011} = \mathcal{I}_{01} : \langle \psi_3 \rangle^{\infty}, \\ \mathcal{I}_{100} = \mathcal{I}_{10} + \langle \psi_3 \rangle, ~~ \mathcal{I}_{101} = \mathcal{I}_{10} : \langle \psi_3 \rangle^{\infty}, ~~ \mathcal{I}_{110} = \mathcal{I}_{11} + \langle \psi_3 \rangle, ~~ \mathcal{I}_{111} = \mathcal{I}_{11} : \langle \psi_3 \rangle^{\infty}. \end{align*} $$

It follows from Proposition 3.1 that

$$ \begin{align*} \mathbb{V}_{\mathbb{R}}(\mathcal{J}) &= \mathbb{V}_{\mathbb{R}}(\mathcal{I}_0) \cup \mathbb{V}_{\mathbb{R}}(\mathcal{I}_1) \\ &= \left(\mathbb{V}_{\mathbb{R}}(\mathcal{I}_{00}) \cup \mathbb{V}_{\mathbb{R}}(\mathcal{I}_{01})\right) \cup \left(\mathbb{V}_{\mathbb{R}}(\mathcal{I}_{10}) \cup \mathbb{V}_{\mathbb{R}}(\mathcal{I}_{11})\right) \\ &= \left( \left(\mathbb{V}_{\mathbb{R}}(\mathcal{I}_{000}) \cup \mathbb{V}_{\mathbb{R}}(\mathcal{I}_{001})\right) \cup \left(\mathbb{V}_{\mathbb{R}}(\mathcal{I}_{010}) \cup \mathbb{V}_{\mathbb{R}}(\mathcal{I}_{011})\right) \right) \\ &\quad \cup \left( \left(\mathbb{V}_{\mathbb{R}}(\mathcal{I}_{100}) \cup \mathbb{V}_{\mathbb{R}}(\mathcal{I}_{101})\right) \cup \left(\mathbb{V}_{\mathbb{R}}(\mathcal{I}_{110}) \cup \mathbb{V}_{\mathbb{R}}(\mathcal{I}_{111})\right) \right) = \bigcup_{b \in \{ 0,1\}^3} \mathbb{V}_{\mathbb{R}}(\mathcal{I}_{b}). \end{align*} $$

The affine varieties $\mathbb {V}_{\mathbb {R}}(\mathcal {I}_{000}), \mathbb {V}_{\mathbb {R}}(\mathcal {I}_{100}), \mathbb {V}_{\mathbb {R}}(\mathcal {I}_{010}), \mathbb {V}_{\mathbb {R}}(\mathcal {I}_{001})$ are empty sets. On the other hand,

$$ \begin{align*} \mathbb{V}_{\mathbb{R}}(\mathcal{I}_{011}) &= \left\{ \left(0, \frac{3}{4}, \frac{8}{9}, \pm 1, \pm\frac{1}{2}, \pm\frac{1}{3}\right) \right\}, & \mathbb{V}_{\mathbb{R}}(\mathcal{I}_{101}) &= \left\{ \left(\frac{3}{4}, 0, \frac{5}{9}, \pm\frac{1}{2}, \pm1, \pm\frac{2}{3}\right) \right\}, \\ \mathbb{V}_{\mathbb{R}}(\mathcal{I}_{110}) &= \left\{ \left(\frac{8}{9}, \frac{5}{9}, 0, \pm\frac{1}{3}, \pm\frac{2}{3}, \pm 1\right) \right\}, & \mathbb{V}_{\mathbb{R}}(\mathcal{I}_{111}) &= \left\{ \left(1, 1, 1, 0, 0, 0\right) \right\} \end{align*} $$

are non-empty proper subsets of $\mathbb {V}_{\mathbb {R}}(\mathcal {J})$ . Thus, the solution space of (12) is divided into the four affine varieties: $\mathbb {V}_{\mathbb {R}}(\mathcal {I}_{011})$ , $\mathbb {V}_{\mathbb {R}}(\mathcal {I}_{101})$ , $\mathbb {V}_{\mathbb {R}}(\mathcal {I}_{110})$ , and $\mathbb {V}_{\mathbb {R}}(\mathcal {I}_{111})$ .

3.1.3 Gröbner bases

Example 2 showed that the solution space can be successfully decomposed by constructing sums and saturations. To obtain some simple sub-problems of the system of algebraic equations (6) and (9) for maximum likelihood factor analysis, in addition, we need to compute generator sets of the sum $\mathcal {J} + \mathcal {K}$ and the saturation $\mathcal {J}:\mathcal {K}^{\infty }$ . The sum $\mathcal {J} + \mathcal {K}$ is generated by $f_1, \ldots , f_r, h_1, \ldots , h_s$ , that is, $\mathcal {J} + \mathcal {K} = \langle f_1,\ldots , f_r, h_1,\ldots ,h_s \rangle $ . A generator set of the saturation $\mathcal {J}:\mathcal {K}^{\infty }$ is obtained by using a Gröbner basis. The definition of the Gröbner basis requires a monomial ordering. Therefore, we define monomial orderings before giving the definition of Gröbner bases.

Let us denote $\boldsymbol {z}^{\boldsymbol {a}} = z_1^{a_1} \ldots z_m^{a_m}$ for $\boldsymbol {a}=(a_1,\ldots ,a_m) \in \mathbb {Z}_{\geq 0}^m$ , which is called a monomial, where $\mathbb {Z}_{\geq 0}^m = \{ (\alpha _1,\ldots ,\alpha _m) \in \mathbb {Z}^m : \alpha _1, \ldots , \alpha _m \geq 0\}$ . Let $M(\boldsymbol {z}) = \{ \boldsymbol {z}^{\boldsymbol {a}} : \boldsymbol {a} \in \mathbb {Z}_{\geq 0}^m \}$ . The definition of monomial orderings is the following (Cox et al., Reference Cox, Little and O’Shea2015, Section 2.2, Definition 1).

Definition 1. A monomial ordering $\succ $ on $\mathbb {R}[\boldsymbol {z}]$ is a relation on the set $M(\boldsymbol {z})$ satisfying that

1. 1. $\succ $ is a total ordering on $M(\boldsymbol {z})$ , that is, $\succ $ satisfies that

   1. (a) $s \succeq s$ for any $s \in M(\boldsymbol {z})$ ,

   2. (b) $s \succeq t$ and $t \succeq s$ implies $s = t$ for any $s,t \in M(\boldsymbol {z})$ ,

   3. (c) $s \succeq t$ and $t \succeq u$ implies $s \succeq u$ for any $s,t,u \in M(\boldsymbol {z})$ , and

   4. (d) $s \succ t$ or $s = t$ or $s \prec t$ for any $s,t \in M(\boldsymbol {z})$ ,

2. 2. $\boldsymbol {z}^{\boldsymbol {a}} \succ \boldsymbol {z}^{\boldsymbol {b}} \Longrightarrow \boldsymbol {z}^{\boldsymbol {a}} \boldsymbol {z}^{\boldsymbol {c}} \succ \boldsymbol {z}^{\boldsymbol {b}} \boldsymbol {z}^{\boldsymbol {c}}$ for any $\boldsymbol {c} \in \mathbb {Z}_{\geq 0}^m$ , and

3. 3. any nonempty subset of $M(\boldsymbol {z})$ has a smallest element under $\succ $ .

Definition 2. Fixing a monomial ordering on $\mathbb {R}[\boldsymbol {z}]$ and given a polynomial $f = \sum _{\boldsymbol {a} \in \mathbb {Z}_{\geq 0}^m} c_{\boldsymbol {a}} \boldsymbol {z}^{\boldsymbol {a}}$ for $c_{\boldsymbol {a}} \in \mathbb {R}$ , we call the monomial $\mathrm {LM}(f) = \max (\boldsymbol {z}^{\boldsymbol {a}} : c_{\boldsymbol {a}} \neq 0)$ the leading monomial of f.

We give two examples of monomial orderings. The first one is called a lexicographic order (or the lex order for short) (Cox et al., Reference Cox, Little and O’Shea2015, Section 2.2, Definition 3). The second one is called a graded reverse lexicographic order (or the grevlex order for short) (Cox et al., Reference Cox, Little and O’Shea2015, Section 2.2, Definition 6). These two monomial orderings satisfy the conditions of Definition 1 (Cox et al., Reference Cox, Little and O’Shea2015, Proposition 4 and Exercise 5, Section 2.2).

Definition 3 (Lexicographic Order; Lex).

Let $\boldsymbol {a}, \boldsymbol {b} \in \mathbb {Z}_{\geq 0}^m$ . We say $\boldsymbol {z}^{\boldsymbol {a}} \succ _{\mathrm {lex}} \boldsymbol {z}^{\boldsymbol {b}}$ if the leftmost nonzero entry of the vector difference $\boldsymbol {a} - \boldsymbol {b}$ is positive. This ordering is called a lexicographic order (or the lex order for short).

Although a Gröbner basis with respect to the lex order defined above transforms a given system of algebraic equations into an equivalent system of algebraic equations that is easier to solve, it is widely known in computational algebra that such computations tend to be heavy in general. Therefore, we next define the graded lexicographic order, a monomial order under which Gröbner basis computation is generally known to be more efficient. Of course, a Gröbner basis with respect to the graded reverse lexicographic order also provides an equivalent transformation.

Definition 4 (Graded Reverse Lexicographic Order; Grevlex).

Let $\boldsymbol {a} = (a_1,\ldots , a_m) \in \mathbb {Z}_{\geq 0}^m$ and $\boldsymbol {b} = (b_1,\ldots , b_m) \in \mathbb {Z}_{\geq 0}^m$ . We say $\boldsymbol {z}^{\boldsymbol {a}} \succ _{\mathrm {grevlex}} \boldsymbol {z}^{\boldsymbol {b}}$ if $|\boldsymbol {a}| = \sum _{i=1}^{m} a_i> |\boldsymbol {b}| = \sum _{i=1}^{m} b_i$ or the rightmost nonzero entry of the vector difference $\boldsymbol {a} - \boldsymbol {b}$ is negative. This ordering is called a graded reverse lexicographic order (or the grevlex order for short).

As shown in the following example, the leading monomial of a given polynomial is determined depending on a fixed monomial ordering.

Example 3 (Lexicographic Order; Lex).

When we fix the lex order, $\mathrm {LM}(3 z_1^2 + 2 z_1 z_2^3) = z_1^2$ . In addition, $\mathrm {LM}(f_4) = \psi _1$ , $\mathrm {LM}(f_5) = \psi _2$ , and $\mathrm {LM}(f_6) = \psi _3$ , where $f_4,f_5,f_6 \in \mathbb {R}[\boldsymbol {\psi },\boldsymbol {\lambda }]$ are presented in Example 1.

Example 4 (Graded Reverse Lexicographic Order; Grevlex).

When we fix the grevlex order, $\mathrm {LM}(3 z_1^2 + 2 z_1 z_2^3) = z_1 z_2^3$ . In addition, $\mathrm {LM}(f_4) = \lambda _{11}^3$ , $\mathrm {LM}(f_5) = \lambda _{11}^2 \lambda _{21}$ , and $\mathrm {LM}(f_6) = \lambda _{11}^2 \lambda _{31}$ , where $f_4,f_5,f_6 \in \mathbb {R}[\boldsymbol {\psi },\boldsymbol {\lambda }]$ are presented in Example 1.

Now, we define Gröbner bases (Cox et al., Reference Cox, Little and O’Shea2015, Section 2.5, Definition 5).

Definition 5. Fix a monomial ordering on $\mathbb {R}[\boldsymbol {z}]$ . A finite subset $G = \{ g_1, \ldots , g_t \}$ of an ideal $\mathcal {J} \subseteq \mathbb {R}[\boldsymbol {z}]$ is said to be a Gröbner basis of $\mathcal {J}$ if we have

$$ \begin{align*} \langle \mathrm{LM}(g_1), \ldots, \mathrm{LM}(g_t) \rangle = \langle \mathrm{ \{ LM}(f) : f \in \mathcal{J} \} \rangle. \end{align*} $$

Note that every polynomial ideal has a Gröbner basis (Cox et al., Reference Cox, Little and O’Shea2015, Corollary 6, Section 2.5), and that any Gröbner basis of $\mathcal {J}$ generates the ideal $\mathcal {J}$ .

For the case $\mathcal {K} = \langle h \rangle \subseteq \mathbb {R}[\boldsymbol {z}]$ , a generator set of the saturation $\mathcal {J}:\mathcal {K}^{\infty }$ is obtained by using a Gröbner basis as follows (Cox et al., Reference Cox, Little and O’Shea2015, Theorem 14(ii), Section 4.4):

1. 1. let $\tilde {\mathcal {J}} = \langle f_1, \ldots , f_r, 1 - y h \rangle \subseteq \mathbb {R}[y, \boldsymbol {z}]$ ,

2. 2. compute a Gröbner basis of $\tilde {\mathcal {J}}$ , say $\tilde {G}$ , with respect to the lex order,

3. 3. $G = \tilde {G} \cap \mathbb {R}[\boldsymbol {z}]$ is a Gröbner basis of the saturation $\mathcal {J}:\mathcal {K}^{\infty }$ .

Example 5. Let us reconsider Example 1, which is related to the system of algebraic equations (6) and (9). Gröbner bases of the sums and saturations described above give us simplified sub-problems. For $b \in \{0,1\}^3$ , a Gröbner basis $G_b$ of $\mathcal {I}_{b}$ with respect to the lex order is given by

$$ \begin{align*} G_{000} &= \left\{ 1 \right\},~~G_{001} = \left\{ 1 \right\},~~G_{010} = \left\{ 1 \right\},~~G_{100} = \left\{ 1 \right\}, \\ G_{011} &= \left\{ \psi_1, \psi_2-\frac{3}{4}, \psi_3-\frac{8}{9}, \lambda_{11}-3 \lambda_{31}, \lambda_{21} - \frac{3}{2} \lambda_{31}, \lambda_{31}^2-\frac{1}{9} \right\}, \\ G_{101} &= \left\{ \psi_1-\frac{3}{4}, \psi_2, \psi_3-\frac{5}{9}, \lambda_{11}-\frac{3}{4} \lambda_{31}, \lambda_{21} - \frac{3}{2} \lambda_{31}, \lambda_{31}^2-\frac{4}{9} \right\}, \\ G_{110} &= \left\{ \psi_1-\frac{8}{9}, \psi_2-\frac{5}{9}, \psi_3, \lambda_{11}-\frac{1}{3} \lambda_{31}, \lambda_{21} - \frac{2}{3} \lambda_{31}, \lambda_{31}^2-1 \right\}, \\ G_{111} &= \left\{ \psi_1+\frac{81}{25} \lambda_{31}^2-1, \psi_2-1, \psi_3+\lambda_{31}^2-1, \lambda_{11}+\frac{9}{5} \lambda_{31}, \lambda_{21}, \lambda_{31}^3 + \frac{5}{27} \lambda_{31} \right\}. \end{align*} $$

So the Gröbner bases $G_{011}, G_{101}, G_{110}, G_{111}$ divide (12) into the following simple sub-problems:

$$ \begin{align*} \left\{\begin{array}{l} \psi_1=0, \\ \psi_2=\frac{3}{4}, \\ \psi_3=\frac{8}{9}, \\ \lambda_{11}=3 \lambda_{31}, \\ \lambda_{21}=\frac{3}{2} \lambda_{31}, \\ 0 = \lambda_{31}^2-\frac{1}{9}, \end{array}\right. ~~ \left\{\begin{array}{l} \psi_1=\frac{3}{4}, \\ \psi_2=0, \\ \psi_3=\frac{5}{9}, \\ \lambda_{11}=\frac{3}{4} \lambda_{31}, \\ \lambda_{21} = \frac{3}{2} \lambda_{31}, \\ 0 = \lambda_{31}^2-\frac{4}{9}, \end{array}\right. ~~ \left\{\begin{array}{l} \psi_1=\frac{8}{9}, \\ \psi_2=\frac{5}{9}, \\ \psi_3=0, \\ \lambda_{11}=\frac{1}{3} \lambda_{31}, \\ \lambda_{21} = \frac{2}{3} \lambda_{31}, \\ 0=\lambda_{31}^2-1, \end{array}\right. ~~ \left\{\begin{array}{l} \psi_1=1-\frac{81}{25} \lambda_{31}^2, \\ \psi_2=1, \\ \psi_3=1-\lambda_{31}^2, \\ \lambda_{11}=-\frac{9}{5} \lambda_{31}, \\ \lambda_{21}=0, \\ 0 = \lambda_{31}^3 + \frac{5}{27} \lambda_{31}. \end{array}\right. \end{align*} $$

Gröbner bases have a good property that plays an important role in solving a system of algebraic equations of the form (11). The following proposition shows that Gröbner bases transform (11) into easily solvable problems, without being subject to chance. The following proposition is called the finiteness theorem (Cox et al., Reference Cox, Little and O’Shea2015, Theorem 6, Section 5.3).

Proposition 3.2. Fix a monomial ordering on $\mathbb {R}[\boldsymbol {z}]$ and let G be a Gröbner basis of an ideal $\mathcal {J} \subseteq \mathbb {R}[\boldsymbol {z}]$ . Consider the following four statements:

1. 1. for each $i = 1, \ldots , m$ , there exists $t_i \geq 0$ such that $z_i^{t_i} \in \langle \mathrm {LM}(f) : f \in \mathcal {J} \rangle $ ,

2. 2. for each $i = 1, \ldots , m$ , there exists $t_i \geq 0$ and $g \in G$ such that $\mathrm {LM}(g) = z_i^{t_i}$ ,

3. 3. for each $i = 1, \ldots , m$ , there exists $t_i \geq 0$ and $g \in G$ such that $\mathrm {LM}(g) = z_i^{t_i}$ and $g \in \mathbb {R}[z_i, \ldots , z_m]$ , if we fix the lex order,

4. 4. the affine variety $\mathbb {V}_{W}(\mathcal {J})$ is a finite set.

Then, statements 1–3 are equivalent and they all imply statement 4. An ideal $\mathcal {J}$ satisfying statement 1, 2, or 3 is called a zero-dimensional ideal. Otherwise, it is called a non-zero-dimensional ideal. Furthermore, if $W = \mathbb {C}$ , then statements 1–4 are all equivalent.

The above proposition shows that if we can obtain a Gröbner basis, we can determine whether the affine variety is a finite or infinite set, and if the affine variety is a finite set, the Gröbner basis with respect to the lex order gives a triangulated representation (statement 3 of Proposition 3.2). This representation produces a much simpler system of algebraic equations than the original one; for example, when we need to solve linear equations, its Gröbner basis provides a triangular matrix, corresponding to the triangularization of linear equations. In general, the computation of a Gröbner basis with respect to the grevlex order is more efficient than that with respect to the lex order. Also, if $\mathcal {J}$ is zero-dimensional, we have an efficient algorithm which converts a Gröbner basis from one monomial ordering to another, which is called the FGLM algorithm (Faugère et al., Reference Faugère, Gianni, Lazard and Mora1993). Hence, in our algebraic approach to maximum likelihood factor analysis, we first compute Gröbner bases with respect to the grevlex order, and then convert them to Gröbner bases with respect to the lex order.

Remark 2. Primary ideal decompositions can more completely give us simplified sub-problems of a system of algebraic equations. However, unfortunately, the primary ideal decompositions require heavy computational loads in general. Thus, we employ Proposition 3.1 instead of primary ideal decompositions.

3.2 An algorithm to find all solutions

We provide an algorithm to find all solutions to Eq. (4). To get the solutions, we solve the following equation:

(13) $$ \begin{align} \left\{S-\Lambda\Lambda^{\top}-\mathrm{diag}\left(S-\Lambda\Lambda^{\top}\right)\right\}S^{-1}\Lambda =0, \end{align} $$

which is derived by substituting (6) into (9). Note that, according to Theorem 2.1, the system of algebraic equations (13) is a necessary condition for the likelihood equation (4). Since the solution space of (13) contains that of (4), we can find all solutions to (4) by solving (13). Eq. (13) is much easier to solve than Eq. (4); this is because Eq. (4) includes the inverse matrix $\Sigma ^{-1}$ that consists of rational functions, while Eq. (13) does not include any inverse matrix.

As described in Example 2 in Section 3.1, (13) can have a large number of, and sometimes infinite, solutions when $\Psi $ is singular. Therefore, as in Example 2, we use Proposition 3.1 to get some sub-problems of Eq. (13).

Considering a rotational indeterminacy, we suppose that the upper triangular elements of $\Lambda $ are fixed by zero. Let us consider the following ideal in $\mathbb {R}[\boldsymbol {\lambda }] = \mathbb {R}[\lambda _{ij} : 1 \leq i \leq p, 1 \leq j \leq \min (i,k)]$ :

$$ \begin{align*} \mathcal{J} = \left\langle \left(\{S-\Lambda\Lambda^{\top}-\mathrm{diag}(S-\Lambda\Lambda^{\top})\}S^{-1}\Lambda\right)_{ij} : 1 \leq i \leq p, 1 \leq j \leq k \right\rangle. \end{align*} $$

We have to choose the ideals of the form $\mathcal {K} = \langle h \rangle \subseteq \mathbb {R}[\boldsymbol {\lambda }]$ well to get sub-problems of (13).

Recall that (6) is equivalent to the lower part of (4) when $\Psi $ is not singular, and is a necessary condition of (4) when $\Psi $ is singular, as shown in Theorem 2.1. Focusing on (6), we consider $\psi _1, \ldots , \psi _p$ as polynomials defined by

$$ \begin{align*} \psi_{i} = s_{ii} - \sum_{j=1}^{\min(i,k)} \lambda_{ij}^2 \in \mathbb{R}[\boldsymbol{\lambda}] ~~ (i=1,\ldots,p). \end{align*} $$

For each $\mathcal {K}_{i} = \langle \psi _i \rangle \subseteq \mathbb {R}[\boldsymbol {\lambda }]$ , we construct sums and saturations as in our algorithm.

Now, we propose our algorithm in Algorithm 1. In Step 3 of Algorithm 1, the sum $\mathcal {I}_0 = \mathcal {J} + \mathcal {K}_{1}$ and the saturation $\mathcal {I}_1 = \mathcal {J}:\mathcal {K}_{1}^{\infty }$ are constructed, where $\mathcal {J}, \mathcal {K}_1, \ldots , \mathcal {K}_p$ are given in Steps 1–3 of Algorithm 1. Proposition 3.1 implies that

$$ \begin{align*} \mathbb{V}_{W}(\mathcal{J}) = \mathbb{V}_{W}(\mathcal{I}_0) \cup \mathbb{V}_{W}(\mathcal{I}_1). \end{align*} $$

Second, we construct the following sums and saturations as in $i=1$ of Steps 4–8 of Algorithm 1:

$$ \begin{align*} \mathcal{I}_{00} = \mathcal{I}_{0} + \mathcal{K}_{2}, ~~ \mathcal{I}_{01} = \mathcal{I}_0:\mathcal{K}_{2}^{\infty}, ~~ \mathcal{I}_{10} = \mathcal{I}_{1} + \mathcal{K}_{2}, ~~ \mathcal{I}_{11} = \mathcal{I}_1:\mathcal{K}_{2}^{\infty}. \end{align*} $$

Since Proposition 3.1 implies

$$ \begin{align*} \mathbb{V}_{W}(\mathcal{I}_0) = \mathbb{V}_{W}(\mathcal{I}_{00}) \cup \mathbb{V}_{W}(\mathcal{I}_{01}), ~~ \mathbb{V}_{W}(\mathcal{I}_1) = \mathbb{V}_{W}(\mathcal{I}_{10}) \cup \mathbb{V}_{W}(\mathcal{I}_{11}), \end{align*} $$

we have

$$ \begin{align*} \mathbb{V}_{W}(\mathcal{J}) = \mathbb{V}_{W}(\mathcal{I}_{00}) \cup \mathbb{V}_{W}(\mathcal{I}_{01}) \cup \mathbb{V}_{W}(\mathcal{I}_{10}) \cup \mathbb{V}_{W}(\mathcal{I}_{11}). \end{align*} $$

In the same way, for each $b \in \{ 0, 1\}^{i}$ , $2 \leq i \leq p - 1$ , we sequentially construct

$$ \begin{align*} \mathcal{I}_{b0} = \mathcal{I}_{b} + \mathcal{K}_{i+1}, ~~ \mathcal{I}_{b1} = \mathcal{I}_b:\mathcal{K}_{i+1}^{\infty} \end{align*} $$

as in each $2 \leq i \leq p-1$ of Steps 4–8 of Algorithm 1. By Proposition 3.1, we have

$$ \begin{align*} \mathbb{V}_{W}(\mathcal{J}) = \bigcup_{b \in \{ 0, 1\}^{p}} \mathbb{V}_{W}(\mathcal{I}_{b}). \end{align*} $$

Finally, we compute all solutions to (13) by using the following steps for each $b \in \{ 0, 1\}^{p}$ as in Steps 9–17 of Algorithm 1:

1. 1. compute a Gröbner basis $\tilde {G}_b = \{ \tilde {g}_1, \ldots , \tilde {g}_s\}$ of $\mathcal {I}_{b}$ with respect to the grevlex order as in Step 10 of Algorithm 1,

2. 2. if $\mathcal {I}_{b}$ is zero-dimensional, compute a Gröbner basis $G_b = \{ g_1, \ldots , g_{r}\}$ of $\mathcal {I}_b$ with respect to the lex order by converting $\tilde {G}_b$ to a Gröbner basis with respect to the lex order as in Step 12 of Algorithm 1, and then solve the equation $g_1 = \cdots = g_r = 0$ as in Step 13 of Algorithm 1,

3. 3. else, compute a sample point from each connected component defined by $\tilde {g}_1 = \cdots = \tilde {g}_s = 0$ as in Step 15 of Algorithm 1.

Remark 4. The ideals $\mathcal {I}_b$ may have multiple roots. Therefore, in practice, Step 10 computes a Gröbner basis of the radical ideals for $\mathcal {I}_{b}$ . The radical ideal $\sqrt {\mathcal {I}_b}$ is defined by

$$ \begin{align*} \sqrt{\mathcal{I}_b} = \{ p : \text{there exists }t \in \mathbb{Z}_{\geq 0}\text{ such that }p^t \in \mathcal{I}_b \}. \end{align*} $$

In general, radical ideals have only solutions with multiplicity 1 (i.e., simple roots). So our algebraic approach can treat very simple sub-problems, and does not generate multiple solutions.

3.3 An algorithm to identify a solution pattern

We provide an algorithm to categorize the MLEs into three patterns: “Proper solution,” “Improper solution,” and “No solutions.” The algorithm is detailed in Algorithm 2. First, we find all solutions to Eq. (4) with Algorithm 1. Then, we get a solution ( $\hat {\Lambda }_{t^*}, \hat {\Psi }_{t^*}$ ) that minimizes the discrepancy function in (2) as in Step 1 of Algorithm 2. When the Hessian matrix $H_q$ at the solution ( $\hat {\Lambda }_{t^*}, \hat {\Psi }_{t^*}$ ) is positive definite, there exists an MLE that minimizes (2), where the Hessian matrix $H_q$ is defined by the following:

$$ \begin{align*} H_q = \begin{pmatrix} \frac{\partial^2 q}{\partial \lambda_{11} \partial \lambda_{11}} & \cdots & \frac{\partial^2 q}{\partial \lambda_{11} \partial \lambda_{1m}} & \cdots & \frac{\partial^2 q}{\partial \lambda_{11} \partial \lambda_{p1}} & \cdots & \frac{\partial^2 q}{\partial \lambda_{11} \partial \lambda_{pm}} & \frac{\partial^2 q}{\partial \lambda_{11} \partial \psi_1} & \cdots & \frac{\partial^2 q}{\partial \lambda_{11} \partial \psi_p} \\ \vdots & \cdots & \vdots & \cdots & \vdots & \cdots & \vdots & \vdots & \cdots & \vdots \\ \frac{\partial^2 q}{\partial \lambda_{1m} \partial \lambda_{11}} & \cdots & \frac{\partial^2 q}{\partial \lambda_{1m} \partial \lambda_{1m}} & \cdots& \frac{\partial^2 q}{\partial \lambda_{1m} \partial \lambda_{p1}} & \cdots & \frac{\partial^2 q}{\partial \lambda_{1m} \partial \lambda_{pm}} & \frac{\partial^2 q}{\partial \lambda_{1m} \partial \psi_1} & \cdots & \frac{\partial^2 q}{\partial \lambda_{1m} \partial \psi_p} \\ \vdots & \cdots & \vdots & \cdots& \vdots & \cdots & \vdots & \vdots & \cdots & \vdots \\ \frac{\partial^2 q}{\partial \lambda_{p1} \partial \lambda_{11}} & \cdots & \frac{\partial^2 q}{\partial \lambda_{p1} \partial \lambda_{1m}} & \cdots& \frac{\partial^2 q}{\partial \lambda_{p1} \partial \lambda_{p1}} & \cdots & \frac{\partial^2 q}{\partial \lambda_{p1} \partial \lambda_{pm}} & \frac{\partial^2 q}{\partial \lambda_{p1} \partial \psi_1} & \cdots & \frac{\partial^2 q}{\partial \lambda_{p1} \partial \psi_p} \\ \vdots & \cdots & \vdots & \cdots& \vdots & \cdots & \vdots & \vdots & \cdots & \vdots \\ \frac{\partial^2 q}{\partial \lambda_{pm} \partial \lambda_{11}} & \cdots & \frac{\partial^2 q}{\partial \lambda_{pm} \partial \lambda_{pm}} & \cdots& \frac{\partial^2 q}{\partial \lambda_{pm} \partial \lambda_{p1}} & \cdots & \frac{\partial^2 q}{\partial \lambda_{pm} \partial \lambda_{pm}} & \frac{\partial^2 q}{\partial \lambda_{pm} \partial \psi_1} & \cdots & \frac{\partial^2 q}{\partial \lambda_{pm} \partial \psi_p} \\ \frac{\partial^2 q}{\partial \psi_{1} \partial \lambda_{11}} & \cdots & \frac{\partial^2 q}{\partial \psi_{1} \partial \lambda_{1m}} & \cdots& \frac{\partial^2 q}{\partial \psi_{1} \partial \lambda_{p1}} & \cdots & \frac{\partial^2 q}{\partial \psi_{1} \partial \lambda_{pm}} & \frac{\partial^2 q}{\partial \psi_{1} \partial \psi_{1}} & \cdots & \frac{\partial^2 q}{\partial \psi_{1} \partial \psi_p} \\ \vdots & \cdots & \vdots & \cdots& \vdots & \cdots & \vdots & \vdots & \cdots & \vdots \\ \frac{\partial^2 q}{\partial \psi_{p} \partial \lambda_{11}} & \cdots & \frac{\partial^2 q}{\partial \psi_{p} \partial \lambda_{1m}} & \cdots& \frac{\partial^2 q}{\partial \psi_{p} \partial \lambda_{p1}} & \cdots & \frac{\partial^2 q}{\partial \psi_{p} \partial \lambda_{pm}} & \frac{\partial^2 q}{\partial \psi_{p} \partial \psi_{1}} & \cdots & \frac{\partial^2 q}{\partial \psi_{p} \partial \psi_p} \end{pmatrix}.\end{align*} $$

We are able to test whether $H_q$ at the solution ( $\hat {\Lambda }_{t^*}, \hat {\Psi }_{t^*}$ ) is positive definite or not by using its observed Fisher information matrix that is defined as $H_q/2$ . So Step 2 of Algorithm 2 uses its observed Fisher information matrix for testing whether $H_q$ at the solution ( $\hat {\Lambda }_{t^*}, \hat {\Psi }_{t^*}$ ) is positive definite or not. If it is positive definite, the estimate is categorized as either proper or improper according to the sign of the unique variances, as in Step 4 or 6 of Algorithm 2. When all unique variances are positive, the estimate is referred to as a “proper solution,” as in Step 4 of Algorithm 2. Meanwhile, when at least one of the unique variances is not positive, the estimate is referred to as an “improper solution,” as in Step 6 of Algorithm 2. When the Hessian matrix at the solution is not positive definite, there are no optimal MLEs. In this case, the estimate is referred to as a “no solutions,” as in Step 9 of Algorithm 2. With Algorithm 2, we can investigate the tendency of the solution pattern through Monte Carlo simulations, which will be described in the following section.

4.1 Investigation of MLEs

We investigate the MLEs obtained by computational algebra and factanal. Table 1 shows the solution pattern of the MLE computed using Algorithm 2. For S3, we obtain proper solutions in all ten simulations, implying that the results for S3 are stabler than S1 and S2. Meanwhile, several improper solutions are obtained, or no solutions are found for S1 and S2. In particular, S1 tends to produce “no solutions” more frequently than S2.

Table 1 Solution pattern obtained by Algorithm 2: proper solution (P), improper solution (I), and no solution (NA)

With factanal function, however, it would be difficult to identify the solution pattern as in Table 1 because the factanal provides estimates under a restriction that all unique variances are larger than some threshold, such as 0.005. Moreover, the observed Fisher information matrices of the MLEs computed with factanal function are positive definite for all simulated datasets; therefore, it is remarkably challenging to distinguish “improper solution” and “no solutions” using factanal function.

To study the solution pattern in more detail, we investigate the number of solutions obtained by computational algebra, as shown in Table 2. When counting the number of solutions, the identifiability with respect to the sign of column vectors in a loading matrix is addressed. We categorize the algebraic solutions based on the signs of the unique variances:

• • For (a)–(b), all unique variances are positive.

• • For (c)–(d), all unique variances are non-zero but some are negative.

• • For (e)–(h), some unique variances are zero.

When some unique variances are zero, the solution does not always satisfy the original normal equation in (4). Therefore, we check whether the solution satisfies (4) in such cases numerically. Meanwhile, when all unique variances are non-zero, the solutions to (6) and (9) are equivalent to those of (4). For this reason, we do not verify (4) for cases where all unique variances are non-zero.

Table 2 The number of solutions obtained by computational algebra

Note: When counting the number of solutions, the identifiability with respect to the sign of column vectors in a loading matrix is addressed.

Moreover, interestingly, in our simulation, a very complicated solution space is computed when all unique variances are zero. In such cases, it is not possible to evaluate the value of the log-likelihood function. Therefore, we do not analyze these cases in this simulation.

For all cases, the number of solutions that satisfy the conditions (6) and (9) is relatively large, but most of them do not have a positive-definite observed Fisher information matrix. Thus, it would be essential to check if the observed Fisher information matrix is positive definite or not. Another noteworthy point is that although we get a number of solutions with zero unique variances using (6) and (9), about 98.5 % of them are not the solution of (4) (see (e) and (f)), implying that most of the solutions to (6) and (9) with zero unique variances are inappropriate. Furthermore, in some simulations, such solution spaces consist of a finite number of points and subspaces containing an infinite number of points. We note that, for subspaces consisting of an infinite number of points, (e) and (f) in Table 2 show the number of sample points from the connected components computed by our algorithm.

Here, we denote that the jth dataset on simulation model Si is Si-j ( $i=1,2,3; j=1,\dots 10$ ). Further empirical comparison for each case in Table 2 is presented below.

4.2 Investigation of the best estimate obtained by computational algebra and factanal

We further compare the characteristics of the loading matrix obtained by computational algebra and factanal. Figure 3 shows the heatmaps of the estimated factor loadings along with the true ones. The loading matrices of computational algebra are unavailable when the optimal solution does not exist (i.e., S1-2, S1-3, S1-4, S1-5, and S2-4).

Figure 3 Heatmap of the factor loadings obtained by computational algebra and factanal function.

We first investigate whether the estimates can approximate the true loading matrix. For S1, the results for the heatmap of the loading matrix in Figure 3 suggest that both computational algebra and factanal function can well approximate the true loading matrix. In particular, even if no solution is found with the computational algebra, the factanal function can approximate the true loadings. This result is surprising to us because the true loading matrix has an identifiability issue that has been considered as a severe problem (Kano, Reference Kano, Rizzi, Vichi and Bock1998). In contrast to our expectation, the estimates of S2 are much more unstable than those of S1. Thus, the small communalities rather than the identification problem would make the estimation difficult in our simulation settings. For S3, we obtain proper and stable estimates for all settings.

We compare the computational algebra’s estimates and factanal’s ones. In many cases, these estimates are essentially identical. In some cases where the computational algebra provides negative unique variances (i.e., S1-7, S2-3, S2-6, and S2-7), we find more or less difference between the two methods. This is because the computational algebra provides negative unique variances while the factanal provides 0.005 due to the restriction of the parameter space. However, the overall interpretation of these two estimates is essentially identical for S1-7, S2-6, and S2-7. An exception is S2-3; in this case, the results are entirely different from each other. Such a difference may be caused by the fact that the factanal function uses only one initial value. We prepare 100 initial values from $U(0,1)$ , perform the factanal function, and select a solution that maximizes the likelihood function; as a result, the factanal’s estimate is essentially identical to computational algebra’s one. The result suggests that performing factanal function with multiple initial values would be essential when we obtain the improper solutions.

4.3 A numerical algorithm to identify the solution pattern

Because the computational algebra remains a challenge in heavy computation, it would be difficult to apply it to a large model. Still, our empirical result is helpful for constructing a method to determine whether the estimate exists or not with a numerical approach. Specifically, we obtain the estimate by Jennrich and Robinson’s (Reference Jennrich and Robinson1969) algorithm with multiple initial values. Subsequently, we employ the same procedure as Algorithm 2 based on the numerical approach. The algorithm for identifying the solution pattern with the numerical algorithm is detailed in Algorithm 3.

With Algorithm 3, it is possible to further investigate the tendency of solution pattern with a number of simulation iterations. Table 3 shows the distribution of solution pattern with Algorithm 3 over 100 simulation runs. The results indicate that only 30% of S1 provides proper solutions, while S2 has more than double the percentage of proper solutions compared to S1. Thus, S1 tends to provide more improper solutions than S2. Interestingly, we get several “improper solutions” and “no solutions” for both S1 and S2, suggesting that various solution patterns are obtained, whether they arise from identifiability issue (S1) or sampling fluctuation (S2). In other words, it would be difficult to identify the cause of improper solutions by using the solution pattern. However, the distribution of the solution pattern appears to differ among S1–S3. Therefore, it would be essential to find the distribution of solution pattern with a number of simulations (e.g., bootstrap method) for identifying the cause of improper solutions.

Table 3 Solution pattern obtained by Algorithm 3 over 100 simulation runs

5.1 Algebraic approach

In this section, we analyze the actual data studied in Davis (Reference Davis1944). The data have the following sample correlation matrix:

$$ \begin{align*} S = \begin{pmatrix} 1 & 0.72 & 0.41 & 0.28 & 0.52 & 0.71 & 0.68 & 0.51 & 0.68 \\ & 1 & 0.34 & 0.36 & 0.53 & 0.71 & 0.68 & 0.52 & 0.68 \\ & & 1 & 0.16 & 0.34 & 0.43 & 0.42 & 0.28 & 0.41 \\ & & & 1 & 0.30 & 0.36 & 0.35 & 0.29 & 0.36 \\ & & & & 1 & 0.64 & 0.55 & 0.45 & 0.55 \\ & & & & & 1 & 0.76 & 0.57 & 0.76 \\ & & & & & & 1 & 0.59 & 0.68 \\ & & & & & & & 1 & 0.58 \\ & & & & & & & & 1 \\ \end{pmatrix}. \end{align*} $$

As reported in Jöreskog (Reference Jöreskog1967), by applying the maximum likelihood factor analysis model with two common factors to this data, we obtain an improper solution. It is not easy for Algorithm 1 to handle data with nine observed variables and two common factors from a realistic computation time/memory point of view. Therefore, we select five observed variables out of nine, and perform analysis based on these selected variables. In other words, we select five rows and five columns from the sample correlation matrix above. Note that such a variable selection procedure would not typically be performed by an analyst in practice. We reduce the number of variables only to enable the algebraic approach within a realistic computation time and to investigate whether the resulting solution patterns are consistent with those observed in the simulation study in Section 4.

By applying a variable selection algorithm with group lasso (Hirose & Konishi, Reference Hirose and Konishi2012) to the sample correlation matrix above, rows/columns 1, 2, 6, 7, 9 are selected. For this input, that is, for the sample correlation sub-matrix, say $S_{\{ 1,2,6,7,9\}\{ 1,2,6,7,9\}}$ , the factanal function in R outputs a proper solution. Changing the third row and column from 6 to 3, we obtain the sample correlation sub-matrix $S_{\{ 1,2,3,7,9\}\{ 1,2,3,7,9\}}$ . For $S_{\{ 1,2,3,7,9\}\{ 1,2,3,7,9\}}$ , the factanal function produces an improper solution. Moreover, we select rows/columns 1, 2, 3, 4, 5. For $S_{\{ 1,2,3,4,5\}\{ 1,2,3,4,5\}}$ , the factanal function provides an improper solution as well. In this section, we analyze the sample correlation submatrixes $S_{\{ 1,2,6,7,9\}\{ 1,2,6,7,9\}}$ , $S_{\{ 1,2,3,7,9\}\{ 1,2,3,7,9\}}$ and $S_{\{ 1,2,3,4,5\}\{ 1,2,3,4,5\}}$ by using Algorithm 2.

The MLEs are computed based on Algorithm 1 from the sub-matrices expressed as real numbers to the second decimal place. The computation time for the case of the second decimal place was about several months. Table 4 shows the solution patterns obtained by Algorithm 2 for $S_{\{ 1,2,6,7,9\}\{ 1,2,6,7,9\}}$ , $S_{\{ 1,2,3,7,9\}\{ 1,2,3,7,9\}}$ , and $S_{\{ 1,2,3,4,5\}\{ 1,2,3,4,5\}}$ .

Table 4 Solution pattern obtained by Algorithm 2: proper solution (P), improper solution (I), and no solution (NA)

For the selected subsets of rows and columns, a proper solution (P) is obtained for $\{1,2,6,7,9\}$ , an improper solution (I) for $\{1,2,3,7,9\}$ , and no solution (NA) for $\{1,2,3,4,5\}$ . As in the investigation in Section 4.1, proper solutions, improper solutions, and no solutions are obtained from algebraic computations for the actual data.

We also investigate the performance of the numerical algorithms. Specifically, we apply factanal and the method of Jennrich & Robinson (Reference Jennrich and Robinson1969) to $S_{\{ 1,2,6,7,9\}\{ 1,2,6,7,9\}}$ , $S_{\{ 1,2,3,7,9\}\{ 1,2,3,7,9\}}$ , and $S_{\{ 1,2,3,4,5\}\{ 1,2,3,4,5\}}$ . Several important findings are summarized as follows:

• • For $S_{\{ 1,2,6,7,9\}\{ 1,2,6,7,9\}}$ , factanal and computational algebra provide almost identical results with proper solutions.

• • For $S_{\{ 1,2,3,7,9\}\{ 1,2,3,7,9\}}$ , factanal and computational algebra yield similar estimates for most elements of the loading matrix. An exception is the estimate of $\lambda _{22}$ due to the negative unique variance ( $\hat {\psi }_2 = -0.212$ ) with computational algebra, compared to $\hat {\psi }_2=0.005$ with factanal.

• • For $S_{\{ 1,2,3,4,5\}\{ 1,2,3,4,5\}}$ , the algorithm based on Jennrich & Robinson (Reference Jennrich and Robinson1969) produces a very large negative value of unique variance ( $\hat {\psi }_1 = -11642.5$ ).

These results indicate that the numerical findings are largely consistent with those obtained with Monte Carlo simulation.

We perform Algorithm 3 to numerically investigate the solution pattern for such findings. The results are identical to those obtained by algebraic computation shown in Table 4. As Algorithm 3 completes the analysis in just seconds with a standard laptop, the algorithm would be useful when sufficient computational resources are unavailable. Furthermore, Algorithm 3 can be applied to larger models. We apply Algorithm 3 to the original covariance matrix, S; in this case, no optimal solution is obtained.

Table 5 shows the number of solutions to (6) and (9), which closely resembles the results in Table 2 from the Monte Carlo simulation in Section 4.1. When a proper solution is obtained, we get a unique stationary point of the log-likelihood function, with a positive definite observed Fisher information matrix and positive unique variances. When an improper solution is obtained, a stationary point of the log-likelihood function is obtained, with a positive definite observed Fisher information matrix and some negative unique variances. In addition, it can be seen that “NA” does not have these stationary points. Interestingly, most of the solutions to (6) and (9) with zero unique variances do not satisfy (4), as is the case with Monte Carlo simulations in Section 4. Thus, it is confirmed that the results of the actual data analysis in Table 5 show patterns similar to those of the Monte Carlo simulation presented in Table 2.

Table 5 The number of solutions obtained by computational algebra

Note: When counting the number of solutions, the identifiability with respect to the sign of column vectors in a loading matrix is addressed.

The overall analysis of the actual data provides results similar to those from the Monte Carlo simulations. Further investigation would involve resampling methods, such as bootstrap, but this is beyond the scope of our research. We would like to take this as a future research topic.

5.2 Numerical computation

The previous section uses a low-dimensional data to apply algebraic computation. In this section, we apply Algorithm 3, a numerical algorithm to investigate the behavior of improper solutions, to five commonly-used benchmark datasets. It is known that each dataset exhibits improper solutions under certain conditions: the Davis data in the previous example ( $n=421$ and $p=9$ , improper at $m=2$ ) (Davis, Reference Davis1944); Harman’s 24 mental tests ( $n=145$ and $p=24$ , improper at $m=6$ ) (Harman, Reference Harman1976); Maxwell’s normal children dataset ( $n=810$ and $p=10$ , improper at $m=4$ ) (Maxwell, Reference Maxwell1961); the Bechtoldt dataset based on the first half-sample of Thurstone and Thurstone’s (Reference Thurstone and Thurstone1941) 17-variable test battery ( $n=212$ and $p=17$ , improper at $m=6$ ) (Bechtoldt, Reference Bechtoldt1961); and Kendall’s job applicant ratings dataset ( $n=48$ and $p=15$ , improper at $m=4$ ) (Kendall, Reference Kendall1975). For each dataset, we conduct the following three analyses:

1. (i) We apply Algorithm 3 to the original data or its correlation matrix to investigate the solution patterns of improper solutions.

2. (ii) We investigate the distribution of solution patterns by simply applying the parametric bootstrap with $B=100$ resamples. The MLEs with varimax rotation of the loading matrix and unique variances are treated as the true values, and datasets with the same sample sizes are generated from a multivariate normal distribution using these parameters.

3. (iii) Same as (ii), but with a different loading matrix. We select the column vector with the smallest sum of squares and modify it so that the two elements with the largest absolute values are retained and all others are set to zero. As a result, the column vector with the smallest sum of squares has only two nonzero entries. Therefore, the loading matrix is not identifiable because the condition of theorem 5.6 in Anderson & Rubin (Reference Anderson and Rubin1956) is not satisfied. With this modified loading matrix, the unique variances are computed such that the diagonal elements of the covariance matrix of the observed variables are equal to one; that is, $\Psi = \mathrm {diag}(I - LL^\top )$ . The bootstrap is then performed using this loading matrix and the corresponding unique variances.

Table 6 shows the results for the above-mentioned three analyses. When we apply Algorithm 3 to the original data, the first four datasets yield “no solutions,” while Kendall’s dataset exhibits an improper solution with a negative unique variance.

Table 6 Distribution of solution patterns for five datasets

Note: (i) Original data; (ii) bootstrap (estimated parameters); (iii) bootstrap (modified loading matrices). P: proper solution, I: improper solution, and NA: no solution.

Next, in (ii), we employ the parametric bootstrap. Based on the results from the original datasets, we expect a large number of cases with no solutions for the first four datasets, while improper solutions are expected for Kendall’s dataset. However, we actually observe that the first two datasets (Davis and Harman) yield mostly proper solutions, and the next two (Maxwell and Thurstone) produce mostly improper solutions. For Kendall’s dataset, we observe many improper solutions, but the proportion is still only 51%. Therefore, the distribution is quite dispersed rather than concentrated, and does not always coincide with the solution pattern of the original dataset.

In (iii), the bootstrap is conducted using loading matrices that are not identifiable. Therefore, we expect most of the results for the solution patterns to be no solutions. However, even when the true loading matrix is not identifiable, we often observe negative unique variances (especially with the Davis data), and sometimes observe positive unique variances as well.

These bootstrap results are consistent with Table 1 from the simulation study. The results suggest that it would be difficult to infer the true loading matrix based solely on the observed solution patterns. One possible reason is that the loss (discrepancy) function is highly complex. Because of this complexity, even when the true loading matrix suffers from identifiability issues, it is still possible to obtain either a proper or improper solution. As a result, the boundaries between proper solutions, improper solutions, and no solutions are ambiguous, and the outcome can easily change depending on sampling variability. While this kind of consideration may have been discussed at an intuitive level, our algebraic approach (Algorithm 2) and its numerical mimic (Algorithm 3) provide further formal support.