Consistent Estimation of Low-Dimensional Latent Structure
in High-Dimensional Data

Low-dimensional latent variable models are often used to capture systematic structure in the conditional mean space of high-dimensional random variables (rv’s). This has been a popular strategy in high-dimensional probabilistic modeling and data analysis, and it serves as an attractive strategy for dimension reduction and recovering latent structure. Examples include factor analysis (Bartholomew et al., 2011), probabilistic principal components analysis (PCA) (Tipping and Bishop, 1999), non-negative matrix factorization (Lee and Seung, 1999), asymptotic PCA (Leek, 2011), latent Dirichlet allocation (LDA) (Blei et al., 2003), and exponential family distribution extensions of PCA (Collins et al., 2001).

Let $\mathbf{Y}_{k\times n}=\left(y_{ij}\right)$ be an observed data matrix of $k$ variables (one variable per row), each with $n$ observations, whose entries are rv’s such that

where $\mathbb{E}$ is the expectation operator, $\mathbf{M}_{r\times n}$ is a matrix of $r$ latent variables, and $\boldsymbol{\Phi}_{k\times r}$ is a matrix of coefficients relating the latent variables to the observed variables. Furthermore, the dimensions are such as $k\gg n\geq r$. In model (1), the conditional mean $\theta_{ij}$ of $y_{ij}$ for a fixed $i$ only depends on the $j$th column $\mathbf{m}^{j}$ of $\mathbf{M}$, and each row of the conditional mean matrix $\boldsymbol{\Theta}$ lies in the row space $\Pi_{\mathbf{M}}$ of $\mathbf{M}$. The latent structure in $\boldsymbol{\Theta}$ is therefore induced by $\Pi_{\mathbf{M}}$.

The above model is a general form of several highly used models. This includes instances of factor analysis (Bartholomew et al., 2011), probabilistic PCA (Tipping and Bishop, 1999), mixed membership clustering in population genetics (Pritchard et al., 2000; Alexander et al., 2009) which is closely related to LDA, and non-negative matrix factorization (Lee and Seung, 1999). Whereas the specialized models are often focused on the probabilistic interpretation of the columns of $\mathbf{M}$, we are instead here interested in its row space, $\Pi_{\mathbf{M}}$. This row space is sufficient for: (i) characterizing systematic patterns of variation in the data $\mathbf{Y}$, which can be used for exploratory data analysis or dimension reduction; (ii) accounting for the latent variables in downstream modeling of $\mathbf{Y}$ (Leek and Storey, 2007, 2008) that requires adjustment for these variables; (iii) potentially identifying suitable initial values or geometric constraints for algorithms that estimate probabilistically constrained versions of $\mathbf{M}$; (iv) or recovering $\mathbf{M}$ itself if additional geometric properties are known (e.g., as in Arora et al. (2013)). Furthermore, many of the above models make assumptions about the probability distribution of $\mathbf{M}$ that may be untrue on a given data set. One can compare our estimate of $\Pi_{\mathbf{M}}$ (which makes minimal assumptions) to the space induced by the model-based estimates of $\mathbf{M}$ to gauge the accuracy of the model assumptions and fit. Therefore, we focus on estimation of the latent variable space $\Pi_{\mathbf{M}}$. Estimation of $\boldsymbol{\Phi}$ may be tractable, but we do not focus on that here.

Leek (2011) and Anandkumar et al. (2012, 2015) have carried out work that is complementary to that presented here. They both study moment estimators of linear latent variable models applied to high-dimensional data. We explain how our work is related to Leek (2011) in Section 4.3 and how it is related to Anandkumar et al. (2012, 2015) in Section A. The strategies employed in these papers have ties to what we do here; however, they each consider different probabilistic models with theory that does not directly apply to the models we study.

