Distributed Learning for Principle Eigenspaces without Moment Constraints

Principal component analysis (PCA) is one of the most important statistical tool for dimension reduction, which extracts latent principal factors that contribute to the most variation of the data. For the classical setting with fixed dimension, the consistency and asymptotic normality of empirical principal components have been raised since Anderson (1963). In the last decade, high-dimensional PCA gradually attracted the attention of statisticians, see for example Onatski (2012), Wang & Fan (2017), Kong et al. (2021), Bao et al. (2022). The existing work on high-dimensional PCA typically assume the Gaussian/sub-Gaussian tail property of the underlying distribution, which is really an idealization of the complex random real world. Heavy-tailed data are common in research fields such as financial engineering and biomedical imaging (Jing et al., 2012; Kong et al., 2015; Fan et al., 2018; Li et al., 2022). Thus it is desperately needed to find a robust way to do principal component analysis for heavy-tailed data. Elliptical family provides one a proper way to capture the different tail behaviour of various common distributions such as Gaussian and $t$-distribution, and has been widely studied in various high-dimensional statistical problems, see for example, Han & Liu (2014); Yu et al. (2019); Hu et al. (2019); Chen et al. (2021). Han & Liu (2018) presented a robust alternative to PCA, called Elliptical Component Analysis (ECA), for analyzing high dimensional, elliptically distributed data, in which multivariate Kendall’s tau matrix plays a central role. In essence, for elliptical distributions, the eigenspace of the population Kendall’s tau matrix coincides with that of the scatter matrix The multivariate Kendall’s tau is first introduced in Choi & Marden (1998) for testing dependence, and is later adopted for estimating covariance matrix and principal components (Taskinen et al., 2012). Properties of ECA in low dimension were considered by Hallin et al. (2010), Hallin et al. (2014), and some recent works related to high-dimensional ECA include but not limited to Feng & Liu (2017); Chen et al. (2021); He et al. (2022). Noticeably, He et al. (2022) extended ECA to factor analysis under the framework of elliptical factor model.

With rapid developments of information and technology, the modern datasets exhibit the characteristic of extremely large-scale and we are now embracing the big-data era. Efficient statistical inference algorithms on such enormous dataset is unprecedentedly desirable. Distributed computing provides an effective way to deal with large-scale datasets. In addition to computation efficiency, distributed computing is also relatively robust to possible failures in the subsevers. There are lots of other reasons for establishing rigorous statistical theories on distributed computing, such as privacy protection, data ownerships and limitation of data storage. Two patterns of data segmentation were considered in the literature over the past few years, which are horizontal and vertical. “Horizontal” reserves all the features in each subserver, while data are scattered to different machines. Conversely, “Vertical” means that features are divided into several parts, each storage has full access of data but part of the features, which is common in signal processing and sensor networks. Some representative works include but not limited to Qu et al. (2002), Zhang et al. (2012), Fan et al. (2019) and Fan et al. (2021). In particular, Fan et al. (2019) proposed a distributed PCA algorithm: first each submachine computes the top $K$ eigenvectors and transmits them to the central machine; then the central machine aggregates the information from all the submachines and conducts a PCA based on the aggregated information.

In the current work, we consider a robust alternative to the existing distributed PCA algorithm. We adopt the horizontal division and propose a distributed algorithm for estimating the principle eigenspaces without any moment constraint on the underlying distribution. In detail, assume that we have $m$ subservers under the elliptical family framework, we first extract the leading $K$ orthonormal eigenvectors (as columns of $\widehat{\bm{V}}_{K}^{(i)}$) of the local sample multivariate Kendall’s tau matrix with observations on the $i$-th subserver, for $i=1,\ldots,m$. We next transport these eigenspace estimators $\{span(\widehat{\bm{V}}_{K}^{(i)})\}$ to the central server, where $\{span(\widehat{\bm{V}}_{K}^{(i)})\}$ denotes the linear space spanned by the columns of $\widehat{\bm{V}}_{K}^{(i)}$. In essence, the transported local eigenspace estimators $\{span(\widehat{\bm{V}}_{K}^{(i)})\}$ can be viewed as points (representatives of equivalent class) in the Grassmann manifold, see Figure 1 for better illustration. We then find the “center” ( the red point $span(\widetilde{\bm{V}}_{K})$ in Figure 1) of these points by minimizing the projection metrics on Grassmann manifolds in the central server, which is analogous to the physical notion of barycenter. Mathematically, the “center” $span(\widetilde{\bm{V}}_{K})$ satisfies

