Distributed Stochastic Algorithms for High-rate Streaming Principal Component Analysis
Authors: Haroon Raja
Waheed Bajwa
Venue: Accepted by TMLR ()

Abstract: This paper considers the problem of estimating the principal eigenvector of a covariance matrix from independent and identically distributed data samples in streaming settings. The streaming rate of data in many contemporary applications can be high enough that a single processor cannot finish an iteration of existing methods for eigenvector estimation before a new sample arrives. This paper formulates and analyzes a distributed variant of the classical Krasulina's method (D-Krasulina) that can keep up with the high streaming rate of data by distributing the computational load across multiple processing nodes. The analysis improves upon the one in (Balsubramani et al., 2013) for the original Krasulina's method and shows that---under appropriate conditions---D-Krasulina converges to the principal eigenvector in an order-wise optimal manner; i.e., after receiving $M$ samples across all nodes, its estimation error can be $O(1/M)$. In order to reduce the network communication overhead, the paper also develops and analyzes a mini-batch extension of D-Krasulina, which is termed DM-Krasulina. The analysis of DM-Krasulina shows that it can also achieve order-optimal estimation error rates under appropriate conditions, even when some samples have to be discarded within the network due to communication latency. Finally, experiments are performed over synthetic and real-world data to validate the convergence behaviors of D-Krasulina and DM-Krasulina in high-rate streaming settings.

Published in Transactions on Machine Learning Research (10/2022)
Distributed Stochastic Algorithms for High-rate Streaming Principal Component Analysis
Haroon Raja
hamon..najo@rutters.ell
Deportment of Electrical and Computer Engineering -New
Rulgers University New Brunswick, Piscalaway, NJ 08854 USA
Wahced U. Bajwa
waheed. bajuan &rutgers edu
Deportment Electrical and Computer Engineering Department of Statistics
Rutgers University-Nor Brunswick, Piscataway, NJ 08854 USA
Reviewed on OpenReview: /ttps://openreeee net/ forum? id= pp2
Abstract
This paper considers the problem of estimating the principal eigenvector of a covariance matrix from independent and identically distributed data samples in streaming settings. The streaming rate of data in many contemporary applications can be high enough that a single processor cannot finish an iteration of existing methods for eigenvector estimation before a new sample arrives. This paper formulates and analyzes a distributed variant of the classical Krasulina's method (D-Krasulina) that can keep up with the high streaming rate of data by distributing the computational load across multiple processing nodes. The analysis improves upon the one in (Balsubramani et al. 2013) for the original Krasulina's method and shows that under appropriate conditions D-Krasulina converges to the prin- cipal cigenvector in an order-wise optimal manner; i.c., after receiving M samples across all nodes, its estimation error can be (1/MM). In order to reduce the network communica- tion overhead, the paper also develops and analyzes a mini-batch extension of D-Karsulina, which is termed DM-Krasuling
The analysis of DM-Krasulina shows that it can also achieve order-optimal estimation error rates under appropriate conditions, even when some samples have to be discarded within the network due to communication latency. Finally, experiments are performed over synthetic and real-world data to validate the convergence behaviors of D-Krasulina and DM-Krasulina in high-rate streaming settings.
1 Introduction
Dimensionality reduction and feature learning methods such as principal component analysis (PCA), sparse PCA, independent component analysis, and autoencoder form an important component of any machine learn- ing pipeline. For data lying in a d-dimensional space, such methods try to find the A < d variables seriabllesfeaturrs features that are most relevant for solving an application-speccfic task (c.g., classification, regression, estimation, data compression, etc.). The focus of this work is on PCA, where the objective to compute k-features that capture most of the variance in data. The proliferation of big data (both in terms of dimensionality and number of samples) has resulted in an increased interest. in developing new algorithms for PCA due to the fact that classical numerical solutions (e.g., power iteration and Lanczos method (Golub & Van Loan 2012)) ) for computing eigenvectors of symmetric matrices do not scale well with high dimensionality and large sample sizes. The main interest in this regard has been on developing algorithms that are cheap in terms of both memory and computational requirements as a function of dimensionality and number of data samples.
In addition to high dimensionality and large number of samples, another defining characteristic of modern data is their streaming nature in many applications; examples of such applications include the internet-of- things, high-frequency trading, meteorology, video surveillance, autonomous vehicles, social media analytics,

---

Published in Transactions on Machine Learning Research (10/2022)
etc. Several stochastic methods have been developed in the literature to solve the PCA problem in streaming settings (Krasulina 1969 Oja & Karhunen 1985 Sanger 1989 Warmuth & Kuzmin 2007 Zhang & Balzano, 2016). These methods operate under the implicit assumption that the data arrival rate is slow enough so that each sample can be processed before the arrival of the next one. But this may not be true for many modern applications involving high-rate streaming data. To overcome this obstacle corresponding to high- rate streaming data, this paper proposes and analyzes distributed and distributed, mini-batch variants of the classical Krasulina's 's method (Krasulina, 1969). Before providing details of the proposed methods and their relationship to prior work, provide brief overview of the streaming PCA problem.
1.1 Principal Component Analysis (PCA) from Streaming Data
For data lying in Rd PCA learns k-dimensional subspace with maximum data variance. Let x E be a random vector that is drawn from some unknown distribution Pa with zero mean and E covariance matrix. For the constraint set V == {V e Rdxk; VTV== optimization problem:
= 1}, we can pose PCA as the following constrained
Q' = arggmax Ep, { (T(VTXXTV))
VEV
where Tr(.) denotes the trace operator. The solution for the statistical risk marimization problem [ is
the matrix Q" with top k eigenvectons of E. In practice, however, ( cannot be solved in its current form since Ps is unknown. But if we have T data samples, {x+]E= drawn independently from Pz, then we can accumulate these data samples to calculate the sample covariance matrix as:
Σ
A
(2)
where A, == X,x," Instead of solving
I we can now solve an empirical risk mazimization problem
Q = arg max (((VTTV) = arg max
Σ T((VTA.V).
(3)
VEV
VEV
In principle, we can solve by computing the singular value decomposition (SVD) of sample covariance Ar. But this is computationally intensive task that requires (dd) multiplications and that has mem- ory overhead of O(f). In contrast, the goal in high-dimensional PCA problems is often to have O(d*k) computational complexity and O(dk) memory complexity (Li et al. 2016).
iteration and Lanczos method (Golub & Van Loan 2012 Chapter 8). Although these methods improve overall computational complexity of PCA to O(d)k), they still have memory requirements on the order of Odd). In addition, these batch methods that require computing the sample covariance matrix Ar, which results in (ddT) multiplication operations. Further, in streaming settings where the goal is real-time decision making from data, it is infeasible to compute Är. Because of these reasons, stochastic approximation methods such as Krasulina's method [Kassulina] 1969) and Oja's rule OOja & Karhunen 1985) are often favored for the PCA problem. Both these are simple and extremely efficient algorithms, achieving (dd) computational and memory complexity per iteration, for computing the principal eigenvector (i.e., k = 1) covariance matrix in streaming settings. Recent years in particular have seen an increased popularity of these algorithms and we will discuss these recent advances in Section
More efficient (and hence popular) approaches for PCA use methods such as the power
Both Oja's rule and Krasulina's method share many similarities. In this paper, we focus on Krasulina's method with the understanding that our findings can be mapped to Oja's rule through some tedious but straightforward calculations. Using t for algorithmic iteration, Krasulina's method estimates the top eigen- vector by processing one data sample in each itcration as follows
VT xxX, VI-IVe-1
Vi Vt-1 + 7o Xixi V
(4)
V
II contrast, the iterate of Oja's rule given by Vs Ve-1 t (xexF - vi-ve1---
2

---

Published in Transactions on Machine Learning Research (10/2022)
where "Y denotes the step size. Going forward, we will be using A, in place of xxx;" in expressions such as (4) for notational compactness. In practice, however, one should neither explicitly store A, nor explicitly use it for calculation purposes.
Note that one can interpret Krasulina's method as a solution to an optimization problem. Using Courant- Fischer Minimax Theorem (Golub & Van Loan 2012 Theorem 8.1.2), the top eigenvector computation (i.e., 1-PCA, which is the version of (T) can be posed as the following optimization problem:
v Avv
q1 == arg min f(v) arg min vER*
(5)
vERd
1v
In addition, the gradient of the function S(v) defined in 5 is:
(v AvVv
Vf(v) 1
A.
(6)
1
Looking at 4 we see that 4 is very similar to applying stochastic gradient descent (SGD) to the nonconvex problem 5 with the only difference being the scaling factor of 1/1lv||3. Nonetheless, since 5 is nonconvex problem and we are interested in global convergence behavior of Krasulina's method, existing tools for analysis of the standard SGD problem (Bottoul 2010 Recht et al. 2011 Dekel et al. 2012 Reddi et al.) (016b[c) do not lend themselves to the fastest convergence rates for Krasulina's 's method.
Despite its nonconvexity, however, (5) has a somewhat benign optimization landscape and a whole host of algorithmic techniques and analytical tools have been developed for such structured nonconvex problems recent years that guarantce fast convergence a global solution. In this paper, we leverage some of these recent developments to guarantee near-optimal global convergence of two variants of Krasulina's method in the case of high-rate streaming data.
Before proceeding further, it is worth noting that while Krasulina's method primarily focuses on the 1-PCA problem, it can be used to solve the k-PCA problem. But such an indirect approach, which involves repeated use of the Krasulina's method & times, can be inefficient in terms of sample complexity (Allen AhllnnZZuu & Li, 2017a, Section 1). We leave investigation of a ncar-optimal direct method for the k-PCA problem involving high-rate streaming data for future work.
1.2 Our Contributions
In this paper, we propose and analyze two distributed variants of Krasulina's 's method for estimating the top eigenvector of a covariance matrix from fast streaming, independent and identically distributed (i..d.) data samples. Our theoretical analysis, well numerical experiments ncar-optimality of the proposed algorithms. In particular, our analysis shows that the proposed algorithms can achieve the optimal convergence rate of O(1/M) for 1-PCA after processing a total of O(M) data samples (see (Jain et al.) 2016. Theorem 1.1) and (Allen- Zhu & Li 2017a Theorem 6)). In terms of details, following are our key contributions:
synthetic and real data, establish
1. Our first contribution corresponds to the scenario in which there is mismatch of N E 14 > between the data streaming rate and the processing capability of a single processor, i.e., one iteration of Krasulina's 's method on one processor takes as long as N data arrival epochs. Our solution to this problem, which avoids discarding samples, involves splitting the data stream into N parallel streams that are then input to N interconnected processors. Note that this splitting effectively reduces the streaming rate at each processor by a factor of N. We then propose and analyze a distributed variant of Krasulina's 's method-tenee
D D-Krrsulina-tte that solves the 1-PCA problem for this distributed setup consisting of N processing nodes. Our analysis substantially improves the one in (Balsubramani et al., 2013) for Krasulina's method and shows that D-Krasulina can result in an improved convergence rate of (1/NN)) after iterations (Theorem [] as opposed to the ((1/l) rate for the classical Krasulina's method at any one of the nodes seen in isolation. Establishing this result involves novel analysis of Krasulina's method that brings out the dependence of its convergence rate on the variance of the sample covariance matrix; this analysis coupled with a variance reduction argument leads to the convergence rate of (1/NNt) for D-Krasulina under appropriate conditions.
3

---

Published in Transactions on Machine Learning Research (10/2022)
2. Mini-batching of data samples has long been promoted strategy in stochastic methods to reduce the wall-clock time. Too large of a mini-batch, however, can have an adverse effect on the algorithmic performance; see, e.g. (Shamir & Srebro, 2014, Sec. VIII). One of the challenges in mini-batched stochastic methods, therefore, is characterizing the baniibbach size that lcads to near- optimal convergence rales terms of the number of processed samples. In (Agarwal & Duchi
2011
Cotter et al., 2011 Dekel et al.] 2012 Shamir & Srebro 2014 Ruder 2016 Golmant et al.
2018
Goyal et al. 2017], for example, the authors have focused on this challenge for the case of
mini-batch SGD for convex and nonconvex problems. In the case of nonconvex problems, however, the guarantees only hold for convergence to first-order stationary points. In contrast, our second contribution is providing comprehensive understanding of the global convergence behavior of mini- batched Krasulina's 's method. In fact, our analysis of D-Krasulina is equivalent to that of a mini- batch (centralized) Krasulina's ' method that uses a mini-batch of N samples in each iteration. This analysis, therefore, guarantecs near-optimal convergence rate with arbitrarily high probability for an appropriately mini-batchedee Krasulina's method in a centralized setting. This is in contrast to (Jain et al.) 2016 Yang et al.] 2018), where the focus is on Oja's rule and the probability of success is upper bounded by 3/4 for a single algorithmic run F In addition, in the case of high- rate streaming data that requires splitting the data stream into N parallel ones, we characterize the DM-Krasulina in terms of the batch size, This involves specifying the conditions under which mini- batches of size B/N per node can lead to near-optimal convergence rate of Q(1BBt) after t iterations of DM-Krasulina (Theorem 2). An implication of this analysis is that for a fixed (network-wide sample budget of T samples, DM-Krasulina can achieve O(1/T) rate after t == T/B iterations provided the (network-wide) mini-batch size B satisfies B = O(T ) for some constant Co > 2 (Corollary
global convergence behavior of a mini-batch generalization of --Krasulina-temma
3. Our next contribution is an extended analysis of DM-Krasulina that concerns the scenario where (computational and/or communication) resource constraints translate into individual nodes still receiving more data samples than they can process in one iteration of DM-Krasulina. This resource- constrained setting necessitates DM-Krasulina dropping μ E Z, samples across the network in each iteration. Our analysis in this setting shows that such loss of samples need not result in sub-optimal performance. In particular, DM-Krasulina can still achieve near-optimal convergence rate a function of the number of samples arriving in the network for both infinite-sample and finite-sample regimes-as as long as μ = O(B) (Corollary 2
4. We provide numerical results involving both synthetic and real-world data to establish the usefulness of the proposed algorithms, validate our theoretical analysis, and understand the impact of the number of dropped samples per iteration of DM-Krasulina on the convergence rate. These results in particular corroborate our findings that increasing the mini-batch size improves the performance of DM-Krasulina up to a certain point, after which the convergence rale starts to decrease. Since the focus of this work is on systems theory, the reported results do not focus on some of the large- scale system implementation issues such unexpected processor failures, high network coordination costs, etc. Such large-scale implementation issues, while relevant from practical perspective, are beyond the scope of this paper and provide interesting research directions for future work.
1.3 Related Work
Solving the PCA problem elficiently in a number of settings has been an active area of research for decades. [Krasulina]
1969
Oja & Karhunen 1985) are among the earliest and most popular methods to solve PCA
in streaming data settings. Several variants of these methods have been proposed over the years, includ-
ing (Bin Yang 1995 Chatterjee
2005 Doukopoulos & Moustakides 2008). Like earlier developments in
"Note that it is possible to improve the probability success for Oja's rule, as derived in (Jain et al. 2016 Yang et al. 2018 to 1, d E (0, 1), by running (llg((1/)) instances of the algorithm and combining their outcomes. But such strategy, which also adds to the computational and storage overhead because of the 'combine' step, impractical for the streaming settings considered in this paper. Additional differences between the results of this paper pertaining the centralized PCA problem for streaming data and that of (Jain et al. 2016 Yang al. 2018) discussed in Section 1.3 It is, however, important to point out that this work and (Jain et al. 2016 Yang al. (2018) are analyxing two related, yet different, algorithmic approaches that are based Krasulina's as method and Oja's 8 rule, respectively, and their results are therefore complementary in nature.
4

---

Published in Transactions on Machine Learning Research (10/2022)
stochastic approximation methods (Robbins & Monro 1951 such variants were typically shown to converge asymptotically. Convergence rate analyses for stochastic optimization in finite-sample settings (Shapiro & Homem-de Mello 2000 Linderoth et al. 2006) paved the way for non-asymptotic convergence analysis of different variants of the stochastic PCA problem, which is fundamentally a nonconvex optimization problem. Within the context of such works, the results that are the most relevant to the algorithmic strategy devised in this paper can be found in (Jain et al. 2016) and (Yang et al.] 2018). The authors in these two papers provide variance-dependent convergence guarantees for Oja's rule in the finite-sample regime, thereby mak- ing their results also translatable to the algorithmic framework being considered in this paper for high-rate streaming PCA. However, as noted earlier, the results derived in (Jain et al. 2016: Yang et al. 2018) only hold with probability 3/4, which is in contrast to the high-probability results of this paper. And while one could increase the probability of success in (Jain el al. 2016 Yang et al. 2018 through multiple algorithmic runs, this is not a feasible strategy in the streaming settings.
In order to complement the pioncoring results of (Jain et al. 2016 Yang ct al. 2018 in streaming settings, we shift our focus away from Oja's rule as the base algorithm and develop a different proof strategy that substantially extends and generalizes the analysis in (Balsubramani et al.] 2013) for Krasulina's method. The analysis in (Balsubramani et al.] 2013 assumes that the Eg norm of the data samples is bounded by a positive constant, but it does not take into account the variance of the sample covariance matrices. Such an analysis leads to convergence results that are independent of the variance and hence are unable to capture any improvement in the convergence rate due to mini-batching and/or distributed implementations for computing the top cigenvector of covariance matrix. We overcome this limitation of the analysis in [Balsubramani et al. 2013) by developing a proof in this work that explicitly accounts for the variance of the sample covariance matrices. Note that this extension/ extension/geenerlizatiin generalization of the analysis in (Balsubramani et al., 2013 despile the intuitive nature of our final set of results is nontrivial task. There are in particular two main technical challenges that are addressed in this paper: (i) Using concentration of measure results that allow for incorporation of the variance within the analysis, as opposed to the Hoeffding inequality (Boucheron et al. 2013) utilized in (Balsubramani et al. 2013); and (ii Utilizing the variance-dependdnt concentration guarantees within two terms in Lemma 1 one of which depends on the norm bound on data samples and the other on the variance careful decoupling of these two quantities being critical to obtain convergence results in which the dominant term depends only on the variance. Finally, another aspect that distinguishes
this paper from prior works such as (Balsubranani et al. 2013 Jain et al., 2016 Yang et al.) 2018) is that it provides a formal framework for studying the communications and computation tradeoffs involved in solving the 1-PCA problem in distributed streaming settings. This framework is described in detail in Section 3.2 with theoretical characterization of the proposed framework in terms of the communications and computation costs described in Section 4.2
Because of the vastness of literature on (stochastic) PCA, this work is also tangentially or directly related to number of additional prior works. We review some of these works in the following under the umbrellas of different problem setups, with the understanding that the resulting lists of works are necessarily incom- plete. Much of our discussion in the following focuses on solving the PCA problem in (fast) streaming and distributed data settings, which is the main theme in this paper.
Sketching for PCA. Sketching methods have long been studied in the literature for solving problems involving matrix computations; see (Woodruff) 2014) for a review of such methods. The main idea behind these methods is to compress daia using either randomized or deterministic sketches and then perform
computations the resulting low-dimensional data. While sketching has been used a tool to solve the
PCA problem in an efficient manner (sce, c.g., (Warmuth & Kuzmin 2007 Halko et al. 2011 Liberty 2013:
Leng et al. 2015 Karnin & Liberty 2015)), the resulting methods cannot be used to ezactly solve in the fast streaming settings of this paper.
Online PCA. The PCA problem has also been extensively studied in online settings. While such settings also involve streaming data, the main goal in online PCA is to minimize the cumulative subspace estimation error over the entire time horizon of the algorithm. The online PCA framework, therefore, is especially useful
in situations where either the underlying subspace changes over time or there is some adversarial noise in
the sampling process. Some of the recent works in this direction include (Garber et al. 2015) Allen- Zhu &
2017b
Garber 2018
Marinov et al. 2018 Kotlowski & Neu 2019).
5

---

Published in Transactions on Machine Learning Research (10/2022)
Stochastic convex optimization for PCA. One approach towards solving (3) in streaming settings is to relax the PCA problem to convex optimization problem and then use SGD to solve the resulting stochastic convex optimization problem (Arora et al. 2013 Garber & Hazan 2015 Nie et al.) 2016). The benefit of this approach is that now one can rely on rich literature for solving stochastic convex problems using SGD. But, the tradeoff is that one now needs to store an iterate of dimension Rdxdd as opposed to an ilerate of dimension Rd when we solve the PCA problem in its original nonconvex form. Due to these high memory requirements of Odd)), we limit ourselves to solving PCA in the nonconvex form.
Streaming PCA and nonconvex optimization. The PCA problem in the presence of streaming data can also be tackled as an explicit constrained nonconvex optimization program (Zhang & Balzano 2016 De Sa et al. 2015). In (Zhang & Balzano, 2016), for instance, the problem is solved as an optimization program over the Grassmannian manifold. The resulting analysis, however, relies on the availability of a good initial guess. In contrast, the authors in (De Sa et al. 2015) analyze the use of the SGD for solving certain nonconvex problems that include PCA. The resulting approach, however, requires the step size to be significantly small constant for eventual convergence (e.g.; 10-12 12 for the Netflix Prize dataset); this translates into slower convergence in practice.
Classical stochastic approximation methods for PCA. Recent years have seen an increased interest in understanding the global convergence behavior of classical stochastic approximation met such as Krasulina's ' method (Krasulina 1969) and Oja's rule {Oja & Karhunen 1985) for the PCA problem in non-
asymptotic settings (Allen-Zhu & Li 2017a; Chatterjee [005] HHrddt & Price, 2014; Shamir 2015; 2016
Jain et al.) 2016) Li et al. 2016 Tang 2019 Henriksen & Ward 2019) Amidd & Warmuth 2019). Some of these works, such as (Shamir 2015 and (Shamir) 2016), use variance reduction techniques to speed-up the algorithmic convergence. Such works, however, require multiple passes over the data, which makes them ill-suited for fast streaming settings. The analysis in (Shamir) 2015) and (Shamir 2016) also requires an initialization close to the true subspace, which is somewhat unlikely in practice. Among other works, the authors in (Allen-Zhu & Li 2017a) provide cigengap-free convergence guarantecs for Oja's rule. Since the results in this work do not lake into account the variance of data samples, they do not generalize to mini-batche b batch/distributed distributed streaming settings. As stated earlier, the authors in (Jain et al.] 2016) do provide variance-dependent convergence guarantees for Oja's rule, which makes this work the most relevant to ours. In particular, the authors in YYang et al. 2018 have extended the initial analysis in (Jain et al., 2016
to mini-batch settings. But limitation of the analysis in (Jain et al., 2016 Yang et al. 2018) is that it guarantees convergence of Oja's rule only up to a probability of 3/4. And while the probability of success can be increased to 1-8 8 by running and combining the outcomes of O(log(1//)))) (1)5) instances of Oja's rule, such an approach has two major drawbacks in streaming settings. First, since new data samples arrive continuously in a streaming setting, multiple runs of an algorithm in this case can only be achieved through multiple replicas of the processing system. Such strategy, therefore, leads to substantial increase in system costs. Second, the outcomes of the multiple runs need to be appropriately combined. In (Jain et al., 2016), it is suggested this be done by computing the geometric median of the multiple outcomes, which requires solving an additional optimization problem. This then adds to the computational and storage overhead for the PCA problem. We conclude by remarking on two key distinctions bet ween our results and those in (Jain) et al. 2016 Yang et al., 2018). First, the arbitrarily high probability of success in our analysis requires the initial step size in Krasulina's method to decrease with an incrcase in the probability; we refer the reader to the discussion following Theorem [ in this paper for further details on this point. Second, our convergence guarantees have the flavor of 'convergence in mean' as opposed to the 'convergence in probability' nature of the results in (Jain et al.) |2016] Yang et al. 2018). A straightforward application of Markov's inequality, however, leads to results that are directly comparable to the ones in (Jain et al. 2016 Yang et al. 2018 we refer the reader to Corollary 3 in Appendix D as an illustrative example of this.
2016 Boutsidis et al. 2016; Garber et al.] [017 De Sa et al.] 2018) have focused on the PCA problem in distributed settings. Among these works, the main focus in (Balcan et al., 2016 Boutsidis et al. 2016 Garber et al.,
Distributed PCA and streaming data. Several recent works such as (Balcan et al.
2017 is on improving the communications efficiency. This is accomplished in (Balcan et al. 2016 Boutsidis
et al. 2016) by sketching the local iterates and communicating the resulting compressed iterates to central
server in each iteration. In contrast, (Garber et al. 2017 provides a batch solution in which every node
6

---

Published in Transactions on Machine Learning Research (10/2022)
in the network first computes the top eigenvector of its local (batch) covariance matrix and then, a last step of the algorithm, all the local eigenvector estimates are summed up at central server to provide an eigenvector estimate for the global covariance matrix. In contrast to these works, our focus in this paper is on establishing that distributed (mini-batch) variants of stochastic approximation methods such as Oja's rule and Krasulina's 's method can lead to improved convergence rates, as a function of the number of samples, for the PCA problem in fast streaming settings. In this regard, our work is more closely related to (De Sa iet al., 2018), where the authors use the momentum method to accelerate convergence of power method and further extend their work to stochastic settings. However, the approach of (De Sa et al.) 2018) relies a variance reduction technique that requires a pass over the complete dataset every once in a while; this is impractical in streaming settings. In addition, theoretical guarantees in (De Sa et al. 2018 are based on the assumption of a "good" initialization; further, an implicit assumption in (De Sa el al. 2018 is that inter-node communications is fast enough that there are no communication delays.
Decentralized PCA. Another important practical sctting under which the PCA problem has been studied is when the data are distributed across an interconnected set. of nodes that. form complete graph but not fully connected one. We refer the reader to (Blondel et al. 2005 Aysal et al.
2009
Khan et al. 2009
Dimakis et al. 2010) for a general understanding of this decentralized sctting. Some contributions to the
PCA problem in this setting are (Kempe et al. 2008 Li et al. 2011 Korada et al.) 2011 Wu et al.
2017 2018
Gang et al. 2021 Gang & Bajwa 2021 2022) Most of these contributions consider batch data settings and their extensions to the streaming setting are not evident. And while contributions such as (Li et al. 2011
do consider the streaming data setting, the convergence rates established in such works are asymptotic nature and the theoretical analyses cannot be used to derive conditions for linear speed-up in convergence rates in the mini-batch sctting.
Connections to low-rank matrix approximation. The low-rank metriz approzimation problem (Frieze et al. 2004 Clarkson & Woodruff 2017), which involves computing a low-rank approximation of given matrix, is closcly related to the PCA problem. The overarching goal in both problems is the same: find a subspace that best upprozinates the data samples / given malria. In the case of PCA, however, the focus is fundamentally on finding an orthogonal basis for the subspace that is precisely given by the top eigenvectors of the data covariance matrix. Notwithstanding this difference between the two problems, (Yun et al., 2015; Tropp et al. 2019) are two works within the low-rank matrix approximation literature that are the most related to this paper. The setting in (Yun et al., 2015) corresponds to a large-scale but fixed matrix whose columns are presented to a low-rank approximation algorithm in a streaming manner. This is in contrast to the streaming setting of this work that is akin to having a matrix with infinitely many columns. In addition, the algorithm being studied in Yun et al. 2015) requires computing the top principal component directions of a random submatrix of the larger matrix in a batch setting as its first step, which is again departure from the streaming data setting of this work. Next, the mathematical model in (Tropp et al. 2019) corresponds to matrix that is given by the sum of a long sequence of low-rank and sparse matrices that are presented to a low-rank approximation algorithm in a streaming manner. This summation aspect of the data model in (Tropp et al. 2019) does not coincide with the mathematical model of the streaming data samples in this work. Finally, and in stark costr to this paper, neither of these works are concerned with the interplay between the streaming data rate, data processing rate, and the number of interconnected processors.
Connections to stochastic nonconvex optimization. Recent years have also seen an increased focus understanding (variants of) SGD for general (lypically unconstrained) stochastic nonconvex optimization problems. Among such works, some have focused on classical SGD (Ge et al. 2015 Hazan et al. 2016
2017 Li ct
al. 2021 Mokhtari et al. 2017), some have studied varianco-reduction variants of SGD Reddi
et al. 2016a b. Qureshi et al. 2021). and some have investigated accelerated variants of stochastic nonconvex
optimization [AllenZhhu] 2018bla) a). In particular, works such as (Reddi et al., 2016a Allen-Zhu & Hazan 2016) are directly relevant to this paper since these works also use mini-batches to reduce sample variance and improve on SGD convergence rates. While (implicit, through the distributed framework, and explicit) mini-batching is one of the key ingredients of our work also, this paper differs from such related works because of its ability to prove convergence to global optimum of the 1-PCA problem. In contrast, aforementioned works only provide guarantees for convergence to first-order stationary points of (typically unconstrained) stochastic nonconvex optimization problems.