We show that both the row space $\Pi_{\mathbf{M}}$ of $\mathbf{M}$ (in Section 2) and the row rank of $\mathbf{M}$ (in Section 3) can be consistently estimated using information from a suitably adjusted $n\times n$ matrix $k^{-1}\mathbf{Y}^{T}\mathbf{Y}$. In Section 4, we specialize these general results to $y_{ij}$’s that, conditional on $\mathbf{M}$, come from exponential family distributions. In particular, we explicitly construct a nonparametric, consistent estimator of the row space $\Pi_{\mathbf{M}}$ of $\mathbf{M}$ for $y_{ij}$ rv’s that follow the natural exponential family (NEF) with quadratic variance function (QVF) using information only up to their second moments, and the estimator is computationally straightforward to implement. In Section 5, we extend the results of previous sections to the case where $\boldsymbol{\Phi}$ is random. A simulation study is conducted in Section 6 to check and confirm our theoretical findings, and we apply the estimators to a genomics data set in Section 7. Finally, we end the article with a discussion in Section 8, collect in Section B all technical proofs, and present the the full set of results from simulation studies in Section C and Section D.

The big picture strategy we take is summarized as follows. Carrying out a decomposition related to Leek (2011), we first expand $k^{-1}\mathbf{Y}^{T}\mathbf{Y}$ into four components:

We then show the following about each of the components as $k\rightarrow\infty$, under the assumptions given in detail below:

• •

  $(i)$: this term may be estimated arbitrarily well as $k\rightarrow\infty$ by a diagonal matrix defined and studied below;

• •

  $(ii)$ and $(iii)$: these terms converge to the zero matrix;

• •

  $(iv)$: this term converges to a symmetric matrix with $r$ positive eigenvalues whose leading $r$ eigenvectors span the row space of $\mathbf{M}$.

Once the convergence of these terms is rigorously established, the strategy we take is to form an estimate of term $(i)$, denoted by $\mathbf{\hat{D}}_{k}$, and show that the space spanned by the leading $r$ eigenvectors of $k^{-1}\mathbf{Y}^{T}\mathbf{Y}-\mathbf{\hat{D}}_{k}$ converges to the row space of $\mathbf{M}$ as $k\rightarrow\infty$. We then also provide a framework to estimate the dimension $r$ of the row space of $\mathbf{M}$ and incorporate it into this estimation framework.

The strategy we take to consistently estimate the dimension of the latent variable space $r$ is to carefully scale the ordered eigenvalues of $\mathbf{R}_{k}=k^{-1}\mathbf{Y}^{T}\mathbf{Y}-\mathbf{D}_{k}$ and identify the index of the eigenvalue whose magnitude separates the magnitudes of these eigenvalues into two particular groups when $k$ is large. Recall that, by Theorem 2 and Corollary 3, the difference between the vector of descendingly ordered eigenvalues of $\mathbf{R}_{k}$ and that of those of $\mathbf{H}=\mathbf{M}^{T}\mathbf{WM}$ converges to zero as $k\rightarrow\infty$. However, since ${\rm rank}\left(\mathbf{H}\right)=r$ and $\mathbf{W}>0$, we know that the $r$ largest eigenvalues of $\mathbf{H}$ are positive but the rest are zero. This means that the $r$ largest eigenvalues of $\mathbf{R}_{k}$ are all strictly positive as $k\rightarrow\infty$, while the smallest $n-r$ eigenvalues of $\mathbf{R}_{k}$ converge to $0$ as $k\rightarrow\infty$. Depending on the speed of convergence of the $n-r$ smallest eigenvalues of $\mathbf{R}_{k}$ to $0$, if we suitably scale the eigenvalues of $\mathbf{R}_{k}$, then the scaled, ordered eigenvalues will eventually separate into two groups, those with very large magnitudes and the rest very small. The index of the scaled, ordered eigenvalues for which such a separation happens is then a consistent estimator of $r$. If we replace $\mathbf{R}_{k}$ with an estimator $\hat{\mathbf{R}}_{k}=k^{-1}\mathbf{Y}^{T}\mathbf{Y}-\hat{\mathbf{D}}_{k}$ that satisfies a certain level of accuracy detailed below, then the previous reasoning applied to the eigenvalues of $\hat{\mathbf{R}}_{k}$ will also give a consistent estimator of $r$.

