A Linearly Convergent Algorithm for Distributed Principal Component Analysis

The modern era of machine learning involves leveraging massive amounts of high-dimensional data, which can have large computational and storage costs. To combat the complexities arising because of the high dimensions of data, dimensionality reduction and feature learning techniques play a pivotal and necessary role in information processing. The most common and widely used technique for this task is Principal Component Analysis (PCA) [2] which, in the simplest of terms, transforms data into uncorrelated features that aid conversion of data from a high-dimensional space to a low-dimensional space while retaining maximum information. Simultaneously, the enormity of the amount of available data makes it difficult to manage it at a single location. There are multiple and an increasing number of scenarios where data is distributed across different locations, either due to storage constraints or by its inherent nature like in the Internet-of-Things [3]. This aspect of the modern-world data have led researchers to explore distributed algorithms, which can process information across different locations/machines [4]. These aforementioned issues have motivated us to study and develop algorithms for distributed PCA that are efficient in terms of computations and communications among multiple machines, and that can also be proven to converge at a fast rate.

When the data is available at a single location, one of the goals of PCA is to find a $K$-dimensional subspace, given by the column space of a matrix $\mathbf{X}\in\mathbb{R}^{d\times K}$, such that the zero-mean data samples $\mathbf{y}\in\mathbb{R}^{d}\leavevmode\nobreak\ (d\gg K)$ retain maximum information when projected onto $\mathbf{X}$. In other words, when reconstructed as $\mathbf{X}\mathbf{X}^{\mathrm{T}}\mathbf{y}$ (subject to $\mathbf{X}^{\mathrm{T}}\mathbf{X}=\mathbf{I}$), the data samples should have minimum reconstruction error. It can be shown that this minimal error solution is given by the projection of data onto the subspace spanned by the eigenvectors of data covariance matrix. This implies that for dimensionality reduction, learning any basis of that subspace is sufficient. This is referred to as the subspace learning problem. But while simple dimension reduction does not necessarily need uncorrelated features, most downstream machine learning tasks like classification, pattern matching, regression, etc., work more efficiently when the data features are uncorrelated. In the case of image coding, e.g., PCA is known as the Karhunen–Loeve transform [5], wherein images are compressed by decorrelating neighbouring pixels. With this goal in mind, one needs to aim to find the specific directions that not only have maximum variance, but that also lead to uncorrelated features when data is projected onto those directions. These specific directions are given by the eigenvectors (also called the principal components) of the data covariance matrix, and not just any set of orthogonal basis vectors spanning the same space. Mathematically, along with minimum reconstruction error, the other goal of PCA is to ensure the condition that the off-diagonal entries of $\mathbb{E}[\mathbf{X}^{\mathrm{T}}\mathbf{y}\mathbf{y}^{\mathrm{T}}\mathbf{X}]$ are zero (i.e., the data gets decorrelated), while finding the eigenspace of the covariance matrix $\mathbb{E}[\mathbf{y}\mathbf{y}^{\mathrm{T}}]$.

As explained above, the true and complete purpose of PCA is served when the search for the optimal solution ends with the specific set of eigenvectors of the data covariance matrix, and not just with the subspace it spans. It is known that getting the principal components from any other basis of the subspace would only require performing singular value decomposition (SVD) of the obtained subspace. Although true, the SVD operation has a high computational complexity, which makes it an expensive step for big data. The traditional solutions for PCA were developed to overcome the cost of SVD and hence reverting back to it defeats the whole purpose.

Thus, even though the problem of dimensionality reduction of data has many optimal solutions (corresponding to all the sets of basis vectors spanning the $K$-dimensional space), our goal is to find only the ones that give the eigenvectors as the basis. In terms of optimization geometry of the PCA problem in which one tries to minimize the mean-squared reconstruction error under an orthogonality constraint, it is a non-convex strict-saddle function. In a strict-saddle function, all the stationary points except the local minima are strict saddles wherein the Hessians have at least one negative eigenvalue that helps in escaping these saddle points [6, 7]. Also, in the case of PCA the local minima are the same as the global minima. These geometric aspects make PCA, despite being non-convex, a “nice and solvable” problem whose optimal solution can be reached efficiently. However, note that the set of global minima contains, along with the set of eigenvectors as basis, all other possible bases that are rotated with respect to the eigenvectors. And our goal is not to find just any of the global minima but to look into a very particular subset of it, where the basis is not rotated.

