{'fraction_non_alphanumeric': 0.04582496752096374, 'fraction_numerical': 0.01901499940947207, 'mean_word_length': 4.073217495506291, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 25, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A two-way factor model for high-dimensional matrix data

In this article, we introduce a two-way factor model for a high-dimensional data matrix and study the properties of the maximum likelihood estimation (MLE). The proposed model assumes separable effects of row and column attributes and captures the correlation across rows and columns with low-dimensional hidden factors. The model inherits the dimension-reduction feature of classical factor models but introduces a new framework with separable row and column factors, representing the covariance or correlation structure in the data matrix. We propose a block alternating, maximizing strategy to compute the MLE of factor loadings as well as other model parameters. We discuss model identifiability, obtain consistency and the asymptotic distribution for the MLE as the numbers of rows and columns in the data matrix increase. One interesting phenomenon that we learned from our analysis is that the variance of the estimates in the two-way factor model depends on the distance of variances of row factors and column factors in a way that is not expected in classical factor analysis. We further demonstrate the performance of the proposed method through simulation and real data analysis.


Introduction
The factor model, as a classical model in multivariate statistics, has been widely used in the undertaking of high-dimensional data analysis in a variety of scientific areas including finance, psychology, biology. By introducing latent variables (known as factors), a factor model assumes that observed variables are independent from each other after being introduced to small numbers of latent factors and, as a result, provides a simplified but meaningful framework to summarize the effects of latent factors as well as covariance structures between observed variables. Conventional factor-model-based methods focus mainly on analyzing vector-valued data, in which the observable attributes are converted into a vector, and the relevant observations are considered to be independent or approximately independent samples (Anderson and Rubin (1956); Anderson and Amemiya (1988), et al.). In recent studies, Bai (2003), , Fan et al (2008) and Fan et al (2011) have further generalized these methods under a high-dimensional framework and obtained a series of remarkable theoretical results. Nevertheless, these methods are faced with difficulties when applied to matrix data analysis. This is mainly because (a) there are often no replicates for an observation of matrix-based variables; (b) relationships between the attributes of rows and between those of columns should be considered separately; (c) vectorization procedure in matrix-based observations usually ignores distinguishing information about row and column, and thus cannot separate the effects of attributes along rows and columns. In terms of statistical methods that have been developed to address these issues, Gupta and Nagar (2000), Werner et al (2008), Leng and Tang (2012) and Ding and Cook (2018) focused mainly on a matrix data set with replicates. Tsai et al (2016) considered doubly constrained factor models for a data matrix. Though common factors could be interpreted by measuring effect of row, column and interaction according to their models, the doubly factor models indeed deal with vector-valued data if the known constrained matrices are absorbed by common factors or factor loadings. Wang et al (2016), Wang et al (2019) and Chen et al (2019) proposed factor-model-based methods for high-dimensional matrix-valued data with replicates, in which a low-sized matrix is used to represent the structures of hidden factors. However, their methods do not make an explicit separation between the row and the column effect and therefore can only provide evaluations of their joint behaviors. Faced with the reality of no replicates of matrix data, Zhou (2014) and Michael et al (2018) proposed two methods which are derived from a general matrix variate distribution, but they additionally require sparsity assumptions about the covariance matrix in order to obtain the parameter estimations as well as the large sample properties. Moreover, their ideas are essentially not factor domain.
In this work, we extend the idea of factor models to the analysis of high-dimensional matrix data, such as X p×q , with no replicates. To the best of our knowledge, this is the first attempt to apply factor analysis to studying complex correlation structures of matrix-wise variables with a single observation. Our interest stems from an environmental study in which volume readings of 14 chemicals (columns), including SO 2 , CO, collected from 338 cities (rows), are reported and the aim is to discover any patterns exhibited by cities and pollutants that could provide a systematic explanation of the status of the air, especially when extreme air pollution occurs. This motivates us to consider a comprehensive method that can simultaneously model possible pollution patterns through latent factors, decompose the information from the data matrix by rows and columns, and evaluate the effect of row and column factors based on the observed data matrix.
The core thinking behind our model is that the behaviors of each entry or variable in the data matrix X are affected by two groups of latent factors. One group summarizes the effects of the row attributes, while the other summarizes the effects of column attributes. The complex correlated relationships among the entries can then simply be decomposed and explained by the latent row and column factors. More specifically, we can assume that the data matrix X can be regarded as a sum of two unobservable matrices where U is the 'row effect' describing matrix and is made up of p independent row vectors, while V is the 'column effect' describing matrix and is made up of q independent column vectors. We further assume that the relationships between the variables in each row vector of U and those in each column vector of V can be respectively explained by row hidden factors F and column hidden factors E, U i. = LF i + η i , i = 1, ..., p, V .j = ΛE j + ξ j , j = 1, ..., q. As a result, we call this model two-way factor model (2wFM), because both row and column effects are described by hidden factors. Figure 1 (a) is a Bayesian network representation for 2wFM. It can be seen that 2wFM is essentially a generalization of classical factor models in analyzing matrix data. It inherits the dimension-reduction idea from classical factor models and can directly distinguish between row and column effects by introducing a hidden group of factors F and E. Figure 1 (b) further illustrates the differences between 2wFM and classical factor models on one-way scale. In the 2wFM framework, hidden factors affect the observable variables with a specific selection.
Our contributions to the factor analysis of high-dimensional matrix-valued data in this work exist in two parts. In the first we achieve maximum likelihood estimation (MLE) of all parameters in the settings of 2wFM. The specific structure in the covariance matrix of vec(X), denoted by Σ X , brings difficulties in terms of getting an analytical expression of the log-likelihood function. Here, we generalize a conclusion introduced by Miller (1981) on calculating the inverse of a special kind of matrix with 'a nonsingular matrix plus a singular matrix' form and obtain the exact form of Σ −1 X as well as |Σ X | under a group of identification conditions. Due to the entanglements that exist between factor loadings and the variance parameters of random factors and noises, we further implement a block alternating maximizing strategy to get the MLE for each parameter. The proposed algorithm includes alternatively updating factor loadings and the variance parameters.
The second part of our contributions refers to studying the theoretical properties of the MLE. Under general conditions, we derive the consistency properties as well as the asymptotic distribution, i.e., the central limit theorem, for each estimator. As can be seen in Subsection 2.5, factor loadings of row and column factors are combined with each other in the estimating equations, which makes it difficult to study each of them separately. Moreover, without replication information from the original data set, the convergence rates of many terms that constituted by the random factors and the estimated factor loadings cannot be directly identified. This fact has motivated us to undertake an in-depth study on the likelihood function and MLE. The final results present a phenomenon in which the variance of the estimates under the two-way factor model depends on the distance of variances of row factors and column factors in a way that could not have been expected in the classical factor analysis. The asymptotic variance of row factor loadings is a trade-off between its own variance and the variance of column factors (and vice versa). The distance between the variances of row and column factors has a heavy influence on the estimations of factor loadings, a very small difference may result in large fluctuations in the factor loading estimations. On the other hand, a small positive lower-bounded distance would lead to an effective asymptotic variance of both row and column estimated factor loadings.
The rest of the paper is organized as follows. In Section 2, we mathematically describe and explain our model, show the analytical expression of the likelihood function, analyze its basic structures, and present our block alternating algorithm so as to obtain the MLE for each parameter. In Section 3, we systematically discuss the main theoretical results for the estimators. Results from simulations and real data analysis are given in Section 4. All proof is displayed in the Appendix and supplementary material. Throughout the paper, A p×q denotes a p × q matrix. In particular, 0 q×r represents a q × r zero matrix, and 0 r is a zero vector with r entries. 1 r denotes a r dimensional vector with each entry being equal to 1, 1 r×c is a r × c matrix with each entry being equal to 1. e (m) k represents a k dimensional vector where the mth entry is 1 and other entries are 0. I k = e (1) k e (2) k · · · e (k) k denotes the k dimensional identity matrix. δ ij is a δ function. ||.|| 2 denotes the squared Euclidean distance of a matrix. For simplicity, if A is a r dimensional vector and there exists a δ > 0 such that each entry of A can be controlled by p δ as p → ∞, then A is written as A = O p (p δ )1 r . If A is a r × c matrix, we write A as O p (p δ )1 r×c .

The two-way factor model
Throughout this paper, we focus on covariance matrix analysis for high-dimensional matrix-valued data and ignore the mean-related parameter inference. We assume that the number of row and column factors are already known and do not need to be estimated. We first introduce our model in Subsection 2.1; we then illustrate the MLE problem and the analytical expression of the log-likelihood function in Subsection 2.2 and 2.3. Computation details for obtaining the MLEs are discussed in Subsection 2.4. Subsection 2.5 concerns on the estimating equations of MLE.

The two-way factor model and its identification conditions
In summary, our model can be directly written as Em , for k = 1, ..., r, m = 1, ..., c.
Here r and c are the numbers of row and column factors, L ≡ (L 1 · · · L r ) ∈ R q×r and Λ ≡ (Λ 1 · · · Λ c ) ∈ R p×c are the factor loadings for F and E, and L j ∈ R q , j = 1, . . . , r, Λ i ∈ R p , i = 1, . . . , c. To make the effects of each factor distinguishable, we make a common assumption, which has also appeared in several other works (Anderson and Rubin (1956), ), that the components of F and E have different variances, respectively, with an ordered relationship in (MC5). Similar to classical factor model, 2wFM assumes that the covariance structure between the entries of the data matrix X can be effectively described by the linear combinations of several hidden factors. When the components of L or Λ are all equal to zero, the above model would be reduced to the classical factor model.
where A = LΨ F L T is a q × q symmetric matrix, and B = ΛΨ E Λ T is a p × p symmetric matrix. It is natural to compare 2wFM to the well-known normal matrix variate model (MVM) where Φ records the covariance of each column in X, and Ω records the covariance of each row in X, as mentioned in Gupta and Nagar (2000), Adhikari (2007), Allen and Tibshirani (2012) and Zhou (2014). 2wFM uses a different approach in decomposing the covariance matrix of X. Not that provided by some kind of matrix decomposition, such as direct production in MVM or Singular Value Decomposition (SVD) as proposed in (Wang et al (2016)), 2wFM provides a direct separation to the row and column effects by assuming an additive structure existing between the row and column factors and the random noises. This setting induces the covariance matrix Σ X consisting of (the sum of) three parts. The first part is a full-ranked diagonal matrix corresponding to the noise structure, and the second and third parts are Kronecker products of an identity matrix with a low-rank matrix representing the sparse effect induced by row and column latent factors. This result can be further considered as a generalization to the covariance structure obtained under classical factor model settings.
Without any further restrictions, parameters in (2.1) cannot be directly identified. This is a common problem in classical factor models that has been studied carefully by Anderson and Rubin (1956) and . In the context of our model, since there are hidden factors for both row and column, this problem should be reconsidered. In the following proposition, we present the identification conditions which should be held throughout this work.
Identification Conditions. Let θ = (L, Λ, Ψ F , Ψ E , σ 2 ) be all the parameters. For model (2.1), there are two conditions (IC1) L T L qσ 2 = I r , Λ T Λ pσ 2 = I c ; (IC2) Two factor loadings, such as L and L , are considered as equivalent if and only if there exists a diagonal matrix D with diagonal entries equal to 1 or −1 such that L = L D.
(IC1) and (IC2) are two restrictions that are often imposed in factor analysis. Proposition 1 states that by ignoring the effect of two directions on a line in (IC2), it is enough for someone to impose restrictions on row and column factor loadings and the variance of random errors, and no more restrictions need to be considered for the parameter identification problem in model (2.1). In the next section, we will show that under (IC1) and (IC2), model (2.1) not only has concise form on likelihood function but also provides a convenient way to update the parameters in the iterative process of obtaining the MLE.

The MLE problem
For the model in (2.1), the log-likelihood function can be written as ln (θ; X) ∝ − ln |Σ X | − vec(X) T Σ −1 X vec(X). (2. 2) The maximum likelihood estimation for θ can then be defined aŝ (2.3) is a difficult optimization problem because (a) the number of parameters contained in θ is proportional to p and q and (b) the specific structure of Σ X makes the log-likelihood function in (2.3) difficult to obtain. Although several methods have been suggested for when this situation occurs (Sheena and Gupta, 2003;Vandenberghe and Boyd, 2004;Wainwright et al, 2006;Yuan and Lin, 2007;Yuan, 2009), they are only designed to handle problems, and most of them essentially follow indirect ways to handle the relevant optimization problems. Additionally, they often require more assumptions, such as sparsity, on the structure of Σ X in order to obtain the theoretical properties (Sheena and Gupta, 2003;Wainwright et al, 2006;Friedman et al, 2007;Bickel and Levina, 2008;Won et al, 2013;Dahl et al, 2008). In this work, we instead seek a more direct way to solve (2.3). We strive to obtain the closed form of Σ −1 X and |Σ X |, then provide an analytical expression of (2.2). As can be seen in Subsection 2.4, this results in a direct strategy for solving the problem (2.3).

The analytical expression of the log-likelihood function
A and B in Σ X can be written as To get Σ −1 X and |Σ X |, we generalize Miller's results (Miller (1981)) in the following proposition. Proposition 2. (a) Let W = G ⊗ I N + I M ⊗ E, where E is an N -dimensional semipositive definite matrix of rank r, and G is an M -dimensional positive definite matrix.
i and e i are respectively the ith eigenvalue and eigenvector of E with e T i e j = δ ij , δ ij being the delta function, is equal to 1 if i = j, otherwise, equal to 0.
In the case of there being single row and column factors, i.e., r = 1 and c = 1, Σ X can be simplified as By (a) in Proposition 2, Σ −1 X can then be written as . . , σ 2 ) are diagonal matrices with entries in the principal diagonal equal to the eigenvalues of A and B + σ 2 I p .
According to the property of the Kronecker product, Σ X can then be decomposed as As a result, |Σ X | can be written as When there are multiple row and column factors in 2wFM, we have Proposition 3. In general case, i.e., r ≥ 1, c ≥ 1, with (IC1) and (IC2), According to Proposition 3, the log-likelihood function (2.2) can be written as

