Sparsest Factor Analysis of Gene Expression Data
Kohei Adachi∗ and Nickolay T. Trendafilov†
Abstract — We consider factor analysis (FA) of an n-observations × p-genes data matrix with p >> n. For suchhigh-dimensional data, FA produces factor loading matrix with too many rows to be easily interpreted. To dealwith this difficulty, we consider an FA procedure subject to the loading matrix being sparsest in that each rowof the matrix has only one nonzero loading. As such a sparsest FA procedure, a least squares
matrix-decomposition approach has been presented by the authors, but it is not feasible for high-dimensionalcase. In this paper, a maximum-likelihood sparsest FA procedure is proposed which is feasible for the data withp >> n. Its key point is using the EM algorithm in which n > p is indispensable. The proposed method is applied
for a gene expression data set.
Keyword: High-dimensional data; The sparsest loadings; EM algorithms

1

Introduction

In factor analysis (FA), the variations of p observed variables are assumed to be explained by the twotypes of unobserved factors called common factors and unique factors, with their uncorrelated mutually.
Here, the number of the common factors, which we denote as m, are supposed to be less than that ofobserved variables, i.e., m < p, and those factors serve for explaining the variations of all p variables. On
the other hand, the number of the unique factors equals that of observed variables (p) and the factors havethe one-to-one relationship to p variables: each of unique factors explains specifically the variation of thecorresponding variable.
Using x (p × 1) for a p-variate vector with its expectation E[x] equaling the zero vector 0p, the FAmodel is expressed as x = Λf + e. Here, f is an m × 1 common factor vector with m < p, e is a p × 1 uniquefactor vector, and Λ is a p × m factor loading matrix describing the relationships between the p variables
and m factors. The p × p covariance matrix Ψ for e is assumed to be a diagonal matrix, whose diagonal
elements are called unique variances, as they stand for the amount of variances of variables accounted forby unique factors. In FA, Λ, Ψ, and factor correlation matrix Φ are estimated from the n-observations ×
p-variables data matrix X whose rows are the realizations of x′.
In this paper, we focus on gene-expression data matrices X with p >> n, i.e., much more variables (i.e.,genes) than observations. For such a case, the resulting p × m loading matrix Λ has too many rows so thatit cannot be interpreted easily. A strategy for avoiding the difficulty is constraining Λ to be sparse, i.e., itto have a number of zero elements. Among the sparse Λ, the sparsest one is the most interpretable. SuchΛ = [λ1, … , λp]′ = (λij) can be formulated as
λi being filled with zeros except a single element with unknown location.

(1)

An FA procedure for obtaining Λ with (1) has been proposed by [3], in which a least squares approachis taken with the matrix-decomposition formulation of FA ([1], [7]). However, this procedure is not
feasible for the data with p > n. On the other hand, a maximum likelihood (ML) FA procedure withEM-algorithm is feasible for such data, as shown in [2]. In this paper, we thus propose a modified MLFA
procedure in which an EM algorithm is used for obtaining sparsest loadings. The modification is detailedin Section 3 after the EM algorithm for the standard FA is introduced in the next section.
2EM Algorithm for FA
The normality of f and e is assumed as f ~ Nm(0m, Φ) and e ~ Np(0p, Ψ), with 0m the m × 1 zero vector.

∗Graduate Schools of Human Sciences, Osaka University, 1-2 Yamadaoka, Suita, Osaka 565-0871, Japan, E-mail:adachi@hus.osaka-u.ac.jp, Tel: +81-6-6879-4040.
† Department of Mathematics and Statistics, Open University, Walton Hall, Milton Keynes MK7 6AA, UK.

Using it in the FA model, we have x ~ Np(0m, Σ) with Σ = ΛΦΛ′+Ψ, which leads to the log-likelihoodf(Λ,Ψ,Φ) = −log| ΛΦΛ′+Ψ| − trS(ΛΦΛ′+Ψ)−1,(2)where S the p × p sample covariance matrix. For maximizing (2) over Λ,Ψ, and Φ, E- and M-steps areiterated in the EM-algorithm.
The E-step requires us to obtaing(Λ,Ψ,Φ) = −Σjlogψi − tr(S+trΛQΛ′−2trCΛ′)Ψ−1 − log|Φ| − trQΦ−1 .

(3)Here, ψi denotes the jth diagonal element of Ψ, Q = A′SA + U, and C = SA, with A and U obtained as A= Σ−1ΛΦ and U = Im − A′ΛΨ−1Φ using Σ, Λ, and Φ in the previous round of iteration.
The M-step requires to obtain Λ,Ψ, and Φ that increase (3). The diagonal Ψ maximizing (3) ispresented in [6]. For maximizing (3) over correlation matrix Φ, we can reparameterize it as Φ = T′T withdiag(T′T) = Im. It allows the task to be reformulate as the maximization as minimizingf(T) = log|T′T| + tr(T′T)−1Q + Constant(4)

over T subject to diag(T′T) = Im with Constant the term irrelevant to T. The constrained minimization of(4) can be attained using the GP algorithm [5], in which T is iteratively updated as Tnew =VTdiag(VTVT′)−1/2. Here, VT = T − α∂f(T)/∂T with α a positive value leading to f(T) ≥ f(Tnew) and∂f(T)/∂T = 2T′TT′−1 − 2T′−1QT−1T′−1 . The remaining task is minimization of (3) over Λ. As detailedin the next section, the modification with (1) is incorporated in the task.
3Sparsest Modification
Let us consider maximizing (3) over Λ subject to (1). We can rewrite (3) asf(T) = const − tr(ΛQΛ′−2CΛ′)Ψ−1 = const − ∑ip=1 g i (λi ) /ψ i ,(5)with ψi > 0, const being an expression irrelevant to Λ, andgi(λi) = λi′Qλi −2Σjcijλij = ΣjΣkqjkλijλik −2Σjcijλij(6)