and it can be shown that $\widetilde{\bm{V}}_{K}$ is exactly the leading $K$ eigenvectors of the average projection matrix $\widetilde{\bm{\Sigma}}=1/m\sum_{i=1}^{m}\widehat{\bm{V}}_{K}^{(i)}\widehat{\bm{V}}_{K}^{(i)\top}$.

The contributions of the current work include the following aspects: firstly we for the first time propose a robust alternative to distributed principal PCA. The proposed distributed algorithm is easy to implement and avoids huge transportation costs between central processor and subservers; secondly we investigate the bias and variance for the robust distributed estimator and derive its convergence rate which depends on the effective rank of the scatter matrix, eigengap, and the number of submachines. We show that the distributed estimator performs as if we have full access of whole data; thirdly, we extend the distributed algorithm to the elliptical factor model in a horizontal partition regime, which is the first distributed algorithm for robust factor analysis as far as we know.

The rest of this paper is organized as follows. Section 2 includes some notations and preliminary results on elliptical family and multivariate Kendall’tau matrix. In Section 3, we present the details on our robust distributed algorithm for estimating principal eigenspace. In Section 4, we investigate the theoretical properties of the robust distributed estimator. In Section 5, we extend the distributed algorithm to the distributed elliptical factor analysis. Simulation results are given in Section 6. We also applied the distributed algorithm to analyze a large-scale real macroeconomic dataset. Further discussions and future work direction are left in Section 7.

We first introduce some notations used throughout the article. Vectors and matrices will be written in bold symbol, scalars are written by regular letters. For two sequences $\{a_{n}\}_{n=1}^{\infty},\{b_{n}\}_{n=1}^{\infty}$, $a_{n}\lesssim b_{n}$ if there is a universal constant $C>0$, such that $a_{n}\leq Cb_{n}$, for which we also adopt another notation $a_{n}=O(b_{n})$, while $a_{n}=o(b_{n})$ means $a_{n}/b_{n}\rightarrow 0$, as $n\rightarrow\infty$. Index set $\{1,2,\ldots,p\}$ is simply denoted by $[p]$. For vectors, we define $\bm{e}_{i}$ to be the unit vector with 1 in the $i$-th component, $\|\cdot\|_{p}$ is the $\ell_{p}$-norm. For random variable $X$, Orlicz $\psi_{1}$ norm $\|X\|_{\psi_{1}}$ is defined as $\mathop{sup}_{p\geq 1}(\mathbb{E}\lvert X\rvert^{p})^{1/p}/p$. If random vectors $\bm{X},\bm{Y}$ share the same distribution, we write $\bm{X}\overset{d}{=}\bm{Y}$. For matrix $\bm{A}\in\mathbb{R}^{p\times p}$, its transpose is $\bm{A}^{\top}$; $\lambda_{i}(\bm{A})$ is the $i$-th largest eigenvalue of $\bm{A}$; $\|\bm{A}\|_{2}$ is the spectral norm and $\|\bm{A}\|_{F}$ is the Frobenius norm. Given a fixed $K\in[p]$, the eigen gap of $\bm{A}$ is $\Delta(\bm{A})=\lambda_{K}(\bm{A})-\lambda_{K+1}(\bm{A})$, and $r=r(\bm{A}):=Tr(\bm{A})/\lambda_{1}(\bm{A})$ is the effective rank of $\bm{A}$. The space spanned by the columns of $\bm{A}$ is denoted by $span(\bm{A})$. For two matrices with orthogonal columns $\bm{A}\in\mathbb{R}^{p\times n_{1}},\bm{B}\in\mathbb{R}^{p\times n_{2}}$($n_{1},n_{2}\leq p$), the spatial distance between $span(\bm{A})$ and $span(\bm{B})$ is $\rho(\bm{A},\bm{B})=\|\bm{A}\bm{A}^{\top}-\bm{B}\bm{B}^{\top}\|_{F}$.

Elliptical family is a large distribution family containing common distributions such as Gaussian and $t$-distribution, which is widely used for modeling heavy-tailed data in finance and macroeconomics. To begin with, we recall its definition and some of its nice properties. For further details on the elliptical distribution, see for example Hult & Lindskog (2002).