A very popular tool that has been used to learn features of data, and hence compress it, is autoencoders. It was shown in [8] that the globally optimum weights of an autoencoder for minimum reconstruction error are the principal components of the covariance matrix of the input data. In [9], Oja described how using the Hebbian rule for updating the weights of a linear neural network would extract the first principal component from the input data. Several other Hebbian-based rules like Rubner’s model, APEX model [10], Generalized Hebbian Algorithm (GHA) [11], etc., were proposed to extend this idea of training a neural network for finding the first eigenvector to extract the first $K$ principal components (eigenvectors) of the input covariance matrix. Given the parallelization potential, a feedforward linear neural network-based solution for PCA seems to be very attractive.

The other aspect of modern day data is, as mentioned earlier, its massive size. Collating the huge amount of raw data is usually prohibitive due to communications overhead and/or privacy concerns. These reasons have encouraged researchers over the last couple of decades to develop algorithms that can solve various problems for non-collocated data. The algorithms developed to deal with such scenarios can be broadly classified into two categories: (1) the setups where a central entity/server is required to co-ordinate among the various data centers to yield the final result, and (2) the setup where the data is scattered over an arbitrary network of interconnected data centers with no availability of a central co-ordinating node. The authors in [12, 3] talk about these different setups and the algorithms developed for each of them in more detail. The second scenario is more generic and usually algorithms developed for such setups can be easily modified to be applied to the first scenario. The terms distributed and decentralized are used interchangeably for both the setups in the literature. In this paper, we focus on the latter scenario of arbitrarily connected networks and henceforth call it distributed setup. Hence, here our goal is to solve the PCA problem in the distributed manner when data is scattered over a network of interconnected nodes such that all the nodes in the network eventually agree with each other and converge to the true principal components of the distributed data.

Principal Component Analysis (PCA) aims at finding the basis of a low-dimensional space that can decorrelate the features of data points and also retain maximum information. More formally, for a random vector $\mathbf{y}\in\mathbb{R}^{d}$ with $\mathbb{E}\begin{bmatrix}\mathbf{y}\end{bmatrix}=\mathbf{0}$, PCA involves finding the top-$K$ eigenvectors of the covariance matrix $\boldsymbol{\Sigma}:=\mathbb{E}\begin{bmatrix}\mathbf{y}\mathbf{y}^{\mathrm{T}}\end{bmatrix}$. The zero mean assumption is taken here without loss of generality as the mean can be subtracted in case data is not centered. Mathematically, PCA can be formulated as

The constraint $\Big{(}\mathbb{E}\begin{bmatrix}\mathbf{X}^{\mathrm{T}}\mathbf{y}\mathbf{y}^{\mathrm{T}}\mathbf{X}\end{bmatrix}\Big{)}_{lq}=0,\forall l\neq q$, ensures that $\mathbf{X}$ decorrelates the features of $\mathbf{y}$. Now, $\mathbb{E}\begin{bmatrix}\mathbf{X}^{\mathrm{T}}\mathbf{y}\mathbf{y}^{\mathrm{T}}\mathbf{X}\end{bmatrix}=\mathbf{X}^{\mathrm{T}}\mathbb{E}\begin{bmatrix}\mathbf{y}\mathbf{y}^{\mathrm{T}}\end{bmatrix}\mathbf{X}=\mathbf{X}^{\mathrm{T}}\boldsymbol{\Sigma}\mathbf{X}$ and it is straightforward to see that this quantity is diagonal only if $\mathbf{X}$ contains the eigenvectors of $\boldsymbol{\Sigma}$. This explains why the search for a solution of PCA ends with the eigenvectors and not the subspace spanned by them. In practice, we do not have access to $\boldsymbol{\Sigma}$ and so a covariance matrix estimated from the samples of $\mathbf{y}$ is used instead. Specifically, for a dataset with $N$ samples $\{\mathbf{y}_{l}\}_{l=1}^{N}$, or equivalently, for a data matrix $\mathbf{Y}:=\begin{bmatrix}\mathbf{y}_{1},\mathbf{y}_{2},\dots,\mathbf{y}_{N}\end{bmatrix}$, the sample covariance matrix can be written as $\mathbf{C}=\frac{1}{N}\mathbf{Y}\mathbf{Y}^{\mathrm{T}}$ such that $\boldsymbol{\Sigma}:=\mathbb{E}\begin{bmatrix}\mathbf{C}\end{bmatrix}$. The true solution for PCA is then obtained by finding the eigenvectors of the covariance matrix $\mathbf{C}$, which are also the left singular vectors of the data matrix $\mathbf{Y}$. The empirical form of (1) is thus