To find the scaling sequence for the magnitudes of the eigenvalues of $\mathbf{R}_{k}$ or equivalently the speed of convergence of the $n-r$ smallest eigenvalues of $\mathbf{R}_{k}$ to $0$, we define the $k\times n$ matrix $\mathbf{E}$ with entry $\mathbf{E}(i,j)=e_{ij}=y_{ij}-\theta_{ij}$, and study as a whole the random part of $k^{-1}\mathbf{Y}^{T}\mathbf{Y}$ defined by

Note that the terms $(i)$, $(ii)$ and $(iii)$ in $\mathbf{F}$ correspond to those from equation (3). We will show that $\mathbf{F}$ configured as a $2n^{2}$ vector possesses asymptotic Normality after centering and scaling as $k\rightarrow\infty$. This then reveals that the scaling sequence for the eigenvalues of $\mathbf{R}_{k}$ should be no smaller than being proportional to $k^{-1/2}$.

Let $\boldsymbol{\Delta}=\left(\delta_{ij}\right)$, $\boldsymbol{\Delta}_{i}={\rm diag}\left\{\delta_{i1},...,\delta_{in}\right\}$ and define

for $1\leq i\leq k$, where $\mathbf{e}_{i}$ and $\mathbf{m}^{j}$ are respectively the $i$th row of $\mathbf{E}$ and the $j$th column of $\mathbf{M}$. We have:

With the concentration property of $\mathbf{F}$ established by Proposition 4, we will be able to explore the possibility of consistently estimating $r$ by studying the magnitudes of eigenvalues $\left\{\hat{\alpha}_{k,i}\right\}_{i=1}^{n}$ of

where each $\hat{\delta}_{kj}$ is an estimate of $\bar{\delta}_{kj}$ for $1\leq j\leq n$, and the $\hat{\alpha}_{k,i}$’s are ordered into $\hat{\alpha}_{k,i}\geq\hat{\alpha}_{k,i+1}$ for $1\leq i<n-1$.

Theorem 5 shows that the speed of convergence, $\left\{\tau_{k}\right\}_{k\geq 1}$, of the eigenvalues of $\hat{\mathbf{R}}_{k}$ (and those of $\mathbf{R}_{k}$) is essentially is determined by those related to $k^{-1}\boldsymbol{\Phi}^{T}\boldsymbol{\Phi}$ and $\mathbf{D}_{k}$; see equations (13)–(15). Further, it reveals that, when $\left\{\tau_{k}\right\}_{k\geq 1}$ is known, with probability approaching to $1$ as $k\rightarrow\infty$, the scaled, ordered eigenvalues $\left\{\tilde{\tau}_{k}^{-1}\hat{\alpha}_{k,i}\right\}_{i=1}^{n}$ eventually separates into two groups: $\tilde{\tau}_{k}^{-1}\hat{\alpha}_{k,i}$ for $1\leq i\leq r$ all lie above $\tilde{c}$ but the rest all lie below $\tilde{c}$ for any chosen $\tilde{c}>0$. In other words, for a chosen $\tilde{c}>0$, the number of $\tilde{\tau}_{k}^{-1}\hat{\alpha}_{k,i}>\tilde{c}$ when $k$ is large is very likely equal to $r$. However, in practice $\left\{\varepsilon_{k}\right\}_{k\geq 1}$ and $\left\{c_{k}\right\}_{k\geq 1}$ are unknown, and even if they are known or can be estimated, unfortunately the hidden constants in (14) are unknown (even when $k=\infty$). Further, when $k$ is finite, the $n-r$ smallest eigenvalues of $\mathbf{\hat{R}}_{k}$ may not yet be identically zero, the hidden constants may have a large impact on estimating the scaling sequence $\tilde{\tau}_{k}$, and rates slightly different than $\tilde{\tau}_{k}$ may have to be chosen to balance the effects of the hidden constants; see in Section 6 a brief discussion on the estimation of $\tilde{\tau}_{k}$ and the effects of choosing $\tilde{\tau}_{k}$ on estimating $r$ for finite $k$.