In the following we give two equivalent definitions, one by its characteristic function, the other by its stochastic representation. A random vector $\bm{X}=(X_{1},\ldots,X_{p})^{\top}$ belongs to the elliptical family, denoted by $\bm{X}\sim ED_{p}(g;\bm{\mu},\bm{\Sigma})$, if its characteristic function has the form:

where $g(\cdot)$ is a proper function defined on $[0,\infty)$, $\bm{\mu}\in\mathbb{R}^{p}$ is the location parameter, $\bm{\Sigma}$ is called scatter matrix with $rank(\bm{\Sigma})=q$. Equivalently,

where $\bm{U}$ is a random vector uniformly distributed on the unit sphere $S^{q-1}$ in $\mathbb{R}^{q}$, $\xi\geq 0$ is a scalar random variable indepent of $\bm{U}$, $\bm{A}\in\mathbb{R}^{p\times q}$ is a deterministic matrix satisfying $\bm{AA}^{\top}=\bm{\Sigma}$. Thus we can also denote as $\bm{X}\sim ED_{p}(\bm{\mu},\bm{\Sigma},\xi)$ and $\xi$ and $g(\cdot)$ are mutually determined. Elliptically distributed random vectors have the same nice properties as Gaussian random vectors, such as linear combination of elliptically distributed random vectors also follows elliptical distribution and their marginal distributions are also elliptical.

For Gaussian distribution, scatter matrix is the covariance matrix up to a constant. For non-Gaussian distribution, even with infinite variance (such as Cauchy distribution), scatter matrix still measures the dispersion of a random vector. To estimate the principal eigenspace of elliptically distributed data, simply performing PCA to the sample covariance matrix is not satisfactory due to the possible inexistence of the population covariance matrix. To address the problem, we turn to a tool, the spatial Kendall’s tau matrix. For $\bm{X}\sim ED(\bm{\mu},\bm{\Sigma},\xi)$ and its independent copy $\widetilde{\bm{X}}$, the population spatial Kendall’s tau matrix of $\bm{X}$ is defined as:

The spatial Kendall’s tau matrix was first introduced in Choi & Marden (1998) and has been used for a lot of statistical problems such as covariance matrix estimation and principal eigenspace estimation, see Visuri et al. (2000), Fan et al. (2018) and Han & Liu (2018). A critical property of spatial Kendall’s tau matrix is that for elliptical vector $\bm{X}$, it shares the same eigenvectors of scatter matrix with the same ordering of eigenvalues, see Han & Liu (2018) for detailed proof. Suppose $\left\{\bm{X}_{t}\right\}_{t=1}^{n}$ are $n$ independent data points from $\bm{X}$, the sample spatial Kendall’s tau matrix is naturally defined as a U-statistic:

In this section, a new distributed algorithm for estimating principal eigenspace robustly is introduced in detail, which generalizes the distributed PCA by Fan et al. (2019) to the elliptical family setting.

Suppose there are $N$ samples of $\bm{X}\sim ED_{p}(\bm{\mu},\bm{\Sigma},\xi)$ in total, and we are interested in recovering the eigenspace spanned by the leading $K$ eigenvectors of the scatter matrix $\bm{\Sigma}$. The $N$ samples spread across $m$ machines, and on the $l$-th machine, there are $n=N/m$ samples. For $l\in[m]$, let $\{\bm{X}_{t}^{(l)}\}_{t=1}^{n}$ denote the samples stored on the $l$-th machine. First calculate the local sample spatial Kendall’ tau matrix on the $l$-th server as

We further compute the leading $K$ eigenvectors of sample Kendall’s tau matrix $\widehat{\bm{K}}^{(l)}$ on the $l$-th server and denote $\widehat{\bm{V}}_{K}^{(i)}$ as the matrix with columns being these eigenvectors. We next transport these $\widehat{\bm{V}}_{K}^{(i)}$ to the central server. The communication cost of the proposed distributed algorithm is of order $O(mKp)$. In contrast, to share all the data or entire spatial Kendall’tau matrix, the communication cost will be of order $O(mp\min(n,p))$. In most cases $K=o(min(n,p))$, the distributed algorithm requires much less communication cost than naive data aggregation.