In the literature, however, PCA is usually posed with a ‘relaxed’ orthogonality constraint of $\mathbf{X}^{T}\mathbf{X}=\mathbf{I}$ instead of $\Big{(}\mathbf{X}^{\mathrm{T}}\mathbf{Y}\mathbf{Y}^{\mathrm{T}}\mathbf{X}\Big{)}_{lq}=0,\forall l\neq q,$ as follows:

The optimization formulation in (3) with this constraint will only lead to a subspace spanned by the eigenvectors of $\mathbf{C}$ as the solution, thus actually making it a Principal Subspace Analysis (PSA) formulation. In other words, although the formulation (3) gives a solution on the Stiefel manifold, the actual PCA formulation (2) requires the solution to be within a very specific subset of that manifold that corresponds to the eigenvectors of $\mathbf{C}$. The accuracy of the solutions given by the PCA and PSA formulations will be the same when measured in terms of the principal angles between the subspace estimates and the true subspace spanned by the eigenvectors of the covariance matrix. Specifically, if $\mathbf{X}=\begin{bmatrix}\mathbf{x}_{1},\cdots,\mathbf{x}_{K}\end{bmatrix}$ is an estimate of the basis of the space spanned by the eigenvectors $\mathbf{Q}=\begin{bmatrix}\mathbf{q}_{1},\cdots,\mathbf{q}_{K}\end{bmatrix}$, then the principal angles between $\mathbf{Q}$ and $\mathbf{X}$ given by either (2) or (3) will be the same. But a more suitable measure of accuracy for any PCA solution should be the angles between $\mathbf{x}_{i}$ and $\mathbf{q}_{i}$ for all $i=1,\cdots K$, which motivates us to judge the efficacy of any solution with respect to this metric instead of the principal subspace angles.

In the distributed setup considered in this paper, we consider a network of $M$ nodes such that the undirected graph, $\mathcal{G}:=(\mathcal{V},\mathcal{E})$, describing the network is connected. Here $\mathcal{V}=\{1,2,\dots,M\}$ is the set of nodes and $\mathcal{E}$ is the set of edges, i.e., $(i,j)\in\mathcal{E}$ if there is a direct path between $i$ and $j$. The set of neighbors for any node $i$ is denoted by $\mathcal{N}_{i}$. Under the setup of samples being distributed over the $M$ nodes, let us assume that the $i^{th}$ node has a data matrix $\mathbf{Y}_{i}$ containing $N_{i}$ samples such that $N=\sum_{i=1}^{M}N_{i}$. Thus each node has access to only a local covariance matrix $\mathbf{C}_{i}=\frac{1}{N_{i}}\mathbf{Y}_{i}\mathbf{Y}_{i}^{\mathrm{T}}$ instead of the global covariance matrix but one can see that $N\mathbf{C}=\sum_{i=1}^{M}N_{i}\mathbf{C}_{i}$. In this setting, a straightforward approach might be that each node finds its own solution independent of the data at all the other nodes. While this might seem viable, this approach will have major drawbacks. Recall that the sample covariance $\mathbf{C}$ approximates the population covariance $\boldsymbol{\Sigma}$ at a rate of $\mathcal{O}(f(N^{-1}))$, where $f$ is some function (depending on the distribution) of the number of samples $N$ [33]. Since the local data has smaller number of samples than the global data, working with the local covariance matrix $\mathbf{C}_{i}$ alone instead of somehow using the whole data will lead to a larger error in estimation of the eigenvectors. Also, since uniform sampling from the underlying data distribution is not guaranteed in distributed setups, the samples at a node may end up being from a narrow part of the entire distribution, thus being more biased away from the true distribution. This invites the need for the nodes to collaborate amongst themselves in a way that all the data is utilized to find estimates of the eigenvectors at each node while ensuring that all the nodes agree with each other. Thus, for a distributed setting, the PCA problem in (1) can be rewritten here as