the function of the ith row of Λ. It shows that the optimal Λ maximizing (5) subject to (1) can be obtainedby minimizing (6) over λi under (1) for each i. Indeed, considering (1) and letting J(i) denote the locationof the element to be given a nonzero value in λi, we can further rewrite gi(λi) as
2ci2, J (i )ci2, J (i )ci , J (i ) gi(λi) = qJ(i),J(i)λi,J(i) − 2ci,J(i)λi,J(i) = q J (i ), J (i )  λi , J (i ) −≥.
−−q J (i ),J (i )q J (i ),J (i )  q J (i ), J (i )
2(7)The inequality in (7), which follows from the positive-definiteness of Q implying qjj > 0, shows that thelower limit of gi(λi) is −ci,J(i)2/qJ(i),J(i) which is attained for λi,J(i) = −ci,J(i)2/qJ(i),J(i). Further, the optimal J(i) isthe index (from 1 to m) for which the limit −ci,J(i)2/qJ(i),J(i) is minimal. This selection of the optimal nonzeroloading can be expressed as
c 2 / q jj
 0λij =  ij

iff j = arg min1≤k ≤m − cik2 / qkk.
otherwise(8)We can find that the formulas require in this and the last sections only require thenonnegative-definiteness of S: it may not be positive-definite. Thus, the proposed method is feasible todata matrices with p > n.
4Application
We performed the proposed method for the yeast cell cycle data matrix of n = 17 by p = 384 (genes)Phase 1

Phase 2Phase 3Phase 4Phase 5(A) Heat maps of the genes × factors blocks in the loading matirx

Factor
12
3
4

0. 0

0.2

0.4

0.6

0.8(B) Histogram of unique variances

1
1

2
-0.13
1

30.04-0.52
140.710.08
0.06
11.0(C) Factor Correlations

Figure 1: Solution of the proposed factor analysis for gene expression data

presented by [8]. This data matrix, which is publicly available at http://faculty.washington.edu/kayee/pca,
have been first log-transformed and then standardized so that the column averages and variances are zeroand one, respectively. The solution resulting with m set at 4 as in [4] is shown in Figure 1.
The solution of the loading matrix is presented block-wise in Figure 1(A). As the 384 genes are
categorized into five phases of cell cycles [8], the five blocks (genes × components) in (A) correspond to
the five phases. The loadings are considered to be reasonable, as each phase has a specific feature ofloadings: [a] The genes in Phases 1, 2, and 4 are positively loaded by Factors 1, 2, and 3, respectively; [b]Phase 5 are characterized by positive loadings for Factor 4 and negative ones for 2; [c] Phase 3 consists ofthe genes positively loaded by Factor 2 or 3 and by both.
The histogram of the resulting unique variances is shown in Figure 1(B). Their range was [0.03,0.89]: a solution with ψj = 0 was not found, which suggests that the unique variances were reasonablyestimated.
Figure 1(C) shows the resulting factor correlations, where we can find factor 2 is negativelycorrelated with factor 3, while the factor 1 are considerably positively correlated with factor 4.
Those results are summarized as the path diagram in Figure 2.
5

Final Remarks
We proposed a sparsest FA procedure with an EM algorithm, which is feasible for high-dimensionalG5 • • • G5

G1 • • • G1•• ••• •f4

f1

G2 • • • G2

G3 • • • G3 G4 • • • G4•• ••• •f2•• •f3Figure 2: Path diagram for the solution in Figure 1 with f# and G standing forfactors and genes, respectively.

data and provides a factor loading matrix with each row having single nonzero element. The applicationof the procedure to a gene expression data set demonstrated that genes can reasonably be classified into a
few groups characterized by common factors.
Acknowledgment
This work is supported by a grant RPG-2013-211 from The Leverhulme Trust, UK and KAKENHI(26330039; Grant-in-Aid for Scientific Research (C)) from JSPS, Japan.
References[1] Adachi, K.: Some contributions to data-fitting factor analysis with empirical comparisons to covariance-fitting
factor analysis. Journal of the Japanese Society of Computational Statistics, 25(1): 25-38 (2012).
[2] Adachi, K. (2013). Factor analysis with EM algorithm never gives improper solutions when sample covarianceand initial parameter matrices are proper. Psychometrika, 78(2): 380-394.
[3] Adachi, K. and Trendafilov, N.T. (2015). Sparse factor analysis for identifying optimal perfect simple structure.
Proceedings of the 2013 IASC satellite for the ISI WSC and the 8th IASC-ARS conference: pp. 285-290.
[4] Adachi, K. and Trendafilov, N.T. (2015). Sparse principal component analysis subject to prespecified cardinalityof loadings. Computational Statistics, http://link.springer.com/article/10.1007/s00180-015-0608-4 (in press).
[5] Jennrich, R. I. (2002). A simple general method for oblique rotation. Psychometrika, 67(7): 7-20.
[6] Rubin, D. B. and Thayer, D. T. (1982). EM algorithms for ML factor analysis. Psychometrika, 47(1): 69-76.
[7] Unkel, S., Trendafilov, N.T.: Simultaneous parameter estimation in exploratory factor analysis: An expository
review. International Statisitical Review, 78(3): 363-382 (2010).
[8] Yeung, K. Y. and Ruzzo, W. L. (2001). Principal component analysis for clustering gene expression data.
Bioinformatics, 17(9): 763-774.