The transported local estimator $\widehat{\bm{V}}_{K}^{(i)}$ would span an eigenspace $\{span(\widehat{\bm{V}}_{K}^{(i)})\}$, which can be viewed as points in the Grassmann manifold, as illustrated in Figure 1. We then find the “center” of these points by minimizing the projection metrics on Grassmann manifolds in the central server, i.e,

Denote $P_{\bm{V}}=\bm{V}\bm{V}^{\top}$ and $P_{\widehat{\bm{V}}_{K}^{(i)}}=\widehat{\bm{V}}_{K}^{(i)}\widehat{\bm{V}}_{K}^{(i)\top}$, which are the projection matrices of $span(\bm{V})$ and $span(\widehat{\bm{V}}_{K}^{(i)})$ respectively. Then we have:

The first two terms on the right hand side are fixed, so it’s equivalent to maximizing the last term, namely $\mathcal{L}(\bm{V})=\mathop{\mathrm{tr}}\left[\bm{V}^{\top}(\sum_{i=1}^{m}P_{\widehat{\bm{V}}_{K}^{(i)}})\bm{V}\right]$. Thus it’s clear that $\widetilde{\bm{V}}_{K}$ is the leading $K$ eigenvectors of the average projection matrix

At last we take $\widetilde{\bm{V}}_{K}$ as the final distributed estimator of the principal eigenspace. The detailed algorithm is summarized in Algorithm 1.

Fan et al. (2019) proposed a distributed PCA method in which they explain the aggregation procedure from the perspective of semidefinite programming (SDP) problem and minimize the sum of squared losses. In the current work, we provide a rationale of the aggregation procedure from the perspective of the Grassmann manifold, which is quite easy to understand. The viewpoint from the Grassmann manifold is similar to Li et al. (2022), where they proposed manifold PCA for matrix elliptical factor model but not from the perspective of distributed computation.

In this section, we investigate the theoretical property of the distributed estimator $\widetilde{\bm{V}}_{K}$. We focus on the distance of the spaces spanned by the columns of $\widetilde{\bm{V}}_{K}$ and those of $\bm{V}_{K}$. We adopt the metric in Fan et al. (2019) to analyse the statistical error of $span\{\widetilde{\bm{V}}_{K}\}$ as an estimator of $span\{\bm{V}_{K}\}$. At the beginning, we give some assumptions for technical analysis.

In the following analysis, we decompose the error $\rho(\widetilde{\bm{V}}_{K},\bm{V}_{K})$ into two parts: the bias term and the variance term. In detail,

where $\bm{V}_{K}^{*}$ is the $K$ leading eigenvectors of $\bm{\Sigma}^{*}=\mathbb{E}(\widehat{\bm{V}}_{K}^{(i)}\widehat{\bm{V}}_{K}^{(i)^{\top}})$. In the following we present the upper bound of the variance term and the bias term in Theorem 4.1 and Theorem 4.2 respectively.

The measure of concentration is essential in high-dimensional statistical analysis. Fan et al. (2019) exerted the sub-Gaussian assumption to derive a form of concentration inequality. However, the proposed algorithm maps each eigenspace estimator to the Grassmann manifolds where concentration inequalities are easily derived due to the compactness therein, thus the results also hold even for Cauchy distribution without any moment. Next we give the upper bound of the bias term.

By combining the results of the bias term and the variance term, we have the following corollary.

By Corollary 4.1 and the theoretical result in Han & Liu (2018), we conclude that our algorithm is not only communication-efficient but also outputs an estimator which shares the same convergence rate as the ECA estimator with all samples available at one machine when $m$ is not larger than the bounds given in Corollary 4.1. This is also confirmed in the following simulation study.

In this section, we assume that the data points $\bm{X}_{t}$ has a elliptical factor model, i.e., $\left\{\bm{X}_{t}\right\}_{t=1}^{N}\subset\mathbb{R}^{p}$,