Block alternating maximizing strategy for MLE
Here, θ is split into three groups L, Λ, (Ψ F , Ψ E , σ 2 ), and a block alternating maximizing strategy is proposed to calculate their MLE. This design induces to nondecreasing updates for the value of the log-likelihood function (2.4) (Proposition S1 in Subsection A.4 in the supplementary material). More specifically, updating each parameter group proceeds as follows: E , σ 2(0) (the method for choosing the initial values can be referred to in Subsection A.2 in the supplementary material for more detail) and set err 0 = 0.01 and 0 = 0.005.

Given
j L j subject to L T L = I r through the following iterative steps; , else t = t + 1 and repeat steps (2.2)-(2.3).
Remarks. 1. Given Λ, (Ψ F , Ψ E , σ 2 ) and L T L = qσ 2 I r , (2.4) can be written as .., r. When r = 1, (2.5) can be simplified as which is just a quadratic form of L. (2.6) is maximized when When r ≥ 1, (2.5) is the sum of a series of quadratic form as follows where W L 1 , . . . , W L r are positive semidefinite q × q matrices. It is well known that an analytical solution can be obtained for the problem of maximizing g(L) when W L 1 = · · · = W L r . However, orthogonal restrictions on each pair of (L i , L j ) i, j = 1, . . . , r, with distinct structures of W L 1 , . . . , W L r lead to difficulties in obtaining the solutions of maximizing (2.7). Here, two previous works (Bolla et al, 1998;Bolla, 2001) are introduced to overcome these difficulties. Bolla et al (1998) and Bolla (2001) indicate that for any solution of (2.7), there exists a r × r symmetric matrix A L such that this solution must satisfy (W L 1 L 1 · · · W L r L r ) = LA L . Thus, an iterative updating process can be implemented in three steps, (a) assembling obtaining SVD to W L and (c) getting an update of L. Bolla et al (1998) and Bolla (2001) proved that this process can make the value of the object function g converge to a local maximum. 2. Symmetrically, given L, (Ψ F , Ψ E , σ 2 ) and Λ T Λ = pσ 2 I c , (2.4) can be written as .., c. 3. Given L and Λ, updating (Ψ F , Ψ E , σ 2 ). As opposed to (2.5) and (2.8), there are generally no analytical solutions to updating (Ψ F , Ψ E , σ 2 ). Here, the EM method are adopted to update Ψ F , Ψ E and σ 2 in order to definitely be able to obtain a local optimal solution. The details can be referred to Subsection A.3 in the supplementary material.