It is easy to see that $\sum_{i=1}^{M}f_{i}(\mathbf{X})=f(\mathbf{X})$. Also, in a distributed setup, each node $i$ maintains its own copy $\mathbf{X}_{i}$ of the variable $\mathbf{X}$ due to the difference in local information (local data) they carry. Thus, all nodes need to agree with each other to ensure the entire network reaches the same true solution. Hence, the true distributed PCA objective is written as

Note that (1)–(5) are non-convex optimization problems due to the non-convexity of the constraint set. One possible solution to the PCA problem is to instead solve a convex relaxation of the original non-convex function [34, 35]. The issue with these solutions is that they require $O(d^{2})$ memory and computation, which can be prohibitive in high-dimensional settings. In addition, due to $O(d^{2})$ iterate size these solutions are not ideal for distributed settings. Also, these formulations, without any further constraints, will not necessarily give a basis that is the set of dominant eigenvectors. Instead, they might end up giving a rotated basis as explained earlier, thereby not completing the task of decorrelating features. Hence, in this paper we use an algebraic method based on GHA for neural network training, which has $O(dK)$ memory and computation requirements, to solve the distributed PCA problem. Our goal is to converge to the true eigenvectors of the global covariance matrix $\mathbf{C}$ at every node of the network. As noted earlier in Section I-A, distributed variants of the power method exist in the literature [22, 24, 23] that can find the dominant eigenvector but these methods employ two time-scale approaches that involve several consensus averaging rounds for each iteration of the power method. Such two time-scale approaches can be expensive in terms of communications cost. In this paper, we propose a one time-scale method that finds the top $K$ eigenvectors of the global sample covariance matrix $\mathbf{C}$ at each node through local computations and information exchange with neighbors. The proposed method also converges linearly up to a neighborhood of the true solution when the error metric considered is the angle between the estimates and the true eigenvectors.

In [11], Sanger proposed a generalized Hebbian algorithm (GHA) to train a neural network and find the eigenvectors of the input autocorrelation matrix (same as the covariance matrix for zero-mean input). The outputs of such a network, when the weights are given by the eigenvectors, are the uncorrelated features of the input data that allow data reconstruction with minimal error, hence serving the true purpose of PCA. The algorithm was originally developed to tackle the centralized PCA problem in the case of streaming data, where a new data sample $\mathbf{y}_{t},t=1,2,\ldots$, arrives at each time instance.

In this paper we consider a batch setting, but the alignment of GHA with our basic goal of finding the eigenvectors motivates us to leverage it for our distributed setup. The rationale behind the idea of extrapolating the streaming case to a distributed batch setting is simple: since $\mathbb{E}[\mathbf{y}_{t}\mathbf{y}_{t}^{\mathrm{T}}]=\mathbb{E}[\mathbf{Y}_{i}\mathbf{Y}_{i}^{\mathrm{T}}]=\boldsymbol{\Sigma}$, the sample-wise distributed data setting can be seen as a mini-batch variant of the streaming data setting. In the context of neural network training, our approach can be viewed as training a network at each node with a mini-batch of samples in a way that all nodes end up with the same trained network whose weights are given by the eigenvectors of the autocorrelation matrix of the entire batch of samples.

The iterate for the GHA as given in [11] has the following update for the matrix of eigenvectors (i.e., the neural network weight matrix) $\mathbf{X}$ when the $t^{th}$ sample $\mathbf{y}_{t}$ arrives at the input of the neural network:

where $\boldsymbol{\mathcal{U}}:\mathbb{R}^{K\times K}\rightarrow\mathbb{R}^{K\times K}$ is an operator that sets all the elements below the diagonal to zero and $\alpha_{t}$ is the step size. For $K=1$, and defining $\boldsymbol{\Sigma}_{t}=\mathbf{y}_{t}\mathbf{y}_{t}^{\mathrm{T}}$, it was shown in [9] that the term $(\mathbf{X}^{(t)})^{\mathrm{T}}\mathbf{y}_{t}\mathbf{y}_{t}^{\mathrm{T}}\mathbf{X}^{(t)}=(\mathbf{X}^{(t)})^{\mathrm{T}}\boldsymbol{\Sigma}_{t}\mathbf{X}^{(t)}$ is the consequence of a power series approximation of Oja’s rule in lieu of the explicit normalization used in the case of the power method. In the case of $K>1$, $\boldsymbol{\mathcal{U}}\Big{(}(\mathbf{X}^{(t)})^{\mathrm{T}}\mathbf{y}_{t}\mathbf{y}_{t}^{\mathrm{T}}\mathbf{X}^{(t)}\Big{)}=\boldsymbol{\mathcal{U}}\Big{(}(\mathbf{X}^{(t)})^{\mathrm{T}}\boldsymbol{\Sigma}_{t}\mathbf{X}^{(t)}\Big{)}$ helps combine Oja’s algorithm with the Gram–Schmidt orthogonalization step as follows:

Thus, for any $k=1,\ldots,K$, the term involving $\boldsymbol{\mathcal{U}}(\cdot)$ in (6) includes an implicit normalization term $(\mathbf{x}_{k}^{(t)})^{\mathrm{T}}\boldsymbol{\Sigma}_{t}\mathbf{x}_{k}^{(t)}\mathbf{x}_{k}^{(t)}$ as well an orthogonalization term $\sum_{p=1}^{k-1}(\mathbf{x}_{p}^{(t)})^{\mathrm{T}}\boldsymbol{\Sigma}_{t}\mathbf{x}_{k}^{(t)}\mathbf{x}_{p}^{(t)}$, which—analogous to the Gram–Schmidt orthogonalization procedure—forces the estimate $\mathbf{x}_{k}^{(t)}$ to be orthogonal to all the estimates $\mathbf{x}_{p}^{(t)},p=1,\ldots,k-1$. Another important thing to note about the GHA algorithm is that, in order to estimate the dominant $K$ eigenvectors, it only requires the corresponding top $K$ eigenvalues to be distinct (and nonzero). In other words, it does not require the covariance matrix to be non-singular.

In the deterministic setting, where we have the full-batch instead of new samples every instance, this iterate changes to

Here, we term $\boldsymbol{\mathcal{H}}:\mathbb{R}^{d\times d}\times\mathbb{R}^{d\times K}\rightarrow\mathbb{R}^{d\times K},\boldsymbol{\mathcal{H}}(\mathbf{C},\mathbf{X}^{t}):=\Big{(}\mathbf{C}\mathbf{X}^{t}-\mathbf{X}^{t}\boldsymbol{\mathcal{U}}\Big{(}(\mathbf{X}^{t})^{\mathrm{T}}\mathbf{C}\mathbf{X}^{t}\Big{)}\Big{)}$ as the Sanger direction. An iterate similar to (8) has been proven to have global convergence in [17] for some very specific choice of the step sizes that are dependent on the iterate itself. Its straightforward extension to the distributed case is not possible as that would lead to different step sizes at different nodes of the network, making it difficult to talk about its convergence guarantees. Hence, to adapt this iterative method to our distributed setup, we use the typical combine and update strategy used quite richly in the literature for distributed algorithms such as [32, 36, 37, 38]. The main contributions of such works lie in showing that the resulting distributed algorithms achieve consensus (i.e., all nodes will have the same iterate values eventually) and, in addition, the consensus value is the same as the centralized solution. The convergence guarantees for these methods are mainly restricted to convex and strongly convex problems though. Our distributed version of (8) for PCA, which is non-convex, is based on similar principles of combine and update.