where $\bm{X}_{t}=(X_{1t},\ldots,X_{pt})^{\top},\bm{f}_{t}\in\mathbb{R}^{K}$ are the latent factors, $\bm{L}=(\bm{l}_{1},\ldots,\bm{l}_{p})^{\top}$ is the factor loading matrix and $\bm{u}_{t}$’s are the idiosyncratic errors. Note that for the elliptical factor model, the loading space $span(\bm{L})$ rather than $\bm{L}$ is identifiable. We assume that random vectors $(\bm{f}_{t}^{\top},\bm{u}_{t}^{\top})^{\top}$ are generated independently and identically from the elliptical distribution $ED_{K+p}({\bm{0}},\bm{\Sigma}_{0},\xi)$, where $\bm{\Sigma}_{0}=\text{diag}(\mathbf{I}_{K},\bm{\Sigma}_{u})$, see He et al. (2022) for more details on elliptic distribution. By the property of the elliptical family, the scatter matrix of $\bm{X}_{t}$ has the form: $\bm{\Sigma}=\bm{L}\bm{L}^{\top}+\bm{\Sigma}_{u}$. We assume that the loading matrix $\bm{L}$ satisfies the pervasive condition $\bm{L}^{\top}\bm{L}/p^{\alpha}\rightarrow\bm{Q}$ as $p\rightarrow\infty$, where $\bm{Q}$ is a $K\times K$ positive definite matrix, and the eigenvalues of $\bm{\Sigma}_{u}$ are bounded from below and above. Then $\bm{\Sigma}$ satisfies Assumption 4.1. Naturally, we use the space spanned by the $K$ leading eigenvectors of the sample Multivariate Kendall’s tau matrix as the estimator of the loading space $span(\bm{L})$. In the article, we assume that the $N$ observations spread across $m$ machines, and on the $l$-th machine, there are $n=N/m$ observations. For $l\in[m]$, let $\{\bm{X}_{t}^{(l)}\}_{t=1}^{n}$ denote the observations stored on the $l$-th machine.

By the proposed distributed algorithm, we further propose a distributed procedure to estimate the factor loading space $span(L)$ and the factor scores. The algorithm goes as follows: firstly for $l=1,\ldots,m$, compute the $K$ leading eigenvectors $\{\widehat{\bm{v}}_{1}^{(l)},\ldots,\widehat{\bm{v}}_{K}^{(l)}\}$ of the sample Kendall’s tau matrix on the $l$-th server, and let $\widehat{\bm{V}}_{K}^{(l)}=(\widehat{\bm{v}}_{1}^{(l)},\ldots,\widehat{\bm{v}}_{K}^{(l)})_{p\times K}$; then transport $\widehat{\bm{V}}_{K}^{(l)}$ to the central processor; then perform spectral decomposition for the average of projection operators $\widetilde{\bm{\Sigma}}={1}/{m}\sum_{l=1}^{m}\widehat{\bm{V}}_{K}^{(l)}\widehat{\bm{V}}_{K}^{(l)^{\top}}$ to get its $K$ leading eigenvectors $\{\widetilde{\bm{v}}_{1},\ldots,\widetilde{\bm{v}}_{K}\}$, and let $\widetilde{\bm{V}}_{K}=(\widetilde{\bm{v}}_{1},\ldots,\widetilde{\bm{v}}_{K})$. At last $span(\widetilde{\bm{V}}_{K})$ is set as the final estimator of $span(\bm{L})$. If we are also interested in the factor scores on each server, then we simply transport $\widetilde{\bm{V}}_{K}$ back to each subserver, and solve a least squares problem to get the estimators of the factor scores. We summarized the procedure in Algorithm 2.

We proposed a robust distributed algorithm to estimate principal eigenspace for high dimensional, elliptically distributed data, which are stored across multiple machines. The distributed algorithm reduce the transmission cost, and performs much better than distributed PCA when data are extremely heavy-tailed. The algorithm goes as follows: in the first step, we compute the $K$ leading eigenvectors of sample Kendall’s tau matrix on each server and then transport them to the central server; in the second step, we take an average of the projection operators and perform a spectral decomposition on it to get the $K$ leading eigenvectors, which spans the final estimated eigenspace. We derive the convergence rate of the robust distributed estimator. Numerical studies show that the proposed algorithm works comparably with the full sample ECA while performs better than distributed PCA when data are heavy-tailed. We also extend the algorithm to learn the elliptical factor model in a distributed manner, and the theoretical analysis is more challenging and we leave it as future work. It’s also interesting to study the “vertical” division for distributed robust PCA and elliptical factor model.

The authors gratefully acknowledge National Science Foundation of China (12171282, 11801316), National Statistical Scientific Research Key Project (2021LZ09), Young Scholars Program of Shandong University, Project funded by China Postdoctoral Science Foundation (2021M701997).