Estimating equations ofθ
In this section, we study the estimating equations ofθ (Anderson and Rubin, 1956;. Although there are restrictions that (IC1) and (IC2) impose on θ, we can prove (using Lagrange multiplier techniques) that the properties ofθ that satisfyθ can be studied based on the estimating equations with (IC1) and (IC2): This means that the above equations include all the information fromθ in terms of discovering the relationships between each parameter. A detailed expression of each estimating equation is presented in the supplementary material (Subsection A.8).

Asymptotic properties of MLE
In this section, the statistical properties of the MLE in a large sample framework are discussed. It is important to note that the number of parameters in our model, such as parameters contained in factor loadings L and Λ, increases as p and q diverge. This fact broadly exists in high-dimensional data analysis and has been becoming a growing concern in recent theoretical studies. As illustrated by , methods that originate from the Taylor expansion such as the delta method cannot be applied directly to this situation. This is mainly because that the tail terms, which could be ignored in classical fixed or lower dimensional factor models, become the sum of infinity terms and each of them is o p (1) as p and q diverge. Thus, it is difficult to bound the tail terms in proper convergence orders.
Rather than approximating the distance between the true values of parameters θ * and their estimationsθ by using the variations of the linear part of the likelihood function valued in the local area ofθ, as is done in the Taylor expansion, we instead follow the technical route proposed by . Broadly speaking, we look for a direct algebraic decomposition of the log-likelihood function (2.4) such that the rates of each term obtained from the decomposition can be identified through technical analysis. However, due to the complexity of (2.4) as well as there being no replication information from the original data set, this is an extremely difficult task.
Before introducing our theoretical results in detail, we shall first form the following asymptotical conditions: (AC1) p/q → y ∈ (0, ∞); (AC2) L m· 2 = o(p 0.5 ), m = 1, . . . , q and Λ k· 2 = o(p 0.5 ), k = 1, . . . , p; (AC3) There exists a large enough positive constant C, such that Remarks. Conventional studies often impose assumptions on the relationships between the number of variables, p, and the sample size, n, while in matrix-valued data analysis, special attention must be paid to the size of the data set X. (AC1) sets p and q with an equal speed of divergence. We think that this restriction is natural and reasonable because if the diverging order of p, for example, is larger than q, then the parameters in the factor loadings Λ will increase far faster than those in L, and the information brought about by new introduced data may not be enough to distinguish and measure the effects of the row and column factors. On the other hand, inspired by the case in the classical factor model, people usually make assumptions that the order of p is no larger than the order of n Fan et al, 2008). We follow this method and find that this can induce a consistency result for each estimator. (AC2) and (AC3) are similar to (C.1) and (C.2) listed in . (AC2) is a generalized version of (C.1) in which each row entry of factor loadings L and Λ is O(1).
Proposition 4 presents the boundedness property for (Ψ F ,Ψ E ,σ 2 ). This condition is a basic result for deriving the large sample properties of the estimators, although it has often appeared as an underlying assumption in previous work Doz et al, 2012). Along with Theorem 1, a rigorous proof of boundedness has been given here. We compare the marginal values of the log-likelihood function (2.2) when Ψ F , Ψ E and σ 2 individually approach to 0 + and +∞, and to the value when θ locates at its true value θ * . We find that the values of ln (θ)/p and ln (θ)/q at θ * are bounded in probability and will tend to −∞ when Ψ F , Ψ E and σ 2 approach to 0 + or +∞, thus verifying the conclusion in Proposition 4.
Theorem 1 (Consistency). Under the model conditions (MC1)-(MC5), identification conditions (IC1)-(IC2), and asymptotic conditions (AC1)-(AC3), as p, q → ∞, we have (a) in simple case that r = 1, c = 1, Here, the idea of 'average convergence' is taken from  and Doz et al (2012) for the consistency expression of the factor loading estimators (L,Λ) since there will be infinity estimators for L and Λ as p and q diverge. This is a crucial step in the traditional M-method (in ) in proving sup θ∈Θ |R(θ)| = o p (1), where R(θ) is a remainder obtained after decomposing the original likelihood function ln (θ; X) into two parts ln (θ; X) = ln (θ; X) + R(θ), where ln (θ; X) is maximized at the true value θ * ; therefore, if sup θ∈Θ |R(θ)| = o p (1), it means thatθ and θ * can asymptotically be closed in some sense. However, in 2wFM, due to the fact that there are no replicates for X and the temporary lack of the boundedness property ofθ, this conclusion for R(θ) cannot be obtained. Therefore, more investigations into ln (θ; X) shall be undertaken.
The process of getting the consistency results is divided into three steps. In the first step, we focus on (L,Λ). A thorough study on the MLE is performed, which are induced by a degenerate form of (2.3) and (2.4), that is, As an optimization problem first appeared in Lemma 6A, solutions to (3.1), (L (I) ,Λ (I) ), are easy to obtain and have close connections to the real MLE (L,Λ). Relevant conclusions on (L (I) ,Λ (I) ) can be fully generalized to the case of (L,Λ) (Lemma 8). Based on the results regarding (L,Λ) from the first step and the order estimations of some basic terms, in the second step we go back to the original likelihood function (2.2), construct connections between (L,Λ) and (Ψ F ,Ψ E ,σ 2 ), and conclude that (Lemma 7) In the third step, the terms in these relations are connected to their true value θ * (Lemma 9) and are substituted into the estimating equation ofσ 2 , the consistency property of each estimator is then straightforward. Remarks. 1. Here, only the main results obtained from r = 1 and c = 1 are explained. In general case, r ≥ 1, c ≥ 1, all terms can be structured into fixed dimensional matrix form and the proof follows the same course. Lemma 6B, Lemma 7, Lemma 8 and Lemma 9 can be referred to for more details.
2. The key problem during the process of obtaining the consistency results for each estimator is to confirm the orders of the four terms L * TL qσ 2 , The first two can be regarded as a projection of the estimated factor loadings in the space constituted by their true values, their order estimations also being the key problem in one-way high-dimensional factor models . The final two terms can be regarded as a projection of the estimated factor loadings in the space constituted by row and column factors. In the classical high-dimensional factor model framework, only L * TL qσ 2 needs to be dealt with, while in the context of our model, it is necessary to simultaneously estimate the orders of these four terms, which makes our problem more general and more of a challenge. The optimization function in (3.1) is essentially a function of four highly similar terms L * TL qσ * 2 , which therefore obtains the essential structure of (3.1). As a result, the following conclusions can be drawn: This means the space comprised of L * or Λ * , respectively, keeps the most part ofL andΛ, which reveals a fact that we can deal with them separately.
Theorem 2 (Central Limit Theorem). Under the model conditions (MC1)-(MC5), identification conditions (IC1)-(IC2) and asymptotic conditions (AC1)-(AC3), as p, q → ∞, we have (a) in simple case that r = 1 and c = 1, (b) In general case that r ≥ 1, c ≥ 1, Theorem 2 presents the behavior of each estimator in a large sample scenario when both row and column factors exist. Results forσ 2 F andσ 2 E are similar to the results from . The size of σ * 2 F or σ * 2 E is positive correlation with the asymptotic variance ofσ 2 F orσ 2 E . Results onσ 2 in Theorem 2 show thatσ 2 is an asymptotic unbiased estimation for σ * 2 . Different from the results obtained from the classical factor model, the variances of row and column factors have interaction effects on the behaviors ofL andΛ. In detail, the asymptotic variance for each entry ofL when r = 1 and c = 1 can be written as Intuitively, a small σ * 2 indicates a small disturbance from noises, while a large σ * 2 F lets the value of F be distributed across a wide range, thus bringing more useful information on the inference of L. Therefore, a small σ * 2 or a large σ * 2 F will lead to a small asymptotic variance ofL. The expression from (3.2) agrees well with this judgement. Term 1+y∆ (∆−1) 2 in (3.2) represents the degree of impact from the alternative column factor E. Given σ * 2 , σ * 2 F and y, the size of (3.2) relates to the value of ∆, the ratio of σ * 2 F and σ * 2 E . It can be seen that a small distance between σ * 2 F and σ * 2 E will result in a large value of (3.2). Furthermore, (3.2) may be very large when ∆ is close to 1 (Figure 2 shows more details for the case in which σ * 2 = y = σ * 2 F = 1). As a result, it may need an extremely large number of p and q in order to make each entry ofL be consistent with its true value. On the other hand, a large distance between σ * 2 F and σ * 2 E , especially when σ * 2 F = 1, σ * 2 E → 0 + , will make ∆ close to 0, and the size of (3.2) in this situation will coincide with the conclusion derived from the same settings of the classical factor model. Similar conclusions can be obtained forΛ and for (L,Λ) in general case (r ≥ 1, c ≥ 1). Remarks.
1. In order to obtain the results in Theorem 2, we need to deal with estimating equations which report the restrictions thatθ must satisfy, and contain terms with important information on reflecting the limiting behaviors of each parameter inθ. However, these terms often remain together and appear in several equations as a series of cliques, and as a result, their orders cannot be identified directly. This difficulty motivate us to find a group of basic elements that connect these terms. It is done by sufficiently utilizing the relationships that exist in the estimating equations and substituting these relationships into the estimating equations of L and Λ. The whole process is cumbersome with constant decomposing of each term, calculating of the orders of new appeared terms, and keeping of the terms with higher or unknown orders while unifying the terms with known lower orders into one single term represented by their common upper bound. This process continues until there are no more elemental terms to be discovered. Finally, two terms, E T (L−L * ) , come out to the surface and it can be concluded that all the other terms are essentially the functions of these two terms (Lemma 12). We obtain their limiting distributions by solving a group of two linear equations that are induced during the simplifying process of the estimating equations of L and Λ (Lemma 13). Results for other terms are then straightforward.

