Distributed Principal Component Analysis with Limited Communication
Authors: Foivos Alimisis
Peter Davies
Bart Vandereycken
Dan Alistarh
Venue: NeurIPS 2021 Poster (2021)
Keywords: Principal Component Analysis
Leading Eigenvector
Bit Complexity
Riemannian Optimization
Sphere

Abstract: We study efficient distributed algorithms for the fundamental problem of principal component analysis and leading eigenvector computation on the sphere, when the data are randomly distributed among a set of computational nodes. 
We propose a new quantized variant of Riemannian gradient descent to solve this problem, and prove that the algorithm converges with high probability under a set of necessary spherical-convexity properties. We give bounds on the number of bits transmitted by the algorithm under common initialization schemes, and investigate the dependency on the problem dimension in each case.

Distributed Principal Component Analysis with Limited Communication
Foivos Alimisis Department of Mathematics University of Geneva
Peter Davies
Bart Vandereycken
Department of Computer Science University of Surrey
Department of Mathematics University of Geneva
Dan Alistarh
IST Austria & Neural Magic, Inc.
Abstract
We study efficient distributed algorithms for the fundamental problem of principal component analysis and leading eigenvector computation on the sphere, when the data are randomly distributed among a set computational nodes, We propose a new quantized variant of Riemannian gradient descent to solve this problem, and prove that the algorithm converges with high probability under a set of necessary coherrcal-convenee properties. We give bounds the number of bits transmitted by the algorithm under common initialization schemes, and investigate the dependency the problem dimension in each case.
1 Introduction
The need to process ever increasing quantities of data motivated significant research interest into efficient distributed algorithms for standard tasks in optimization and data analysis; see, e.g., [66 27
пл In this context, one the fundamental metries of interest is communication cost, measured as the number of bits sent and received by the processors throughout computation. Communication- efficient variants are known for several classical tasks and algorithms, from basic ones such as mean estimation (28 20 to variants of gradient descent [19 [2] even higher-order optimization 2 15).
By contrast, much less is known concerning the communication cost for classical problems in data analysis, such as Principal Component Analysis (PCA) [13] [5], spectral clustering (23] In this paper, we will focus on the closely-related problem of computing the leading eigenvector the data covariance matrix, i.e., determining the first principal component, in the setting where d-dimensional data comes from an unknown distribution, and is split among n machines.
Given the fundamental nature of PCA, there is a rich body of work on solving sequential variants, both in the Euclidean domain (24 25 266 and, more recently, based on Riemannian optimization. This is because maximizing the Rayleigh quotient can be re-interpreted in terms of hyper-sphere, leading to an unconstrained Riemannian optimization problem
These references have exclusively focused on the sequential setting, where the data can be stored and processed by single node. However, many settings the large volume of data requires processing by multiple nodes, which has inspired significant work reducing the distribution overheads of PCA, e.g. 4 specifically on obtaining efficient algorithms in terms round complexity.
More precisely, an impressive line of work considered the deterministic setting where the input partitioned arbitrarily between machines, for which both round-efficient and communication-effiices solutions are now known. The general approach is to first perform a careful low-rank approximation of the local data and then communicate iteratively reconstruct faithful
35th Conference on Neural Information Processing Systems (NcurIPS 2021).

---

estimate of the covariance matrix. Due to the worst-case setting, the number of communication rounds scales polynomially in the eigengap ó and in the precision E, making them costly in practice.
Thus, on the practical side, it is common to assume a setting with stochastic rather than worst-case data partitioning. Specifically, Garber et al. 11 consider this setting, and note that standard solutions based on distributed Power Method or Lanczos would require number of communication rounds which increases with the sample size, which is impractical at scale. They circumvent this issue by replacing explicit Power Method iterations with a set of convex optimization problems, which can be solved in a round-efficient manner. Huang and Pan [4] proposed round- cond-efficieet solution, by leveraging the connection to Riemannian optimization. Their proposal relies on a procedure which allows the update of node-local variables via a surrogate gradient, whose deviation from the global gradient should be bounded. Thus, communication could only be used to reach consensus on the local variables among the nodes. Unfortunately, upon close inspection, we have noted significant technical gaps in their analysis, which we discuss in detail in the next section.
The line of work in the stochastic setting focused primarily on latency cost, i.e., the number of com- munication rounds for solving this problem. For highly-dimensional data and/or highly-distributed settings, it is known that bandwidth cost, i.e., the number of bits which need be transmitted per round, can be a significant bottleneck, and an important amount of work deals with it for other classic problems, e.g. complexity of distributed PCA.
28 12 15]. By contrast, much less is known concerning the communication
Contribution. In this paper, we provide a new algorithm for distributed leading eigenvector com- putation, which specifically minimizes the total number of bits sent and received by the nodes, Our approach leverages the connection between PCA and Riemannian optimization; specifically. we propose the first distributed variant of Riemannian gradient descent, which supports quantized communication. We prove convergence of this algorithm utilizing set of necessary convexity -type properties, and bound its communication complexity under different initializations.
At the technical level, our result is based on a new sequential variant of Riemannian gradient descent on the sphere. This may be of general interest, and the algorithm can be generalized to Riemannian manifolds of bounded sectional curvatures using more sophisticated geometric bounds. We strengthen the gradient dominance property proved in 31, by proving that our objective is reakly-quasi-conver in a ball of some minimizer and combining that with a known quadratic-growtt condition. We leverage these properties by adapting results of (7) in the Riemannian setting and choosing carefully the learning-rate in order to guarantee convergence for the single-node version.
This algorithm is specifically
to be amenable to distribution: in the tangent spaces of the sphere, using variant of the lattice quantization technique of [8) which we port from the context of standard gradient descent. We provide a non-trivial analysis for the communication complexity of the resulting algorithm under different initializations. We show that, under random initialization, solution can be reached with total communication that is quadratic in the dimension d, and linear in 1 1/6, where d is the gap between the two biggest cigenvalues; when the initialization done to the leading eigenvector of one of the local covariance matrices, this dependency can be reduced to linear in d.
specifically, we port it the
distributed setting, and add communication-comprees
2 Setting and Related Work
Setting. We assume to be given m total samples coming from some distribution D, organized as global m x d data matrix M, which is partitioned row-wise among n processors, with node i being assigned the matrix M,, consisting of m, consecutive rows, such that > m.; 17. As is common,
let A == M = Σi= M M, be the global covariance matrix, and A₁ := .M." M, the local covariance matrix owned by node i. We denote by λ1, A2,... ..., da the eigenvalues of A in descending order and by 01, U2 Va the corresponding eigenvectors. Then, it well-knows that, if the gap 8- = X1 A between the two leading eigenvalues of A is positive, we can approximate the leading eigenvector by solving the following empirical risk minimization problem up to accuracy E:
x² argmin
(-1/0/22
x Aa
argminn 2xr 1= (-x² Au) argmin (-x"Ax), Ar), (1)
where sd-1 is the (d - 1) diimension sphere.
2

---

We define f gd-1 + R. with f(z) -x" Ar and ft sd + R. with fi(z) -x"A,a A,a. Since the inner product is bilinear, we can write the global cost as the sum of the local costs: f():) = E fi()).
In our analysis we use notions from the geometry of the sphere, which we recall in Appendix A
Related Work. As discussed, there has been a significant amount of research efficient variants of PCA and related problems [3] 24 25 26 1 32 [31. Due to space constraints, we focus on related work on communication--ffiienn algorithms. In particular, we discuss the relationship to previous round-efficien algorithms; to our knowledge, this is the first work to specifically focus on the bit complexity of this problem in the setting where data is randomly partitioned. More precisely, previous work on this variant implicitly assumes that algorithms are able to transmit real numbers at unit cost.
The straightforward approach to solve the minimization problem (I) in distributed setting, where the data rows are partitioned, would be to use a distributed version of the Power Method, Riemannian gradient descent (RGD), or the Lanczos algorithm. In order achieve an e-approxination of the minimizer z², the latter two algorithms require ( log(1 /c) rounds, where the O notation hides poly-togarithmic factors in d. Distributed Lanczos and accelerated RGD would improve this by an O(VA//5) /5) factor. However, Garber et al. [1] point out that, when is small, e.g. 5 = (11/ mm) then unfortunately the number of communication rounds would increase with the sample size, which renders these algorithms non-scalable.
Standard distributed convex approaches, e.g. [16 27l. do not directly extend to our setting due to non-convexity and the unit-norm constraint. Garber et al. M proposed a variant of the Power Method called Distributed Shift- Shift-and-Inveri (DSI) which converges in roughly O b Huang and Pan [14 aimed to improve the dependency of the algorithm on & and 8, and proposed an algorithm called algorithm is shown to have round complexity O h log(1/s)), which does not scale with the O
(Vovmloo"" log" log(1/0)) rounds, where b is a bound on the squared 6--norm of the data.
Distributed Riemannian Eigensolver (CEDRE). This
sample size for 2211/m)) and has only logarithmic dependency on the accuracy E.
Technical issues in [1]. Huang and Pan 14] proposed an interesting approach, which could provide the most round-efficien distributed algorithm to date. Despite the fact that we find the main idea of this paper very creative, we have unfortunately identified a significant gap in their analysis, which we now outline.
Specifically, one of their main results, Theorem 3, uses the local gradient dominance shown in 31: yet, the application of this result is invalid, as it is done without knowing in advance that the iterates of the algorithms continue to remain in the ball of initialization. This is compounded by another error on the constant of gradient dominance (Lemma 2), which we believe is caused by typo in the last part of the proof of Theorem 4 in [31. This typo is unfortunately propagated into their proof. More precisely, the objective is indeed gradient dominated, but with constant which vanishes when we approach the equator, in contrast with the 2/5 which is claimed globally. (We reprove this result independently in Proposition 4). Thus, starting from a point lying in some ball of the minimizer can lead to a new point where the objective is more poorly gradient dominated.
This is non-trivial technical problem which, in the case of gradient descent, can be addressed by choice of the learning rate depending on the initialization. Given these issues, we perform new and Normally-corred analysis for gradient descent based on new convexity properties which guarantee convergence directly in terms of iterates, and not just function values. We would like however to note that our focus in this paper is on bit and not round complexity.
3 Step 1: Computing the Leading Eigenvector in One Node
3.1 Convexity-like Properties and Smoothness
We first analyze the case that all the matrices A; are simultaneously known to all machines. This practically means that we can run gradient descent one of the machines (let be the master), since we have all the data stored there. All the proofs for this section are in Appendix B
3

---

Our problem reads
-x_ Ar E8
min
where A M M is the global covariance matrix. I Y A2 > 0, this problem has exactly two global minima: v1 and -t1. Let æ E 3 be an arbitrary point. Then æ can written in the form t Ei= Ogty Fixing the minimizer UI, we have that ball in g around vi is of the form
B = {x E gd-1 ai 2 a} = {a E g (z, v1) 2
for some a. Without loss of generality, we may assume a > 0 (otherwise, consider the ball around -1 to establish convergence -1). In this case, the region as defined in Def. TT is identical to B. and xx' : 2/1
We begin by investigating the convexity properties of the function -z Ax. In particular, prove that this function is weakly-que quasi-conver (see Def. [3} with constant 2a in the ball B,, See also Appendix A for definition of the Riemannian gradient grad and the logarithm log.
Proposition 1. The function f(x) xx Ax. satisfies
2a(f(x) f*) I (grad (z), log,(x"))
for any x Bg with a >
We continue by providing quadratic growth condition for our cost function, which can also be found in (30 Lemma 2] in slightly different form. Here dist is the intrinsic distance in the sphere, that we also define in Appendix
Proposition 2. The function f(x) xx
TT Ar satisfies
[(x) - 5² dist
for any E Ba with > 0.
According to Def. 14 this means that satisfies quadratic growth with constant μ ==
Next, we prove that quadratic growth and weak-quasi-convent imply "strong- woooggg-aakcconvxiit"" property. This similar 7 Proposition 3.1].
Proposition 3. A function f: s + R which satisfies quadratic growth with constant μ > 0 and is quueaakl-quas--convee with respect minimizer x² in ball of this minimizer, satisfies for any x in the same ball the inequality
f(x) - f (grad (x), log, (xg),
11 dist xx)
2
Interestingly, using Proposition we can recover the gradient dominance property proved Theorem 4]. Proposition 4. The function f(x) - Ar satisfies
B1
fx))l 2 5a² (f2((()) ff)
for any x E B, with > 0.
The proof Proposition 4 can be done independently based only on Proposition 3 and 7 Lemma 3.2]. This verifies that our results until now optimal regarding the quality of initialization since gradient dominance proved in 31 is optimal in that sense (i.e. the inequalities in the proof are tight).
Proposition 5. The function f(x) = -T Ax is geodesically 2 on the sphere.
Thus, the smoothness constant of as defined in Def. 12 equals 7 = 2 Similarly, let 73 denote the smoothness constant of Si, which equals twice the largest eigenvalue of A,. In order to estimate in distributed fashion using only the local data matrices A,, we shall use the over-approximation
y max:=1.... %t

---

3.2 Convergence
We now consider Riemannian gradient descent with learning rate 17 starting from point
r00)) Ba:
((++ =
(x)
where exp is the exponential map of the sphere, defined in Appendix
Using Proposition 3 and a proper choice of 7,, we can establish a convergence rate for the instrinsic distance of the iterates the minimizer.
Proposition 6. An iterate of Riemannian gradient descent for f(x) -x? Ax. starting from point a((() Ba and with step-size 1i S a/7 where Z 2x1, produces point x((+1) that satisfies
dist (x: xx) (1 -ajm) dist?( ,," for μ 6/2.
Note that this result implies directly that if our initialization 20)) lies in the ball B., the distance of x(1) to the center x² decreases and thus all subsequent iterates continue being in the initialization ball. This is essential since it guarantees that the convexity-like properties for f continue hold during the whole optimization process. The proof of Proposition is in Appendix B
This result is also valid in more general setting when applied arbitrary f: S R that 20-weakly quasi convex and y-smooth, and that satisfies quadratic growth with constant μ. The Euclidean analogue of this result be found in 7 Lemma 4.4].
4 Step 2: Computing the Leading Eigenvector in Several Nodes
We now present our version of distributed gradient descent for leading eigenvector computation and measure its bit complexity until reaching accuracy in terms of the intrinsic distance of an iterate (F{) from the minimizer x² (x* is the leading eigenvector closest to the initialization point).
Lattice Quantization. For estimating the Riemannian gradient in distributed manner with limited communication, we use a quantization procedure developed in [8. The original quantization scheme involves randomness, but vector that we want to encode. This is similar to the quantization scheme used by 3 and has the following properties.
deterministic version of it, by picking up the closest point the
Proposition 7. 18 3) Demoting by b the number of bits that each machine uses to communicate, there exists a quantization function
Q: Rd-1 Rd- x
Rd11
which, for each > O, consists of encoding function enc R + (0. 13 and a decoding one dec {0, 13 Rd-1 + RR 1- such that, for all x' E Rd-1
dec (enc w. (I), x²) (xx x² w). if [| x' $ y- 1Q(z, w) 2|| S w. if || x'
Ify/w 1, the cost of the quantization procedure in number of bits satisfies b b O((d (1)o,, (2))
Odlog
In the following, the quantization takes place in the tangent space of each iterate Sd-1, 1 which is linearly isomorphic to dd-1 We denote by Q2 the specification of the function Q at T_ The vector inputs of the function Qu are represented in the local coordinate system of the tangent space that the quantization takes place at each step. For decoding at > 0, we use information obtained in the previous step, that we need translate the same tangent space. We do that using parallel transport I, which we define in Appendix AA]
Algorithm
We present now main algorithm, which is inspired by quantized gradient descent firstly designed by 22,. and its similar version in 3
5

---

1. Choose an arbitrary machine to be the master node, let it be in-
2. Choose (0)) E s (we analyze later specific ways to do that).
3. Consider the following parameters
pm, K == 2
0= vo(1 - vo)
Vo
4
V = 0K + Vo, R = >K((E))
where D an over-upproximation for dist(z x²) (the way to practically choose D is discussed later).
We assume that (os(D) S z. otherwise we run the algorithm with 2-
In
Sd=
4. Compute the local Riemannian gradient gi at x (0)) in each node.
5. Encode gr (z!") in each node and decode in the master node using its local information:
(R(O)
= Qu:co fi grad (x) (z) 2,YY,,
2n
6. Sum the decoded vectors in the master node:
grad (
7. Encode the sum in the master and decode in each machine using its local information:
d (0) Q2.00) (,,00) gradfi(>!")). (P(
2y 010
For O:
8. Take gradient step using the exponential map:
(1+1) = (x
with step-size (the step-size is discussed later).
In T T S11:
9. Compute the local Riemannian gradient grad (xx(+1)) at (++11) in each node.
Encode gradf((x(t+1))) (x(++1) in each node and decode in the master node using its (parallelly transported) local information from the previous step:
R(++1)
OR(-1))
Q (211)
(rad)((z((+1)),,
rwitt1) mic)
di
2n
Sum the decoded vectors in the master node:
(t-1)
grad (x((t))
12. Encode the sum in the master and decode in each machine using its local information the previous step after parallel transport:
(1+1) = Quicon (fOO-
pe (+D 200
(0)) (1+
R(+1)
(R(++I)
2
Convergence We first control the convergence of iterates simultaneously with the convergence of quantized gradients.
Note that
1 1- Vo =
1 Vo
V2V1+0 oo
I+o
VE
2 <
2
2
2

---

Lemma 8. If S Coo(D) the previous quantized gradient descent algorithm produces iterates x!
and quantized gradients q(t)) that satisfy
(R()
dist (x!" ,x) D² llu: (1) - - gradff(x
S
2n 11 grad S OR!
The proof is a Riemannian adaptation of the similar one in 22 and 3 and is presented in Appendix C
We recall that since the sphere is positively curved, it provides a landscape easier for optimization.
is quite direct to derive a general Riemannian method for manifolds of bounded curvature using more advanced geometric bounds, however this exceeds the scope of this work which focuses on leading eigenvector computation.
We now move main complexity result.
Theorem 9. Let S cos D! Then, the previous quantized gradient descent algorithm needs most
b b=o -0(mm 0mll
1 log
log D) (9))
cos( (D)ó
cos(D)ōn
bits in total to estimate the leading eigenvector with accuracy measured in intrinsic distance.
The proof is based on the previous Lemma S in order to count the number of steps that the algorithm needs to estimate the minimizer with accuracy and Proposition to count the quantization in each round. We present the proof again in Appendix
Step 3: Dependence on Initialization
5.1 Uniformly random initialization
The cheapest choice initialize quantized gradient descent is point in the sphere chosen uniformly random. According Theorem 4 in [3I], such random point (x))) will lie, with probability at least 1 - p, in a ball Ba (see Def. where
dist(x ,x) arccos
(3)
V
Here, is universal constant (we estimated numerically that e is lower bounded by 1, see Appendix D thus we can use = 1). By choosing the step-size 7)
1) cp
ddy
and the parameter D as
P Vdl
= arccos
we are guaranteed that cost D) and (0) ,x*) S with probability at least - p. Our
analysis above therefore applies (up to probability - p) and the general communication complexity result becomes
nd (D)n
log D
)
-0
nd mynn
n nyynn
log = - 0 nd
D
n log 2
D
log cos D)ón
log
log
17
2²2 dr the number bits satisfies (up to probability 1 - p) the upper bound
Substituting 1?2
ndy p2 pô
d
d2
b =
log
log

---

5.2 Warm Start
reasonable strategy to get more accurate initialization to perform an eigenvalue decomposition to one of the local covariance matrices, for instance Ai, (in the master node to), and compute its leading cigenvector, let it be w'o For simplicity we will assume here that each machine hosts the same number of data points m, ms Then we communicate y to all machines in order to the normalized quantized approximation (0) as initialization. We define:
x (0) Q 0 ,0,,,', v", 111 - v'1l, (000 x²)
2(v2 +
((0)
2 (0)
||
where Q is the lattice quantization scheme in Rd (i.e. we quantize the leading eigenvector of the master node as a vector in I& and then project back to the sphere). The input and output variance in this quantization can be bounded by constants that we can practically estimate (see Appendix DD.
Proposition 10. Assume that our data are i.i.d and sampled from a distribution D bounded lz norm by constant h. Given that the eigengap 8, the number of machines and the total number of data points m satisfy
Vmvva log d p
8 8 Ω
(4)
we have that the previous quantization costs ((nd) many bits and (x:1), x²)) is lower bounded by constant with probability at least p.
Thus, if bound is satisfied, then the communication complexity becomes
T2 8
log 1 - - (nd
b O
nd
log
many bits in total with probability least (notice that this be further simplified using bound This because D) in Theorem 9 is upper bounded by a constant and the communication cost of quantizing ydo does affect the total communication cost. If
estimate the specifie relation in bound 4 (the constant hidden inside !2), then we can compute estimations of the quantization
parameters in the definition of (0)
Condition 4 is quite typical in this literature; see 14] and references therein (beginning of page 2), as we also briefly discussed in the introduction. Notice that vmm appears in the numerator and not the denominator, only because we deal with the sum of local covariance matrices and not the average, thus our eigengap is m times larger than the eigengap the normalized covariance matrix. Denoting by 8' the eigengap of the normalized covariance matrix, bound can be written equivalently as
Vn Vm
8' 2
log
where m4 is the number of data points owned by the master node (and any other machine) and Ω hides logarithmic factors from the lower bound.
6
Experiments
We evaluate our approach experimentally, comparing the proposed method of Riemannian gradient quantization against three other benchmark methods:
Full-precision Riemannian gradient descent: Riemannian gradient descent, as described in Section 3.2 is performed with the vectors communicated full (64-bit) precision.

---

Euclidean gradient difference quantization: the 'naïve' approach to quantizing Riemannian gradient descent. Euclidean gradients are quantized and averaged before being projected to Riemannian gradients and used to take step. To improve performance, rather than quantizing Fuclidean gradients directly, we quantize the difference between the current local gradient and the previous local gradient, at each node. Since these differences are generally smaller than the gradients themselves, we expect this quantization to introduce lower error.
Quantized power iteration:
we also use as benchmark a quantized version of power iteration, a common method for
computation given by the update rule |m Ax(1)) can be computed distributed fashion by communicating and summing the vectors A,x i S n. It these vectors that we quantize.
x(++)
Az(t)
All three of the quantized methods use the same vector quantization routine, for fair comparison.
Figure 1 Convergence results on real datasets
Comparison quartization methods for Remannian gra diant descant
Compa on of quantization methods for Romannian gradient descant Mice Protein Expression, m 1UBU *77 eta -3e-4o *a bitts =4
110 humann Activey, m =24075 d -00 eta 11-00 n -5 bits
prokiin Sermaniie G3
-10
Suckkeeng Qusnibbr ternton
Sweone groder unitration
-00
-10
-12
Epech
Epoch
Comparsos of quauniizatiion mathods Riemannian gradient descant
Comparison of quuantization methods Rlemannian gradiont descent Libras
0 x199
Spambase, =4001 1e-09 n bits
Rlenaminn GD
Ful prid Rarnara G3
Euskddan gradien areenco quartaaee
Récsamian grodeni cunnttation
90
Eudicesn gracier itheence quatitation
Thens en grocket connttatton
Quuntibed power beorinn
Quarltadd posertarsion
-40
Coel x Ax
.00
--3-
-80
-4
-100
-5
50 -120
We show convergence results (Figure T for the methods on four real datasets: Human Activity from the MATLAB Statistics and Machine Learning Toolbox, and Mice Protein Expression, Spambase, and Libras Movement from the UCI Machine Learning Repository [9] All results are averages over 10 runs of the cost function -x- Ax (for each method, the iterates 1- 1-ball, so a lower value of -z -zTAx Ax corresponds to being closer to the principal eigenvector).
normalized lie in the
All four datasets display similar behavior: our approach of Riemannian gradient quantization out- performs naive Euclidean gradient quantization, and essentially matches the performance of the full-precision method while communicating only bits per coordinate, for a 16x compression. Power iteration converges slightly faster (as would be expected from the theoretical convergence guarantees), but is much more adversely affected by quantization, and reaches a significantly suboptimal result.
We also tested other compression levels, though we do not present figures here due to space constraints: for fewer bits per coordinate (1-3), the performance differences between the four methods are in the same order but more pronounced; however, their convergence curves become more erratic due to increased variance from quantization. For increased number of bits per coordinate, the performance of the quantized methods reaches that of the full- preccisioioon methods at between 8-10 bits per coordinate depending on input, and the observed effect of quantized power iteration converging to suboptimal result gradually diminishes. Our code is publicly available
https: http:///gtthab..oossS---AALLL///
9

---

7 Discussion and Future Work
We proposed communication efficient algorithm compute the leading eigenvector of covariance data matrix in a distributed fashion using only limited amount of bits. To the best of our knowledge, this is the first algorithm for this problem with bounded bit complexity. We focused on gradient descent since it provides robust convergence in terms of iterates, and can circumvent issues that arise due to a cost function that only has local convexity-like properties; second, it known be amenable to quantization.
We would like briefly to discuss the round complexity of our method, emphasizing inherent round complexity trade-offs due to quantization. Indeed, the round complexity our quantized RGD algo- rithm for accuracy e is (,)²g), where e appears inside a logarithm. Specifically, the initialization 0(1vV)). If we have warm start
(0()) is uniformly random in the sphere (section 5.1), then (section 5.2), then a = O(). Both statements hold with high probability.
For random initialization, the round complexity of distributed power method is o O(3) (see the table at the beginning of page 4 in [II.. This is because the convergence rate of power method is tan(dis(( x)) S (1 \1 /5)* (dnntis(("), x")). The dimension is again present in that rate, because tan(dist(zro.. )) = @(vd), d). but it appears inside a logarithm in the final iteration complexity.
Since we want to insert quantization, the standard way to do so by starting from a rate which is contracting (see e.g. (222. It may in theory be possible to add quantization in a non-contracting algorithm, but we are unaware of such a technique. Thus, we firstly prove that RGD with step-size $ is contracting (based on weak-convexity and quadratic growth) and then we add quantization, paying the price of having an extra O(1/a²) = (dd) in our round complexity. Thus, our round complexity is slightly worse than distributed power method, but this is in some sense inherent due to the introduction of quantization. Despite this, our algorithm can achieve low bit complexity and is more robust in perturbations caused by bit limitations, as our experiments suggest.
We conclude with some interesting points for future work.
First, our choice of step-size is quite conservative since the local convexity properties of the cost function improve when approaching the optimum. Thus, instead of setting the step-size to be 7 a/>, where is a lower bound for (x) xx"), we can re-run the whole analysis with adaptive step-size 1| (z) _x*/ /y based on an improved value for from (a((), (x("),")) This should improve the dependency on the dimension of our complexity results.
Second, the version of Riemannian gradient descent that we have chosen is expensive, since one needs to calculate the exponential map in each step using sin and cos calculations. Practitioners usually use cheaper retraction- versions, and it would be of interest to study how our analysis can be modified if the exponential map substituted by a retraction (i.e., any first-order approximation).
Third, it would be interesting to see whether the ideas and techniques presented in this work can also be used for other data science problems in a distributed fashion through communication--ffiiiaaa Riemannian methods. This could be, for example, PCA for the top k eigenvectors using optimiza- tion on the Stiefel Grassmann manifold (T) and low-rank matrix completion the fixed-rank manifold 22
Acknowledgments and Disclosure of Funding
We would like to thank the anonymous reviewers for helpful comments and suggestions. We also thank Aurelien Lucchi and Antonio Orvieto for fruitful discussions an early stage of this work. FA is partially supported by the SNSF under research project No. 192363 and conducted part of this work while at IST Austria under the European Union's Horizon 2020 research and innovation programme (grant agreement No. 805223 ScaleML). PD partly conducted this work while at IST Austria and was supported by the European Union's 's Horizon 2020 programme under the Marie Skłodowska-Cure grant agreement No. 754411.
10

---

Bibliography
[1] P-A Absil, Robert Mahony, and Rodolphe Sepulchre. Optimization algorithms on matrix manifolds. Princeton University Press, 2009.
[2] Foivos Alimisis, Peter Davies, and Dan Alistarh. Communicaiion-eficiennn distributed optimiza- tion with quantized preconditioners. In International Conference on Machine Learning. PMLR, 2021.
[3] Dan Alistarh and Janne H. Korhonen. Improved communication lower bounds for distributed optimisation. arXiv preprint arXiv:2010.08222. 2020.
[4] Maria-Florina Balcan, Vandana Kanchanapally, Yingyu Liang, and David Woodruff. Improved distributed principal component analysis. arXiv preprint 1oXiv:1408.5823, 2014.
[5] Christopher M Bishop. Puttern recognition and machine learning. Springer, 2006.
[6] Christos Boutsidis, David Woodruff, and Peilin Zhong. Optimal principal component analysis in distributed and streaming models. In Proceedings of the forty-eighth ammual ACM symposium on Theory of Computing, pages 236-249, 2016.
[7] Jingjing Bu and Mehran Mesbahi. A note on nesterov's 's accelerated method in nonconvex optimization: weak estimate sequence approach, 2020.
[8] Peter Davies, Vijaykrishna Gurunathan, Niusha Moshrefi, Salch Ashkboos, and Dan Alistarh. New bounds for distributed mean estimation and variance reduction. In International Conference on Learning Representations, 2021.
[9] Dheeru Dua and Casey Graff. UCI machine learning repository, 2017.
[10] Jianqing Fan, Dong Wang. Kaizheng Wang, and Ziwei Zhu. Distributed estimation of principal eigenspaces. Annals of statistics, 47(6))3009, 2019.
[11] Dan Garber, Ohad Shamir, and Nathan Srebro. Communication-effitii algorithms for dis- tributed stochastic principal component analysis. In International Conference Machine Learning, pages 1202-1212. PMLR, 2017.
[12] Samuel Horváth, Dmitry Kovalev, Konstantin Mishchenko, Sebastian Stich, and Peter Richtárik. Stochastic distributed learning with gradient quantization and variance reduction. arXiv preprint arXiv:1904.05115, 2019.
[13] Harold Hotelling. Analysis a complex of statistical variables into principal components. Journal of educational psychology, 24(6):417, 1933.
[14] Long-Kai Huang and Sinno Pan. Communication-efiin distributed PCA by Riemannian optimization. International Conference on Machine Learning, pages 4465-4474. PMLR, 2020.
(15) Rustem Islamov, Xun Qian, and Peter Richtárik. Distributed second order methods with fast rates and compressed communication. arXiv preprint arXiv:2102.07158, 2021.
[16] Martin Jaggi, Virginia Smith, Martin Takač, Jonathan Terhorst, Sanjay Krishnan, Thomas Hofmann, and Michael I Jordan. Communication-efficienn distributed dual coordinate ascent. arXiv preprint arXiv: 1409.1489. 2014.
[17] Ravi Kannan, Santosh Vempala, and David Woodruff. Principal component analysis and higher correlations for distributed data. In Conference Learning Theory, pages 1040-1057. PMLR, 2014.
[18] Hamed Karimi, Julie Nutini, and Mark Schmidt. Linear convergence of gradient and proximal- gradient methods under the polyak-lojasiewicz condition, 2020.
[19] Jakub Konečný, H Brendan McMahan, Felix X Yu, Peter Richtárik, Ananda Theertha Suresh, and Dave Bacon. Federated learning: Strategies for improving communication efficiency. arXiv preprint
2016.
11

---

[20] Jakub Konečný and Peter Richtárik. Randomized distributed mean estimation: Accuracy vs. communication. Frontiers in Applied Mathematics and Statisties, 4:62, 2018.
[21] Li. Concise Formulas for the Area and Volume Hyperspherical Cap. Asian Journal of Mathemalics & Statistics, 4(1):66-70, 2011.
[22] Sindri Magnússon, Hossein Shakki-Chadikollei, and Na Li. On maintaining linear convergence of distributed learning and optimization under limited communication. IEEE Transactions on Signal Processing, 08:6101-6116 2020.
[23] Andrew Ng, Michael Jordan, and Yair Weiss. On spectral clustering: Analysis and algorithm. Advances in neural information processing systems, 14:849-856, 2001.
[24] Erkki Oja and Juha Karhunen. On stochastic approximation of the eigenvectors and eigenvalues of the expectation of random matrix. Journal of mathematical analysis and applications, 106(1):665 84, 1985.
[25] Ohad Shamir. A stochastic pea and svd algorithm with exponential convergence rate. In International Conference on Machine Learning, pages 144-152. PMLR, 2015.
[26] Ohad Shamir. Fast stochastic algorithms for svd and pea: Convergence properties and convexity.
In International Conference on Machine Learning, pages 248-256. PMLR, 2016.
[27] Ohad Shamir, Nati Srebro, and Tong Zhang. Communicatiion-ffiiiienn distributed optimization using an approximate newton-type method. In International conference on machine learning, pages 1000-1008. PMLR, 2014.
[28] Ananda Theertha Suresh, X Yu Felix, Sanjiv Kumar, and H Brendan McMahan. Distributed mean estimation with limited communication. In International Conference on Machine Learn- ing, pages 3329-3337, PMLR, 2017.
[29] Bart Vandereycken. Low-rank matrix completion by riemannian optimization. SIAM Journal on Optimization, 23(2):1214-1236 2013.
[30] Zhiqiang Xu, Xin Cao, and Xin Gao. Convergence analysis of gradient descent for eigenvector computation. International Joint Conferences on Artificial Intelligence, 2018.
[31] Hongyi Zhang, Sashank J Reddi, and Suvrit Sra. Riemannian svrg: Fast stochastic optimization on riemannian manifolds. arXiv preprint 17Xiv:1605.07147, 2016.
[32] Hongyi Zhang and Suvrit Sra. First-order methods for geodesically convex optimization. In Conference Learning Theory, pages 1617-1638. PMLR, 2016.
12


Reviews:
[{"review_id": "upCZdAXLFqw", "reviewer": "Reviewer_m9Bn", "summary": "The paper considers the problem of distributed computation of the leading eigenvector of a data covariance matrix. To solve this problem, a distributed variant of Riemannian gradient descent is proposed. Convergence is proven and an upper bound is given on the number of bits that need to be transmitted to estimate the leading eigenvector. Finally, the dependency of the method on its initialization is analyzed.\n", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "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": "", "main_review": "Originality:\n\nThe paper presents a novel algorithm for the computation of the leading eigenvector of a covariance matrix, as it is done in principal component analysis. The paper fills a gap in the literature since compared to previous approaches, the proposed technique is specifically designed to minimize the bandwidth cost, i.e. the number of bits transmitted between the nodes. The main contribution of the paper is the theoretical analysis of the approach, in particular the convergence is proven and a bound on its communication complexity is given.\n\nQuality: \n\nThe derivation of the algorithm is technically sound, and the authors show theoretically that the approach is indeed communication efficient. However, there is no experimental evaluation of the approach. For a paper proposing a novel algorithm, I would expect at least some experiments that show its usefulness in practice. Moreover, while the approach is motivated by the application to principal component analysis, the paper only focuses on computing the first eigenvector, while in practice typically multiple eigenvectors are required.\n\nClarity: \n\nThe paper is clearly written and well-organized. However, as stated above, a section containing the experimental validation of the approach as well as details on the implementation are missing.\n\nSignificance:\n\nDue to the ever increasing amounts of data, coming up with novel distributed algorithms is a worthy goal. As such, the paper makes an important contribution to the field, in particular since it is specifically designed to be efficient in terms of communication complexity. Moreover, the distributed variant of Riemannian gradient descent may turn out to be useful in different contexts. On the other hand, the paper appears to be incomplete as it lacks any sort of experimental evaluation. Moreover, for the approach to be adopted for PCA it should be extended to compute multiple eigenvectors of the covariance matrix.\n\n\nPost Rebuttal:\n\nMy main concern was the lack of experimental evaluation, which was addressed by the authors in their rebuttal. I have therefore increased my score.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The authors discuss limitations of the work in the last section of the paper. Since this is a methodological contribution, there is no direct societal impact.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "5"}, "extra_text": {}}, {"review_id": "fglIylyJ54", "reviewer": "Reviewer_F3cp", "summary": "This paper studies the problem of finding the leading eigenvector in a distributed setting, where the $i$-th worker holds $m_i$ pieces of local data and the goal is to find the leading eigenvector of the global covariance data matrix while minimizing the total communication cost. Unlike most of the previous results that largely focus on round complexity, this paper considers bit complexity. \n\nThe proposed algorithm is a combination of the Riemannian gradient descent over an $\\ell_2$ sphere and the vector quantization method based on Davies et al. (2021). The main technical contribution of this paper is showing a \"strong-weak-convexity\" result of the quadratic function $-x^TAx$ over a Riemannian manifold, which resolves a technical issue in Huang et al. (2020).  ", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "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": "", "main_review": "Though this paper fixes a technical issue in the previous work Huang et al. (2020),  I do feel the novelty is rather limited since both the Riemannian gradient descent and the vector quantization technique are well-studied before. In addition, I have the following concerns/suggestions:\n\n1. The comparison with prior works can be improved. For instance, although previous works mainly focus on the round complexity, one can apply standard vector quantization techniques to convert the round complexity to the bit complexity. The authors should compare their proposed scheme with this baseline. \n\n\n2. Following the previous point, to study the round complexity of the proposed scheme, one can check in how many iterations the Riemannian terminates. However, in the Algorithm section (Page 6) the authors did not specify the stopping criteria.\n\n\n3. The final bit complexity involves a term $\\gamma$ (which is defined as $n\\max \\lambda_i$). To better compared with other results, the authors can show how it scales as a function of $n, b$ under the i.i.d. setting (and perhaps assume $m_i = m$).\n\n\n4. An experimental result to verify the proposed scheme would highly strengthen the paper.\n\n\n5. The presentation can be improved. For instance, in the main body, several important definitions that appear in the appendix are referenced. Some of them should be moved to the main paper. Also, the algorithm section on Page 6 can be made more concise.\n\n\n6. A minor typo in line 77: $A = \\sum_{i}^n M_i^TM_i$ instead of  $\\sum_{i}^d M_i^TM_i$.\n\n\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The authors have listed the limitations as part of future works. This work does not present any foreseeable societal consequences due to its theoretical nature.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "6"}, "extra_text": {}}, {"review_id": "0LiZyFyqWH", "reviewer": "Reviewer_YLwt", "summary": "This paper proposed a communication efficient distributed algorithm to compute the leading eigenvector of a covariance data matrix where the rows of the matrix are distributed on multiple nodes. They presented a distributed Riemannian gradient descent to compute the leading eigenvector and avoid projection onto constraint operation. To prove the convergence of Riemannian gradient descent, they first derive the weakly-quasi-convex property of the objective.  A quantization scheme is also proposed to save communication costs. Finally, they proved the convergence of the algorithm, and the communication complexities for different initialization schemes.\n", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "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": "", "main_review": "Distributed eigenvector estimation is important in many distributed ML algorithm's implementation, and communication cost is also a pain point in distributed learning. This paper, for the first time, tackles these two aspects together. Using Riemannian gradient descent to optimize manifold-constrained problems is a popular idea, and so is the quantized idea to save communication cost. The authors are the first to combine these ideas and apply them to distributed PCA. The structure of this paper is clear and easy to follow. I believe a reader without prior knowledge can follow this paper. The paper is well-organized and the presentation of the theorem/lemma is clear. \n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "1. For the quantization algorithm, it will be great if the authors can give the explicit mapping rule in the main paper, since the readers may not be familiar with lattice quantization. \n\n2. I would like the authors to discuss the relation between this work and federated PCA [Grammenos, Andreas, et al.19], where they also consider data matrix is distributed on different nodes.\n \n\n3. I would like to see some experiments on the effectiveness of this algorithm, so that we can see how the quantization idea improves the communication cost compared to the non-quantized version.\n", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "2"}, "extra_text": {}}]

Author rebuttal:

Meta-review: All reviewers support acceptance for the contributions on communication-efficient distributed PCA and the theoretical aspects of the proposed approach. I also recommend acceptance. However, please make sure to incorporate experimental results in the final version to validate the theoretical findings. A discussion on how to compute multiple eigenvectors would also strengthen the paper.
Decision: Accept (Poster)