Before we apply our theory to $y_{ij}$’s that follow specific parametric distributions, we pause to comment on the key results obtained so far. Theorem 2 and Corollary 3 together ensure that asymptotically as $k\rightarrow\infty$ we can span the row space $\Pi_{\mathbf{M}}$ of $\mathbf{M}$ by the $r$ leading eigenvectors of $\mathbf{R}_{k}=k^{-1}\mathbf{Y}^{T}\mathbf{Y}-\mathbf{D}_{k}$ (see the definition of $\mathbf{D}_{k}$ in (5)). However, the conclusions in these results are based on a known $r$ and $\mathbf{D}_{k}$. In contrast, Proposition 4 and Theorem 5 retain similar assertions to those of Theorem 2 and Corollary 3 by replacing the unknown $\mathbf{D}_{k}$ by its consistent estimate $\hat{\mathbf{D}}_{k}$, and they enable us to construct a consistent estimate $\hat{r}$ of $r$ such that the linear space spanned by the leading $\hat{r}$ eigenvectors of $\hat{\mathbf{R}}_{k}=k^{-1}\mathbf{Y}^{T}\mathbf{Y}-\hat{\mathbf{D}}_{k}$ consistently estimates $\Pi_{\mathbf{M}}$ as $k\rightarrow\infty$. This reveals that to consistently estimate $\Pi_{\mathbf{M}}$ in a fully data driven approach using our theory, it is crucial to develop a consistent estimate $\hat{\mathbf{D}}_{k}$ of $\mathbf{D}_{k}$.

We specialize the general results obtained in Section 2 and Section 3 for model (1) to the case when $y_{ij}$ follow the single parameter exponential family probability density function (pdf) given by

where $\eta$ is the canonical link function, and $g$ and $h$ are known functions such that $f(y;\theta)$ is a proper pdf. The values that $\theta$ can take are $\Omega=\left\{\theta\in\mathbb{R}:\int f\left(y;\theta\right)=1\right\}$. The following corollary says that our general results on consistent estimation of the latent space $\Pi_{\mathbf{M}}$ and its dimension $r$ hold when the link function $\eta$ is bounded on the closure of the parameter space $\Omega_{0}\subseteq\Omega$.

We remark that the boundedness of $\eta$ on the closure of $\Omega_{0}$ is not restrictive, in that in practice either $\Omega_{0}$ is bounded in Euclidean norm or $\eta$ is bounded in supremum norm. The proof of Corollary 6 is a simple observation that $g$ is analytic in $\eta$ when $\theta\in\Omega_{0}$ (see, e.g. Letac and Mora, 1990) and that the derivative of $g$ in $\eta$ of any order is bounded when $\eta$ is bounded; the proof is thus omitted.

We now discuss the case where $\mathbf{\Phi}$ is random but then conditioned. Assume that $\boldsymbol{\Phi}$ is a random matrix with entries $\phi_{ij}$ defined on the probability space $\left(\Omega_{1},\mathcal{F}_{1},\mathbb{P}_{1}\right)$, and that the entries $y_{ij}$ of $\mathbf{Y}$, conditional on $\mathbf{\Phi}$ and $\mathbf{M}$, are defined on the probability space $\left(\Omega,\mathcal{F},\mathbb{P}\right)$. Rather than model (1), consider the following model:

Suppose assumption A4) holds (see (13)) and:

A1^′)