Specifically, the node $i$ at iteration $t$ carries a local copy $\mathbf{X}_{i}^{(t)}$ of the estimate of the eigenvectors of the global covariance matrix $\mathbf{C}$. In the combine step, each node $i$ exchanges the iterate values with its immediate neighbors $j\in\mathcal{N}_{i}$, where $\mathcal{N}_{i}$ denotes the neighborhood of node $i$, and then takes a weighted sum of the iterates received along with its local iterate. Then this sum is updated independently at all nodes using their respective local information. Since node $i$ in the network only has access to its local sample covariance $\mathbf{C}_{i}$, the update is in the form of a local Sanger’s direction given as

The details of the proposed distributed PCA algorithm, called the Distributed Sanger’s Algorithm (DSA), are given in Algorithm 1. The weight matrix $\mathbf{W}=[w_{ij}]$ in this algorithm is a doubly stochastic matrix conforming to the network topology [25] in the sense that for $i\neq j$, $w_{ij}\neq 0$ when $(i,j)\in\mathcal{E}$ and $w_{ij}=0$ otherwise. Also, $\forall i,w_{ii}\neq 0$, i.e., there is a self loop at each node. Note that connectivity of the network, as discussed in Section II, is a necessary condition for convergence of DSA. The connectivity assumption, in turn, ensures the Markov chain underlying the graph $\mathcal{G}$ is aperiodic and irreducible, which implies that the second-largest (in magnitude) eigenvalue of $\mathbf{W}$, $\beta=\max\{|\lambda_{2}(\mathbf{W})|,|\lambda_{M}(\mathbf{W})|\}$, is strictly less than $1$. While DSA shares algorithmic similarities with first-order distributed optimization methods [32, 39] in which the combine-and-update strategy is used, our challenge is characterizing its convergence behavior due to the non-convex and constrained nature of the distributed PCA problem. To this end, we first provide a general result in Section IV where we prove the convergence of a modified form of GHA. Then we utilize that result, along with some linear algebraic tools and additional lemmas provided in the appendices, to characterize the dynamics of the distributed setup in Section V and prove the convergence of the proposed algorithm.

Let $\mathbf{X}^{(t)}=\begin{bmatrix}\mathbf{x}_{1}^{(t)}&\mathbf{x}_{2}^{(t)}&\cdots&\mathbf{x}_{K}^{(t)}\end{bmatrix}\in\mathbb{R}^{d\times K}$, $K\leq d$, be an estimate of the $K$-dimensional subspace spanned by the eigenvectors of the covariance matrix $\mathbf{C}$ after $t$ iterations and $\mathbf{q}_{l},l=1,\ldots,d$, be the eigenvectors of $\mathbf{C}$ with corresponding eigenvalues $\lambda_{l}$. On expanding (8) using (7), it is clear that the GHA update equation for estimation of the $k^{th}$ eigenvector using a constant step size $\alpha$ is as follows:

We now slightly modify (10) by replacing $\mathbf{x}_{p}^{(t)}$ for $p<k$ by the true eigenvectors $\mathbf{q}_{p}$. We term the resulting update equation modified GHA and note that this is not an algorithm in the true sense of the term as it cannot be implemented because of its dependence on the true eigenvectors $\mathbf{q}_{p}$. The sole purpose of this modified GHA is to help in our ultimate goal of providing convergence guarantee for the DSA algorithm. The update equation of the modified GHA for “estimation” of the $k^{th}$ eigenvector of $\mathbf{C}$, $k=1,\ldots,K$, has the form

Note that similar to the original GHA, this modified GHA assumes that $\mathbf{C}$ has $K$ distinct eigenvalues, i.e., $\lambda_{1}>\lambda_{2}>\ldots>\lambda_{K}>\lambda_{K+1}\geq\cdots\geq\lambda_{d}\geq 0$. Now, since $\mathbf{q}_{l},l=1,\ldots,d$, are the eigenvectors of a real symmetric matrix, they form a basis for $\mathbb{R}^{d}$ and can be used for expansion of any $\mathbf{x}_{k}^{(t)}$ as

where $z_{k,l}^{(t)}$ is the coefficient corresponding to the eigenvector $\mathbf{q}_{l}$ in the expansion of $\mathbf{x}_{k}^{(t)}$. Multiplying both sides of (11) by $\mathbf{q}_{l}^{T}$ and using the fact that $\mathbf{q}_{l}^{T}\mathbf{q}_{l^{\prime}}=0$ for $l\neq l^{\prime}$, we get