---

Published in Transactions on Machine Learning Research (10/2022)
1.4 Notational Convention and Paper Organization
We use lower-case (a), bold-faced lower-case (a), and bold-faced upper-case (A) letters represent scalars, vectors, and matrices, respectively. Given a scalar and vector a, [a] denotes the smallest integer greater than or equal to a, while ||a]|2 denotes the f2-norm of a. Given a matrix A, || denotes its spectral norm and || denotes its Frobenius norm. In addition, assuming A E Rdxd to be a positive semi-definite matrix, Ax(A) denotes its i-th largest eigenvalue, i.e,, || == (A) 2 (2 2 ... 2 >a( 2 0. Whenever obvious from the context, we drop A from (A) for notational compactness. Finally, E{-} denotes the expectation operator, where the underlying probability space (Ω, F, P) is either implicit from the context is explicilly pointed out in the body.
The rest of this paper is organized as follows. We first provide a formal description of the problem and the system model in Section 2 The two proposed variants of Krasulina's method that can be used to solve the 1- PCA problem in fast streaming settings are then presented in Section 3 In Section ! we provide theoretical guarantces for the proposed algorithms, while proofs / outlines of the proofs of the main theoretical results are provided in Section 5
Finally, numerical results using both synthetic and real-world data are presented in Section B while appendices used for detailed proofs of some of the theoretical results.
2 Problem Formulation and System Model
Our goal is to use some variants of Krasulina's method (cf. 44) in order to obtain an estimate of the top eigenvector of covariance matrix from independent and identically distributed (i.i.d)) data samples that are fast streaming into system. The algorithms proposed in this regard and their convergence analysis rely the following sets of assumptions concerning the data and the system.
2.1
Data Model
We consider a streaming data setting where a new data sample x Id independently drawn from an unknown distribution Pa arrives at system each sampling time instance t'. We assume a uniform data arrival rate of R2 samples per second and, without loss of generality, take the data arrival index t' 2 1 to be an integer. We also make the following assumptions concerning our data, which aid in our convergence analysis.
[A1] (Zero-mean, norm-bounded samples) Without loss of generality, the data samples have zero mean, iee, Ep, {x²}= = In addition, the data samples are almost surely bounded in norm, i.e., ||xella S where let the bound 2 without loss of generality.
[A2] (Spectral gop of the covariance matrix) The largest cigenvalue of Σ = Ep, (xxxx) is strictly greater than the second largest eigenvalue, i.e., A1 (E) > ( > A
Σ 2 ... 2 A Σ 2 0.
Note that both Assumptions A1 and A2 are standard in the literature for convergence analysis of Kra- sulina's method and Oja's rule (ef. (Balsubramani et al. 2013 Oja & Karhunen, 1985 Allen-Zhu & Li 2017a: Jain et al. ).16)).
We also associate with each data sample Xy rank-one random matrix Ay == xxx,, which is trivial unbiased estimate of the population covariance matrix Σ. We then define the variance of this unbiased estimate as follows.
Definition (Variance of sample covariance matrix). We define the variance of the sample covoriance matria Ar == x as follows:
0² = Ep,
- =1
Note that all moments of the probability distribution P, exist by virtue of the norm boundedness of x, (cf. Assumption A1]]. The variance ²² of the sample covariance matrix Ay as defined above, therefore, exists and is finite.
8

---

Published in Transactions on Machine Learning Research (10/2022)
X1; -
Splitter
NN, x
X2, X2+NN
X2: X N
NN; X2N,
X1, N
X1, XN+1
(a)
(b)
The distributed PCA problem, which involves distributod processing of data over a network of N processors, can arise in two conlexts. (a) A data splitter can split data stream into N parallel streams, for each processor in the network. In relation to the original data stream, this effectively reduces the data arrival rate for each parallel stream by a factor of N. (b) Data can be inherently distributed, as in the Internet-of- Things systems, and can arrive at N different processing nodes as N separate data streams.
Figure 1:
The two algorithms proposed in this paper, namely, D-Krasulina and DM-Krasulina, are initialized with a random vector Vo E I2d that is randomly generated over the unit sphere R with respect to the uniform (Haar) measure. All analysis in this paper is with respect to the natural probability space FD,F,P) given by the stochastic process (Vo, X1, x2,...) and filtered versions of this probability space.
2.2 System Model
Let Rp denote the number of data samples that single processing node in the system process in second using ileration of the form 44 The focus of this paper is on the high-rate streaming setting, which corresponds to the setup in which the data arrival rate R, is strictly greater than the data processing rate Rg A naive approach to deal with this computation streaming mismatch is to discard (per second) a fraction a := R Rp of samples in the system. Such an approach, however, leads to an equivalent reduction the convergence rate by a. We pursue an alternative to this approach in the paper that involves the simultaneous use of N [a]] interconnected processors, each individually capable of processing Rg samples per second, within the system. In particular, we advocate the of such network of N processors in the following two manners to achieve near-optimal convergence rates (as function of the number of samples arriving at the system) for estimates of the top cigenvector of Σ in high-rate streaming settings.
2.2.1 Distributed Processing Over a Network of Processors
We assume the fast data stream terminates into a data splitter, which splits the original stream with data rate R, samples per second into N parallel streams, each with data rate R,/N samples per second, that are then synchronously input to the interconnected network of N processors; see Figure a) for schematic, rendering of such splitting. In order to simplify notation in this setting, we reindex the data samples associated with the i-th processor / dala ream the following as (x,,)ee)) where the reindexing map (i, () + is simply defined as t' = i N...
Note that the parameter Rp, among other factors, function of the problem dimension this dependence being suppressed the notation for case exposition.

---

Published in Transactions on Machine Learning Research (10/2022)
We also assume the network of processors implements some message passing protocol that allows it to com- N pule sums of locally stored vectors, i.e., Ei=1 a; for the set of local vectors {a.) within the network. This could, for instance, be accomplished using either Reduce or AllReduce primitives within most message passing implementations. We let (N) denote the number of these primitive (sum) operations that the message passing protocol can carry out per second in the network of N processors. Note that, in addition to the number of nodes N, this parameter also depends upon the problem dimension d, message passing implementation, network topology, and inter-node communications bandwidth, all of which are being ab- stracted here through (NV). We will also be suppressing the dependence of the parameter Re on N within the notation beyond this section for ease of exposition.
Data splitting among this network of N processors effectively slows down the data streaming rate at each processing node by a factor of N. It is under this system model that we present a distributed variant of Krasulina's ' method, termed D-Krasulina, in Section B.T that involves two key operations after each round of splitting of the N samples: (i) per-processor computation of the form (4) which requires 1/ Rp seconds for completion, and (ii) nelwork-wide vector-sum operation, which requires an additional 1/R seconds for completion. Under such algorithmic framework, can express the effective network-wide data processing rate, denoted by Re, in terms of the per-node data processing rate Rg and the network-wide sum rate Re (N) by noticing that the effective time it takes to process one sample within the network is given by
1 seconds/sample.
Te NR +
RN(N)
Here, the first algorithmic operation gives rise to the first term in the above expression, since it takes 1/R, seconds for computation of the form 4 for N samples, and the wide vector-sun operation gives rise to the second term. It then follows that the effective data processing rale is R. = 1/T, = (NN)
(N)
In the high-rate streaming setting, D-Krasulina can only operate as long as the streaming rate R, does not exceed the effective processing rate Re, ie., R. I Rg. This gives rise to the following condition on the number of processors N that enables our algorithmic framework to cope with the high-rate streaming data:
N (Re(N)- - R.) ) 2 (N)
N (gR Rp (R.(N) R.))
(Ro(N) - R.) > 0.
(7)
Rp
Our developments in the following will be based on the assumption that the condition on N in 7 is met. The main analytical challenge for us is understanding the scenarios under which distributed processing over network of N processors using D- Krassulina still yields near-optimal performance; we address this challenge
in Section 4.1 by improving on the analysis in (Balsubramani et al., 2013) for Krasulina's 'ss method.
We conclude by remarking on the nature of the lower bound in 7 on the number of processors N, given that the bound itself is a function of N through its dependence on the parameter (NV that is expected to decrease with an increase in N. The pergh-pefforaccc computing community has access to copious amounts of trace data for different parallel computing environments that allows one to relate the parameter Re (N) for given dimension d to the mumber of processors N; see, e.g., (Kavulya et al., 2010). One such relationship could be, for instance, that (N) x I/N" for some E (0, 1]. The final bound on the number of processors N can then be obtained by plugging such a relationship into the provided bound
Remark 1. It is struightforword to see that our developments in this paper are also applicuble to the selling in which dato naturally arrives in a distributed manner at N different nodes, as in Figure (b). In addition, our analysis of D-Krasulina. is equivalent to that of mini--atch Krasulina's 's method running on powerful- enough single processor that uses a mini- mini-batch of N samples in each iteration.
As illustrative toy example, let consider scenario for fixed d in which R. 10 sec-1, Rp 50 sec-1, and the parameter Re(N) scales as Re(N)= = Cs sec for a constant CR that depends on the message passing implementation, network topology, and inter-node communications bandwidth. In "slow" network, corresponding to Cn < 10 there exists N that satisfies the condition 7 In "faster" network, corresponding to CR 1 necessary condition on the number of processors is N 10-3CR. But while this necessary coudition is satisfied for 1 N < 80 for CR high as 8 10 there still does exist any N that satisfies ( In the case of even faster networks, however, we are able find feasible values of N for running D-Krasulina without discarding any samples; e.g., any 28 < N < 72 satisfies the condition for CR 10a.
N
10

---

Published in Transactions on Machine Learning Research (10/2022)
2.2.2 Distributed Processing Coupled with Mini-batching
Mini-batching in (centralized) stochastic methods, as discussed in Section 1.2 helps reduce the wall-clock time by reducing the number of read operations per iteration. Mini-batching of samples in distributed settings has the added advantage of reduction in the average number of primitive (sum) operations per processed sample, which further reduces the wall-clock time. It is in this vein that we put forth mini-batched variant of D- D-Krasulina, which is termed DM-Krasulina, in Section 3.2
Similar to the case of D-Krasulina (cf. Figure D there are several equivalent system models that can benefit from the DM-Krasulina framework. In keeping with our theme of fast streaming data, as well as for the sake of concreteness, we assume the system buffers (i.e., mini-batcchss) B == bN 2 [R,/Ra, samples of the incoming data stream every B/R, R., seconds for some parameter b E Z4. This network- widd mini-batch of B samples is then split into N parallel (local) mini-batches, each comprising b B/N samples, which are then synchronously imput to the interconnected network of N processors at rate of R,/N samples per second and collaboratively processed by DM-Krasulina. In each iteration i of DM-Krasulina, therefore, the network processes total of B 2 N samples, as opposed to N samples for D-Krasulina. In order to simplify notation in this mini-batched distributed setting, we reindex the data samples in the mini-batch associated with the i-th processor in iteration of AM-Krasulina as {x.j,))f where the reindexing map (..j.t) + is defined as t' = j + (i - 1)6 (t - B)B.
Tt is straightforward to see from the discussion surrounding (7) that the DM-Krasulina framework can process all data samples arriving at the system as long as N Z BR.RR R,[bH.-K,]] and (bR. - Ra) > 0. However, when this condition is violated due to faster streaming rate Rg, slower processing rate Rp, slower summation rate Re, or any combination thereof, it becomes necessary for DM-Krasulina to discard μ == 6R, samples at the splitter per iteration. The main analytical challenges for DM-Krasulina are, therefore, twofold: first, assuming lr 0, characterize the mini-batch size B that leads to near-optimal convergence rates for DM-Krasulina in terms of the total number of samples arriving at the system; second, when discarding of samples becomes necessary, characterize the interplay between B and that allows DM-Krasulina to still achieve (order-wise) near-optimal convergence rates. We address both these challenges in Section 4.2
+ N.. B R.
3 Proposed Distributed Stochastic Algorithms
We now formally describe the two stochastic algorithms, termed D-Krasulina and DM-Krasulina, that can be used to solve the 1-PCA problem from high-rate streaming data under the two setups described in Section 2.2.1 and Soction 2.2.2 respectively.
3.1 Distributed Krasulina's Method (D-Krasulina) for High-rate Streaming Data
Recall from the discussion Section 2.2. that each node in the network receives data sample X in iteration of the distributed implementation, which comprises N processing nodes. Unlike the centralized Krasulina's method (cf. 4), therefore, any distributed variant of Krasulina's method needs to process N samples in every iteration t. Using Ai, as shorthand for xi,x one natural extension of (4] that processes N samples in each iteration is as follows:
Σ A
1
A vc-1 + 7
Vi Vt-1 + 7t
|
One natural question here is whether S can be computed within our distributed framework. The an- R.R. swer to this is the affirmative under the assumption N 2 R,(R.-R,)) R,): with the implementation (termed D-Krasulina) formally described in Algorithm
Notice that unlike classical Krasulina's method, which processes a total of t samples after t iterations, D-Krasulina processes a total of Nł samples after iterations in order to provide an estimate v, of the
Here, and the following, the implicit assumption that the quantity R. R. (resp., dR. R.) is strictly positive for D-arssulnna (resp., DM- Krasulina).
11

