Online Distributed Estimation of Principal Eigenspaces

TL;DR

This paper introduces an online distributed PCA algorithm that efficiently estimates principal eigenspaces, demonstrating faster computation and maintained accuracy across tasks like low-rank approximation and clustering.

Contribution

It presents a novel online distributed PCA method with proven convergence rates, linking performance to network size, data rank, and spectral gaps.

Findings

Significant computational speed-up over existing methods

Maintains high accuracy in eigenspace estimation

Effective for low-rank approximation and clustering tasks

Abstract

Principal components analysis (PCA) is a widely used dimension reduction technique with an extensive range of applications. In this paper, an online distributed algorithm is proposed for recovering the principal eigenspaces. We further establish its rate of convergence and show how it relates to the number of nodes employed in the distributed computation, the effective rank of the data matrix under consideration, and the gap in the spectrum of the underlying population covariance matrix. The proposed algorithm is illustrated on low-rank approximation and $\boldsymbol{k}$-means clustering tasks. The numerical results show a substantial computational speed-up vis-a-vis standard distributed PCA algorithms, without compromising learning accuracy.

Equations

$${\Vert \Delta ({\tilde{\mathbf{V}}}_K(T),{\tilde{\mathbf{W}}}_K)\Vert }_{{\psi }_1}\lesssim \kappa \sqrt{\frac{Kr}{Tmn}}.$$

$${\Vert \Vert \tilde{\mathbf{\Sigma }}(T)-\mathbb{E}\tilde{\mathbf{\Sigma }}(T){\Vert }_F\Vert }_{{\psi }_1}\lesssim \frac{\alpha }{\sqrt{T}}.$$

Full text

**Index Terms— ** distributed estimation, principal components analysis, local and global aggregations, streaming data.

1 Introduction

Principal components analysis constitutes a fundamental technique for reducing the dimension of the data in numerous areas. It extracts a $K$ dimensional (affine) subspace that accounts for most of the variation in the data, thus revealing its low-dimensional structure. Further, it can serve as a preprocessing step to reduce the data dimension in various machine learning tasks including regression analysis, $\boldsymbol{k}$-means clustering [1], Non-negative Matrix Factorization (NMF) [2], and Latent Dirichlet Allocation (LDA) models [3].

In many applications, large data sets are distributed across different locations. This leads to some new challenges in both the design and analysis of the learning procedures. Indeed, different considerations need to be carefully addressed including local data collection, privacy/sharing issues, reliability, and communication costs. To address these issues, one needs to design and analyze distributed estimation procedures. A typical method computes low-rank approximations for subsets of the full data, and aggregates the resulting estimates at a central node [4]. A number of such techniques have been presented in the literature, including distributed sketching [5, 6, 7, 8, 9] and estimation [10, 11] methods. Recently, Fan et al. [11] studied a distributed estimation algorithm, wherein each local node computes the top $K$ eigenvectors for its own subset of the data, and then transmits them to a central node. Then, the central node runs an eigen-decomposition algorithm on the matrix comprising of all the eigenvectors received from the local nodes. Further, the estimation accuracy of the above method compared to PCA of the entire data is provided in the aforementioned reference.

In contrast to previous work on distributed estimation of principal eigenspaces, we are interested in an online scheme where the data becomes available in a sequential manner. We propose an online distributed algorithm consisting of two aggregation steps. In the first step, the incoming batch of the data is further split into subsets and allocated to different computation nodes. Each node computes a low-rank approximation of the corresponding subset, followed by a local aggregation step to obtain an estimate of the principal eigenspaces of the current batch. Then, the final estimate is obtained by passing the intermediate batch results to a fusion center that aggregates across all batches (global aggregation). Hence, the key contribution of this work is to introduce a fast distributed PCA algorithm, capable of providing high quality results for streaming data.

The remainder of the paper is organized as follows. Section 2 introduces the problem of estimating the principal eigenspaces. Then, we discuss algorithms for distributed, as well as online and distributed estimation of principal eigenspaces in Section 3. Section 4 establishes theoretical properties, while Section 5 illustrates the methodology on selected real data sets.

2 Problem Formulation

