Factor analysis  is a multivariate statistical technique that characterizes dependence among variables using a small number of latent factors. Suppose that we have a sample from the
-variate Gaussian distribution
with mean vectorand a covariance matrix . We assume that where is a matrix of rank
that describes the amount of variance shared among thecoordinates and is a diagonal matrix with positive diagonal entries representing the unique variance specific to each coordinate. Factor analysis of Gaussian data for was first formalized by  with efficient maximum likelihood (ML) estimation methods proposed by [3, 4, 5, 1] and others. These methods however do not apply to datasets with that occur in applications such as the processing of sequencing data [6, 7, 8], analyzing the transcription factor activity profiles of gene regulatory networks using massive gene expression datasets , portfolio analysis in stock return . In such cases, the available computer memory may be inadequate to store the sample covariance matrix or to make multiple copies of the dataset needed during the computation.
The expectation-maximization (EM) approach of  can be applied to datasets with but is computationally slow. So, here we develop a profile likelihood method for high-dimensional Gaussian data. Our method allows us to compute the gradient of the profile likelihood function at negligible additional computational cost and to check first-order optimality, guaranteeing high accuracy. We develop a fast sophisticated computational framework called FAD (Factor Analysis of Data) to compute ML estimates of and . Our framework is implemented in an R  package called fad.
The remainder of this paper is organized as follows. Section II describes the factor model for Gaussian data and an ML solution using the EM algorithm, and then proposes the profile likelihood and FAD. The performance of FAD relative to EM is evaluated in Section III. Section IV applies our methodology on a functional magnetic resonance imaging (fMRI) dataset related to suicidal behavior . Section V concludes with some discussion. An online supplement, with sections, tables and figures referenced here with the prefix “S”, is available.
Ii-a Background and Preliminaries
Let be the data matrix with as its th row. Then, in the setup of Section I, the ML method profiles out using the sample mean vector and then maximizes the log-likelihood,
where , , where is the vector of 1s. The matrix is almost surely singular and has rank when . The factor model (1) is not identifiable because the matrices and
give rise to the same likelihood for any orthogonal matrixSo, additional constraints [1, 5] are imposed.
Ii-A1 EM Algorithms for parameter estimation
The EM algorithm [11, 14] exploits the structure of the factor covariance matrix by assuming -variate standard normal latent factors and writing the factor model as where ’s are i.i.d and ’s are independent of ’s. The EM algorithm is easily developed, with analytical expressions for both the expectation (E-step) and maximization (M-step) steps that can be speedily computed (see Section S1.1).
Although EM algorithms are guaranteed to increase the likelihood at each iteration and converge to a local maximum, they are well-known for their slow convergence. Further, these iterations run in a -dimensional space that can be slow for very large . Accelerated variants [15, 16] show superior performance in low-dimensional problems but come with additional computational overhead that dominates the gain in rate of convergence in high dimensions. EM algorithms also compromise on numerical accuracy by not checking for first-order optimality to enhance speed. So, we next develop a fast and accurate method for EFA that is applicable in high dimensions.
Ii-B Profile likelihood for parameter estimation
We start with the common and computationally useful identifiability restriction on that constrains to be diagonal with decreasing diagonal entries. This scale-invariant constraint is completely determined except for changes in sign in the columns of . Under this constraint, can be profiled out for a given as described in the following
Suppose that is a given p.d. diagonal matrix. Suppose that the largest singular values of
largest singular values ofare and the corresponding -dimensional right-singular vectors are the columns of Then the function is maximized at where is a diagonal matrix with th diagonal entry as The profile log-likelihood equals,
where is a constant that depends only on and but not on Furthermore, the gradient of is given by:
See Section S1.2. ∎
The profile log-likelihood in (2) depends on only through the largest singular values of So, in order to compute and we need to only compute the largest singular values of and the right singular vectors. For , as is usually the case, these largest singular values and singular vectors can be computed very fast using Lanczos algorithm.
Further, is optimized over a -dimensional space. In the absence of constraints on and analytically tractable optima, this is the largest possible reduction on the dimension of the parameter space for feeding into a numerical optimization method. Further constraints on (e.g. when all diagonals are equal) can be easily incorporated in the method if needed. Overall, this reduction in dimensionality further accelerates convergence. Also, is available in closed form that enables us to check first-order optimality and to ensure high accuracy.
Finally, is expressed in terms of However, ML estimators are scale-invariant, so we can estimate and using the correlation matrix and scale back to . A particular advantage of using the sample correlation matrix is that needs to be optimized over a fixed bounded rectangle that does not depend on the data and is conceivably numerically robust.
Ii-C Matrix–free computations
Ii-C1 A Lanczos algorithm for calculating partial singular values and vectors
In order to compute the profile likelihood and its gradient, we need the largest singular values and right singular vectors of We use the Lanczos algorithm [17, 18] with reorthogonalization and implicit restart. Suppose that and that is any random vector with Let and For let reorthogonalize and set Also within this loop, if update reorthogonalize , and set
Next, consider the bidiagonal matrix with diagonal entries with th entry for and all other entries as 0. Now suppose that are the singular values of and that ’s and ’s are the corresponding right and left singular vectors, which can be computed via a Sturm sequencing algorithm . Also, let and . Then it can be shown that for all and where is the last entry of Because and is approximately the largest singular value of the algorithm stops if holds for where is some prespecified tolerance level, and and are accurate approximations of the largest singular values and corresponding right singular vectors of .
Convergence of the reorthogonalized Lanczos algorithm often suffers from numerical instability that slows down convergence. To resolve this instability,  proposed restarting the Lanczos algorithm, but instead of starting from scratch, they initialized with the first singular vectors. To that end, let and reset Then for , let , and reset and Define and For compute and This yields a matrix with entries and for and for and for , and all other entries 0. The matrix
is not bidiagonal but is still small-dimensioned matrix whose singular value decomposition can be calculated very fast. Convergence of the Lanczos algorithm can be checked as before. This restart step is repeated until all thelargest singular values converge.
The only way that enters this algorithm is through matrix-vector products of the forms and , both of which can be computed without explicitly storing Overall, this algorithm yields the largest singular values and vectors in computational cost using only additional memory, resulting in substantial gains over the traditional methods [3, 4]
. These traditional methods require a full eigenvalue decomposition ofthat is of computational complexity and requires memory storage space.
Having described a scalable algorithm for computing and , we detail a numerical algorithm for computing the ML estimators.
Ii-C2 Numerical optimization of the profile log-likelihood
On the correlation scale, ’s lie between 0 and 1. Under this box constraint, the factanal function in R and factoran function in MATLAB® employ the limited-memory Broyden-Fletcher-Goldfarb-Shanno quasi-Newton algorithm  with box-constraints (L-BFGS-B) to obtain the ML estimator of However, in high dimensions, the advantages of the L-BFGS-B algorithm are particularly prominent. Because Newton methods require the search direction , where is the Hessian matrix of , they are computationally prohibitive in high dimensions in terms of storage and numerical complexity. The quasi-Newton BFGS replaces the computation of the exact search direction by an iterative approximation using the already computed values of and The limited-memory implementation, moreover, uses only the last few (typically less than 10) values of and instead of using all the past values. Overall, L-BFGS-B reduces the storage cost from to and the computational complexity from to . Interested readers are referred to [20, 21] for more details on the L-BFGS-B algorithm.
The L-BFGS-B algorithm requires both and to be computed at each iteration. Because is available as a byproduct while computing (see Section Sections II-B and II-C1), we modify the implementation to jointly compute both quantities with a single call to the Lanczos algorithm at each L-BFGS-B iteration. In comparison to R’s default implementation (factanal) that separately calls and in its optimization routines, this tweak halves the computation time.
Iii Performance evaluations
Iii-a Experimental setup
The performance of FAD was compared to EM using 100 simulated datasets with true or 5 and . For each setting, we simulated the diagonal entries of to be from and entries of and from . We also evaluated performance with to match the settings of our application in Section IV: the true were set to be the ML estimates from that dataset.
For the EM algorithm, was initialized as the first principal components (PCs) of the scaled data matrix computed via the the Lanczos algorithm while was started at . FAD requires only to be initialized, which was done in the same way as the EM. We stopped FAD when the relative increase in was below 100 and where is the machine tolerance, which in our case was approximately . The EM algorithm was terminated if the relative change in was less than and , or if the number of iterations reached 5000. Therefore, FAD and EM had comparable stopping criteria. For each simulated dataset, we fit models with factors and chose the number of factors by the Bayesian Information Criterion (BIC): , where is the maximum log-likelihood value with factors. All experiments were done using R  on a workstation with Intel E5-2640 v3 CPU clocked @2.60 GHz and 64GB RAM.
Because BIC always correctly picked , we evaluated model fit for each method in terms of and where and and are the correlation matrices corresponding to and .
Iii-B1 CPU time
Iii-B2 Parameter estimation and model fit
The final values from the FAD and EM are identical. Figure S5 shows that the two algorithms also yield essentially identical values of and of under the best fitted models.
Iii-C Additional experiments in high-noise scenarios
We conclude this section by evaluating performance in situations where ostensibly, weak factors are hardly distinguished from high noise by SVD methods and where EM may be preferable . We applied FAD and EM to the simulation setup of 
: Here, the uniquenesses were sampled from three inverse Gamma distributions with unit means and variances of, and , and . Figure S9 shows that our algorithm was substantially faster while having similar accuracy as EM.
Iv Suicide ideation study
We applied EFA to data  from an fMRI study conducted while 20 words connoting negative affects were shown to 9 suicide attempters, 8 suicide non-attempter ideators and 17 non-ideator control subjects. For each subject-word combination,  provided voxel-wise per cent changes in activation relative to the baseline in image volumes. Restricting attention to the 24547 in-brain voxels yields datasets for the attempters, ideators and controls of sizes . We assumed each dataset to be a random sample from the multinormal distribution. Our interest was in determining if the variation in the per cent relative change in activation for each subject type can be explained by a few latent factors and whether there are differences in these factors between the three groups of subjects.
For each dataset, we performed EFA with factors and using both FAD and EM. Table I demonstrates the computational benefits of using FAD over EM. We also used BIC to decide on the optimal () and obtained 2-factor models for both suicide attempters and ideators, and a 4-factor model for the control subjects. Figure 12 provides voxel-wise displays of the factor
loadings, obtained using the quartimax criterion , for each type of subject. All the factor loadings are non-negligible only in voxels around the ventral attention network (VAN) which represents one of two sensory orienting systems that reorient attention towards notable stimuli and is closely related to involuntary actions . However, there are differences between the factor loadings in each group and also among them.
For instance, for the suicide attempters, each factor is a contrast between different areas of the VAN, but the contrasts themselves differ between the two factors. The first factor for the ideators is a weighted mean of the voxels while the second factor is a contrast of the values at the VAN voxels. For the controls, the first three factors are different contrasts of the values at different voxels while the fourth factor is more or less a mean of the values at these voxels. Further, the factor loadings in the control group are more attenuated than for either the suicide attempters or ideators. While a detailed analysis of our results is outside the purview of this paper, we note that EFA has provided us with distinct factor loadings that potentially explains the variation in suicide attempters, non-attempter ideators and controls. However, our analysis assumed that the image volumes are independent and Gaussian: further approaches relaxing these assumptions may be appropriate.
In this paper, we propose a new ML-based EFA method called FAD using a sophisticated computational framework that achieves both high accuracy in parameter estimation and fast convergence via matrix-free algorithms. We implement a Lanczos method for computing partial singular values and vectors and a limited-memory quasi-Newton method for ML estimation. This implementation alleviates the computational limitations of current state-of-the-art algorithms and is capable of EFA for datasets with . Although not demonstrated in this paper, FAD is also well-suited for distributed computing systems because it only uses the data matrix for computing matrix-vector products. FAD paves the way to develop fast methods for mixtures of factor analyzers and factor models for non-Gaussian data in high-dimensional clustering and classification problems.
T. W. Anderson, An Introduction to Multivariate Statistical Analysis
, ser. Wiley Series in Probability and Statistics. Wiley, 2003. [Online]. Available:https://books.google.com/books?id=Cmm9QgAACAAJ
-  D. N. Lawley, “The estimation of factor loadings by the method of maximum likelihood,” Proceedings of the Royal Society of Edinburgh, vol. 60, pp. 64–82, 01 1940.
-  K. G. Jöreskog, “Some contributions to maximum likelihood factor analysis,” Psychometrika, vol. 32, pp. 443–482, 12 1967.
-  D. N. Lawley and A. E. Maxwell, “Factor analysis as a statistical method,” Journal of the Royal Statistical Society. Series D (The Statistician), vol. 12, pp. 209–229, 1962.
-  K. V. Mardia, J. T. Kent, and J. M. Bibby, Multivariate analysis. Elsevier Amsterdam, 2006.
-  J. T. Leek and J. D. Storey, “Capturing heterogeneity in gene expression studies by surrogate variable analysis,” PLoS Genetics, vol. 3, no. 9, p. e161, 2007.
-  J. T. Leek, “Svaseq: Removing batch effects and other unwanted noise from sequencing data,” Nucleic Acids Research, vol. 42, no. 21, p. e161, 2014.
-  F. Buettner, N. Pratanwanich, D. J. McCarthy, J. C. Marioni, and O. Stegle, “f-sclvm: Scalable and versatile factor analysis for single-cell rna-seq,” Genome biology, vol. 18, no. 1, p. 212, 2017.
-  I. Pournara and L. Wernisch, “Factor analysis for gene regulatory networks and transcription factor activity profiles,” BMC Bioinformatics, vol. 8, p. 61, 02 2007.
-  C. T. Ng, C. Y. Yau, and N. H. Chan, “Likelihood inferences for high-dimensional factor analysis of time series with applications in finance,” Journal of Computational and Graphical Statistics, vol. 24, pp. 866–884, 11 2014.
-  D. B. Rubin and D. T. Thayer, “Em algorithms for ml factor analysis,” Psychometrika, vol. 47, pp. 69–76, 02 1982.
-  R Core Team, R: A Language and Environment for Statistical Computing, R Foundation for Statistical Computing, Vienna, Austria, 2019. [Online]. Available: https://www.R-project.org/
M. A. Just, L. A. Pan, V. L. Cherkassky, D. L. McMakin, C. B. Cha, M. K. Nock, and D. P. Brent, “Machine learning of neural representations of suicide and emotion concepts identifies suicidal youth,”Nature Human Behaviour, vol. 1, pp. 911–919, 2017.
-  G. McLachlan and T. Krishnan, The EM Algorithm and Extensions, ser. Wiley Series in Probability and Statistics. Wiley, 1996. [Online]. Available: https://books.google.com/books?id=iRSWQgAACAAJ
-  C. Liu and D. B. Rubin, “Maximum likelihood estimation of factor analysis using the ecme algorithm with complete and incomplete data,” Statistica Sinica, vol. 8, pp. 729–747, 08 2002.
-  R. Varadhan and C. Roland, “Simple and globally convergent methods for accelerating the convergence of any em algorithm,” Scandinavian Journal of Statistics, vol. 35, pp. 335–353, 06 2008.
-  J. Baglama and L. Reichel, “Augmented implicitly restarted lanczos bidiagonalization methods,” SIAM Journal on Scientific Computing, vol. 27, no. 1, pp. 19–42, 2005.
-  S. Dutta and D. Mondal, “An h-likelihood method for spatial mixed linear model based on intrinsic autoregressions,” Journal of the Royal Statistical Society: Series B (Statistical Methodology), vol. 77, pp. 699–726, 09 2014.
J. H. Wilkinson, “The calculation of the eigenvectors of codiagonal matrices,”Computer Journal, vol. 1, pp. 90–96, 02 1958.
-  R. H. Byrd, J. N. P. Lu, and C. Zhu, “A limited memory algorithm for bound constrained optimization,” SIAM Journal on Scientific Computing, vol. 16, pp. 1190–1208, 1995.
-  C. Zhu, R. H. Byrd, P. Lu, and J. Nocedal, “Algorithm 778: L-bfgs-b: Fortran subroutines for large-scale bound-constrained optimization,” ACM Transactions on Mathematical Software, vol. 23, pp. 550–560, 1994.
-  G. E. Schwarz, “Estimating the dimension of a model,” The Annals of Statistics, vol. 6, pp. 461–464, 03 1978.
-  A. B. Owen and J. Wang, “Bi-cross-validation for factor analysis,” Statistical Science, vol. 31, pp. 119–139, 03 2015.
-  A. B. Costello and J. Osborne, “Best practices in exploratory factor analysis: Four recommendations for getting the most from your analysis,” Practical Assessment, Research & Evaluation, vol. 10, pp. 1–9, 01 2005.
-  S. Vossel, J. J. Geng, and G. R. Fink, “Dorsal and ventral attention systems: Distinct neural circuits but collaborative roles,” The Neuroscientist, vol. 20, no. 2, pp. 150–159, 2014.
S1 Supplementary materials for Methodology
s1.1 The EM algorithm for factor analysis on Gaussian data
The complete data log-likelihood function is
where is a constant that does not depend on the parameters.
s1.11 E-Step computations
Since the ML estimate of is , at the current estimates and , the expected complete log-likelihood or so called function is given by
Since . Then,
s1.12 M-Step computations
The parameters and are obtained by maximizing following equation S2. Specifically, given , and , the maximizer is given by
where . By Woodbury matrix identity (Henderson and Searle, 1981), , so S5 can be simplified as
Next, given , and , the ML estimate of is given by
Substitute with S5, we get
s1.2 Proof of Lemma 1
From equation 1, the ML estimates of and are obtained by solving the score equations
From , we have
Suppose that and that the diagonal elements in are in decreasing order with . Let with and containing the largest eigenvalues . The corresponding eigenvectors forms columns of matrix so that . Then, if , S10 shows that
The square roots of are the largest singular values of and columns in are then the corresponding right-singular vectors. Hence, conditional on , is maximized at , where is a diagonal matrix with elements .
From the construction of and , we have and hence, .
S2 Additional results for simulation studies
s2.1 Average CPU time
s2.2 Evaluations in terms of and
s2.3 Performance of FAD compared to EM for high-noise models
Henderson, H. V. and Searle, S. R. (1981), “On deriving the inverse of a sum of matrices,” SIAM Review, 23, 53–60.