The Multiple Roots Phenomenon in Maximum Likelihood Estimation for Factor Analysis

TL;DR

This paper investigates the multiple roots problem in maximum likelihood estimation for factor analysis, revealing the existence of uncountably many solutions and characterizing their structure in simple cases.

Contribution

It provides a rigorous analysis of the multiple roots phenomenon in factor analysis MLE, including proofs of the solution set structure in basic scenarios.

Findings

Likelihood equations have uncountably many solutions.

Solutions form a one-dimensional real curve in simple cases.

Multiple roots cause computational and inferential challenges.

Abstract

Multiple root estimation problems in statistical inference arise in many contexts in the literature. In the context of maximum likelihood estimation, the existence of multiple roots causes uncertainty in the computation of maximum likelihood estimators using hill-climbing algorithms, and consequent difficulties in the resulting statistical inference. In this paper, we study the multiple roots phenomenon in maximum likelihood estimation for factor analysis. We prove that the corresponding likelihood equations have uncountably many feasible solutions even in the simplest cases. For the case in which the observed data are two-dimensional and the unobserved factor scores are one-dimensional, we prove that the solutions to the likelihood equations form a one-dimensional real curve.

Equations

$$\ell ({\tau }^2,\beta ,R)=-\frac{1}{2}nf({\tau }^2,\beta ,R)$$

$$f({\tau }^2,\beta ,R)=\log \det ({\tau }^2+{\beta }^{\prime }R\beta )+\mathrm{tr}{C}_{yy}({\tau }^2+{\beta }^{\prime }R\beta {)}^{-1}.$$

$$\frac{\partial}{\partial t}\log\det(A(t))=\operatorname{tr}\Big{[}A(t)^{-1}\frac{\partial}{\partial t}A(t)\Big{]}$$

$$\frac{\partial}{\partial t}A(t)^{-1}=-A(t)^{-1}\Big{[}\frac{\partial}{\partial t}A(t)\Big{]}A(t)^{-1}.$$

$$\frac{\partial }{\partial {\tau }_k}({\tau }^2+{\beta }^{\prime }R\beta )=2{\tau }_k{E}_{kk},$$

$$\frac{\partial }{\partial {\tau }_k}\log \det ({\tau }^2+{\beta }^{\prime }R\beta )$$

$$=2{\tau }_k\det ({\tau }^2+{\beta }^{\prime }R\beta {)}^{-1}{\mathrm{Cof}}_{kk}({\tau }^2+{\beta }^{\prime }R\beta ),$$

$$\frac{\partial }{\partial {\tau }_k}\mathrm{tr}{C}_{yy}({\tau }^2+{\beta }^{\prime }R\beta {)}^{-1}$$

$$=-2{\tau }_k\mathrm{tr}{C}_{yy}({\tau }^2+{\beta }^{\prime }R\beta {)}^{-1}{E}_{kk}({\tau }^2+{\beta }^{\prime }R\beta {)}^{-1}.$$

$$(\tau^{2}+\beta^{\prime}R\beta)^{-1}=\det(\tau^{2}+\beta^{\prime}R\beta)^{-1}\Big{(}\operatorname{Cof}_{ij}(\tau^{2}+\beta^{\prime}R\beta)\Big{)}$$