The objective is to learn the $K$ dimensional principal eigenspace of the unknown population covariance matrix $\boldsymbol{\Sigma}$ of the data. Technically, we want to estimate the subspace spanned by the eigenvectors corresponding to the top $K$ eigenvalues of $\boldsymbol{\Sigma}\in\mathbb{R}^{d\times d}$, for a given $K\in[d]\equiv\left\{1,\cdots,d\right\}$. The data consists of streaming i.i.d samples $\{\mathbf{X}_{i}\left(t\right)\}_{i=1}^{N}\subseteq\mathbb{R}^{d}$, being observed over the time interval $1\leq t\leq T$. Hence, the positive semidefinite matrix $\boldsymbol{\Sigma}$ corresponds to the population covariance matrix of the centered samples; ${\mathbb{E}}\mathbf{X}_{1}\left(t\right)=\mathbf{0}$, ${\mathbb{E}}\left[\mathbf{X}_{1}\left(t\right)\mathbf{X}_{1}^{\top}\left(t\right)\right]=\boldsymbol{\Sigma}$.

We further assume that the samples are generated from a sub-Gaussian distribution; $\|\|\mathbf{X}_{1}\left(t\right)\|_{2}\|_{\psi_{2}}<\infty$, where $\|\cdot\|_{\psi_{i}}$ denotes the Orlicz norm: $\|X\|_{\psi_{i}}=\sup_{p\geq 1}({\mathbb{E}}|X|^{p})^{1/p}p^{-1/i}$. To proceed, consider the eigenvalue decomposition $\boldsymbol{\Sigma}=\mathbf{V}\boldsymbol{\Lambda}\mathbf{V}^{\top}$, where $\boldsymbol{\Lambda}={\rm diag}(\lambda_{1},\lambda_{2},\cdots,\lambda_{d})$ with $\lambda_{1}\geq\lambda_{2}\geq\cdots\geq\lambda_{d}$, and $\mathbf{V}=(\mathbf{v}_{1},\cdots,\mathbf{v}_{d})\in{\cal{O}}_{d\times d}$ (the set of $d\times d$ orthonormal matrices). Further, let $\mathbf{V}_{K}=(\mathbf{v}_{1},\cdots,\mathbf{v}_{K})$ denote the orthonormal matrix comprising of the $K$ eigenvectors corresponding to the largest $K$ eigenvalues $\lambda_{j},j=1,\cdots,K$. Hence, we want to identify the subspace $\text{Col}(\mathbf{V}_{K})$, utilizing the streaming data $\{\mathbf{X}_{i}\left(t\right)\}_{i=1}^{N}$.

In order to ensure that $\text{Col}(\mathbf{V}_{K})$ is identifiable, we assume that there is a gap between the eigenvalues; $\delta:=\lambda_{K}-\lambda_{K+1}>0$. Throughout the paper, let $\kappa:=\lambda_{1}/\delta$ and $r=r(\boldsymbol{\Sigma}):={\rm Tr}(\boldsymbol{\Sigma})/\lambda_{1}$ denote the condition number and the effective rank of $\boldsymbol{\Sigma}$, respectively.

3 Algorithm

In the offline (i.e. non-streaming) setting where $T=1$, the standard procedure for estimating $\text{Col}(\mathbf{V}_{K})$ is to use the $K$ leading eigenvectors of the empirical covariance matrix***For $T=1$, we drop the trivial time indicator $(1)$. $\widehat{\boldsymbol{\Sigma}}=N^{-1}\sum_{i=1}^{N}\mathbf{X}_{i}\mathbf{X}_{i}^{\top}$. Let $\widehat{}\boldsymbol{\Sigma}=\widehat{}\mathbf{V}\widehat{}\boldsymbol{\Lambda}{\widehat{}\mathbf{V}}^{\top}$ be the eigenvalue decomposition of $\widehat{}\boldsymbol{\Sigma}$, analogous to that of $\boldsymbol{\Sigma}$ defined in the previous section. Then, $\text{Col}(\widehat{}\mathbf{V}_{K})$ is being used to estimate $\text{Col}(\mathbf{V}_{K})$.