Identifying the orders ofL
is also a key step in getting the limiting distributions of eachθ. As opposed to the results of  which show that the order ofL is a lower one compared to the order ofσ 2 , we find that this term has the same order asσ 2 − σ * 2 in our model settings.

Synthetic Data
In this section, we implement four groups of simulations. In the first and second simulation, we evaluate the accuracy of our estimators under single and multiple row and column factors, respectively. In the third simulation, we study the precision ofL andΛ as the ratio of σ 2 F and σ 2 E varies across a wide range. In the fourth simulation, we set F and E to follow Chi-square distributions and to test the robustness performance of our method.
In the accuracy study, we first set r = 1, c = 1, p = {50, 100, 200, 500, 1000}, q = {50, 100, 200, 500, 1000} and Ψ F = {1.5, 2, 4, 8}. For each pair of p, q and each Ψ F , we set Ψ E = 1, σ 2 = 0.01 (the same below) and draw 1000 samples (the same below) from the model (2.1). The elements of the factor loadings are first independently sampled from uniform distribution U [0, 1] and then normalized by the squared norms in order to make the resulting factor loadings satisfy the identification condition (IC1). People could refer to Subsection A.1 in supplementary material for the sampling process when all parameters are given. MLE are achieved according to the algorithm in Subsection 2.4. Table 1 Table S1 and S2 (in Subsection A.5). Since the number of parameters in L and Λ are proportional to q and p, we fit the estimated factor loadingsL andΛ to their true values and calculate the average R 2 as a measurement of the accuracies(the same below). ForΨ F ,Ψ E andσ 2 , we list the average value of each estimator with their mean absolute error (MAE) and mean square error (MSE) (the same below). In Table 1 and 2, there are two numbers in each grid. The first is the value of the MAE and the second one in brackets is the value of the MSE. The magnitudes of the MAE and MSE forσ 2 are respectively 10 −4 and 10 −7 (the same below). It can be seen that the precisions ofL andΛ are, respectively, closely related to the size of q and p. For each pair of p and q,L andΛ become more precise as Ψ F varied from 1.5 to 8.σ 2 converges to its true value more quickly than Ψ F andΨ E , whereas the value changes of Ψ F and Ψ E have little impact on the results forσ 2 . As p and q diverge, all estimators tend to approach their true values, and these results are consistent with the conclusions from Theorem 1. We also provide Q-Q plots in order to give more details on the behaviors of each estimator. Figure 3 provides Q-Q plots to give more details on the behaviors ofΨ E . A noticeable improvement can be found as p and q increase. We also obtain similar results for other estimators (see Figure S2, S3, S4 and S5 in Subsection A.5 in the supplementary material for more details). In the second simulation, we focus on the situation with multiple row and column factors. At this time, r is set to 2, c is set to 3, and Ψ F and Ψ E are set to diag(10, 8) and diag(6, 4, 2), respectively. All other parameters, p, q and σ 2 are set to the same values used in the first simulation. The results in Tables S3 and S4 (Subsection A.6) show that the number of factors, regardless of the rows or columns, have inverse effects on the precision ofL andΛ, and this phenomenon coincides with the conclusion for the asymptotic variance ofL andΛ in Theorem 2. On the other hand, this effect can be omitted if r and c are not too large compared to the sizes of p and q.
In the third simulation, we set r = 1, c = 1, Ψ E = 4, p = q ∈ {200, 500} and let ∆ = Ψ F Ψ E = {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1, 1.1, 1.2, 1.3, 1.5, 1.7, 2, 3, 5, 7}. For each pair (∆, p), we derive the MLE, obtain the average estimated accuracy ofL and Λ, then present the two accuracy curves in Figure 4. The curves indicate that the distance of the values of σ 2 F and σ 2 E has a serious impact on the precision ofL andΛ. With the same p or q, the asymptotic variances ofL andΛ become less as the ratio of σ 2 F and σ 2 E beyond 1 and as p or q increases, the accuracy ofL andΛ are more close to 1. Those facts confirm the conclusions from Theorem 2. Another interesting phenomenon, which also appears in the first simulation, is that, when ∆ < 1, the accuracy curves ofΛ are superior to these ofL (and vice versa).
To evaluate the robustness of our method, in the fourth simulation we set r = 1, c = 1, p = {200, 300, 400}, q = {200, 300, 400} and assume that factors E and F are generated from a Chi-square distribution with both means equal to 0 and Ψ E = 1 and Ψ F = {1.5, 2, 4}, respectively. The relative results are collected in Table S5 and S6 (see Subsection A.7 in the supplementary material). Our method performs well for each parameter and can estimate parameters with high precision even when common factors both F and E with non-normal distribution.