$1\leq{\rm rank}\left(\mathbf{M}\right)=r<n$ and $n$ is finite; conditional on $\mathbf{M}$ and $\boldsymbol{\Phi}$, $\left\{y_{ij}\right\}_{i,j}$ are jointly independent with variance $\mathbb{V}\left[y_{ij}|\mathbf{M},\boldsymbol{\Phi}\right]=\delta_{ij}$ such that $\sup\limits_{k}\max\limits_{i,j}\mathbb{E}\left[y_{ij}^{8}|\mathbf{M},\boldsymbol{\Phi}\right]\leq C$ (which implies $\sup\limits_{k}\max\limits_{i,j}\mathbb{E}\left[y_{ij}^{4}|\mathbf{M},\boldsymbol{\Phi}\right]\leq C$), where $\mathbb{E}$ and $\mathbb{V}$ are the expectation and variance wrt $\mathbb{P}$.

A2^′)

Either A2) holds $\mathbb{P}_{1}$-a.s., i.e., $\mathbb{P}_{1}\left(\sup_{k\geq 1}\max_{1\leq i\leq k}\left\|\boldsymbol{\phi}_{i}\right\|\leq C\right)=1$ and

$$\mathbb{P}_{1}\left(\left\|k^{-1}\boldsymbol{\Phi}^{T}\boldsymbol{\Phi}-\mathbf{W}\right\|=c_{k},\exists c_{k}\rightarrow 0,\exists\mathbf{W}_{r\times r}>0\right)=1,$$

(27)

or A2) holds in probability $\mathbb{P}_{1}$, i.e., as $k\rightarrow\infty$, $\mathbb{P}_{1}\left(\sup_{k\geq 1}\max_{1\leq i\leq k}\left\|\boldsymbol{\phi}_{i}\right\|\leq C\right)\rightarrow 1$ and

$$\mathbb{P}_{1}\left(\left\|k^{-1}\boldsymbol{\Phi}^{T}\boldsymbol{\Phi}-\mathbf{W}\right\|=c_{k},\exists c_{k}\rightarrow 0,\exists\mathbf{W}_{r\times r}>0\right)\rightarrow 1\text{.}$$

(28)

A3^′)

Conditional on $\mathbf{M}$, $k^{-1}\sum\limits_{i=1}^{k}\mathbb{E}\left[e_{ij}^{4}|\mathbf{M}\right]$, $k^{-1}\sum\limits_{i=1}^{k}\boldsymbol{\phi}_{i}\mathbb{E}\left[e_{ij}^{3}|\mathbf{M}\right]$ and $k^{-1}\sum_{i=1}^{k}\boldsymbol{\phi}_{i}\mathbf{m}^{j}\boldsymbol{\phi}_{i}\mathbf{m}^{l}\boldsymbol{\Delta}_{i}$ for any $1\leq j\leq l\leq n$ are all convergent $\mathbb{P}_{1}$-a.s. as $k\rightarrow\infty$.

Note that assumptions A1^′), A2^′) and A3^′) are the probabilistic versions of assumptions A1), A2 and A3) that also account for the randomness of $\mathbf{\Phi}$. Recall assumption A2) when $\boldsymbol{\Phi}$ is deterministic, i.e., for some $\mathbf{W}_{r\times r}>0$ and some non-negative sequence $c_{k}\rightarrow 0$, $\left\|k^{-1}\boldsymbol{\Phi}^{T}\boldsymbol{\Phi}-\mathbf{W}\right\|=c_{k}$. Assumption A2) implies that (see also Lemma 1):

These two uniform boundedness results in equation (29) are then used to show the a.s. convergence in Theorem 2 which induces the a.s. convergence in Corollary 3, the validity of Lindeberg’s condition as (34) that induces Proposition 4, convergence in probability in equations (14) and (16) that induces Theorem 5, Corollary 6, and Theorem 9.

Let $\tilde{N}_{1}\in\mathcal{F}_{1}$ be such that (27) does not hold, then $\mathbb{P}_{1}\left(\tilde{N}_{1}\right)=0$. On the other hand, if (28) holds, then for any positive constants $C$ and $\tilde{\varepsilon}$ there exists a $k_{0}=k_{0}\left(C,\tilde{\varepsilon}\right)$ such that the set

satisfies