Next, we briefly discuss the Distributed PCA algorithm (Dpca) [11]. In order to estimate the principal eigenspace $\text{Col}(\mathbf{V}_{K})$ in a distributed manner, the data is split across $m$ nodes, each possessing ${n}$ samples; i.e. $N=m{n}$. For each $\ell\in[m]$, the projection matrix $\widehat{\mathbf{V}}_{K}^{(\ell)}$ corresponding to the top $K$ empirical eigenvectors of node $\ell$ are formed. Then, $\bar{\boldsymbol{\Sigma}}=m^{-1}\sum\limits_{i=1}^{m}\widehat{\mathbf{V}}_{K}^{(\ell)}\widehat{\mathbf{V}}_{K}^{(\ell)^{\top}}$ is calculated, followed by its eigenvalue decomposition to obtain the $K$ leading eigenvectors of $\bar{\boldsymbol{\Sigma}}$, denoted by $\bar{\mathbf{V}}_{K}=(\bar{\mathbf{v}}^{(\ell)}_{1},\cdots,\bar{\mathbf{v}}^{(\ell)}_{K})$.

Clearly, the communication cost of Dpca is $O(NKd/{n})$. Indeed, every node needs to send the $d\times K$ matrix $\widehat{\mathbf{V}}^{(\ell)}_{K}$ to the central aggregation node. Thus, for large data sets, in addition to possible storage issues, there is a risk of communication congestion. These issues are further compounded in the presence of streaming data.

To address the above issues, we propose Algorithm 1 (Odpca) for estimating the principal eigenspaces in an online and distributed fashion. Let $N=Tmn$ denote the total number of samples, and let $\{\mathbf{X}_{i}^{(\ell)}(t)\}_{i\in[n]}$ be the batch of data that becomes available to node $\ell\in[m]$, at time $t\in[T]$. In the first round of aggregation all local nodes contribute to obtain the approximation $\bar{\boldsymbol{\Sigma}}(t)$ at time $t$. Then, the final estimate across all data batches, $\widetilde{\mathbf{V}}_{K}(T)$, is calculated by the fusion center. The communication cost for the local processors is divided by the number of time steps; so it is $O(NKd/\left(nT\right))$.

In order to analyze Odpca, first define $\Delta(\mathbf{U},\mathbf{V}):=\|\mathbf{U}\mathbf{U}^{\top}-\mathbf{V}\mathbf{V}^{\top}\|_{F}$, where $\|\cdot\|_{F}$ denotes the Frobenius norm, and $\mathbf{U},\mathbf{V}$ are arbitrary projection matrices with the same number of rows. The metric $\Delta$ is well-defined, and is commonly used in the literature to reflect the distance between subspaces, as well as the projections on them [12, 11]. Next, consider the online problem [13] $\min\limits_{\mathbf{U}\in{\cal O}_{d\times K}}H(\mathbf{U}),$ where $H(\mathbf{U}):=(Tm)^{-1}\sum_{t=1}^{T}\sum_{\ell=1}^{m}\Delta^{2}(\mathbf{U},\widehat{\mathbf{V}}_{K}^{(\ell)}(t))$. Intuitively, a minimizer $\mathbf{U}^{*}$ is a “center” of $\{\widehat{\mathbf{V}}^{(\ell)}_{K}(t)\}_{\ell\in[m],t\in[T]}$, which is by the following result a projection matrix on the principal eigenspace of $\boldsymbol{\Omega}:=(Tm)^{-1}\sum\limits_{t=1}^{T}\sum\limits_{\ell=1}^{m}\widehat{\mathbf{V}}_{K}^{(\ell)}(t){\widehat{\mathbf{V}}_{K}^{(\ell)}(t)}^{\top}$.

Lemma 1**.**

The matrix $\mathbf{U}^{*}\in{\cal O}_{d\times K}$ consists of the top $K$ eigenvectors of $\boldsymbol{\Omega}$, if and only if ${\mathbf{U}}^{*}\in{\rm argmin}_{\mathbf{U}\in{\cal O}_{d\times K}}H(\mathbf{U})$.

Proof.

Letting $\widehat{\mathbf{P}}^{(\ell)}(t)=\widehat{\mathbf{V}}_{K}^{(\ell)}(t)\widehat{\mathbf{V}}_{K}^{(\ell)\top}(t)$, we have $H(\mathbf{U})=\|\boldsymbol{\Omega}\|_{F}^{2}+\|\mathbf{U}\mathbf{U}^{\top}\|_{F}^{2}+(Tm)^{-1}\sum_{t=1}^{T}\sum_{\ell=1}^{m}\|\boldsymbol{\Omega}-\widehat{\mathbf{P}}^{(\ell)}(t)\|_{F}^{2}-2{\rm Tr}(\mathbf{U}\mathbf{U}^{\top}\boldsymbol{\Omega})$.