Real Examples
In the context of city air quality assessment, the standard method for measuring air quality is to calculate the air quality index (AQI) according to the volumes of several monitored pollutants, such as sulfur dioxide (SO 2 ) and nitrogen dioxide (N O 2 ). However, AQI only reports the maximum readings (linear transformed) of all the pollutants and does not consider the geographical relationships between cities. Here, our method is applied in order to give a new explanation of the air quality for each city.
Data are obtained from China National Environmental Monitoring Center website (http://www.cnemc.cn/). We selected a typical data set X containing 338 cities and 14 air pollutant indices. Figure 5(A) presents a heatmap showing the whole data set which reveals possible correlations existing between cities as well as between pollutants. We centralize X by columns, set r = 1 and c = 1 and apply our method to estimate F i , i = 1, . . . , 338 for each city and E j , j = 1, . . . , 14 for each pollutant. Figure 5(B) compares the AQIs withF i , i = 1, . . . , 338 in a scatterplot, which shows a high coincidence with the relative R 2 being as large as 0.895. This result indicates that F i , i = 1, . . . , 338, can be regarded as measurements of the air qualities. In order to investigate the impact of geographical factors, we apply K-means and classifyF i , i = 1, . . . , 338, into 6 clusters. The relevant results are collected in Figure 5(C), in which cities in the same cluster are labeled with the same color. Several areas, such as those located around Shandong and Henan province, with less-favorable air qualities can be identified intuitively. Figure  5(C) further indicates that our method can integrate geographic location information and can thus provide better results than AQI does. For example, Xi'an and Xianyang are two very close cities and are divided into different air quality levels by the AQIbased method, while our method considers them to be at the same level.
As well as illustrating the air quality level for each city, our method further provides an effective evaluation for each pollutant. Here, a new quantity, R 2 , is calculated, for each pollutant, which is the R 2 of the AQI regressed on the readings of 14 pollutants. Obviously, pollutants with large R 2 are probably important components of bad air quality. Figure 5(D) comparesÊ j with R 2 j , j = 1, . . . , 14. All pollutants are intuitively classified into 3 clusters. The six pollutants with the largest absolute values ofÊ j have the smallest R 2 j , whereas two pollutants (PM2.5 and PM10) with the smallestÊ j (near zero) have the largest R 2 j . This means that E j , j = 1, . . . , 14 can be considered as measurements that quantify the degree of each pollutant in the composition of the air pollution.

DISCUSSION
In this paper, we have developed a model-based method for high-dimensional matrix data analysis. Our model, called 2wFM, extends the application scope of the classical factor model from traditional vector-valued data to a general kind of matrixvalued data, in which there exists specific correlation structures between the attributes of rows and those of columns. We construct the identification conditions for 2wFM, derive an explicit expression to the likelihood function, and achieve maximum likelihood estimation for each parameter. Under general conditions, we study and obtain a series of theoretical studies on the resulted estimators, including consistency properties as well as asymptotic distributions. Our results provide detailed discussions on the relationships between the large sample behaviors of estimators and the statistical properties of hidden factors. Simulation studies further confirm our theoretical results and show that our method can efficiently estimate parameters with high precision.
Results on air quality data indicate that our method can extract and synthesize information from both pollutants and geographical locations, and can thus provide a more comprehensive evaluation of the level of air pollution than the traditional AQI-based methods. It would be of interest to develop a method that determines the number of row and column factors based on some statistical criteria such as BIC. Generalizations beyond the assumption of normal distributions for factors and noises may also be a considerable problem.

SUPPLEMENTARY MATERIALS
Supplementary material for "A two-way factor model for a high-dimensional data matrix". The simulation results, rigorous proof of Theorem 1 and Theorem 2 in simple case, as well as the relative lemmas with detailed proof that are useful in proving Theorem 1-2 are all provided in the supplementary material.

APPENDIX
Proof of Proposition 1. It suffices to prove that if Σ X = Σ X holds, then θ = θ . Recall that Σ X = I p ⊗ A + B ⊗ I q + σ 2 I p ⊗ I q , that is, where Σ X = (Σ ij ) p×p , Σ 11 , . . . , Σ pp are q × q matrices, and B = (b km ) p×p . Since Σ X = Σ X , we have Σ ii = Σ ii , i = 1, . . . , p, that is, Note that A + (b ii + σ 2 )I q is 'a low-rank plus a diagonal matrix' decomposition to Σ ii . It is well known that when min{p, q} > max{r, c}, this decomposition is unique under the identification conditions (IC1) and (IC2) (Anderson and Rubin, 1956). This leads to: (5.1) Because Σ X = Σ X , with (5.1), we further have which implies B + σ 2 I p = B + σ 2 I p . Note B + σ 2 I p is also 'a low-rank plus a diagonal matrix' decomposition to a known matrix. This leads to B = B , σ 2 = σ 2 . With (IC1), (IC2) and the above results, Λ = Λ , Ψ E = Ψ E , L = L and Ψ F = Ψ F are immediate consequences of the identification conclusions of the classical factor model.
Proof of Theorem 1.

Proof of Theorem 2.
Through Lemma 14, we have Following the technique proof in Lemma 13, we obtain Λ T XL pqσ 4 = O p (p −1 )1 c×r .
Recalling the estimating equation of L in Subsection A.8, the mth (m ∈ {1, . . . , q}) row of which can be written as With the expression ofL T X T XL pq 2σ6 in Lemma 9 and the above results, the final two terms of (5.21) can be simplified aŝ From the conclusions of (L − L * ) T E qσ 2 in Lemma 13, we have where Σ L = σ * 2 Ψ * −1 Similar results forΛ k· , k = 1, . . . , p can be obtained using the same techniques For the asymptotic distributions ofΨ F andΨ E , the estimating equation of σ 2 F j , j = 1, . . . , r, in Subsection A.8, is With Lemma 10 and Lemma 11, we havê (5.24) then reduces toΨ Similarly, we have