---

Published in Transactions on Machine Learning Research (10/2022)
Algorithm Distributed Krasulina's Method (D-Krasulina)
Input: Incoming data streams at N processors, expressed as xx, l.d. Pz (=1,122,
and step-size sequence
{ne E ER+1EE, Initialize: All processors initialize with vo E I&d randomly generated over the unit sphere
1: for = 2. (In Parallel) Processor receives data sample x,, and updates Eit locally as follows:
Vi (I,...,NN, Eit x V-11 Vi- V Vi-
V 12
N E N ii=1 Ei, in the network using a distributed vector-sum subroutine 4: Update eigenvector estimate in the network as follows: V. f- Ve-1 + 71ht 5: end for Return: An estimate of the eigenvector q of associated with (Z)
Compute & +
top eigenvector q" of Σ. Another natural question, therefore, whether the estimate Vr returned by D-Krasulina can converge to q" at the near-optimal rate of O (1/# of processed samples). Convergence analysis of D-Krasulina in Section 4 establishes that the answer to this is also in the affirmative under appropriate conditions that are specified in Theorem An important interpretation of this result is that our proposed distributed implementation of Krasulina's method can lead to linear speed-up as a function of the number of processing nodes N in the network.
3.2 Mini-batched D-Krasulina (DM-Krasulina) for High-rate Streaming Data
The distributed, mitcheetehed setup described in Section 2.2.2 entails each node receiving mini-batch of b = B/N data samples, {x.j,)}=ı in each iteration t, for total of B bN samples across the network in every iteration. Similar to (8). these B samples can in principle be processed by the following variant of the original Krasulina's 's itcration:
Ve Vi-1 +
B i=1
Aij.Vi-1-
||
1
11
Vt-1
T
B
A A.j. v-
(9)
where Aig.t is a shorthand for xI Practical computation of within our distributed framework, however, requires consideration of two different scenarios.
bR.Rt R,(OK.-K.) This enables incor- poration of every sample arriving at the system into the eigenvector estimate.
Scenario 1:
The mini-batched distributed framework satisfies 2
Scenario 2:
6R. This The mini-batched distributed framework leads to the condition N < (bRl
necessitates discarding of = 6R. NR. B samples per iteration in the system. Stated dif- ferently, the system reccives B + /1 samples per iteration in this scenario, but only B samples per iteration are incorporated into the eigenvector estimate.
R2
We now formally describe the algorithm (termed (MM-Kraullna) that implements under both these scenarios Algorithm
Speaking strictly in terms of implementation, the mini-batched setup of DM-Krasulina allows one relax RR
the condition N 2 Rp(Re-R. R.) associated with D-Krasulina to either N 2 (bR R.j which still incorporates
R.R.
all samples into the eigenvector estimate, N (bl Rt,, which involves discarding of 0 samples
6R.
12

---

Published in Transactions on Machine Learning Research (10/2022)
Algorithm 2 Distributed Mini-batch Krasulina's Method (DM-Krasulina)
Input: Incoming streams of mini-batches { {x.,je P=) ) at N processors, size of the network-wide mini batch == bN, and a step-size sequence {7e E TR+ Initialize: All processors initialize with Vo E IRd randomly generated over the unit sphere 1: for = 2.2.... do
1.d. Ps
(In Parallel) Vi E N). E..
for I....,b do
(In Parallel) Processor i receives data sample Xijit and updates E.t locally follows:
E NJ, Ei, + Eit + Xi.a
V
5:
end for
6:
Compute & + B Ei= Ese in the network using a distributed vector-sum subroutine
Update eigenvector estimate the network as follows: Vt f- Vt-1 + 7
if N bR.R, then [,,bbl,-R,))
9:
The system (e.g., data /buute)r receives (B + 4) additional samples during execution of Steps
7 out of which μ E Z4 samples are discarded
10:
end if
11:
end for
Return: An estimate v, of the eigenvector q" of E associated with A1 (Z)
per algorithmic iteration. While this makes DM-Krasulina particularly attractive for systems with slower communication links, the major analytical hurdle here is understanding the interplay between the different problem parameters that still allows OM-Krasulina to achieve near-optimal convergence rates in terms of the number of samples received at the system. We tease out this interplay as part of the convergence analysis of DM-Krasulina in Section
3.3 A Note on the Processing of Non-centered and Non-i.i.d. Data
Both D-Krasulina and DM-Krasuling have been developed under the assumptions of zero-mean (i.e., centered) and i..dd data samples. In this section, we discuss one possible approach to handling non-centered data using the two algorithms and also provide a rationale for the applicability of D-Krasulina and DM-Krasulina in the face of any shifts the data distribution.
In the case of non-centered data, one simple strategy that works for D-Krasulina and MM-Krraulina is maintain (network-wide) rumning average of the non-centered data samples, and then use it center the data samples at each processor before applying Step 2 (resp., Step 4) in Algorithm [ (resp., Algorithm 2 While such modification requires an extension of the convergence analysis presented in the next section, this can be accomplished in a manner similar to the analytical extension in (Zhou & Bai 2021 for the centralized Oja's rule with non-centered data.
Next, while the forthcoming convergence analysis for D-Krasulina and DM-Krasulina has been provided under the assumption of i.i.d. data samples, the two algorithms are expected remain effective in the non-i.id. data setting. This is because D-Krasulina and DM-Krasulina first essentially compute new gradient-like quantity using the latest batch of data samples at each time t (cf. Step 2 in Algorithm [ and Step 4 in Algorithm 2), and then update their respective cigenvector estimates using this quantity. In particular, any shifts in the data distribution can be tracked by the two algorithms because of such an update rule. It because of this reason that algorithms such as Oja's rule and Krasulina's method are also often employed for the problem of subspace tracking (soe, o.g., (Bin Yang 1995 Chatterjec 2005 Doukopoulos & Moustakides 22008). Since providing formal analysis of such tracking capabilities of D-Krsullina and DM-Krasulina for non-i.i.d. data is beyond the scope of this paper, leave it for future work.
13

---

Published in Transactions on Machine Learning Research (10/2022)
4 Convergence Analysis of D-Krasulina and DM-Krasulina
Our convergence analysis of D-Krasulina and DM-Krasulina is based on understanding the rate at which the so-called potential function Vr of these methods converges to zero as a function of the number of algorithmic iterations t. Formally, this potential function V is defined as follows.
Definition 2 (Potential function). Let q' be the eigenvector of Σ associated with X (Σ) and let Ve be an estimate of q" returned by an iterative algorithm in iteration t. Then the quality of the estimate Vt can be measured in terms of the potential function Vz Vs + [0, 1]
that is defoned as
V = 1 (v (v
(10)
1
Notice that V, is a measure of estimation error, which approaches as v, converges to any scalar multiple of q². This measure, which essentially computes sine squared of the angle between q' and ve, is frequently used in the literature to evaluate the performance of PCA algorithms. In particular, when one initializes an algorithm with a random vector Vo uniformly distributed over the unit sphere in I2 then it can be shown that [[ < - 1/d (Balsubramani et al. 2013). While this is statement in expectation for ( = O, our analysis relies on establishing such a statement in probability for any 2 0 for both D-Krasulina and DM-Krasulina. Spocifically, we show in Theorem 3 that sup:20 V. S O(1/d) with high probability as long as 7/- - cL + 1) for any constant and large-enough constant L.
All probabilistic analysis in the following uses a filtration (Fi)c2o of sub o-algebras of F on the sample space 22 where the o-algebra Fr captures the progress of the itcrates of the two proposed stochastic algorithms up to iteration t. Mathematically, let us define the sample covariance matrix A, as A; A = 1 B 1 - Σ Ej Ai.j, for D-Krasulina and DM-Krasulina, respectively. In order to simplify notation and unily some of the analysis of D-Krasulina and DM-Krasulina, we will be resorting to the use of random matrices A,, opposed to Xi, and Xij,t, in the following. We then have the following definition of a-algebras in the filtration.
:= N Em= A.. and
Definition 3 (o-algebra Fi). The o-algebra F. F sample space Ω for both D-Krosulina and OM-Krasulina is defined the o-algebra generaled by the Vo.AA,.... ,o.AA,,.,AA.. Ar), i.e., Fr = o oVvo, Ai,..,,A)).
random variables
In addition to the filtration (Fi)eza, the forthcoming analysis also uses a sequence of nested sample spaces that is defined as follows.
Definition 4 (Nested sample spaces). Let (tu, (o), (t1, (1), (t2, e2).... (tJ,Es) be a sequence of pairs such that 0 = to < h < </2 <... < i. and 6n > 6 > 62 > > 6J > 0 for any non-negative integer J. We then define a sequence (0))Ezz, of nested sample spaces such that Q 3 22 5 O2 5 each Si is
and
22 = { w QUER:voSies
sup Vi(w) S 1-s Gik
(11)
Here, w denotes an outcome within the sample space Ω and (((w) is the (random) potential function after the 1-th ileration of D- Krasuulina/ / DM Krasulina that is being eaplicilly writlen as a function of the outcomes w in the sample space.
In words, the sample space 22, corresponds to that subset of the original sample space for which the error Vt in all iterations E (ly...,, 1-1 1} is below 1-ej, where j E JJ.. ,. In the following, we use the notation Ee(.} and (-).) to denote conditional expectation and conditional probability, respectively, with respect to 2/.
An immediate implication of Definition 4 that, for appropriate choices of cj's, it allows us to focus on those subsets of the original sample space that ensure convergence of iterates of the proposed algorithms to the top eigenvector q' at the desired rates. The main challenge here is establishing that such subsets have high probability measure, i.e., P (Nesof)) 2 1 - for any 8 > 0. We obtain such a result in Theorem in the following. We are now ready to state our main results for D-Krasulina and DM-Krasulina.
14

---

Published in Transactions on Machine Learning Research (10/2022)
4.1 Convergence of D-Krasulina (Algorithm 1
The first main result of this paper shows that D-Krasulina results in linear speed-up in convergence rale function of the number of processing nodes, ie., the potential function for D-Krasulina converges to 0 at a rate of O1/Ntt. Since the system receives a total of Nt samples at the end of t iterations of D-Krasulina, this result establishes that D-Krasulina is order-wise near-optimal in terms of sample complexity for the streaming PCA problem. The key to proving this result characterizing the convergence behavior of D-Krasulina terms of variance of the sample covariance matrix A, = 인ेr 2i=1 Ai, that is implicitly computed within D-Krasulina. We denote this variance NN, which has the form
N
NN Ep.
X
Σ
(12)
Tt is straightforward to see from Definition and 12 that o = o2/N. This reduction variance of the sample covariance matrix within D-Krasulina essentially enables the lincar speed-up convergence. In terms of specifics, we have the following convergence result for D-Krasulina.
Theorem 1. Fir any 8 (0, 1) and pick = co/2(Ai X2) for any co > 2. Next, define
L = 64edr4 max (1,
In 8¹ L = 512e2 max (1,
4
In 8'
(13)
8²
84
pick any L 2 L1t and choose the sizp-size sequence ii == c/ Then, as long as Assumptions
A1
and [A2] hold, we have for D- --rasullia that there evists a sequence (N)ee of nested sample spaces such
that (Tssoli) 2 1 and
L+1 +
of 2)). (.
50 +C
E. (11)SCO
(14)
L
++++1
where C1 and C2 are constants defined as
1/4ed 4ed
Tim
2 X2)/L
207R/
C1 = 2 82
and C2 = (co 2)
Remark 2. While we can obtain a similar result for the case of co S 2, that result does not lead to ony convergence speed-up. up. In particular, the convergence rate in that case becomes O(t which matches the in (Balsubramani et al. 2013).
Discussion. A proof of Theorem which is influenced by the proof technique employed in (Balsubramani
et al.) 2013), is provided in Scction 5 Here, discuss some of the implications of this result, especially in relation to (Balsubramani et al. 2013 and (Jain et al. 2016). The different problem parameters affecting the performance of stochastic methods for streaming PCA include: (i) dimensionality of the ambient space, d, (ii) eigengap of the population covariance matrix, (Ai-A2), - (iii) upper bound on norm of the received data samples, r, and (év) variance of the sample covariance matrix, σ² and/or d Theorem characterizes the dependence of D-Krasulina on all these parameters and significantly improves on the related result provided in [Balsubramani al. 2013).
First, Theorem 1 establishes D-Krasulina can achieve the convergence rate 0 = O with high probability (cf. (14)). This is in stark contrast to the result in (Balsubramani et al. 2013 which is independent of variance of the sample covariance matrix, thereby only guarantecing convergence rate of /tr4/t) for D-Krasulina and its variants. This ability of variants of Krasulina's methods to achieve faster convergence through variance reduction is arguably one the most important aspects of our analysis. Second, in comparison with (Balsubramani et al. 2013), Theorem also results in an improved lower bound on choice of L by splitting it into two quantities, viz., L1 and L2 (cf. (13 ). This improved bound allows larger step sizes, which also results in faster convergence. In terms of specifics, L in the theorem is the
15

---

Published in Transactions on Machine Learning Research (10/2022)
order of (/(r' which is an improvement over 12( / bound of (Balsubramani et al. 2013 On the other hand, while T22 has same dependence on á and d as (Balsubramani et al. 2013) it depends on oj instead of r² and, therefore, it reduces with an increase in N. Third, the improved lower bound on L also allows for an improved dependence on the dimensionality d of the problem. Specifically, for large enough t and N, the dependence on d in (14g is due to the higher-order (first) term and is of the order O(dzaz+ as opposod to (d +co) ) for [Balubbramani et al., 2013). It worth noting here, however, that this is still loser than the result in (Jain et al. 2016) that has only (d) dependence on d in higher-order error terms.
Fourth, in terms of the eigengap, our analysis has optimal dependence of -A2)) which also matches the dependence in (Balsubramani et al. 2013) and (Jain et al.) 2016). It is important to note here, however, that knowledge of the eigengap (λı - A2) is not necessary to run Oja's rule, Krasulina's method, D-Krasulina, or any of the related stochastic methods in a practical setting. Specifically, it can be seen from Theorem [ that the eigengap only needed to set the step size in D-Krasulina for the optimal convergence rate. In practice, however, step sizes of the form ě/t work well for D-Krasulina and the related methods, and a simple yet highly effective strategy for setting the step size in these methods is to estimate the parameter é by running multiple instances of the method during warm-up phase. Such an approach is akin to approximating several problem-related parameters using single parameter é, and is the one we have followed for the numerical experiments discussed in Section 0
Finally, we compare the recommended step-size sequence "/t = cL(LL + t) in Theorem | to the ones in (Balsubramani et al. 2013) and (Jain et al. 2016). Since the step sizes in these two prior works also take the form % c/(L+ + +), all three works are equivalent to each other in terms of scaling of the step size as function of t. But in terms of the initial step size, and assuming small enough 5 in Theorem 1 i is the largest for (Jain et al., 2016), second-largest for this work, and the smallest for (Balsubramani et al. 2013 In relation to our work, this difference in the initial step size in the case of (Balsubramani et al. 2013] is due to the improved lower bound on L in Theorem attributable to the fact that the paramcter L is independent of 5 in that work. Stated differently, we are able to vary the probability of success - in this work by making the parameter be a function of 6, with the caveat being that the initial step size "1 gets smaller as ô decreases. In contrast, a fixed L in (Jain et all, 2016) can be thought of as one of the reasons the probability of success is fixed at 3/4 in that work. We conclude by noting that this dependence of the performance of D-Kasuulina on different problem parameters is further highlighted through numerical experiments in Section G
In the case of (Jain et al. 2016),, this difference is
Remark 3. While Theorem is for (a distributed variant of) Krasulina's 's method, Oja's rule can also be analyzed using similar techniques; see, e.g., the discussion in (alsubrumanii et al., 2013).
Remark 4. Recall from the discussion in Section 1.1 that an iteration of Krasulina's 's method is similar to that for SGD applied to the optimization problem ( A natural question to ask then is whether D-Krasulina can be "accelerated" in much the same way SGD can be accelerated by adding a momentum term to its iteration expression. The authors in (D Sa et al. 2018), however, argue that naively applying momertum to Oja's 's rule or the power ileralion, both of which are closely related to Krasulina's 's method, resulls in worst performance since this increases the effect of the noise within the iterates. And while the noise within the iterates can be controlled through variance reduction techniques, as done in De Sa et al., 20188) to accelerate the power iteration for eigenvector computation, such techniques typically require multiple data passes and are therefore not suited for the selling in which dala samples continuoushy stream inlo the system.
4.2 Convergence of DM-Krasulina (Algorithm 2
The convergence analysis of DM-Krasulina follows from slight modifications of the proof of Theorem [ for D-Krasulina. The final set of results, which covers the two scenarios of zero data loss (μ 1u=00 0) and some data loss (μ > 0) in each iteration, is characterized in terms of variance of the (mini- batched sample covariance A, : B & Σ/= Ai,j, associated with DM-Krasulina. We denote this variance as of, which is given by
2 0
Ep.
Σ
(15)
j.
16

---

Published in Transactions on Machine Learning Research (10/2022)
It is once again straightforward to that o = oZ/B. We now split discussion of the convergence of DM-Krasulina according to the two scenarios discussed in Section 3.2
4.2.1 Scenario 1 1-DM-Krasulina with no data loss:
BR. R. Rp (blte R.)
p=0
Analytically, this scenario is similar to D-Krasulina, with the only difference being that we are now incorpo- rating an average of B sample covariances x.. in the estimate in each iteration (as opposed to N sample covariances for D-Krasulina). We therefore have the following generalization of Theorem in this scenario.
Theorem 2. Let the parameters and constants be as specified in Theorem is now defined 12 = 512e" d' max(1.c²) In 4.. Then, as long as Assumptions
except that the parameter L
A1 and A hold, we have
for DM-Krsulina that TP (Nisoli) 2 1 5 and
Ex ( S Ci
( ++++1 (1+))
L+1 +C2
(1
B I+L+1
(16)
The proof of this theorem can be obtained from that of Theorem by replacing I/N and o in there with 1 1/B and o. respectively. Similar to the case of D-Krasulina, this theorem establishes that DM-Krasulina can also achieve linear speed-up in convergence as a function of the network-wide mini-batch size B with very high probability, i.c., E, { (v)) = ( O(o2/Bt)
Our discussions of D-Krasulina and DM-Krasulin have so far been focused on the infinite-smmple regime, in which the number of algorithmic iterations t for both algorithms can grow unbounded. We now focus on the implications of our results for the finite-sample rogime, in which a final estimate is produced at the end of arrival lotal T > 1 samples P This finile-sampl regime leads to an interesting interplay between (resp., B) and the total number of samples T for linear speed-up of D-Krasulina (resp., DM-Krasulina). We describe this interplay in the following for DM-Krasulina; the corresponding result for D-Krasulina follows by simply replacing B with N in this result.
Corollary 1. Let the parameters and constants be specified in Theorem Next, pick parameters (T4 I2) such that T. 2 T, and T/2 2 Tn/oks, and define the final number of algorithmic ileralions for MM-KKasulina as To = T/B. Then, as long as Assumptions
A1 and A2
hold and the network-wide
mini-batch size salisfies < T co we have that IP (NOCISTORi)
2 1
8 and
I 2/22
0
ca/2
C20²
ETa (Mra) coC1
+
(17)
T
T
Proof. Substituting & TB in 16 and using simple upper bounds yield
Ern ( Tn I C
L+
L 1+1)+
o +
201
L
+
+C2
BB
TB
TB
TB
Next, substituting L = L3 + L in this expression gives us
L 4 + coCi
( L2
Er, (1TT,
) S Ci
+ C2
(18)
TB
TB
TB
Since 02 =02/B and TB T/B, 18 reduces to the following expression:
ca/2
c//2
BL
a Lb
Cza²
Er, (Vra < caCi
+
+
T
T
T
2
The proof now follows from the assumption that B
<
cg
6 implicit assumption here that large enough that it precludes the batch PCA algorithm.
17

---

Published in Transactions on Machine Learning Research (10/2022)
Discussion. Corollary dictates that linear convergence speed-up for DM-Krasulina (resp., D- D-Krasulina) occurs in the finite-sample regime provided the network-wide mini-batch size. B (resp., number of processing nodes N) scales sublinearly with the total number of samples T'. In particular, the proposed algorithms achiove the best (order-wise) convergence rate of O(1TT) for appropriate choices of system paramcters. We also corroborate this theoretical finding with numerical experiments involving synthetic and real-world data in Section
4.2.2 Scenario 2 2-DM-irasullina with data loss: N
bR. R. R. (bR.. R.)
p 20
The statement of Theorem 2 for DM-Krasulina the lossless setting immediately carries over to the resource- constrained setting that causes loss of μ (> 0) samples per iteration. The implication of this result is that DM-Krasulina achieve convergence rate of O(1/Bt) in the infinite-sample regime after receiving total of (B + (1)t samples. Therefore, it trivially follows that DM-Krasulina can achieve order-wise near-optimal convergence rate in the infinite-sample regime as long as O(B).
We now turn our attention to understanding the interplay between μ4, B, and the total number of samples T arriving at the system for the resource-oonsrriined finite-sample setting for DM-Krasulina To this end, we have the following generalization of Corollary !
Corollary 2. Let the parameters and constants as specified in Corollary & and defone the fonal remmber of algorithmic ilerations for DM-Krasulina as T == T/(B + p). Then, as long as Assumplions
A1
and
A2 hold, we have that IP (Nez0l)) 2 11-60 and
ca/2
20/2
{ Vry S coCi
(B (B+1)T1 T
(B (B+ 10++))22T)
C20² +
Eyg
+ coc
+
(19)
BT
BT
Proof. The proof of this corollary follows from replacing TB with T in 18 and subsequently substituting the values of Th and of in there.
Discussion. Recall that since the distributed framework receives a total of T samples, is desirable to achieve convergence rale of O(/TT). It can be seen from Corollary 2 that the first and the third terms (19 are the ones that dictate whether DM-Krasulina can achieve the (order-wise) optimal rate of O(L/T). To this end, the first term in 19 imposes the condition (B+ S T1-2/00 i.c., the total number of samples received at the system (both processed and discarded) per iteration must scale sublinearly with the final number of samples T. In addition, the third term in 19 imposes the condition μ = O(B), i.e., the number of samples discarded by the system in each iteration must scale no faster than the number of samples processed by the system in each iteration. Once these two conditions are satisfied, Corollary 2 guaaaantces near-optimal convergence for DM-Krasulina.
5 Proof of the Main Result
The main result of this paper is given by Theorem [ which can then be applied to any algorithm that (implicitly or explicitly) involves an iteration of the form (8). We develop a proof of this result in this section, which-similar to the approach taken in (Balsubramani et al. 2013) for the analysis of Krasulina's 's method-consists of characterizing the behavior D-Krasulina three different algorithmic epochs. The result, concerning the initial epoch is described in terms of Theorem 3 in the following, the behavior of the intermediate epoch, which comprises multiple sub-epochs, is described through Theorem 4 while the behavior of D-Krasulina in the final epoch is captured through formal proof of Theorem | at the end of this section.
Before proceeding, recall that our result requires the existence of a sequence (R)rez, of nested sample spaces that are defined in terms a sequence of pains (to = 0, eo), (f1, (ty, Our analysis of the initial epoch involves showing that for the step size ?r chosen as in Theorem [ the error for all t 2 0 will be less than (1 - cu) with high probability for some constant ep- We then define the remaining ej's Ej Deo.j 1....,d J, where J is defined as the smallest integer satisfying e, 2 1/2. Our analysis
18

---

Published in Transactions on Machine Learning Research (10/2022)
in the intermediate epoch then focuses establishing lower bounds on the number of iterations t, for which D-Krasulina is guaranteed to have the error less than 1 1-6 ej with high probability. Stated differently, the intermediate epoch characterizes the sub-epochs {1 + tz-1,tz} during which the error is guaranteed to decrease from (1 c--1 to (1 - -cy) with high probability.
5.1 Initial Epoch
Our goal for the initial epoch is to show that if we pick the step size appropriately, i.e. we set L to be large enough (ef. (I3)), then the error, Pe, will not exceed a certain value with high probability. This is formally stated in the following result.
Theorem 3. Fir any 1), define e (0, as c: 82 /8e, and let
4 8d max 8dr max(1,
4
L >
In
In 62
(20)
g
0
Then, if Assumptions A1
and
A2
hold and we choose step size to be 7e = CL +
we have
P sup vr 21 - 2) s 22c 2
(21)
120
In order prove Theorem g we need several helping lemmas that are stated the following. We only provide lemma statements in this section and move the proofs Appendix A We start by writing the recursion of error metric V, in the following lemma.
Lemma 1. Defining scalar random variable
== 2%
(22)
Vz-1 12
we get the following recursion:
(1 # Ein Air > x/% - Et, and
(1) v. S V + 472
(8) V. < 1 + 72 - 23.
Proof. See Appendix A.1
Part (i) of this lemma will be used to analyze the algorithm in the final epoch for proof of Theorem while Part (it) will be used to prove Theorem for this initial epoch and Theorem 4 for the intermediate phase.
Next we will bound the moment generating function of V, conditioned on Fi-1 (Definition 3 For this, we need an upper bound on conditional variance of 2t, which is given below. Lemma 2. The conditional variance of the rundom variable 2, given by
B{(a, - -(Zz)) S 16
(23)
Proof. See Appendix A.2
Using this upper bound on conditional variance of can now upper bound the conditional moment generating function of V, In order to simplify notation, much of our discussion in the following will revolve around the moment generating function with parameter & E $ == (d/4e, (2/eo) In(4/))}. Note, however, that similar results can be derived for any pasluede-valued parameter E R.
Lemma The conditional moment generuling function of V. for s upper bounded as
Ele S exp sxxp
186111 -sE(x + sf 4 222003
(24)
19

---

Published in Transactions on Machine Learning Research (10/2022)
Proof. See Appendix A.3
Note that this result is similar to (Balsubramani et al. 2013 Lemma 2.3) with the difference being that the last term here is sample variance, TN as opposed to upper bound on input ||xe|la S r in (Balsubramani ett al. 2013, Lemma 2.3). This difference prompts changes in next steps of the analysis of D-Krasulina and it also enables us to characterize improvements in convergence rate of Krasulina's method using iterations of the form (8)
We
ready to prove the statement of Theorem 3 which is based on Lemma
and
Proof of Theorem We start by constructing supermartingale from sequence of errors V, First, restrict- ing ourselves to S, we define quantities
: 22,, C ===syoo
Te == (Bi+ (t), and M; == exp + ( ST:).
1st
Now, taking expectation of M, conditioned on the filtration Fi-1 we get
(a) E(exp (sW.)) exp (STr) S exp (1 + +
E[M,[Fi->> =
+ she + STL)
exp (sW, + s7;-1) = M
Here, (a) is due to Lemma and using the fact that B(2 2 0 Balsubramani et al. 2013 Theorem 2.1).
These calculations show that the sequence {Mr} forms a supermartingale. Using sequence Mc. we can now
use Doob's martingale incquality (Durrctt 2010 pg. 231) to show that V, will be bounded away from 1 with high probability. Specifically, for any A E (0, 1), we have
P
sup V 120
A)
< TP
sup V, + Te 2 A) 120
P supexp (sV, + 87.) e
120
MM ea
= sup M, N ) (20
exp (-o(A - -T)) (00).
Substituting A 1 - c/d and using (Balsubramani et al.
2013
Lemma 2.5) to bound VV get
P sup v, 2 1 - 2) s exp (- (e/d) - -o))ee 120
d
(25)
2s'
Next need to bound Σ1>0 Bi and Elso Gt.
First we get
c
A²2
Σ B =
2ir =
ri
S
(26)
(1
L)2
L
t>0
t>0
150
150
Again using a similar procedure we get
N c
ΣG
(27)
L
t>0
Combining
26 and 27
along with the definition of Ti at the beginning, we gel
To S 2² rd + so (0+
(28)
I
Now using the lower bound on L, we get To s c/d for s = d/4e as shown in Proposition in Appendix
D
Substituting this in (25 we get
P sup v 2 1 - 9) S exp (-s(1 c/d - c/d])e" V
d = exp (2se/d)
d
10
2s
2s
Finally, substituting 8 = d/4e, we get P( (suppoo V 2 1 - 5) s /2e.
20

---

Published in Transactions on Machine Learning Research (10/2022)
5.2 Intermediate Epoch
In Theorem 3 we have shown that if we choose L such that it satisfies the lower bound given in Theorem 3 then we have error P, greater than 1 --60 (here, = 6 /Sed) with probability 8. Next, our aim is to show that if we perform enough iterations {, of D-Krasulina then for any t 2 ł the error in the iterate will be bounded by V, I 1/2 with high probability.
order to prove this, we divide our analysis into different sub-epochs that are indexed by j E {1,.... J}. Starting from 6 we provide a lower bound on the number of iterations ij such that we progressively increase Ej in each sub-epoch until we reach eJ.
Theorem 4. Fix any (0,1) and pick e = co c0/1((1--x)) for any co > 2. Next, let the number of processing nodes N 1, the parameter L 8r max(1,e²) In 3 s0 max(1.c3 In 1 and the step size 7t := cL(+t). Finally, select schedule co), (41, (...... (1,,6)) such that the following conditions are satisfied:
%
[C1] co Sed' 29 S cjt1 s 26, for 0 oSi J, and CJ S 1 and
52 3
(C2)
(f + 2 5/10 (14 (t ((++++1) for jsj sjjj C.
Then P (Nisoli) 2
In order to prove this theorem, we need Lemmas Lemma 4. For t > tj, the moment generating function of W for s E S conditioned on 22 satisfies
7 which are stated as follows.
Eg Erfo e < exp
(s(We-()
s (
cocj)
(t L) + sc203
car
(t + L)2 )
t ++
Proof. See Appendix B.1
Lemma For > is and E S, we have
ty ++
+ (00 )
1 ( =+))
Ex(e"Vi s exp s(1 c()-
(29)
1+L+I
Proof. See Appendix B.2
Using Lemma our next result deals with a specific value of i, namely, l. lj+1-
Lemma Suppose Conditions
Cl
C2 are satisfied. Then for 0 S j < J and 8 E S, we get
Bazss (e (54)) I eex s(1 (41) :-1- se; + (sc2pt + + 2000)) ( 1 -z)
Proof. See Appendix B.3
Lemma 7. Suppose Conditions
C1
C2
are satisfied. Then picking any 0 < 8 < 1, we have
sup V, > e1) s 8 2 121,
>
Pt,
1=1
Proof. See Appendix B.4
Proof. (Proof of Theorem Using results from Lemma and Theorem and applying union bound, we get the statement of Theorem
21

---

Published in Transactions on Machine Learning Research (10/2022)
5.3 Final Epoch
Now that we have shown that V. S 1/2 with probability 1 1 8 for all L 2 LJ, characterize in the final epoch how Vr decreases further as function of algorithmic iterations. The following result captures the rate at which er decreases during this final epoch.
Lemma 8. For any t > ts and e ca/(>1 - >2), the (conditional) expected error in Ve is given by
4c203 (t L)²
E& - 1
2(t 2(t+ L) E ( 2(t + ()1 X2)2
Proof. See Appendix
We ready to prove our main result, which is given by Theorem
Proof. (Proof of Theorem 7 Recall the definitions of the sub-epochs corresponding to the pairs (4i,,) that
satisfy the two conditions in Theorem 4 Following the same procedure in the proof of (Balsubramani et al.
2013 Theorem 1.1), notice that J log2 (1/(2e)) (since EJ 2cs 22 + 2 eee = 1/2co) and therefore Condition
C2
implies
5J (L + 1)
5/(co 22)
4ed 82
5/(co 1n
tu +++++= (L + 1)
exp
= (L [
(30)
co
200
Defining ai :- / - X2)², 02 co/2, b = 4c²0% and using Lemma 8 for t >> is we have
Ed(w,} < 1 +
al
a2 E. E + (t+ L)
b
(t + L)
+ L
Now using Proposition from Appendix
with cp
2, we get
t ++1
1 L++
exp
1
+ L+1
a1
E (11r}
+
b
1+
exp
01
02 -1
t +
L
+
t, + L +
t+1+1
(a) 1
L+1
4ed
(eg
a1
exp
2 + L+1
82
t + L+1
b
2
d1
+
exp
exp
02
tJ +
L+1
tJ + L +
t+L+
1 L.+1 1
/eed
2 + L L+1
82
exp
(L + 1+)ddd///8))) (22)/5222)))
8
2 01
1
+ co 2
exp
(L+ + ()ed (2)dd/22266/2m21)
(t (++
Here, the inequality in (a) is due to 330 and we have also used the fact that (1 + x²)= S (ar) for x
In addition, since (4ed/ (5/2221m2)
1, we get
B.(,,) < + 11
L+1
+ 4ed
a1
SC"OX
(11 + 2
1
exp
exp
82
L+11
CD -2
(t + 1) L+11
$ 1 L+1 1
4ed
Sc 2 (
c 01/L
+
2 / 111
82
Co 2
(1+ L+1)
= =a
L+1
C2
N
2
(31)
+++
++1+11
This completes the proof of the theorem.
22

---

Published in Transactions on Machine Learning Research (10/2022)
10
10
10
10 2
1g" v,""
ig w.
4 B
10
B 10
10%
- B 100
B 500
/ 10
1000
/ 100
10
2000
P 200
10
102
10
10
100
10
10
# samples (T) received
# of samples (T) received
(a) Impact of the mini- batch size the convergence of DM- Ma--Krauulins for the resourceful regime. Note that the plot is effectively Krasulina's method.
(b) Performance of DM-Krasulina &R. R, regime (i.e., N < R,(RRe-B Ra) R.y.. which causes loss of μ sam- ples per iteration; here, (N, B) = (10, 100).
resource-cunsrrineed
Figure 2: Convergence behavior of DM-Krasulina for the case of synthetic data under two scenarios: (a) No data loss (μ = 0) and (b) loss of p > 0 samples per algorithmic iteration.
6 Numerical Results
In this section, we utilize numerical experiments to validate the theoretical findings of this work in terms of the ability of implicit/explicit mini-batched variants of the original Krasulina's method (Krasulina) 1969)
to estimate the top eigenvector of covariance matrix from (fast) streaming data. Instead of repeating the same set of experiments for the original Krasulina's 's method, D-Krasulina, and DM-Krasulina, we present our results that are parameterized by the network-wide mini-batch size B E {1) b Z that appears in DM-Krasulina. This is because B trivially corresponds to the original Krasulina's 's iterations, while N corresponds to iterations that characterize D-Krasulina.
Our goals for the numerical experiments are threefold: (i) showing the impact of immplici/expplicee mini- batching the convergence rate of DM-Krasulina, (di) establishing robustness of DM-Krasulina against the &R.R. loss of μ 0 samples per ileralion for the case when N < RR and (iii) experimental validation for scaling of the convergence findings, namely, eigengap (A1 - -2), dimensionality (d), and upper bound on input samples ([|xe llz S r). In the following, we report results of experiments on both synthetic and real-world data to highlight these points. Since the main purpose is to corroborate the scaling behaviors within the main results and not to investigate additional system-related issues concerned with large-scale implementations the real-world datasets are chosen to facilitate their processing on low-cost compute machines.
in terms of the problem parameters as predicted by our theoretical
6.1 Experiments on Synthetic Data
In the following experiments generate T 10 samples from some probability distribution (specified for each experiment later) and for each experiment we perform 200 Monte-Carlo trials. In all the experiments c/t. We performed experiments with multiple values of
in the following we use step size of the form and here we are reporting the results for the value e which achieves the best convergence rate. Further details about each experiment are provided in the following sections.
6.1.1 Impact of mini-batch size on the performance of DM-Krasulina
For covariance matrix Σ g with Å₁ and eigengap Aι >2 0.2, we generate T 1 samples from (0, Σ) distribution. The first set of experiments here deals with the resourceful regime, i.e., N bR_
(blt, - with
of sizes B {1, 10, 100, 500, 1000, 20003. Note that these values of B can be
23

---

Published in Transactions on Machine Learning Research (10/2022)
10
10
10
i9 v.
10"
A1 À 0.1
10
10.8
A1 b 0.2
- d
A1 A 0.3
d
10
8- A1 A2 0.4
15
À A 0.5
10%
10
100
202
10
10
100
10
10
10
# of samples (7) received
# of samples (2) received
(a)
(b)
Understanding the impact of (a) eigengap (A1 -22) and (b) dimensionality on the convergence behavior of DM-Krasulina, corresponding to B 1000 and μ 0.
Figure 3:
factored into any positive integers and N as long as the condition N R, that governed by the application scenario and the physical system is satisfied. It is, therefore, unnecessary to specify b and N for these experiments, whose results are shown in Figure ((a) These results are obtained for step-size parameter с E {70, 80, 80, 90, 110, 100), which are the values of c resulting in the best convergence rate. As predicted by Corollary [1 we can see that after T/B iterations of DM-Krasulina, the error VI/B is on the order of OLL/T) for B E 10, 100, 500, 1000}, while for = 2000, the error VryB not optimal anymore.
RR.
Next, demonstrate the performance of DM-Krasulina for resource constrained settings, i.e., R,j. which causes the algorithm to discard == [BR R2
N < Rppblle RR
+ NR. B samples per itera- Rh
tion. Using the same data generation sctup as before, we run DM-Krasulina for network of 10 nodes (N 10) with network-wide mini-batch of size B 100 (i.e., b 10). We consider different mismatch factors between streaming, processing, and communication rates in this experiment, which result in the number of samples being discarded μ {0, 10, 100, 200}. The results are plotted in Figure 2(b) which shows that the error (T)(B+)) for 10 is comparable to that for f 0, but the error for μ 200 is an order of magnitude worse than the nominal error.
6.1.2 Impact of the eigengap on the performance of DM-Krasulina
For this set of experiments, we again generate data in RS from normal distribution (/ E), where the covariance matrix Σ has the largest cigenvalue X 1. We then vary the remaining eigenvalues to ensure eigengap that takes values from the set {0.1, 0.2, 0.3, 0.4, 0.5). The corresponding values of that give the best convergence rate for cach unique eigengap satisfy E [180, 110, 90, 70, 60). The final results for these experiments are plotted in Figure 3(a) for the case of B = 1000 and l 0. These results establish that the final gap in error after observing 10 data samples indeed on the order 0(1/(A1 - -A2)?), as suggested by the thcorctical analysis.
6.1.3 Impact of dimensionality on the performance of DM-Krasulina
For this set of experiments, is varied such that d E {5, 10, 15, 20). In addition, we fix the largest cigenvalue of Σ to be \ = and its eigengap to be 0.2. The values of corresponding to each unique value of that provide the best convergence rate in these experiments satisfy {110, 110, 100, 100]; contrary to our theoretical analysis, this soems suggest that the optimal step-size sequence does not have strong dependence on d, at least for small values of d. We also plot the potential function for each d a function of the number of received samples in Figure 3(b) for the case of B 1000 and a 0. Once again, observe little dependence of
generate data in R from normal distribution (V (O Σ) whose dimensionality
24

---

Published in Transactions on Machine Learning Research (10/2022)
100
10% 2
[at v.
10
1
2
10.8
100
102
10
10
# of samples (T) received
Figure 4: Performance of DM-Krasulina for varying upper bound on the norm of the strcaming data.
the performance of DM-Krasulina on d. Both these observations suggest that our theoretical analysis is not tight in terms of its dependence on dimensionality of the streaming data.
6.1.4 Impact of upper bound on the performance of DM-Krasulina
In order to understand the impact of the upper bound ||xe I12 I r on the convergence behavior of DM-Krasulina, we generate Xy E R as X, Cu,- with uy E R having independent entries drawn from uniform distribution u(-a, a) and C chosen ensure an eigengap of 0.2 for the covariance matrix. As vary the value of a within the set f1, 10), generate four different datascts of T 10 samples for which the resulting
(1.45, 2.9, 4.5, 14.5). The values of e that provide best convergence for these values of r satisfy e E {8, 2, 1, 0.08). The final set of results are displayed in Figure ! for B = 1 and 0. It can be seen from this figure that changing r does not affect the convergence behavior of DM-Krasulina. This behavior can be explained by noticing that the parameter r appears in our convergence results terms of a lower bound on (cf. (13) and within the non-dominant term in the error bound. The dependence of on the parameter
already being reflected here in our choice of the step-size parameter e that results in the best convergence result. In addition, hypothesize that the non-dominant error term in our experiments, compared to the dominant one, is significantly small that it masks the dependence of the final error on
6.2 Experiments on Real-world Datasets
In this section, evaluate the performance of DM-Krasulina on two real-world datasets, namely, the MNIST
dataset [LeCun 1998) and the Higgs dataset (Baldi et al. 2014 The MNIST dataset corresponds to d = 784 and has a total of T 6 x 10d samples, while the Higgs dataset is d = 28 dimensional and comprises 1.1 x 10 samples. It is worth noting here that since it is straightforward to store all the samples in these datasets at a single machine, one can always solve the 1-PCA problem for these datasets without resorting to the utilization of a distributed streaming framework. Nonetheless, is still possible to utilize these dataset a simulated distributed strcaming sctting order to highlight the agreement between the scaling behavior predicted by our theoretical results and the scaling behavior observed using real-world datasets; this is indeed the purpose of the following sets of experiments.
Our first set of experiments is for the MNIST datasct, in which we use the step size ? = c/t with e {0.6, 0.9, 1.1, 1.5, 1.6) for nelwork-wide mini-batch sizes B E (1, 10, 100, 300, 1000) in the resourceful regime (μ = 0). The results, which are averaged over 200 random initializations and random shuffling of data, are given in Figure 5(a) It can be seen from this figure that the final error relatively stays the same as increases from 100, but it starts getting affected significantly as the network-win mini-batch size is further increased to B 300 and B 1000. Our second set of experiments for the MNIST dataset corresponds to
25

---

Published in Transactions on Machine Learning Research (10/2022)
10
10%
10
10
ig v.y0
10
- B =1
10
+ 0
10.3
B 10
10
9 B 100
300
1000
10
10
10
10
10
10
100
10
10
10
# of samples (T) received
# of samples (T) received
(a) MNIST Data (nr O): Impact of network- wide mini-batch size on the convergence behavior DM-Krasulina for the resourceful regime.
100): Convergence be- havior of DM-Krasulina resoure-constrained regime, which causes of samples per iteration.
(b) MNIST Data (N
Figure 5: Performance of DM-Krasulina for the MNIST dataset under two scenarios: (a) No data loss (μ 0) and (b) loss of y samples per algorithmic iteration.
10%
10
10
10
(6) vv
10
10
4 F
100
10000
1000
20000
2000
1
10
10
2
10
80
108
10
102
10
10
10
# of samples (I)) received
of samples (I') received
(a) Higgs Data (μ 0): Impact of network-wiid mini-batch size B on the convergence behavior of DM-Krer OM-Krasulina for the resourceful rogime.
(b) Higgs Data (N = 10; B = 1000): Convergence behavior
of DM-Krasulina in a resonree-onsttraneed regime, which causes loss of samples iteration.
Figure 6: Performance of DM-K MraKuraullna for the Higgs dataset under two scenarios: (a) No data loss (μ = 0) and (b) loss of y >> 0 samples per algorithmic iteration.
the resource- consurriee regime with (N, = (10, 100) and step-size parameter E {0.6, 0.9, 1.1, 1.6] for the number of discarded samples / E (0, 10, 20, 40, 1003. The results, averaged over 200 trials and given in Figure 5(b] show that the system can tolcrate loss of some data samples per itcration without significant increase in the final error; the increase in error, however, becomes noticeable as μ approaches B. Both these observations are in line with the insights of our theoretical analysis.
We now turn our attention to the Higgs datasct. Our results for this dataset, averaged over 200 trials and using e 0.07, for the resourceful and resource-consstrained settings are given in Figure [(a] and Figure 6(b) respectively. In the former setting, corresponding to e {1, 102 10a, 10 10*), }, we once again see that the error relatively stays the same for values of B that are significantly smaller than T; in particular, since T for the Higgs dataset is larger than for the MNIST dataset, can accommodate a larger value of B without significant loss in performance. In the latter resource-constrained setting, corresponding to N 10,
26

---

Published in Transactions on Machine Learning Research (10/2022)
B 1000 and μ {0, 10, 100, 1000, 2000), similarly observe that small (relative to B) values of μ do not impact the performance of DM-Krasulina in significant manner. Once again, these results corroborate our research findings.
7 Conclusion
In this paper, we studied the problem of estimating the principal cigenvector of a covariance matrix from independent and identically distributed data samples. Our particular focus in here was developing and an- alyzing two variants, termed D-Krasulina and DM-Krasulina, of a classical stochastic algorithm that can estimate the top cigenvector a ncar-optimal fashion from fast streaming data that overwhelms the pro- cessing capabilities of single processor. Unlike the classical algorithm that must discard data samples in high-rate streaming settings, and thus sacrifice the convergence rate, the proposed algorithms manage the high-rate streaming data by trading off processing capabilities with computational resources and commu- nications infrastructure. Specifically, both D-Krasulina and DM-Krasulina virtually slow down the rate of streaming data by spreading the processing of data samples across of a network of processing nodes. addi- tion, DM-Krasulina can overcome slower communication links and/or or lack of sufficient number processing nodes through a network-wide mini-batching strategy, coupled with discarding of a small number of data samples per iteration.
Our thcoretical analysis, which fundamentally required characterization of the error incurred by the pro- posed algorithms function of the variance of the sample covariance matrix, substantially improved the variance-agnostie analysis in (Balsubramani et al., 2013) and established the conditions under which near- optimal convergence rate is achievable in the fast streaming setting, even when some data samples need to discarded due to lack of sufficient computational and/or or communication resources. We also carried out numerical experiments on both synthetic and real-world data to validate our theoretical findings.
In terms of future work, extension of our algorithmic and analytical framework for estimation of the principal subspace comprising multiple eigenvectors remains open problem. In addition, tightening our theoretical analysis to better elucidate the role of dimensionality of data in the performance of the proposed algorithmic framework is an interesting problem. Finally, investigation of additional practical issues (e.g., processor failures, variable compute costs, and network coordination costs) concerning processing of data in largo-scale systems provides another avenue for future research.
Funding Acknowledgements
This work has been supported part by the National Science Foundation under Awards CCF-1907658, OAC-1940074, and CNS-2148104, and by the Army Research Office under Awards W91-NF-17-1-0546 and W911NF-21-1-0001.
References
Alekh Agarwal and John Duchi. Distributed delayed stochastic optimization. In Proc. Advunces Neural Inform. Process. Syst. (NewrIPS), pp. 873-888 2011.
Zeyuan Allen-Zhu. How to make the gradients small stochastically: Even faster convex and nonconvex SGD. Proc. Advances Neural Inform. Process. Syst. (NewrIPS), pp. 1157 1157-1167, 2018a.
Zeyuan Allen-Zhu. Natasha 2: Faster non-convex optimization than SGD. In Proc. Advances Neural Inform Process. Syst. (NeurIPS), pp. 2680 2691, 2018b.
Zeyuan Allen-Zhu and Elad Hazan. Variance reduction for faster non-convex optimization. In Intl. Conf. Mach. Learning (ICML), pp. 699 707, 2016.
Zeyuan Allen-Zhu and Yuanzhi Li. First efficient convergence for strcaming k-PCA: global, gap-free, and near-optimal rate. In Proc. IEEE 58th Annu. Symp. Found. Comput. Sci. (FOCS), pp. 487-492. IEEE, 2017a.
27

---

Published in Transactions on Machine Learning Research (10/2022)
Zeyuan Allen-Zhu and Yuanzhi Li. Follow the compressed leader: Faster online learning of eigenvectors and faster MMWU. In Proc. 34th Inl. Conf. Mach. Learning. (ICML), pp- 116-125. JMLR. org, 2017b.
Ehsan Amid and Manfred K Warmuth. An implicit. form of Krasulina 's k-PCA updale without the orthonor- mality constraint. arXiv preprint arXiv:1909.00 04800, 2019.
Raman Arora, Andy Cotter, and Nati Srebro. Stochastic optimization of PCA with capped MSG. In Proc. Advances Neural Inform. Process. Syst. (NeurIPS), pp. 1815-1823, 2013.
Tuncer Can Aysal, Mehmet Ercan Yildiz, Anand D Sarwate, and Anna Scaglione. Broadcast gossip algo- rithms for consensus. IEEE Trans. Signal Process., 57(7):2748 2761, 2009.
Maria Florina Balcan, Yingyu Liang, Le Song, David Woodruff, and Bo Xic. Communication efficient distributed kernel principal component. analysis. In Proc. 22nd ACM SIGKDD Jul. Conf. Knowledge Discovery and Data Mining, pp. 725-734, ACM, 2016.
Pierre Baldi, Peter Sadowski, and Daniel Whiteson. Searching for exotic particles in high-energy physics with deep learning. Nature communications, 5:4308, 2014.
Akshay Balsubramani, Sanjoy Dasgupta, and Yoav Freund. The fast convergence of incremental PCA. In Proc. Advances Newral Inform. Process. Syst. (NewrIPS), pp. 3174 3182, 2013.
Bin Yang. Projection approximation subspace tracking. 1EEE Trans. Signal Process., 43(1):95 107, 1995.
Vincent D Blondel, Julien M Hendrickx, Alex Olshevsky, and John N Tsitsiklis. Convergence in multiagent coordination, consensus, and flocking. In Proc. 44th IEEE Conf. Decision and Control, pp. 2996 3000. IEEE, 2005.
Léon Bottou. Large-scale machine learning with stochastic gradient descent. In Proc. 19th Intl. Conf. Computational Statistics (COMPSTAT10), pp- 177-186, 2010. doi: 10.1007 10.1007/978-33 000077978-377900 7908----00-2 16.
Stéphane Boucheron, Gábor Lugosi, and Pascal Massart. Concentration inequalities: A nonasymptolic theory of independence. Oxford University Press, 2013.
Christos Boutsidis, David P Woodruff, and Peilin Zhong, Optimal principal component analysis in distributed and streaming models. In Proc. 48th Annu. ACM Symp. Theory Compuling (STOC), pp. 236 249. ACM, 2016.
Chanchal Chatterjee. Adaptive algorithms for first principal eigenvector computation. Neural Networks, 18: 145-159, 2005. ISSN 08936080. doi: 1.1116//..eumee.20004
Kenneth L Clarkson and David P Woodruff. Low-rank approximation and regression in input sparsity time. J. of the ACM, (6):11 45, 2017.
Andrew Cotter, Ohad Shamir, Nati Srebro, and Karthik Sridharan. Better mini-batch algorithms via ac- celerated gradient methods. In Proc. Advances Newral Inform. Process. Syst. (NeurIPS), pp. 1647 1655, 2011.
Christopher De Sa, Kunle Olukotun, and Christopher Ré. Global convergence of stochastic gradient descent for some non-convex matrix problems. In Proc. 32nd Int. Conf. Machine Learning (TCML), pp. 2341, July 2015. URL http://arxivv http//bs//411....
Christopher De Sa, Bryan He, Ioannis Mitliagkas, Christopher Ré, and Peng Xu. Accelerated stochastic power itcration. Proc. Mach. Learning Res., 84:58, 2018.
Ofer Dekel, Ran Gilad-Bachrach, Ohad Shamir, and Lin Xiao. Optimal distributed online prediction using mini-batches. J. Mach. Learning Res. (JMLR), 13(Jam):155220. 2012.
Alexandros G Dimakis, Soummya Kar, José MF Moura, Michael G Rabbat, and Anna Scaglione. Gossip algorithms for distributed signal processing. Proc. IEFE, 98(11):1847-1864, 2010.
28

---

Published in Transactions on Machine Learning Research (10/2022)
Xenofon G. Doukopoulos and George V. Moustakides. Fast and stable subspace tracking. IEEE Trans.
Signal Process., 56(4):1452-1445 1465,, 2008. ISSN 1053587X. doi: (0.1109/SSP2000.999993 2093339
Rick Durrett. Probability; theory and ezamples. Cambridge University Press, 2010.
Alan Frieze, Ravi Kannan, and Santosh Vempala. Fast Monte-Carlo algorithms for finding low-rank approx- imations. J. of the ACM, 51(())1025-1000 2004.
Arpita Gang and Waheed Bajwa. FAST-PCA: A fast and exact algorithm for distributed principal component analysis. Xixiv preprint or Xiv:21100 12979, 2021.
Arpita Gang and Wahood U Bajwa. A linearly convergent algorithm for distributed principal component analysis. Signal Process., 193:108408, 2022.
Arpita Gang, Bingqing Xiang, and Waheed U Bajwa. Distributed principal subspace analysis for partitioned big data: Algorithms, analysis, and implementation. IEEE Trans. Signo Inform. Process. over Netw., 7: 699-715, 2021.
Dan Garber. On the regret minimization of nonconvex online gradient ascent for online PCA. arXiv preprint arXiv:1809.1044 10491, 2018.
Dan Garber and Flad Hazan. Fast and simple PCA via convex optimization. In Proc. 32nd Ind. Conf. Machine Learning (ICML), 2015.
Dan Garber, Elad Hazan, and Tengyu Ma. Online learning of eigenvectors. In Proc. 32nd Int. Conf. Machine Learning. (ICML), pp. 560-568, 2015.
Dan Garber, Ohad Shamir, and Nathan Srebro. Commuiication-effceenn algorithms for distributed stochas- tic principal component analysis. arXiv preprint arXiu: (tl:1772.0869)) 2017.
Rong Ge, Furong Huang, Chi Jin, and Yang Yuan. Escaping from saddle points Dointe-online stochastic gradient for tensor decomposition. In Proc. Conf. Learning Theory (COLT), pp. 797 842, 2015.
Noah Golmant, Nikita Vemuri, Zhewei Yao, Vladimir Feinberg, Amir Gholami, Kai Rothauge, Michael W Mahoney, and Joseph Gonzalez. On the computational inefficiency of large batch sizes for stochastic gradient descent. arXiv preprint arXim:111.112441 2018.
Gene H. Golub and Charles F. Van Loan. Matriz computations. Johns Hopkins University Press, Baltimore, MD, third edition, 2012.
Priya Goyal, Piotr Dollár, Ross Girshick, Pieter Noordhuis, Lukasz Wesolowski, Aapo Kyrola, Andrew Tulloch, Yangqing Jia, and Kaiming He. Accurate, large minibatch SGD: Training ImageNet in hour. arXiv preprint Kxfx:7006.2000 2017.
Nathan Halko, Per-Gunnal Martinsson, and Joel A Tropp. Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions. STAM review, 53(2):217-288, 2011.
Moritz Hardt and Eric Price. The noisy power method: A meta algorithm with applications. In Proc. Advances Newral Inform. Process. Systs (NewrIPS), pp. 2861 2869, 2014.
Elad Hazan, Kfir Yehuda Levy, and Shai Shalev-Shwartz. On graduated optimization for stochastic non- convex problems. In Proc. Int. Conf. Mach. Learn. (ICML), pp. 1843-1841, 2016.
Elad Hazan, Satyen Kale, and Shai Shalev-Shwartz. Near-optimal algorithms for online matrix prediction. SIAM J. Computing, 46(2):744-773, 2017.
Amelia Henriksen and Rachel Ward. AdaOja: Adaptive learning rates for streaming PCA. arXiv preprint arXiv:190..0 121 15, 2019.
29

---

Published in Transactions on Machine Learning Research (10/2022)
Prateck Jain, Chi Jin, Sham M Kakade, Pranceth Netrapalli, and Aaron Sidford. Streaming PCA: Matching matrix Bernstein and near-optimal finite sample guarantees for Oja's ija's algorithm. In Proc. Conj. Learning Theory (COLT), pp. 1147-1111 1164, 2016.
Zohar Karnin and Edo Liberty. Online PCA with spectral bounds. In Proc. Conf. Learning Theory (COLT): pp. 1129 1140, 2015.
Soila Kavulya, Jiaqi Tan, Rajeev Gandhi, and Priya Narasimhan. An analysis of traces from a production MapReduce cluster. In Proc. 10th IEEE/ACM Int. Conf. Cluster, Cloud and Grid Comput., pp. 94 103, 2010.
David Kempe, Frank McSherry, Kempe David, and Frank McSherry. A decentralized algorithm for spectral analysis. J. Comput. and Syst. Sci., 74(1):70-83 83, 2008.
Usman A Khan, Soummya Kar, and José MF Moura. Distributed sensor localization in random environments using minimal number of anchor nodes. IEEE Trans. Signal Process., 57(5):2000- 2016, 2009.
Satish Babu Korada, Andrea Montanari, and Sewoong Oh. Gossip PCA. ACM SIGMETRICS Performance Evaluation Review, (9(1):169-110) 2011.
Wojeiech Kotlowski and Gergely Neu. Bandit principal component analysis. Tn Proc. Conf. Learning Theory. (COLT), 2019.
TP Krasulina. The method of stochastic approximation for the determination of the least eigenvalue of a symmetrical matrix. USSR Comput. Mathematics and Mathematical Physics, 9(6):189 195, 1969.
Yann LeCun. The MNIST database of handwritten digits. http://yanm. lecun. nom//exb//mnist/ /mnist/ 1998.
Cong Leng, Jiaxiang Wu, Jian Cheng, Xiao Bai, and Hanqing Lu. Online sketching hashing, In Proc. IEEE Conf Comput. Vision and Pattern Recognition (CVPR), pp. 2503-2511, 2015.
Chun- Li, Hsuan-Tien Lin, and Chi-Jen Lu. Rivalry of two families of algorithms for memory-restricted strcaming PCA. In Proc. Int. Conf. Artificial Intell. and Statis. (AISTATS), pp. 473-481, 2016.
Lin Li, Anna Scaglione, and Jonathan H Manton. Distributed principal subspace estimation wireless sensor networks. IEEE J. Sel. Topics Signal Process., 5(4):725-738, 2011.
Zhize Li, Hongyan Bao, Xiangliang Zhang, and Peter Richtárik. PAGE: A simple and optimal probabilistic gradient estimator for nonconvex optimization. In Int. Conf. Mach. Learn., pp. 6286-6295, PMLR, 2021.
Edo Liberty. Simple and deterministic matrix sketching. In Proc. 19th ACM SIGKDD Int. Conf. Knowledge discovery and data mining, pp. 581-588. ACM, 2013.
Jelf Linderoth, Alexander Shapiro, and Stephen Wright. The empirical behavior of sampling methods for stochastic programming. Ann. Operations Research, 142(1):215-24) 2441, 2006.
Teodor Vanislavov Marinov, Poorya Mianjy, and Raman Arora. Streaming principal component. analysis noisy settings. In Proc. 35th Int. Conf. Mach. Learning (ICML), pp. 3410-3419, 2018.
Aryan Mokhtari, Alec Koppel, Gesualdo Scutari, and Alejandro Ribeiro. Large-scale nonconvex stochastic optimization by doubly stochastic successive convex approximation. In JEEE Int. Conf. Acoustics, Speech and Signal Process., pp. 4701 4705, 2017.
Jiazhong Nic, Wojciech Kotlowski, and Manfred K Warmuth. Online PCA with optimal regret. J. Mach. Learning Res. (JMLR), 17(173):1 14:1,11-9999 2016.
Erkki Oja and Juha Karhunen. On stochastic approximation of the eigenvectors and eigenvalues of the expeclation of a random matrix. J. Math. Anal. and Applicul., 106(1):68 84, 1985.
Muhammad I Qureshi, Ran Xin, Soummya Kar, and Usman A Khan. Push-SAGA: A decentralized stochastic algorithm with variance reduction over directed graphs. IEFE Control Syst. Letters, 1207, 2021.
30

---

Published in Transactions on Machine Learning Research (10/2022)
Benjamin Recht, Christopher Re, Stephen Wright, and Feng Niu. Hogwild: A lock-free approach to par- allelizing stochastic gradient descent. In Proc. Advances Newral Inform. Process. Systs (NeurIPS), pp. 693 701, 2011.
Sashank J Reddi, Ahmed Hefny, Suvrit Sra, Barnabas Poczos, and Alex Smola. Stochastic variance reduction for nonconvex optimization. In Proc. Int. Conf. Mach. learning (ICML), pp. 314-323, 2016a.
Sashank Reddi, Suvrit Sra, Barnabás Póezos, and Alex Smola. Fast stochastic methods for nonsmooth nonconvex optimization. In Proc. Advances Neural Inform. Process. Syst. (NewrIPS), 2016b.
Sashank J Reddi, Suvrit Sra, Barnabás Póczos, and Alex Smola. Fast incremental method for smooth nonconvex optimization. In IEEE 55th Conf. Decision and Control (CDC), pp. 1971- 1977-197.. IEEE, 2016c.
Herbert Robbins and Sutton Monro. A stochastic approximation method. Ann. Math. Stat., pp. 400 407, 1951.
Sebastian Ruder. An overview of gradient descent optimization algorithms. arXiv preprint arXiv: 1609.04777, 2016.
Terence Sanger. Optimal unsupervised learning single-layer linear feedforward neural network. Neural networks, 2(6):459 473, 1989.
Ohad Shamir. A stochastic PCA algorithm with an exponential convergence rate. In Proc. Int. Conf. Mach. Learning (ICML), 2015.
Ohad Shamir. Fast stochastic algorithms for SVD and PCA: Convergence properties and convexity. In Proc. 33rd Int. Conf. Mach. Learning, 2016.
Ohad Shamir and Nathan Srebro. Distributed stochastic optimization and learning. In 52nd Annu. Allerton Conf. Commun., Control, and Compute., pp. 850-857. IEEE, 2014.
Alexander Shapiro and Tito Homem-de Mello. On the rate of convergence of optimal solutions of Monte Carlo approximations of stochastic programs. SIAM J. oplimization, 11(1):70 86, 2000.
Cheng Tang- Exponentially convergent stochastic k-PCA without variance reduction. arXiv preprinl arXiv:190,.01750, 2019.
Joel A Tropp, Alp Yurtsever, Madeleine Udell, and Volkan Cevher. Streaming low-rank matrix approximation with an application to scientific simulation. SIAM J. Scientific Comput., 41(4):A2430-12463, 2019.
Manfred K Warmuth and Dima Kuzmin. Randomized PCA algorithms with regret bounds that are logarith- mic in the dimension. Proc. Advances Neural Inform. Process. Syst. (NIPS), volume 19, pp. 1481 1488, 2007. ISBN 1956282 doi: 10.1.1.133.8332.
David P Woodruff. Sketching tool for numerical linear algebra. Foundations and TrendsB in Theoretical Computer Science, 10(1-2):1-157, 2014.
Sissi Xiaoxiao Wu, Hoi-To Wai, Anna Seaglione, and Neil A Jacklin. The Power-Oja method for decentralized subspace estimation/teee ttracking. In Proc. 2017 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pp. 3524 3528. TEFE, 2017.
Sissi Xiaoxiao WU, Hoi-To Wai, Lin Li, and Anna Scaglione. A review of distributed algorithms for principal component analysis. Proc. IEEE, 106(8):1321 1340, 2018.
Puyudi Yang, Cho-Jui Hsieh, and Jane-Ling Wang. History PCA: A new algorithm for streaming PCA. arXiv preprint arXiv: Xxtx:1022. 1ee.1802.0050 05447, 2018.
Se- GouYo Yun, Alexandre Proutiere, et al. Fast and memory optimal low-rank matrix approximation. Advances in Newral Information Processing Systems, 28, 2015.
31

---

Published in Transactions on Machine Learning Research (10/2022)
Dejiao Zhang and Laura Balzano. Global convergence of Grassmannian gradient descent algorithm for subspace estimation. In Proc. Inl. Workshop Artificial Intell. and Statist. (AISTATS), pp- 1460-1468. 2016.
Siyun Zhou and Yanqin Bai. Convergence analysis of Oja's 's iteration for solving online PCA with nonzero- mean samples. Science China Mathematics, 64(4):849 868, 2021.
Appendix A Proofs of Lemmas for the Initial Epoch
A.1 Proof of Lemma |
In order to prove Lemma we first need the following result.
Lemma 9. The second moment of the update vector & in D-Krasulina is upper bounded as
1
f [6
E
2
|
Proof. We start by writing E {llE:
[B( in terms of E [1) follows:
E flll (1))))1))
c E + (B((.)) (E:) {( (E E
E (12) B(E HE{:).
Now defining C, :
and rearranging the above equation, get
(11))1)) = Bill, B(5:)112) + Cr.
Next, substituting value of & from 8
wC get
C2
B
HB
(B{}
1
v -v--Ve-Ve-- -
VT vi-1 -----i- Evi-1-
Ev--1
V
v V
v.T
V
V'v
v Ev-1
(32)
Since Σ is a positive semi-definite matrix, we write its eigenvalue decomposition as == E \ where A > X2 2- 2 Aa 2 0 and qı q"), .2...., qa are the eigenvalues and corresponding eigenvectors of E, respectively. It follows that
C2
)2
d
A
Vi-1111
12
mm1
iml
d
(v
(t=19)) )2
A1 (v 9²)2
d
(vi- q.)
2
A
1iv_1112
12
V 112
V 12
||
1=2
i=2
(vF
d
(v.T
x² (v. 7")
X
+ 11
Σ
1|ve_111
|| 11
||ve-1 11
1=2
X2 (v] "T"1"")
1 (v )²)
(vī,q''²
+ +
1
|v-1111
| 11
| 11
Finally, we get from definition of 4,-1 that
C.
< V1-1((1-1 V.-1) -3)
< SV (X + >2)) S 2277v.11
V
This completes the proof of the lemma.
32

---

Published in Transactions on Machine Learning Research (10/2022)
Using Lemma
9
can now prove Lemma in the following.
" Substituting V, from (8, we get
Proof of Lemma
From
10 we have V = livel
Vi-1 + 7 13 ((v--1 + -xE))= (o) 1 12 +
1/
((v=-1 + met) 9²)²)
1
[val2
1iv.1
(6) Ve-1 1 1/ + 22115: 13 ((V,-1 + 7x51) 99)² 1+
1
((Ve-1 + 7(,)) '')²)
| 11
V/_
V 11
| 12
(v 1" )
+ (57 + 2%
9- (E) q*)
1 + 7/
V| 11
V/-21112
(vī-1q' )²
|
(E. )2
(v (v7-19')
=
1|ve-1112
| 1
27
V 1
||
(V 1 ($F 4')
v 7i
2%/e
(33)
|
12
Here (a) and (6) are due to (Balsubramani et al. perpendicular to &, and (b) is true because [ 1 be bounded as
2013 Lemma A.1), where (a) is true because V is
I Iv I12. The second term in the above incquality can
& 3 V:
E E + BBE))
(c)
E (Il.
B(CH)11)
I
+
1/
12
N
N
1
1
E
A
V
A
Vc-
V
2
N
1=1
1
N
1=1
-E
E
A
N
2
N
N
N
1
1
E
A
A Ac.
N
V2-111
N 1=1
1m1
113
-vi-1t+ +
| -v-
N
NV
V1-
E
Aii - Σ 2)
V
1
12
V--1
T
Ait
E v2)
i=1
N Σ
N
2
1
<4
A
+ 1
A -
E
+ 22
(34)
N
1=1
where (c) is due to Lemma g Substituting we prove Part (it) of the lemma by defining V V1-1 1-l1/lle-lll and noting that
34
in
33
completes the proof of Part (i) of Lemma 1 Next,
&& 1
II(1/)) Σ N E 12
. (L/N2)]
E
V.-1 12
[1 12
ING
(1/N2)
(x vc-1 |
NI
1xx, x
xx,113(()xx,1111 - (x
N 1=1
x l
I max |Xs,ell s4¹
(35)
33

---

Published in Transactions on Machine Learning Research (10/2022)
Here, (d) is by using Saunch-Sccwartzz inquality and the last inequality is due to Assumption [A1]. Now substituting this in 33 completes the proof.
A.2 Proof of Lemma 2
We begin by writing
B[(a - [(2))) B
2% 1"(((=1'")(E7
2
2
l| (
Fr-1]
1
2
(v 9*)²
<.- - B(Erar))
Vc- 1
Substituting value of & in this, we get
4n )2)
N
v. 1 EN
T
[[ B(x)))
E
A
N
"""
Vi-1
N
Vi-1 b
i=1
N
VIIN 1 DN A T
2
-
N
A
V
1=1
4n T 92)2
T 1 E A
T )²
V
IE
A N
Vt-1
91
-=1
N
v VILVE
N Di= Ve-1 99
2
- -V E Σ
AueJa +
Since E 1
N Di=1 Arx] is the covariance matrix E, we gel
E ((zx)E((z)) 1}
1
N
G N
Art - E)v.
T 9
I
A
-)vv11
V/ 18
N
Vt-1
i=1
4v7(vT 19
E
(qTT(
Σ
N
Aur ΣjvA1
VI Σ N
i=1
Aii E) qra V
1 11
| 12
87
N
v ( 1=
Ar, - Eivv-1) 9tT VC-1
VC-1
E.
T
A EjvC-1
+
V 11
N 1=1
1ive-111
87'(v2"
q"' TTT
Li=1
N
Ai,
Σ)ve-1
V 2i=1
N
Ait Σ)v(-1
q T V--1 llve_ll2
2
+
1/ve-1112
|/-21112
V 1
9 (
N A - V)ve-1 (N Li=1
2
T
EN
Ait E)vi-1
S {
(36)
V
2
36 due to the fact that q'Tvi- l|ve-1l1z
We can see that both the remaining terms in 366 are Rayleigh quotients of matrix (E- N EN Ai,)) and hence the largest eigenvalue
where the last incquality in
S 1.
of (Σ- - + >N> A.,)) maximizes both the terms. Using this fact we get
N
N Σ
1
$ - B(( - ([(==E((z)))) ((z 16,
Σ
A S 16
A
N 1=1
N
1=1
Using 12
we get E[((( ((2))) < 20 which completes the proof.
34

---

Published in Transactions on Machine Learning Research (10/2022)
A.3
Proof of Lemma
Using Lemma can write the moment generating function of Mr as follows:
[(ex s B exp (sV,-1 sn2r² 82c Fr-1] exp(sV sy2rr)) exp (- 824 Fi-1
exp(sV + s72r* 21 ((-1))) exp (-s(2) ( Fr-1)
(37)
We can bound this using Bennett's inequality (Proposition in Appendix which requires the variance and range the random variable žt. We have already computed the variance 2t in Lemma 2 Next we compute the boundedness of (2t (2+)) follows:
- E{4} $ 2/2 - 2% | & 2 =: h.
(38)
Here, the last inequality is due to Assumption get
A1
Using parameters ofy and h with Bennett's inequality,
B(exp(sV.)11--11 S exp sV F sggrr 222003 ( eoh -1 s/
(39)
(sh)²
For L 2 In + L2, where Ln and L2 are given by 13 we show in Proposition in Appendix
that
1tth 11-sh) < 1 for S. This implies
B(exp(>V:))Fr-13 S exp sV -B F--) + sygr + 1
which completes the proof of the lemma.
Appendix B
Proofs of Lemmas for the Intermediate Epoch
B.1
Proof of Lemma
Using Lemma
3 we have
8 1/ + 7/2 - F( + s7 ))
W. Fo-13
S exp
8 V.-1 - 27p (>1 --x2) Vz- 11
(1--(1-111)) V.-1)
722 + 07/0%)
exp
(6)
co 1o (1 11
(1-V (o-11
c2p4
+ sea
exp
t 1+1
(1 L)
(1 Lyz
(40)
Here, (a) is due (Balsubramani al. 2 (t + L). Finally, for w E 22 have V (w) < 1 Cj. Now taking expectation over Sj we get the desired result.
2013 Lemma A.3) and (b) is by substituting 7r c/(t L)
B.2 Proof of Lemma 5
and (e(s):== = scrrd 2003 Substituting o and (e(s) in Lemma 4 we get
Define 0g := 1 Cot; z4L
(141)2 + (141)²
E, ( W < E, Le (e exp (S{(s) I E, f exp (Ce(s)).
(41)
35

---

Published in Transactions on Machine Learning Research (10/2022)
Note that the second inequality in 41 is due to Baalubraaani et al. 2013 Lemma 2.8). Applying this procedure repeatedly yields
B 0. < +1 exp (sVi, an ... 01,+)) exp (C(+)) ... exp (Ge++1 1 (.........+)))))
S Be,+t exp (sVi, ar................. 0ty+1) }) exp (G)(s) ... exp (Ge+ (+)).
Substituting values of ar and (e()) in the above, we gel
By ( - B exp (V (1 ++I
CoCs
(1
CoC
tj L
exp (s02)* + +22200%)
1
1
+ (1+++112))
(t+ (t+
$ exp (s(1 - -5) exp (-coe) c0cs 7+ +...
ty +1,++11)
1
exp (10c2r + 2²c²088
1
1
+ +
(42)
(t (++L)* L)2
+ + L + 1)3
Here, the last inequality is true because Ve, (w) 5 for w E Phytt and - x - for x S 1. Next we bound the summations
42 as follows:
[+1 dx = In 1+++1
+ L
t L +1
14+1
+
+ L+1
1
1 (1, L 1)2
dx
1
1
(t L)2
(x+ L)2
1 +
t+
Substituting these bounds in
42
we get the desired result.
B.3 Proof of Lemma 6
This lemma uses Lemma and deals with specific value of t -1+11 For t == t+t11 29 gives
Ex,+ - exp s(1 e,)
t, + +++++++1
+ sc²pp + 2008 -
1 + G+++ ==)+z) (43)
Using conditions
C1 and
C2 and the fact that -2x <1 - for 0 S 5 3/44, we get
t, + L +1 1 00es (e (/ - be, < 1 3, < 1 Ej+1 E3
)
(1 (1-
1i+1
L+11
Substituting this in
43 we obtain the desired result.
B.4 Proof of Lemma
Constructing supermartingale sequence M, in the same way as did in Theorem 3 for E S and applying Doob's martingale inequality, we get.
B(M,} Pt sup v 2 1- -e) s Pt, sup Mr 2 ed(2-s)) s (0(1) 126,
121,
Ef exp (sVi, + sTC, }} (0(())
Ef exp (sW,, , } exp (STr,) e
Using Lemma t then results in
Pr, sup V, > 1 e3) 12g
(x(1-G) exp
1
s(1 - ej) - -sEj-1 + +500)
1 + L7 1T,
1
1
+ tsT,
36

---

Published in Transactions on Machine Learning Research (10/2022)
Substituting a bound Tt from Theorem 3 (sce, e.g., the discussion around 28 get
Pt, sup Wr 2 1 - e3) s exp 126
st s2c2 + S2
( 1 4 (,+L) + s s + scoo2/ 11+
- 8 + s s( se sc
1
exp
1--1 +
Substituting 8 = (2/co) In (4/8) and using the lower bound on L, we get (sec Proposition 5 in Appendix for formal verification)
53-1/50
Pt,
sup 12g
v, 2,
12-16, es) S exp
(- =
8Ej-1
8
I 2
2
4
Summing over j completes the proof of the lemma.
Appendix C Proofs for the Final Epoch
Proof of Lemma
8 From Lemma
Part (1), we have
N Σ
V. S v- + 47
A.,
+ 24-
N
2=1
Taking expectation conditioned on Fr-1, get
E [(V.|Fo-1) & 2))))
-E Fr-1)
where the second term is due to Lemma g
Now using upper bound on - [ from Balsubramani
et al. 2013 Lemma A.4), we get the following:
V < 1 + 4nio -74(1 12) v1-1)
Vー- ?² 27x(A1 -2) Vo-1)) + 47,0
Finally, taking expectation over Di substituting " co/(2(t L) de)), and using the facts that 22 is meessuuuee and for L > L,, V S 1/2 and we lie in sample space S2 with probability greater than we obtain
1 1-8 8 (Theorem
4c
V.V. < E,
V V [1
2(1 + ())
2²)2
2(1 + L)
1)2
4c203
1 1
EK +
2(1+ + L) 11)²
L)
(1+ L)²
2
4c²03
ET- ( V (V,-1) +
2(1 (.) )2)²
2(1 + T.)
(1 T.)2
This completes the proof of the lemma.
Proposition 1. Let a1, b > 0 and be some constants. Consider nonnegalive sequence (u+ : 1 > (,)
that satisfies
u (1
01
02
b
(t + L)²
t L
(t+ L)2
Then we have:
L + 1
a1
a1
2-2
W1 <
[+L+]
exp
20
+
(1+ L
+1
exp
1+11
1
02 -1
37

---

Published in Transactions on Machine Learning Research (10/2022)
Proof. Recursive application of the bound on ue gives:
)
ut S
II (1 1=ty+1
Q1
02 uso Σ
II
[1+
a1
J+1
02
(44)
(i+
1)3
i
i L
(i+ L)²
(++
i=ty+1
j=i+1
Using (aalubrramani et al.
2013 Lemma D.1) we can bound the product terms as
II (1+ j=+t
an
exp Σ
02 j+L
L)²
j+ L
++1++11
exp Σ
a1 (j L)²
(45)
++L+1
j=1
Next, we bound the last term here as
(Lr (xr +
exp Σ
G
S exp
dx
eep [++++1
01
S eep
a1
L)2
1.4 L
i L
Substituting this in
45 we get
II
(1+
Q1 L)²
02 j+ L
i ++
a1
cxp i L +1
++ L+++1|
-=++1
Substituting this in 44
we get
+ L +1
a1
Uts + Σ
b
II
(1
@1
02
44 5
t L +1
exp
t + L
(i+ L)²
[j + L)²
j+ L
i=t++
s=
t + L +1 1
02
a1
b
i i++++
42
a1
exp
ues +
exp
1+ 1+1
(
L
(i+ 1=(++1
1+L+1
i+ ++1+1
t 1+1
a1
01
[++L+1)
1+1+1
exp
f + L+1
Ues + +exp
t + L + 1/
(t+ L + 12
L)
L+1
02
exp
a1
Ue, +exp
Q1
b
L 2
Σ
(i + L++
1+1+1
L + +++
(1 +1+ 1
i=1
Again applying Balsubramani al. 2013
Lemma D.1), we get the final result as follows
42 + tu
b
L
+ 2
2 + L + 1)a2
exp
Ue + exp 1 +1
1+L+1
L
(1 L 1)a2
L+1
a2
[++++1
42
1 (++ L+
a1
2
2 b
1+
exp
f + ++1
UtJ
exp
1, + I
[++1
02
This completes the proof of the proposition.
Appendix D Other Auxiliary Results
Proposition 2 (Bennett's Incquality (Boucheron et al. 2013)) ) Consider a zero-mean, bounded random
variable X; IR (ie., [Xd] < h abmost sureby) with variance of. Then for any 8 R, we have
Ex < exp 0, 8 1- sh
(sh)²
38

---

Published in Transactions on Machine Learning Research (10/2022)
Proposition 3. Let h == 2 and {d/Ae, (2(4))1((0)))) It then follows that 8th -th < 1.
[ah)²
long sh S 7/4. Therefore, in order to prove this proposition, it suffices to show that the lower bound on L implics sh S 7/4 for {d/4c, (2/co) In(4/5))
Proof. It straightforward lo that
-h sh (sh)a
We establish this claim as two separate cases for the two values of 8.
Case 1: For s = d/46, /e, substituting the value of h gives us
sh = d dcr²
der²
dcr²
der²
C 2e 8dr4 max(1, In(4//))
S 1
I
S 4
2e 2(t Liee 2LE
21L1
16
Case II: For = (2/co) In In(4/),, we obtain
sh colt + L)
2l1(4/0)cr² 21n(4/6)er2 21n(4/8)er
e0
$ 5
1 7
coL1
cD
8r max(1,,))
This completes the proof of the proposition.
Proposition 4. Assuming L 2 8dr max(1.e²)
In + 8d max(1.e")
In 1 and the porameter = d/4e, we have
2
I (r + 80 S-
Proof. We prove this by proving the following two statements:
sc'o 2
scoo2
S
and
<
L1
2d
L2
We start by proving the first statement:
d'r
c
20 Next, prove the second
L
Sdi max
Zook < d
statement as follows:
I
This completes the proof.
8d (1,2) &
2d
Proposition 5. For 8r max(1,e²)
In 1 +
80 max(1,
In 4 3: we have
(i) (1-1+L) S 4 and
20 (ii 1+ In 1 20.
Proof. We begin by noting that
2
2
202r4
(t +)) S
S
2
c0
$
8
L1
8r
4
Next we prove the second statement as follows:
20
In
4
20
In
<
20002
20
N
In
2c2
In
2 8 s max c2)
-
80 4
co L)
1a 8
cgLe
eo
This completes the proof of the proposition.
Corollary 3 (Restatement of the Main Result (Theorem 1 Convergence in Probability). Fir any 8 e (0, 1)
and pick c = co/2(A1 - Xx2 for any oo > 2. Next, define
L = 44 max(1, d²)
Le o max max (1, 4
In J.
In 81
82
5
39

---

Published in Transactions on Machine Learning Research (10/2022)
pick any L L1 L2, and choose the step-size sequence n := (L+)), Then, long as Assumptions A1 hold, we have for D Krasulina that there caisls sequence (2t)sEz, of nested sample spaces such
and [22]
that P (Tisoll) 2 21-6 8 and
P v. ci ci
P
++
L 1++1)
+ +
+C²
++L+1
oN
s 26,
(46)
where Cj and C2 are constants defined as
8c2
G
1
4ed
2
x2LL
and C₂
28
82
8(co 2)
Proof. The proof follows by first using Markov's inequality to bound the conditional probability
Pr v. 2 ci (
+ + + c (+.+)))
then using the bound on E, [ in Theorem to further bound the conditional probability by 6, and finally removing the conditioning on the nested sample spaces by using the facts that (i) P(A) S P(A|B) + P(B)) for any two events A and B, and (ii) IP (U2>ol4 s 6.
40


Reviews:
[{"review_id": "sOq3vgtg2T", "reviewer": "Reviewer_Ly9o", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "## Strengths\n* The paper is generally well-written and easy to follow\n* The related work is fairly comprehensive, nicely framing the contribution of this paper\n* Theoretical guarantees are rigorous and establish strong performance bounds for the proposed approach\n\n## Weaknesses\n* The system timing model considered in the paper has some limitations. The processing and communication times $R_p$ and $R_c$ implicitly depend on the problem dimension $d$. The communication time $R_c$ also depends very explicitly on $N$, the number of processors aggregating results. Although this is briefly mentioned in Sec 2.2.1 when $R_c$ is introduced, this brings into question the way that it is used to obtain conditions on $N$ relative to $R_s$, $R_p$, $R_c$, and $b$ throughout the paper.\n* A few older relevant related works are missing (see below).\n* Some aspects of how the main theoretical results are presented could be slightly improved/clarified (see below).\n* Setting a schedule to satisfy the conditions of Theorem 4 requires knowing the gap $\\lambda_1 - \\lambda_2$, among other problem-dependent constants, which may not be possible in some applications.\n* The introduction motivates the problem considered by discussing challenges with high-dimensional data, while the experiments only consider problems with dimension at most 784. The paper would be stronger if the method were demonstrated on much higher-dimensional datasets.\n* The experimental section does not include comparisons to any other method from the literature.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper studies a decentralized version of Krasulina's method for streaming PCA (i.e., eigenvector estimation), with and without mini-batching. The paper provides conditions under which D-Krasulina achieves a linear speedup in the number of distributed workers. In the mini-batch setting, conditions are also provided where a linear speedup is achieved. For the extreme case where, even with mini-batching, the system of N workers is not able to process samples  faster than they are arriving, further analysis quantifies the rate of convergence when some samples are dropped, with an assumption that all samples are iid. In this case, DM-Krasulina still achieves near-optimal convergence rates as long as the fraction of dropped samples is not too large. These results are validated with numerical experiments.", "claims_and_evidence_yn": "", "claims_explanation": "", "audience_yn": "", "audience_explanation": "", "requested_changes": "## Requested changes\nThese are the main points I view as critical.\n* Please clarify/justify the timing model used in the paper, given that $R_c$ depends on $N$ and that this is also used to determine a bound on the number of processors $N$.\n* While the related work section covers most relevant related work, it would be good to also discuss the relationship to [Yun et al., \"Fast and memory optimal low-rank matrix approximation,\" NeurIPS 2015.](https://proceedings.neurips.cc/paper/2015/hash/21be9a4bd4f81549a9d1d241981cec3c-Abstract.html)\n* In the statement of Corollary 1, please clarify whether the event occurring with probability $1 - \\delta$ involves the intersection over $0 < t \\le T_B$, or all $t > 0$ (i.e., the probability $1-\\delta$ is still only achieved asymptotically).\n* Setting a schedule to satisfy the conditions of Theorem 4 requires knowing the gap $\\lambda_1 - \\lambda_2$, among other problem-dependent constants. Please add a discussion of how one might approach this in practice when applying the method.\n* One typical characteristic of streaming data settings (in practice) is that the data distribution may not be iid; rather over time, there may be some distribution shift. It would be good to add some discussion about this iid assumption and ways it may be relaxed in future work.\n\n## Additional questions\n* Given the relationship of Krasulina's method to SGD, might one expect to achieve a superior method by incorporating some sort of momentum or acceleration term?", "broader_impact_concerns": "None noted", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "gLAwRP8Af7", "reviewer": "Reviewer_rnHX", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "**Strengths**: The proof seems to be correct and the paper has provided a good discussion for the proof technique used (which is a key contribution). Overall, the paper has the following strengths:\n\n1. It provides a tighter convergence analysis than (Balsubramani et al., 2013) and show that the asymptotic convergence rate of the Krasulina's algorithm is $O(\\sigma_N^2/t)$ rather than $O(r^2/t)$ as in (Balsubramani et al., 2013), where $t$ is the iteration number. As such, the result shows that the distributed Krasulina's algorithm is able to take advantage of parallel processing to speed up computation of PCA since the simple aggregation of samples allows one to reduce variance of the estimate. \n\n2. It analyzed the effects of skipping data in the multi-batch Krasulina method and show that as long as the skipping data rate does not exceed $\\mu = O(B)$, then the optimal asymptotic convergence rate of $O(1/T)$ can be maintained. \n\n**Weaknesses**: There are some weaknesses with the current submission as listed below \n\n1. The algorithms analyzed are standard and can be considered as the direct extension of classical Krasulina method to the distributed setting. Moreover, important aspects to distributed algorithms such as communication efficiency of the distributed algorithms has not been addressed. In this regard, the paper seems to be better positioned as presenting an improved analysis of Krasulina over (Balsubramani et al., 2013).\n\n2. Compared to prior work (Jain et al. 2016) which has a variance dependent bound unlike (Balsubramani et al., 2013), one of the main claims in this paper is that it provides a high probability bound ($1-\\delta$) to arbitrary $\\delta$. However, the claim is slightly unfair: \n\n- (a) as mentioned in (Jain et al. 2016), it is possible to run $O( log(1/\\delta) )$ of their algorithms to obtain convergence w.p. at least $1-\\delta$; \n\n- (b) to ensure convergence for the Krasulina algorithm in this paper with probability at least $1-\\delta$, it requires setting the step size parameter to be $L \\asymp \\ln(1/\\delta)$ - in other words, to ensure convergence with an arbitrarily high probability, it is required to set an arbitrarily small step size. \n\n- Moreover, it should be noted that the derived convergence results are **expected** convergence rate which holds with high probability, while (Jain et al. 2016) derived a high probability convergence rate. ", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper studies the PCA problem with streaming data. It focuses on a setting where the streaming data arrives at a higher rate than the processing power of a single processor. Therefore it utilizes multiple processors that communicate via a server-worker structure or a decentralized general network structure to process the streaming data in parallel. Theoretical analysis of the Krasulina method-based algorithms in the distributed setting are presented, and they show the convergence rate of $O(1/T)$ when the number of samples observed is $T$.", "claims_and_evidence_yn": "", "claims_explanation": "", "audience_yn": "", "audience_explanation": "", "requested_changes": "The paper is in general OK but the authors should address the issues mentioned in \"Weaknesses\" and make corresponding revisions. There are some additional (but perhaps minor) issues:\n\n1. In (10), the random variable $\\Psi(w)$ is not defined. \n2. In Lemma 8, what is the constant $c$?\n3. It seems strange to impose a new definition in Definition 5 for the variance of distributed samples, while it is clear that the variances are related by $\\sigma_N^2 = \\sigma^2/N$. \n4. At the end of p.17, the upper bound should be $M_{0}$ instead of $M_{t_0}$?", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "EM3IVvPWOi", "reviewer": "Reviewer_Z1FS", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "The paper's writing takes great care to discuss extensive related work, explain the computational model and when and why they are appropriate, and to provide intuition for how the analysis fits together and the main ideas.  \n\nThe topic is well-motivated, as dealing with massive data is an important challenge and the streaming model considered is aligned with how such data is dealt with.  Moreover, estimating the top singular vector is an important challenge in data analysis.  \nHence, this paper is appropriate for TMLR, and I recommend acceptance.  \n\n\nGiven that, I have a couple critiques. \n\nThe experimental section does a fine job of illustrating how well the bounds and analysis of parameters in the data assumptions and algorithms align with empirical performance.  However, the data sets are significantly smaller than the sizes needed to motivate the sort of model studied.  \nNevertheless, I think this is fine, since one needs to verify the accuracy, and running on such motivating sizes would require extensive computing cost.  Also the main claims are still well illustrated.  There are caveats that other unforeseen issues may occur when scaling to sizes discussed in the paper, such as node failure or larger network cost than expected.  However, I rate this concern as minor, and do not think the authors need to extend the empirical study in this way.  \n\nSecond, the paper uses the term PCA, and finding the top principal component, to mean finding the vector with maximum variance.  However, it assumes the data has mean zero.  Typically PCA is more general, and does not require the data to have zero mean, and often the process is to center the data before finding the top singular vector.  As such, the problem the paper actually solves is finding the top singular vector, which does not require centering.  \nThat said the paper does not hide this assumption, and this assumed equivalence between PCA and SVD (by assuming the data is centered) is not uncommon in this sort of analysis.  However, I would appreciate if the authors spent a bit more time discussing this distinction, and the challenges in addressing the full case where the data cannot be assumed to be centered.  \n\n ", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper looks at computing the top singular vector from a set of iid vector samples arriving at a very high rate.  In particular, it is concerned with when the rate is so high that directly adapting the updates via Krasulina's (or Oja's) algorithm is not possible.  The main result of the paper is that this computation can be spread across distributed machines, and the updates averaged still achieve the optimal linear convergence rate.  Batch variants and when a limited amount of data must be dropped, even in the distributed setting, are also analyzed with similar results.  \n\nThe main technical idea is revising existing results showing the convergence based on worst case bound on each update to one that achieves the same linear convergence rate using a variance bound.  This analysis can then be adapted to distributed and batch analysis where averaging estimates across sets decreases variance at the desired linear rate.  \n\nThe paper's writing takes great care to discuss extensive related work, explain the computational model and when and why they are appropriate, and to provide intuition for how the analysis fits together and the main ideas.  \n\nThe topic is well-motivated, as dealing with massive data is an important challenge and the streaming model considered is aligned with how such data is dealt with.  Moreover, estimating the top singular vector is an important challenge in data analysis.  \nHence, this paper is appropriate for TMLR, and I recommend acceptance.  \n\n\nGiven that, I have a couple critiques. \n\nThe experimental section does a fine job of illustrating how well the bounds and analysis of parameters in the data assumptions and algorithms align with empirical performance.  However, the data sets are significantly smaller than the sizes needed to motivate the sort of model studied.  \nNevertheless, I think this is fine, since one needs to verify the accuracy, and running on such motivating sizes would require extensive computing cost.  Also the main claims are still well illustrated.  There are caveats that other unforeseen issues may occur when scaling to sizes discussed in the paper, such as node failure or larger network cost than expected.  However, I rate this concern as minor, and do not think the authors need to extend the empirical study in this way.  \n\nSecond, the paper uses the term PCA, and finding the top principal component, to mean finding the vector with maximum variance.  However, it assumes the data has mean zero.  Typically PCA is more general, and does not require the data to have zero mean, and often the process is to center the data before finding the top singular vector.  As such, the problem the paper actually solves is finding the top singular vector, which does not require centering.  \nThat said the paper does not hide this assumption, and this assumed equivalence between PCA and SVD (by assuming the data is centered) is not uncommon in this sort of analysis.  However, I would appreciate if the authors spent a bit more time discussing this distinction, and the challenges in addressing the full case where the data cannot be assumed to be centered.  \n\n ", "claims_and_evidence_yn": "", "claims_explanation": "", "audience_yn": "", "audience_explanation": "", "requested_changes": "Please discuss challenges associated with not assuming the data is centered.  ", "broader_impact_concerns": "This is analysis and efficient algorithm design for a very standard and generic technique.  There are no significant broader impact concerns.  ", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: The writing looks nice. The reviewers comments were all addressed, and I'm pleased to recommend Acceptance.
Decision: Accept (poster)