Then, $\|\mathbf{U}\mathbf{U}^{\top}\|_{F}^{2}$ is constant since $\mathbf{U}$ is orthonormal. Thus, a minimizer of $H(\cdot)$ maximizes ${\rm Tr}(\mathbf{U}\mathbf{U}^{\top}\boldsymbol{\Omega})={\rm Tr}(\mathbf{U}^{\top}\boldsymbol{\Omega}\mathbf{U})$; i.e. is the projection matrix to a principal subspace. ∎

In the next section, we show that Algorithm 1 provides an accurate approximation of $\mathbf{U}^{*}$. Indeed, the intermediate steps of calculating $\bar{\mathbf{V}}_{K}(t)$ lead to significantly faster computations with negligible growth in the statistical error.

4 Theoretical Analysis

Next, we analyze the learning error in terms of the distance $\Delta(\widetilde{\mathbf{V}}_{K}(T),\mathbf{V}_{K})$ defined in the previous section. The main result states that the statistical error rates are the same as using the full sample $\{\mathbf{X}_{i}^{(\ell)}(t)\}_{i\in[n],\ell\in[m],t\in[T]}$, as long as the sample size $n$ of each node $\ell$ at every time $t$ is large enough. Henceforth, we focus on the centralized error $\Delta(\widetilde{\mathbf{V}}_{K}(T),\widetilde{\mathbf{W}}_{K})$, where $\widetilde{\mathbf{V}}_{K}(T)$ is the output of Algorithm 1 (i.e. principal eigenvectors of $\widetilde{\boldsymbol{\Sigma}}(T)$), and $\widetilde{\mathbf{W}}_{K}$ consists of the top $K$ eigenvectors of ${\mathbb{E}}\left[\widetilde{\boldsymbol{\Sigma}}(T)\right]$. A comprehensive analysis of the deterministic bias $\Delta(\widetilde{\mathbf{W}}_{K},\mathbf{V}_{K})$ for both symmetric and asymmetric distributions is established before [11].

The following result addresses the behavior of the statistical error. Note that if $n$ is large enough, the assumptions of Theorem 1 will be satisfied [11].

Theorem 1**.**

Suppose that in Algorithm 1, $n\geq r$, and $\|\boldsymbol{\Omega}^{*}-\mathbf{V}_{K}\mathbf{V}_{K}^{\top}\|_{2}<1/4$, where $\boldsymbol{\Omega}^{*}:={\mathbb{E}}\left[\widehat{\mathbf{V}}_{K}^{(\ell)}(T){\widehat{\mathbf{V}}_{K}^{(\ell)}(T)}^{\top}\right]$. Them, we have†††the notation “$\lesssim$” states that “$\leq$” holds, modulo a universal constant.

[TABLE]

Proof.

Let $\mathbf{V}^{*}_{K}$ contain the top $K$ eigenvectors of $\boldsymbol{\Omega}^{*}$, and denote $\alpha=\kappa\sqrt{Kr(mn)^{-1}}$. Using Theorem 1 in [11], $\|\boldsymbol{\Omega}^{*}-\mathbf{V}_{K}\mathbf{V}_{K}^{\top}\|_{2}<1/4$ implies that $\left\|\Delta(\bar{\mathbf{V}}_{K}(t),\mathbf{V}_{K}^{*})\right\|_{\psi_{1}}\lesssim\alpha$, for all $t\in[T]$. Further, applying Jensen’s inequality to the definition of $\left\|\cdot\right\|_{\psi_{1}}$, we get $\left\|{\mathbb{E}}\left[\bar{\mathbf{V}}_{K}(t)\bar{\mathbf{V}}_{K}(t)^{\top}-\mathbf{V}_{K}^{*}{\mathbf{V}_{K}^{*}}^{\top}\right]\right\|_{F}\lesssim\alpha$. Therefore, by the triangle inequality, $\lambda_{K}\left({\mathbb{E}}\widetilde{\boldsymbol{\Sigma}}(T)\right)-\lambda_{K+1}\left({\mathbb{E}}\widetilde{\boldsymbol{\Sigma}}(T)\right)$ is a positive number bounded away from zero, as long as $mn$ is sufficiently large. Then, by the triangle inequality, $\left\|\lVert\bar{\mathbf{V}}_{K}(t)\bar{\mathbf{V}}_{K}(t)^{\top}-{\mathbb{E}}\left[\bar{\mathbf{V}}_{K}(t)\bar{\mathbf{V}}_{K}(t)^{\top}\right]\rVert_{F}\right\|_{\psi_{1}}\lesssim\alpha$ holds, which according to Lemma 4 in [11] implies that