whenever $k>k_{0}$. Now if (27) holds, then the results on the a.s. convergence and on convergence in probability for $y_{ij}|\mathbf{M}$ remain true with respect to $\mathbb{P}$ when $y_{ij}|\mathbf{M}$ is replaced by $y_{ij}|\left(\boldsymbol{\Phi},\mathbf{M}\right)$, as long as each $\phi_{ij}\left(\omega\right)$ is such that $\omega\in\Omega_{1}\setminus\tilde{N}_{1}$. In contrast, if (28) holds, then the results on the a.s. convergence for $y_{ij}|\mathbf{M}$ reduce to convergence in probability in $\mathbb{P}$. But those on convergence in probability for $y_{ij}|\mathbf{M}$ remain true with respect to $\mathbb{P}$ when $y_{ij}|\mathbf{M}$ is replaced by $y_{ij}|\left(\boldsymbol{\Phi},\mathbf{M}\right)$, as long as each $\phi_{ij}\left(\omega\right)$ is with $\omega\in\Omega_{1}\setminus N_{1}^{\ast}\left(C,\tilde{\varepsilon}\right)$ (where $\tilde{\varepsilon}$ can be chosen to be small). Therefore, the results (except Lemmas 7 and 8) obtained in the previous sections hold with their corresponding probability statements, whenever, for some $\mathbf{W}_{r\times r}>0$,

holds a.s. or in probability in $\mathbb{P}_{1}$, and they allow $\boldsymbol{\Phi}$ to be random (and conditioned) or deterministic. Statistically speaking, an event with very small or zero $\mathbb{P}_{1}$ probability is very unlikely to occur. This means that the practical utility of these results is not affected when (30) holds in the sense of (27) or (28) when $k$ is large.

Finally, we remark that Lemmas 7 and 8 are independent of the assumptions A1), A2) and A3), and of A1^′), A2^′) and A3^′). In other words, for model (26), Lemma 7 remains valid, and so does Lemma 8 when the assumptions involved there are replaced by those on $y_{ij}|\left(\boldsymbol{\Phi},\mathbf{M}\right)$.

We conducted a simulation study to demonstrate the validity of our theoretical results. The quality of $\hat{\Pi}_{\mathbf{M}}$ is measured by

where $\hat{r}$ is an estimate of $r$, $\mathbf{M}_{\mathbf{V}}=\mathbf{\hat{M}M}^{T}$ and $\mathbf{V}_{\mathbf{M}}=\left(\mathbf{MM}^{T}\right)^{-1}\mathbf{M\hat{M}}^{T}$. In fact, $d\left(\mathbf{M},\mathbf{\hat{M}}\right)$ measures the difference between the orthogonal projections (as linear operators) induced by the rows of $\mathbf{\hat{M}}$ and those of $\mathbf{M}$ respectively, and $d\left(\mathbf{M},\mathbf{\hat{M}}\right)=0$ if and only if $\hat{\Pi}_{\mathbf{M}}=\Pi_{\mathbf{M}}$. We propose a rank estimator $\hat{r}$, directly applied to $\mathbf{\hat{R}}_{k}$, based on Theorem 5. Specifically, $\tilde{\tau}_{k}=\tilde{c}_{k}k^{-\eta}$ for some positive $\tilde{c}_{k}$ and $\eta$ is used to scale the eigenvalues of $\mathbf{\hat{R}}_{k}$, for which $\tilde{c}_{k}$ (dependent on $k$) is determined by a strategy of Hallin and Liska (2007) (that was also used by Leek, 2011); for more details, we refer the reader to these two references. We do so since the assertions in Theorem 5 are valid only for $k=\infty$, assumption A3) may not be satisfied in our simulation study, and the unknown constants in equation (14) for the speed of convergence need to be estimated. However, we caution that, if the scaling sequence $\tilde{\tau}_{k}$ in Theorem 5 defined via (15) is not well estimated to capture the speed of convergence of the eigenvalues of $\hat{\mathbf{R}}_{k}$, $\hat{r}$ defined by (17) as an estimate of $r$ can be inaccurate; see our simulation study results.