This gives

It has been shown in [16] that the update equation given by

for $k=1$ converges to $\pm\mathbf{q}_{1}$ at a linear rate for a certain condition on the step size $\alpha$. Specifically, it was proven that $(z_{1,1}^{(t)})^{2}\rightarrow 1\quad\text{and}\quad\sum_{l=2}^{d}(z_{1,l}^{(t)})^{2}\leq b_{1}\rho_{1}^{t}$, where $b_{1}>0$ is some constant and $\rho_{1}=\big{(}\frac{1+\alpha\lambda_{2}}{1+\alpha\lambda_{1}}\big{)}^{2}<1$. Here, we extend the proof to a general $k$ and show that the update equation given in the form of (11) for any $k=1,\ldots,K,K<d$, converges to the $k^{th}$ dominant eigenvector.

With the analysis of the modified GHA in hand, let us proceed to analyze the proposed DSA algorithm. The iterate of DSA at node $i$ for the dominant $K$-dimensional eigenspace estimate ($K\leq d$) is given as

where $\mathbf{X}_{i}^{(t)}=\begin{bmatrix}\mathbf{x}_{i,1}^{(t)}&\mathbf{x}_{i,2}^{(t)}&\cdots&\mathbf{x}_{i,K}^{(t)}\end{bmatrix}\in\mathbb{R}^{d\times K}$ is an estimate of the $K$-dimensional subspace of the global covariance matrix $\mathbf{C}$ at the $i^{th}$ node after $t$ iterations, $\boldsymbol{\mathcal{H}}_{i}(\mathbf{C}_{i},\mathbf{X}_{i}^{(t)})$ is local Sanger’s direction, and $w_{ij}\geq 0$ is a weight that node $i$ assigns to $\mathbf{X}_{j}^{(t)}$ based on the connectivity between nodes $i$ and $j$ as mentioned before. The Sanger’s direction and the update equation for an estimate of the $k^{th}$ eigenvector is thus given as

Now, let the average of $\mathbf{x}_{1,k}^{(t)},\mathbf{x}_{2,k}^{(t)},\ldots,\mathbf{x}_{M,k}^{(t)}$ after $t^{th}$ iteration be denoted as $\bar{\mathbf{x}}_{k}^{(t)}=\frac{1}{M}\sum_{i=1}^{M}\mathbf{x}_{i,k}^{(t)}$ and given by taking average of (18) over all the nodes $i=1,\ldots,M$ as

where $\mathbf{h}_{k}^{(t)}=\frac{1}{M}\sum_{i=1}^{M}(\boldsymbol{\mathcal{H}}_{i}(\mathbf{x}_{i,k}^{(t)})-\boldsymbol{\mathcal{H}}_{i}(\bar{\mathbf{x}}_{k}^{(t)}))$. We present analysis of the DSA algorithm by first proving convergence of the average $\bar{\mathbf{x}}_{k}^{(t)}$ to a neighborhood of the eigenvector $\mathbf{q}_{k}$ of the global covariance matrix $\mathbf{C}$ while using a constant step size. Then with the help of Lemma 8 in Appendix E, which proves that the deviation of the iterates $\mathbf{x}_{i,k}^{(t)}$ at each node from the average $\bar{\mathbf{x}}_{k}^{(t)}$ is upper bounded, we prove that the iterates at each node also converge to a neighborhood of the true solution. It is noteworthy that the analysis of DSA does not require additional constraints on eigenvalues of $\mathbf{C}_{i}$, i.e., similar to GHA, we only require the top $K$ eigenvalues of $\mathbf{C}$ to be distinct and non-zero.

The complete proof of convergence of DSA is done by induction. First, we show the convergence of $\mathbf{x}_{i,1}^{(t)}$ to a $\mathcal{O}(\alpha)$ neighborhood of $\mathbf{q}_{1}$ and then analyze the rest of the eigenvector estimates $\mathbf{x}_{i,k}^{(t)},k=2,\ldots,K,$ by assuming that the higher-order estimates have converged.

Case I for Induction – $k=1$: The iterate for the dominant eigenvector is