[TABLE]

Thus, since the eigengap of ${\mathbb{E}}\widetilde{\boldsymbol{\Sigma}}(T)$ is bounded from below, according to Davis-Kahan’s Theorem (Theorem 2 in [12]) (4.1) implies the desired result. ∎

5 Real-Data Experiments

Next, we compare the performance of Odpca and Dpca [11]. The experiments correspond to solving rank–$K$ approximation and $\boldsymbol{k}$-means clustering tasks. We set $K=10$ for rank–$K$ approximation, and $k=10$ for $\boldsymbol{k}$-means clustering. The data sets being used in this section are as follows‡‡‡Obtained from the UCI Machine Learning Repository: http://archive.ics.uci.edu/ml.: (i) NewsGroups ($N=18774,d={61188}$); (ii) MNIST ($N=7\times 10^{4},d=784$); (iii) BOW NYTimes ($N=3\times 10^{5},d=102660$).

Algorithm 1 is implemented on a network of a star topology, comprising of $(m,T)=(5,10)$ for NewsGroups and MNIST datasets, and $(m,T)=(5,20)$ for BOW NYTimes.

We analyze both the computational cost, and the learning accuracy, under different scenarios. The following distributed procedure on the global data is being used as the baseline for assessing the performance of both Dpca and Odpca. First, each local node computes all $d$ eigenvectors of its own data, and transmits them to the central node. Then, this node aggregates the information collected from the local nodes, and computes the top $K$ eigenvectors.

The Dpca algorithm is implemented as follows. For each $\ell\in[m]$, the $\ell$-th node computes the projection matrix $\widehat{\mathbf{V}}_{K+Z}^{(\ell)}$ corresponding to the top $K+Z$ empirical eigenvectors of its own data, for some integer number $Z$. Then, the central node calculates $\bar{\boldsymbol{\Sigma}}=m^{-1}\sum\limits_{i=1}^{m}\widehat{\mathbf{V}}_{K+Z}^{(\ell)}\widehat{\mathbf{V}}_{K+Z}^{(\ell)^{\top}}$, followed by another round of eigenvalue decomposition to obtain the $K$ leading eigenvectors of $\bar{\boldsymbol{\Sigma}}$. A similar procedure is employed in the first round of Odpca, at every time $t\in[T]$. Intuitively, the goal of transmitting $Z$ additional or fewer eigenvectors is to demonstrate the computational costs of the algorithms. Henceforth, we refer to $K+Z$ as projection dimension.

For the low-rank approximation tasks, we report the estimation error of the solution obtained by Dpca and Odpca, normalized to that of the baseline algorithm described above. Fig. 1 shows the average results over $5$ replications. Note that the horizontal axis reflects the projection dimension $K+Z$. The left column illustrating the computational time suggests that Odpca is remarkably faster than Dpca. Indeed, the running time of Odpca improves over that of Dpca by a factor between 2 to 5. The column on the right-hand-side of Fig. 1 indicates the estimation performance. It is clear that the normalized error of Odpca is comparable to that of Dpca. Hence, the speed-up technique does not lead to any distinguishable sacrifice in the performance of Odpca.

Next, we illustrate the performance of Odpca on the distributed $\boldsymbol{k}$-means clustering. First, Dpca and Odpca compute the top $K$ eigenvectors, and send them to all nodes. Each node of the distributed $\boldsymbol{k}$-means algorithm projects its own data onto the new feature space, using the aforementioned $K$ eigenvectors. We run the algorithm in the work of Balcan et al. [14] over $5$ replications.

Similar to the previous task, we report the clustering cost of the solution obtained by Dpca and Odpca, to that of results provided by running Lloyd’s method on the full data. Fig. 2 shows the results for $\boldsymbol{k}$-means clustering tasks. It can be easily seen that Odpca’s solution is approximately as accurate as Dpca’s one. Further, a large decrease is observed in the running time of Odpca compared to Dpca.

Both figures indicate that the proposed local and global aggregation steps lead to fast computation of the principal eigenspaces, while preserving the accuracy almost the same.