$$({\tau }^2+{\beta }^{\prime }R$$

$$\displaystyle=\det(\tau^{2}+\beta^{\prime}R\beta)^{-2}\Big{(}\operatorname{Cof}_{ij}(\tau^{2}+\beta^{\prime}R\beta)\Big{)}E_{kk}\Big{(}\operatorname{Cof}_{ij}(\tau^{2}+\beta^{\prime}R\beta)\Big{)}$$

$$\displaystyle=\det(\tau^{2}+\beta^{\prime}R\beta)^{-2}\Big{(}\operatorname{Cof}_{ik}(\tau^{2}+\beta^{\prime}R\beta)\cdot\operatorname{Cof}_{kj}(\tau^{2}+\beta^{\prime}R\beta)\Big{)}.$$

$$\begin{gathered} \mathrm{tr}{C}_{yy}({\tau }^2+{\beta }^{\prime }R\beta {)}^{-1}{E}_{kk}({\tau }^2+{\beta }^{\prime }R\beta {)}^{-1} \\ =\det ({\tau }^2+{\beta }^{\prime }R\beta {)}^{-2}\sum _{i,j=1}^p{c}_{ij}{\mathrm{Cof}}_{jk}({\tau }^2+{\beta }^{\prime }R\beta )\cdot {\mathrm{Cof}}_{ki}({\tau }^2+{\beta }^{\prime }R\beta ). \end{gathered}$$

$$\frac{\partial }{\partial {\tau }_k}f({\tau }^2,\beta ,R)$$

$$-2{\tau }_k\det ({\tau }^2+{\beta }^{\prime }R\beta {)}^{-2}\sum _{i,j=1}^p{c}_{ij}{\mathrm{Cof}}_{jk}({\tau }^2+{\beta }^{\prime }R\beta )\cdot {\mathrm{Cof}}_{ki}({\tau }^2+{\beta }^{\prime }R\beta ).$$

$$\frac{\partial }{\partial {\beta }_{kl}}({\beta }^{\prime }R\beta )$$

$$=E_{kl}^{\prime }R\beta +{\beta }^{\prime }R{E}_{kl}.$$

$$E_{kl}^{\prime }R\beta =({\delta }_{ik}{\delta }_{jl}{)}^{\prime }R\beta$$

$$\displaystyle=\Big{(}\sum_{m=1}^{q}\delta_{il}\delta_{mk}(R\beta)_{mj}\Big{)}=\big{(}\delta_{il}(R\beta)_{kj}\big{)},$$

$$\beta^{\prime}RE_{kl}=(E_{kl}^{\prime}R\beta)^{\prime}=\big{(}\delta_{jl}(R\beta)_{ki}\big{)},$$

$$\frac{\partial}{\partial\beta_{kl}}(\beta^{\prime}R\beta)=\big{(}\delta_{il}(R\beta)_{kj}+\delta_{jl}(R\beta)_{ki}\big{)}.$$

$$\frac{\partial }{\partial {\beta }_{kl}}\log$$

$$\displaystyle=\operatorname{tr}\Big{[}(\tau^{2}+\beta^{\prime}R\beta)^{-1}\frac{\partial}{\partial\beta_{kl}}(\tau^{2}+\beta^{\prime}R\beta)\Big{]}$$

$$\displaystyle=\operatorname{tr}\Big{[}(\tau^{2}+\beta^{\prime}R\beta)^{-1}\big{(}\delta_{il}(R\beta)_{kj}+\delta_{jl}(R\beta)_{ki}\big{)}\Big{]}$$

$$=\det ({\tau }^2+{\beta }^{\prime }R\beta {)}^{-1}\cdot \sum _{i,j=1}^p{\mathrm{Cof}}_{ij}({\tau }^2+{\beta }^{\prime }R\beta )\cdot ({\delta }_{jl}(R\beta {)}_{ki}+{\delta }_{il}(R\beta {)}_{kj}).$$

$$\frac{\partial }{\partial {\beta }_{kl}}({\tau }^2+{\beta }^{\prime }R\beta {)}^{-1}$$

$$\displaystyle=-(\tau^{2}+\beta^{\prime}R\beta)^{-1}\big{(}\delta_{il}(R\beta)_{kj}+\delta_{jl}(R\beta)_{ki}\big{)}(\tau^{2}+\beta^{\prime}R\beta)^{-1},$$

$$\frac{\partial}{\partial\beta_{kl}}\operatorname{tr}C_{yy}(\tau^{2}+\beta^{\prime}R\beta)^{-1}=-\operatorname{tr}C_{yy}(\tau^{2}+\beta^{\prime}R\beta)^{-1}\big{(}\delta_{il}(R\beta)_{kj}+\delta_{jl}(R\beta)_{ki}\big{)}(\tau^{2}+\beta^{\prime}R\beta)^{-1}.$$

$$\frac{\partial }{\partial {\beta }_{kl}}f({\tau }^2,\beta ,R)$$

$$\displaystyle\qquad-\operatorname{tr}C_{yy}(\tau^{2}+\beta^{\prime}R\beta)^{-1}\big{(}\delta_{il}(R\beta)_{kj}+\delta_{jl}(R\beta)_{ki}\big{)}(\tau^{2}+\beta^{\prime}R\beta)^{-1}.$$

$$\beta^{\prime}R\beta=\big{(}\sum_{m,n=1}^{q}\beta_{mi}r_{mn}\beta_{nj}\big{)}=\big{(}\sum_{m=1}^{q}\beta_{mi}\beta_{mj}+2\sum_{1\leq m<n\leq q}r_{mn}\beta_{mi}\beta_{nj}\big{)},$$

$$\frac{\partial}{\partial r_{kl}}\beta^{\prime}R\beta=2r_{kl}\big{(}\beta_{ki}\beta_{lj}\big{)}$$

$$\frac{\partial }{\partial r_{kl}}\log \det ({\tau }^2+{\beta }^{\prime }R\beta )$$

$$\displaystyle=2r_{kl}\operatorname{tr}\big{[}(\tau^{2}+\beta^{\prime}R\beta)^{-1}\big{(}\beta_{ki}\beta_{lj}\big{)}\big{]}.$$

$$\frac{\partial }{\partial r_{kl}}({\tau }^2+{\beta }^{\prime }R\beta {)}^{-1}$$

$$\displaystyle=-2r_{kl}(\tau^{2}+\beta^{\prime}R\beta)^{-1}\big{(}\beta_{ki}\beta_{lj}\big{)}(\tau^{2}+\beta^{\prime}R\beta)^{-1},$$

Full text

In this paper, we study the multiple roots phenomenon in maximum likelihood estimation for factor analysis. We prove that the corresponding likelihood equations have uncountably many feasible solutions even in the simplest cases. For the case in which the observed data are two-dimensional and the unobserved factor scores are one-dimensional, we prove that the solutions to the likelihood equations form a one-dimensional real curve.

1 Introduction

Systems of equations which have multiple roots, such as critical equations in estimation problems, arise in many contexts in statistical inference and there is now a large literature on these problems. We refer to Small, et al. [15] for an extensive account of the history and methods for investigating numerous multiple roots problems; Buot, et al. [5, 6] for enumeration of the roots of likelihood equations arising from mixture models or from the Behrens-Fisher problem; Catanese, et al. [7] for algebraic aspects of some of these problems and Gross, et al. [10] for applications to variance component models; Hoferkamp and Peddada [12] for order-restricted inference; and Zwiernik, et al. [17] for aspects of linear Gaussian covariance models.

In the context of estimation by the method of maximum likelihood, the existence of multiple roots for the likelihood equations causes uncertainty in the computation of maximum likelihood estimators using hill-climbing algorithms. Consequently, the resulting statistical inference faces the deeper difficulty that it is necessary to determine which of the multiple roots produces the maximum value of the likelihood function.

It is well-known that maximum likelihood estimators (MLEs) enjoy many appealing asymptotic properties. Subject to various regularity conditions and as the sample size tends to infinity, the MLE of a parameter converges in probability to the true value of the parameter, a property known as asymptotic consistency. MLEs also are asymptotically efficient: The variance of the MLE converges to the minimal value among all consistent estimators. Moreover, the asymptotic distribution of MLEs are known, and this enables standard statistical procedures, such as the construction of confidence regions and tests of hypotheses.

Despite the large-sample advantages of MLEs, there are instances in which large samples are unavailable, e.g., in the statistical analysis of clinical trials data. On the other hand, from a purely inferential perspective, one must bear in mind that, as the sample size increases, the likelihood of rejecting the null hypothesis increases simply as an artifact of said increasing sample size. Consequently, it is important to investigate maximum likelihood estimation problems in which sample sizes are small.

In this paper, we study the multiple roots phenomenon in maximum likelihood estimation for the problem of factor analysis; cf., Anderson [3], Harman [11]. The factor analysis problem is old and venerable, with many applications to the social sciences, and its difficult mathematical nature is well-known; cf., Adachi [1], Anderson and Rubin [2], and Rubin and Thayer [14]. Starting with the formulation of the factor analysis model in [14], we prove that the corresponding likelihood equations have uncountably many feasible solutions even in the simplest cases. For the case in which the observed data are two-dimensional and the unobserved factor scores are one-dimensional, we prove that the solutions to the likelihood equations form a one-dimensional real curve.

2 A model for factor analysis

Rubin and Thayer [14] considered a problem in factor analysis, as follows. Let $Y$ denote a $n\times p$ matrix of observed data. Let $Z$ be a $n\times q$ matrix of unobserved factor scores, where $q<p$, and each row of $Z$ is assumed to be marginally distributed as $N_{q}(0,R)$, a $q$-dimensional multivariate normal distribution with mean vector [math] and $q\times q$ correlation matrix $R=(r_{ij})$, i.e., $r_{jj}=1$ for all $j=1,\ldots,q$. It also is assumed that the rows of the matrix $(Y,Z)$ are mutually independent and identically distributed.

For $i=1,\ldots,n$, denote by $Y_{i}$ and $Z_{i}$ the $i$th rows of $Y$ and $Z$, respectively. We suppose that the conditional distribution of $Y_{i}$, given $Z$ is given by $Y_{i}|Z\sim N_{p}(\alpha+Z_{i}\beta,\tau^{2})$, where $\tau^{2}=\operatorname{diag}(\tau_{1}^{2},\ldots,\tau_{p}^{2})$ is a diagonal matrix and $\beta=(\beta_{ij})$ is a $q\times p$ matrix called the factor-loading matrix. By centering the matrix $Y$ it can be assumed, with no loss of generality that $\alpha=0$. Then the problem of factor analysis is to calculate the maximum likelihood estimators of the parameters $\beta$, $\tau^{2}$, and $R$.

It was noted in [14, p. 75] that the likelihood function for this factor analysis model generally has multiple local maxima. The issue of the number of critical points, that is the number of roots of the likelihood equations, was left open, however, and this is the focus of the present paper.

3 The derivatives of the likelihood function

Following [14], we write the log-likelihood function in the form

[TABLE]

with $C_{yy}=(c_{ij})$ denoting the sample covariance matrix constructed from $Y$, and

[TABLE]

Thus, we wish to minimize the function $f$ with respect to the parameters $(\beta,\tau^{2},R)$.

To calculate the derivatives of $f$ with respect to $\tau$ and $\beta$, we will apply repeatedly the following result provided by Magnus and Neudecker [13, p. 169 (Theorem 1) and p. 171 (Theorem 3)].

Lemma 3.1**.**

Suppose that the square matrix $A(t)$ is a differentiable function of a real parameter $t$. Then

[TABLE]

and

[TABLE]

Denote by $E_{kk}$ the $p\times p$ matrix with entry $1$ in the $(k,k)$th position and zeros elsewhere. It is straightforward to verify that

[TABLE]

and then it follows from formula (3.2) that

[TABLE]

where $\operatorname{Cof}_{ij}(A)$ is the $(i,j)$th cofactor of the matrix $A$.

Next, we apply the formula (3.3) to obtain

[TABLE]

This latter expression seems formidable initially, but it can be rewritten in terms of the entries of $C_{yy}$ and the cofactors of $\tau^{2}+\beta^{\prime}R\beta$, as follows: Since

[TABLE]

then

[TABLE]

Therefore,

[TABLE]

Collecting together the formulas above, we obtain

[TABLE]

Next, we calculate the derivative of $f(\tau^{2},\beta,R)$ with respect to $\beta$. Let $E_{kl}$ be a $q\times p$ matrix with entry $1$ in the $(k,l)$th position and zeros elsewhere; that is, the $(i,j)$th entry of $E_{kl}$ is $\delta_{ik}\delta_{jl}$, where $\delta_{ij}$ denotes Kronecker’s delta. Then

[TABLE]

Since

[TABLE]

then

[TABLE]

and hence

[TABLE]

Therefore, by formula (3.2),

[TABLE]

Next, by formula (3.3),

[TABLE]

Therefore,

[TABLE]

Consequently,

[TABLE]

Finally, we calculate the derivatives with respect to the parameters $r_{kl}$, $k<l$. On writing

[TABLE]

we obtain

[TABLE]

where, for fixed $(k,l)$, $(\beta_{ki}\beta_{lj})$ denotes the matrix with generic $(i,j)$th entry $\beta_{ki}\beta_{lj}$. Consequently,

[TABLE]

Next,

[TABLE]

and therefore

[TABLE]

4 The solutions of the likelihood equations

The likelihood equations derived in the previous section are rational equations whose denominators are $(\det(\tau^{2}+\beta^{\prime}R\beta))^{2}$. In order to work with polynomial equations, we will substitute $\gamma$ for $(\det(\tau^{2}+\beta^{\prime}R\beta))^{-1}$ in the likelihood equations, add to the system of equations the new equation

[TABLE]

and denote by $\mathbf{f}_{p,q,C_{yy}}=0$ the resulting expanded system of equations.

The set $V_{p,q,C_{yy}}\subset\mathbb{R}^{p+pq+q(q-1)/2+1}$ of real solutions to the system $\mathbf{f}_{p,q,C_{yy}}=0$ is equal to the set

[TABLE]

We now set $p=2$ and $q=1$. Let us consider the polynomials in $\mathbf{f}_{2,1,C_{yy}}$ and treat the entries of $C_{yy}=(c_{ij})$ as unknown variables. Furthermore, for simplification, we make the substitution,

[TABLE]

and add to $\mathbf{f}_{2,1,C_{yy}}$ the equations from (4.2). This new set of polynomials defines an ideal $I_{2,1}\subseteq\mathbb{Q}[c_{11},c_{12},c_{22},\tau_{1},\tau_{2},\beta_{11},\beta_{12},\gamma,x_{11},x_{12},x_{22}]$. As we are interested in the solution set associated with these polynomials, we study the variety $V(I_{2,1})\subseteq\mathbb{C}^{11}$ defined to the set of all points $p=(c_{11},c_{12},c_{22},\tau_{1},\tau_{2},\beta_{11},\beta_{12},\gamma,x_{11},x_{12},x_{22})\in\mathbb{C}^{11}$ such that $g(p)=0$ for all polynomials $g\in I_{2,1}$. To compute the dimension of the set $V(I_{2,1})$, we first decompose it into a union of irreducible varieties. The corresponding algebraic object is the intersection of the primary ideals. For an introduction to ideals, varieties, and primary decompositions see [8].

Using Macaulay2 (Grayson and Stillman [9]), we obtain the decomposition $I_{2,1}=J_{1}\cap J_{2}$ of $I_{2,1}$ into an intersection of primary ideals $J_{1}$ and $J_{2}$, where

[TABLE]

Moreover, we obtain from Macaulay2 that the dimensions of the components $V(J_{1})$ and $V(J_{2})$ corresponding to the ideals $J_{1}$ and $J_{2}$ are 4 and 3, respectively. These computations can be confirmed using numerical algebraic geometry software such as PHCpack (Verschelde [16]). Since there are no relations between the variables $c_{11}$, $c_{12}$, and $c_{22}$, in other words, they are free variables, it follows that if we choose $C_{yy}$ generically from the set of $p\times p$ positive definite matrices then the solution set $V_{2,1,C_{yy}}$ contains a curve (i.e., a one-dimensional component) and a set of points (i.e., a zero-dimensional component).

Moreover, the curve in $V_{2,1,C_{yy}}$ can be described parametrically in terms of $t$ for all values of $t$ such that $c_{12}^{2}/c_{22}\leq t^{2}\leq c_{11}$:

[TABLE]

Since $C_{yy}$ is positive definite then the interval $\{t\in\mathbb{R}:\ c_{12}^{2}/c_{22}\leq t^{2}\leq c_{11}\}$ has positive measure.

The isolated points in $V_{2,1,C_{yy}}$ also can be described analytically:

[TABLE]

In summary, we have proved the following result.

Theorem 4.1**.**

*Suppose that $p=2$ and $q=1$. For a generic sample covariance matrix $C_{yy}$, the likelihood equations for the factor analysis model in Section 2 have an infinite number of real solutions. *

For $p\geq 2$ and $q\geq 1$, we can apply the same approach as in the case $(p,q)=(2,1)$ to construct a sample covariance matrix $C_{yy}$ such that $V_{p,q,C_{yy}}$ is positive dimensional.

Theorem 4.2**.**

For all $p\geq 2$ and $q\geq 1$, there exists a covariance matrix $C_{yy}$ such that $V_{p,q,C_{yy}}$ is positive dimensional.

Proof.

Let $C_{yy}$ be a positive definite matrix of the following form:

[TABLE]

The set of all $(\tau,\beta,R)$ such that $r_{kl}=0$ for $1\leq k<l\leq q$ and

[TABLE]

is a one-dimensional subset of $V_{p,q,C_{yy}}$. Indeed, this set is obtained by solving for $\beta_{11},\beta_{12},\tau_{1},\tau_{2}$ in equations (4.3)-(4.4) and then setting $\beta_{kl}=0$ for $(k,l)\neq(1,1),(1,2)$ and $\tau_{k}=\pm\sqrt{c_{kk}}$. ∎

Remark 4.3**.**

We note that the set of solutions in Theorem 4.1 is the set of solutions to the equation

[TABLE]

Since $f(\tau^{2},\beta,R)$ is a function of $\tau^{2}+\beta^{\prime}\beta$ for the case in which $(p,q)=(2,1)$ then the likelihood function is constant on $V_{2,1,C_{yy}}$. Thus, this set of solutions forms a ridge on the likelihood hypersurface. Although it follows immediately from (3.1) that the log-likelihood function attains the same value at every solution to (4.5), one cannot similarly conclude immediately that these solutions are critical points, the reason being that the first partial derivatives of $f(\tau^{2},\beta,R)$ are not functions of $\tau^{2}+\beta^{\prime}\beta$.

While Theorem 4.2 describes some sample covariance matrices that result in an infinite number of solutions to the likelihood equations, we believe that there are more. However, once we increase $p$ and $q$, even modestly, the computations become infeasible with current methodology. One approach would be to use numerical algebraic geometry, but such approaches would require creative manipulation of the system.

Acknowledgments

The research of Gross was supported by the National Science Foundation award DMS-1304167 and DMS-1620109. The research of Richards was partially supported by the National Science Foundation under grant DMS-1309808. The research of Stasi and Petrović was partially supported by Air Force Office for Scientific Research grant FA9550-14-1-0141. We thank Daniel Brake for helping us to visualize the curve in Theorem 4.1 using bertini_real [4].