Jerome Fellus
email: jerome.fellus@ensea.fr

David Picard
email: picard@ensea.fr

Philippe-Henri Gosselin
email: gosselin@ensea.fr

Asynchronous Gossip Principal Components Analysis

Keywords: Distributed Machine Learning, Dimensionality reduction, Gossip protocols

come

Introduction

Dimensionality reduction plays an important role in solving large scale machine learning problems where input data usually consists of a huge number of observations in a high-dimensional space. Classification, regression, or similarity ranking of such raw data often raise computation and storage issues. In practice, the intrinsic dimensionality of observed phenomena is much lower than the extrinsic dimension of the input space. Dimensionality reduction then aims at pro-jecting input data into a lower-dimensional space such that subsequent learning stages keep a maximal accuracy.

Principal Components Analysis (PCA) is a linear approach to dimensionality reduction [START_REF] Pearson | on lines and planes of closest fit to systems of points in space[END_REF][START_REF] Hotelling | Analysis of a complex of statistical variables into principal components[END_REF]. Given a sample matrix X ∈ R D×n made of n observations in R D , PCA finds an orthonormal basis U ⋆ = [u 1 . . . u q ], u k ∈ R D that projects the input sample X into the q-dimensional subspace, q < D, that retains the maximal variance in X. Equivalently, the PCA solution is the linear projection that best conserves the Gram matrix (i.e., the matrix of pairwise inner-products):

U ⋆ = arg min U∈R D×q X ⊤ X -X ⊤ UU ⊤ X 2 F s.t. U ⋆⊤ = U ⋆-1 (1)
The optimal conservation of the inner product makes PCA particularly suited to feed algorithms that solely rely on the inner product instead of the input data [START_REF] Bishop | Pattern recognition and machine learning[END_REF] (e.g., Support Vector Machines). PCA enjoys a closed-form solution, as U ⋆ is made of the q leading eigenvectors of the sample covariance matrix [START_REF] Hotelling | Analysis of a complex of statistical variables into principal components[END_REF]:

C = 1 n XX ⊤ -µµ ⊤ where µ = 1 n X1 is the sample mean (2)
Like most statistical learning tools, PCA was formulated for centralized setups where all data are available at a single location. This assumes that the solution can be computed by a single machine and that all intermediary results fit in the main memory. However, this assumption is unrealistic in most applicative fields that deal with very large sample matrices. For instance, in biomedical, multimedia, or remote sensing applications, D and n often grows up to millions. The sample and covariance matrices then scale in Terabytes. Moreover, the O(D 3 ) complexity of covariance eigendecomposition translates into an exa-flop computation cost. Besides, along with the democratization of connected devices, data tends to originate from an increasing number of distributed sources with reasonable computing capacity, immersed in large unreliable networks without any central coordination. This has led to a number of so-called Distributed PCA algorithms, designed to deal with the spread of input data over multiple networking nodes.

Because computing µ and C would involve the full data X which is unknown to individual nodes, distributed PCA requires dedicated algorithms combining node-local optimization and inter-node communications. Distributed PCA encompasses two main scenarios, depending of the way the entries of X are spread over the networking nodes. Consider a network made of N (strongly-connected) nodes. In a Distributed Samples (DS) scenario, each node i holds a distinct sample X i ∈ R D×n i of n i observations (i.e., the columns of X are distributed). Con-versely, in a Distributed Coordinates (DC) scenario, each node holds all n observations, but only gets a subset X i ∈ R D i ×n of their components (i.e., the rows of X are distributed). On average, each X i is then N times smaller than X, thus n i ≪ D in DS case while D i ≪ n in DC case. No assumption is made on their exact sizes, as they may be very different at each node. In both DS and DC scenarios, a typical objective is to provide all nodes with compressed representations of their locally-hosted observations that account for the contribution of all components. Nodes then have to find a consensus basis U ⋆ that minimizes the PCA objective defined in Equation ( 1) over the complete data. Usually, one also wants an operator that allows projection of future observations into the same output space, but not all approaches are able to provide such operator at a low cost. In spite of similar goals, the DS and DC scenarios have been tackled separately with very different approaches in the distributed PCA litterature [START_REF] Bertrand | Distributed adaptive estimation of covariance matrix eigenvectors in wireless sensor networks with application to distributed pca[END_REF].

In this work, we consider the asynchronous decentralized PCA problem, which specializes distributed PCA by adding the following constraints:

(C1) No sample exchange -Samples cannot be exchanged between nodes, for size, privacy or property reasons.

(C2) Asynchrony -Nodes must never wait for each other.

(C3) Decentralization -All nodes and links must play the same role. Formally, nodes and links must be selected for communication with the same probability and all nodes must run the same procedures.

(C4) Consensus -All nodes must obtain the same orthogonal basis. Node-local solutions must allow projection of future observations. Satisfying these four constraints, a distributed PCA algorithm gains applicability to a wider range of networking situations such as sensor networks, Internetenabled multimedia devices, cloud computing systems, etc, where central coordination or synchronous functioning can be inapplicable.

In this paper, we propose a decentralized and asynchronous algorithm called Asynchronous Gossip Principal Component Analysis (AGPCA) that satisfy all the above constraints. Our algorithm is a revision and extension of previous work presented in [START_REF] Fellus | Dimensionality reduction in decentralized networks by gossip aggregation of principal components analyzers[END_REF]. The original contributions of this paper are the following:

-We give a formal and in-depth description of AGPCA in the DS case, as well as the intuitions leading to the algorithm.

-We propose an extension to the DC case through a dual transcription.

-Provided a low rank property is met on the data, we give a theoretical garantee that AGPCA yields the exact solution of the centralized PCA.

-We present experiments for both DS and DC scenarios, as well as results for various network topologies.

The remaining of this paper is organized as follows: In the next Section, we detail related works on distributed PCA algorithms. Then, we present our AGPCA algorithm for the DS case in Section 3. We present the extension to the DC case in Section 4. In section 5, we theoretically show that the output of AGPCA is identical to a centralized PCA under a low-rank assumption on the data. The last section gathers experimental results both in DS and DC case, before we conclude.

Related Work

In this section, we present existing algorithms for distributed PCA. We first present methods that integrate prior information on the input data. Then, we compare existing algorithms in terms of decentralization and asynchrony. Finally, we discuss the benefits of approaches based on model aggregation over those based on iterative optimization passes over the data.

Prior knowledge about input data

Existing distributed PCA approaches can be first distinguished by their level of prior knowledge about the input data matrix X. Indeed, X can either carry node-local observations independently of their network relationships or integrate properties of the network graph itself. In the latter case, a typical object of interest is the adjacency matrix of the weighted network graph, whose entries correspond to some scalar relationship between nodes. For instance, when these entries represent estimated geographic distances between neighbouring sensors, their absolute geographic position can be recovered by computing the three principal components through distributed PCA [START_REF] Kempe | A decentralized algorithm for spectral analysis[END_REF].

Other methods have considered the case where the data distribution inherits specific characteristics from the network structure, such as statistical dependencies. This happens when, e.g., data is generated by flowing through directed paths along the network structure, thus making data at downstream nodes dependent of data at upstream nodes, but independent from each other. Properly modeling these statistical dependencies through Graphical Models, either undirected (e.g., Decomposable Gaussian Graphical Models [START_REF] Wiesel | Decomposable principal component analysis[END_REF]) or directed (e.g., Directed Gaussian Acyclic Graphical Models [START_REF] Meng | Distributed principal component analysis on networks via directed graphical models[END_REF]), one can benefit from the natural sparsity of the concentration matrix (i.e., the inverse covariance) or its Cholesky factor to estimate the principal subspace with reduced communication costs.

On the contrary, in this work the network topology has no relevance in the desired result, as information is solely carried by the nodes and not by the links. Still, link-related data can be seen as observations relative to one or both of their ends, making methods aimed at node-related data suitable for link-related data.

Decentralization and asynchrony

Another classification criterion separates decentralized approaches from those which assign node-specific roles and asynchronous approaches from those based on synchronous communications.

In [START_REF] Ordonez | Pca for large data sets with parallel data summarization[END_REF], a parallel PCA algorithm is proposed. Sufficient statistics X i 1 and X i X ⊤ i are locally computed at all nodes and transmitted to a master node that performs a global Singular Value Decomposition (SVD) to obtain the PCA result. This approach is only suitable when the master node can handle the O(D 3 ) complexity of the SVD and assumes that D ×D covariance matrices can be exchanged on the network, which is unrealistic in many large scale contexts.

In [START_REF] Bertrand | Distributed adaptive estimation of covariance matrix eigenvectors in wireless sensor networks with application to distributed pca[END_REF], a distributed PCA algorithm for the DC scenario is proposed. Exchanging only q × q matrices, nodes iteratively maximize the variance retained by the reduced basis. Even though the process is decentralized, nodes have to update their estimates synchronously before performing any further computation, thus violating (C2). The whole system performance is then limited by the slowest networking node. Moreover, synchronous updating is hard to sustain in large networks and can result in overwhelming idle phases even when nodes have identical computing resources.

A fully asynchronous and decentralized Power Iteration method is proposed in [START_REF] Korada | Gossip PCA[END_REF] using random matrix sparsifications and a nested Sum-Weight Gossip averaging protocol to reduce communication costs. However, its original formulation only provides the first principal component. Synchronous passes would be required to sequentially obtain the next ones.

Multiple passes over the data vs. model aggregation

An important feature of some methods is to require only one pass over nodelocal data. This is the case of [START_REF] Ordonez | Pca for large data sets with parallel data summarization[END_REF][START_REF] Bruneau | Aggregation of probabilistic PCA mixtures with a variational-bayes technique over parameters[END_REF], contrarily to [START_REF] Kempe | A decentralized algorithm for spectral analysis[END_REF][START_REF] Korada | Gossip PCA[END_REF] which require multiple access to the X i because of their iterative optimization scheme. Temporary access to data is frequently considered in online approaches like [START_REF] Bertrand | Distributed adaptive estimation of covariance matrix eigenvectors in wireless sensor networks with application to distributed pca[END_REF] to handle stream data under short-time ergodicity and stationarity assumptions.

In [START_REF] Bruneau | Aggregation of probabilistic PCA mixtures with a variational-bayes technique over parameters[END_REF], the authors propose an operator to aggregate Mixtures of Probabilistic PCA models (MPPCA, [START_REF] Tipping | Mixtures of probabilistic principal component analyzers[END_REF]) in a maximum-likelihood sense, without resorting to the original data used to train the models. Multiple models can thus be trained in a first phase, and aggregated in a second phase. Used in conjunction with the Gossip optimization framework proposed in [START_REF] Nikseresht | Gossip-based computation of a gaussian mixture model for distributed multimedia indexing[END_REF], one can easily obtain a decentralized and asynchronous PCA algorithm. However, as developped in [START_REF] Nikseresht | Estimation de modèles de mélange probabilistes: une proposition pour un fonctionnement réparti et décentralisé[END_REF], model aggregation give best results when models are further selected through neighbors cross-validation, which seems non-trivial to achieve in an asynchronous fashion.

In [START_REF] Macua | Consensus-based distributed principal component analysis in wireless sensor networks[END_REF], the authors first compute sufficient local statistics and aggregate them by means of a distributed consensus algorithm before locally computing PCA on the aggregated statistics. Aggregation of the local statistics requires exchanging D × D matrix estimates, which can be prohibitive for high-dimensional input spaces. Our approach is akin to [START_REF] Macua | Consensus-based distributed principal component analysis in wireless sensor networks[END_REF], but solves this problem by computing PCA before the distributed aggregation.

Asynchronous Gossip PCA for Distributed Samples scenario

In this section, we introduce our proposed AGPCA algorithm, which is suitable for both DS and DC scenarios. This section details AGPCA for the DS scenario. We present the extension of AGPCA to the DC scenario in the next section.

In a DS scenario, each node hosts a local sample

X i ∈ R D×n i made of n i observations [x 1 , . . . , x n i ] in R D .
In such scenario, AGPCA provides all nodes with the same orthogonal basis U ⋆ ∈ R D×q that spans the q-principal subspace of the full covariance matrix C as defined by Equation [START_REF] Hotelling | Analysis of a complex of statistical variables into principal components[END_REF]. Each node i can then project its local data X i to obtain Y i = U ⋆⊤ X i , where Y i ∈ R q×n i is the compressed representation of X i . Importantly, any future observation x new can be compressed in the same way: y new = U ⋆⊤ x new , even if X i is deleted or was only accessible during training (e.g., streaming data case).

The intuition behind AGPCA is built upon the following 4 main principles:

1. Distributed PCA can be formulated as a distributed averaging of covariance matrices followed by local eigendecompositions. 2. Asynchronous and decentralized averaging of covariance matrices is possible through Sum-Weight Gossip consensus protocols [START_REF] Kempe | Gossip-based computation of aggregate information[END_REF]. 3. By properly defining uniform scaling and sum operators for orthogonal bases, we can extend the Sum-Weight Gossip protocol to exchange only the

Algorithm 1 AGPCA-DS Emission procedure 1: a i ← X i 1 ; G i ← X ⊤ i X i ; w i ← n i 2: (V i , L i ) ← eigendecompose(G i ) 3: U i ← X i V i Λ -1 2 i 4: loop 5: j ← randomNeighbor(i) 6: (a i , L i , w i ) ← 1 2 (a i , L i , w i ) 7:
Send

(a i , U i , L i , w i ) to j 8: end loop Algorithm 2 AGPCA-DS Reception procedure 1: loop 2:
Upon receipt of (a j , U j , L j , w j )

3:

a i ← a i + a j ; w i ← w i + w j 4: Q 0 ← U i 5: for t ∈ [0, max t[ do 6: (Q t+1 , R t+1 ) ← QR(U i L i U ⊤ i Q t +U j L j U ⊤ j Q t ) 7:
end for 8:

U i ← Q max t ; L i ← diag(R max t ) 9:
end loop q leading eigenvectors of the local covariance matrices, and consequently achieve a better network efficiency.

Using the duality between the eigendecomposition of the covariance and

Gram matrices through Singular Value Decomposition (SVD), we can compute the local Gram matrices (n i × n i ) instead of the local covariance matrices (D × D) and consequently reduce the memory usage at each node.

In the remaining of this section, we detail each of these points leading to AG-PCA procedures for the DS scenario to be concurrently run at each node as presented in Algorithm 1 and 21 .

Distributed PCA as a distributed averaging problem

From Equation ( 1), the sample correlation matrix XX ⊤ and the sample sum X1 are sufficient statistics to compute the PCA solution. As observed in [START_REF] Macua | Consensus-based distributed principal component analysis in wireless sensor networks[END_REF][START_REF] Ordonez | Pca for large data sets with parallel data summarization[END_REF], these statistics can be obtained by computing partial statistics X i 1 and X i X ⊤ i at each node i and summing over i. As a result, the full sample mean and covariance matrix can be computed as a distributed average of locally-computed means µ i =

1 n i X i 1 and correlation matrices B i = X i X ⊤ i : µ = 1 n X1 = 1 n N i X i 1 = 1 i n i N i n i µ i (3)
C = 1 n XX ⊤ -µµ ⊤ = 1 n N i X i X ⊤ i -µµ ⊤ = 1 i n i N i B i -µµ ⊤ (4)
Provided these weighted averages are computed in a decentralized and asynchronous fashion, all nodes get the full covariance matrix C. Locally eigendecomposing C at every nodes finally gives the global PCA solution U ⋆ . Remark that X i 1 and X i X ⊤ i are easily computed together in a single pass over X i . Also notice that once local sufficient statistics are computed, input samples can be dropped without concern. Therefore, such a scheme is suitable when streaming data.

Gossip protocols for asynchronous decentralized averaging

To compute Equations (3-4) in a decentralized and asynchronous fashion, AG-PCA uses a Sum-Weight Gossip protocol [START_REF] Kempe | Gossip-based computation of aggregate information[END_REF][START_REF] Iutzeler | Analysis of sum-weight-like algorithms for averaging in wireless sensor networks[END_REF].

Sum-Weight protocols are an asynchronous subclass of Gossip consensus algorithms for distributed averaging. Gossip consensus algorithms proceed by iterative linear combinations of node estimates through randomized communications.

An example of Gossip averaging protocol is Newscast [START_REF] Jelasity | Newscast computing[END_REF]. In Newscast, random node pairs regularly average their estimates until reaching consensus. Assuming each node i holds a local vector estimate v i in a vector space V, at any random time t two random nodes s and r awake, exchange their estimates and update as follows:

v s (t + 1) = 1 2 (v s (t) + v r (t)) v r (t + 1) = 1 2 (v s (t) + v r (t)) (5)
This update rule entails two properties:

• Mass conservation. The sum of the estimates over the network is conserved:

i v i (t) = i v i (0)
• Convergence to consensus. The variance of the estimates tends to zero:

lim t→∞ i v i (t) - j v j (t) 2 = 0
A trivial consequence of these properties is ∀i, lim t→∞ v i (t) = 1 N j v j (0). Unfortunately, Newscast requires synchronous updating of random node pairs, thus violating constraint (C2).

Sum-Weight protocols aims at removing pairwise synchronization by adding a new estimate w i ∈ R called weight, with initial value w i (0) = 1. In contrast to Newscast, all nodes send their current sum and weight to randomly selected neighbors following independent Poisson emission clocks and select targets independently at random without waiting for their answer. The update rules of the sender node s and receiver node r are modified as follows:

v s (t + 1) = 1 2 v s (t) w s (t + 1) = 1 2 w s (t) (6) v r (t + 1) = v r (t) + 1 2 v s (t) w r (t + 1) = w r (t) + 1 2 w s (t) (7)
As first shown in [START_REF] Kempe | Gossip-based computation of aggregate information[END_REF] under synchronous assumptions, and then in [START_REF] Iutzeler | Analysis of sum-weight-like algorithms for averaging in wireless sensor networks[END_REF] for the general case, the quotient of the two estimates converges to the desired average.

∀i, lim t→∞ 1 w i (t) v i (t) = 1 j w j (0) j v j (0) = 1 N j v j (0) (8)
Moreover, convergence to the consensus is exponential provided that the network has a sufficient conductance [START_REF] Shah | Gossip algorithms[END_REF]. In this case, the number of message exchanges required to achieve a given estimation error ε scales logarithmically with the number of nodes and ε.

As weighted averages, µ and C can be estimated using the above-defined Sum-Weight protocol, by defining node-local estimates a i (t), B i (t) and weights w i (t) such that:

a i (0) = X i 1 B i (0) = X i X ⊤ i w i (0) = n i (9)
By applying the Gossip protocol defined by Equations (6-7) to a i (t),B i (t) and w i (t), we get a covariance estimate C i (t):

C i (t) = B i (t) w i (t) - a i (t)a i (t) ⊤ w i (t) 2 (10)
Note that initial estimates C i (0) are the covariance matrices of their corresponding X i . The limit in [START_REF] Meng | Distributed principal component analysis on networks via directed graphical models[END_REF] shows that each a i (t) w i (t) tends to the global mean µ and each C i tends to the global covariance matrix C:

∀i,          lim t→∞ a i (t) w i (t) = i X ⊤ i 1 i n i = i n i µ i i n i = µ lim t→∞ C i (t) = i B i (0) i w i (0) -µµ ⊤ = 1 n i X i X ⊤ i -µµ ⊤ = C (11)
Once each node gets a sufficiently accurate estimate for C, the final PCA result can be locally computed at any node i by eigendecomposition of C i . Remark this strategy, which will be referred to as Late-PCA in the rest of the paper, is used in [START_REF] Macua | Consensus-based distributed principal component analysis in wireless sensor networks[END_REF] in a synchronous consensus framework. Yet it suffers one major drawback: updating estimates B i using Equations (6-7) requires transmission of D × D matrices, which can be too large to exchange on the network. In [START_REF] Fellus | Dimensionality reduction in decentralized networks by gossip aggregation of principal components analyzers[END_REF], we thus proposed to reduce matrices B i by means of local PCA before their transmission, resulting in what we called an Early-PCA scheme.

Gossiping in the compressed domain: Early PCA

There are two main reasons that make the Late PCA approach unwanted, both based on the fact that such a scheme does not benefit from the rank-deficiency of the local covariance matrices implied by n i ≪ D. Firstly, the size of exchanged information is homogeneous to the input statistics (D × D), not to the output result (D × q). Therefore, we vainly exchange information that will be canceled in the end, because dimensionality reduction happens after the distributed averaging phase. Secondly, the distributed nature of the process does not allow any computational advantage for high dimension, since all nodes have to perform the same O(D 3 ) eigendecomposition operation. Consequently, if the eigendecomposition of C is the computational bottleneck due to large values of D, the distributed scheme is unlikely to bring any speed-up gain.

The Early-PCA approach aims at moving the dimension reduction step before the distributed averaging step, in order to take advantage of the rank deficiency of C i . To this purpose, we reformulate the updates rules in Equations (6-7) so as to handle eigenpairs instead of the full matrices B i . This allows us to implicitly compute the average of very large matrices without ever resorting to their explicit form but rather using their factorized expression.

Observe that the senders update rule ( 6) is a simple uniform scaling. Assuming s holds the decomposed form U s L s U ⊤ s of its estimate B s , the update is simply:

U s (t + 1) = U s (t) L s (t + 1) = 1 2 L s (t) (12)
This update leads to the emission procedure of AGPCA to be concurently run at every node of the network, as presented in Algorithm 1.

Adapting the receivers update rule ( 7) is slightly more involving since we need to compute the eigendecomposition of a sum of two matrices given their own eigendecompositions. We propose to rely on the well-known Orthogonal Iterations technique (as in [START_REF] Kempe | A decentralized algorithm for spectral analysis[END_REF]) to fit an orthonormal basis to the principal subspace of an input matrix by iterating QR decompositions. Assume B s (t) and B r (t) are respectively factorized as U s (t)L s (t)U s (t) ⊤ and U r (t)L r (t)U r (t) ⊤ . Starting from a random D ×q basis Q 0 , and denoting by QR(•) the economy QR decomposition, we iteratively compute

Q τ +1 R τ +1 = QR((B s + B r )Q τ ) = QR(U s L s (U ⊤ s Q τ ) + U r L r (U ⊤ r Q τ )) (13)
Q τ tends to become an orthonormal basis for the q-principal subspace of B s (t) + B r (t), with corresponding eigenvalues on the diagonal of R τ . Thus, Q ∞ gives U r (t+1) and the diagonal entries of R ∞ , denoted by diag(R ∞ ), give L r (t+1).

The parentheses in Equation ( 13) are of great importance. Indeed, observe that instead of expanding U s L s U ⊤ s and U r L r U ⊤ r , we first multiply U ⊤ s and U ⊤ r by Q τ , resulting in a q × q matrix. By doing so, we never store any D × D matrix, and the QR decomposition is rather performed on a D × q matrix. Finally, the dominant eigenvectors U ⋆ of C i are obtained by locally computing its eigendecomposition using Orthogonal Iterations (parentheses are also added to highlight the computational gain):

C i = 1 w i (t) (U i L i )(U ⊤ i ) - 1 w i (t) 2 (a i (t))(a i (t) ⊤ ) (14)
This update leads to the reception procedure concurrently run at every node and presented in Algorithm 2.

Since typically q ≪ D, combining Algorithms 1 and 2 allows to define an averaging protocol that exchanges D × q messages while staying equivalent to the original Sum-Weight protocol defined by Equations (6-7), provided that the rank of XX ⊤ is lower than q. This is formally supported in the theoretical analysis presented in Section 5.

Dual relation between covariance and Gramian eigendecompositions

While Early-PCA gets rid of D × D entities in the networking step, we can further avoid D × D matrices in the entire algorithm. This is particularly useful when covariance matrices do not fit in node memory. To this end, we must avoid explicit computation and storage of the initial B i (0) locally computed in Equation ( 9), reminding that we are only interested in their q principal subspace. Thankfully, the eigendecompositions

B i = U i L i U ⊤ i can be obtained from those of the n i × n i local Gram matrices X ⊤ i X i = V i Λ i V ⊤ i ,
through their well-known relation with the Singular Value Decomposition (SVD) of X. Indeed,

B 2 i = X i X ⊤ i X i X ⊤ i = X i V i Λ i V ⊤ i X ⊤ i = (X i V i Λ -1 2 i )Λ 2 i (X i V i Λ -1 2 i ) ⊤ and (X i V i Λ -1 2 i ) ⊤ (X i V i Λ -1 2 i ) = I Then, U i = X i V i Λ -1 2 i and L i = Λ i (15)
We thus compute, store and factorize

X ⊤ i X i (which is n i × n i ) instead of X i X ⊤ i (which is D × D)
and obtain U i and L i with no additional storage cost.

Asynchronous Gossip PCA for Distributed Coordinates scenario

In this section, we present the DC scenario and show that AGPCA can solve it as well. Contrarily to the DS scenario, each node hosts a local sample Z i ∈ R D i ×n made of D i ≪ D dimensions of the same set of n observations. The full data Z corresponds to stacking the Z i in rows:

Z ⊤ = [Z ⊤ 1 , . . . , Z ⊤ N ] (16)
Observations thus lie in the direct sum of the locally observed subspaces i R D i . Let µ = Z1/n denote the data mean and Z = Z -µ1 ⊤ the centered data.

In DC scenarios, it is useless to provide nodes with a D × q orthogonal basis since nodes only hold part of the input vectors and are thus unable to perform Algorithm 3 AGPCA-DC Emission Procedure

1: C i ← Z i Z ⊤ i ; w i ← D i 2: (U i , L i ) ← eigendecompose(C i ) 3: V i ← Z ⊤ i U i L -1 2 i 4: loop 5:
j ← randomNeighbor(i)

6:

(V i , L i , w i ) ← 1 2 (V i , L i , w i ) 7: Send (V i , L i , w i ) to j 8: end loop Algorithm 4 AGPCA-DC Reception Procedure 1: loop 2:
Upon receipt of (V j , L j , w j )

3:

w i ← w i + w j 4: Q 0 ← V i 5:
for t ∈ [0, max t[ do 6:

(Q t+1 , R t+1 ) ← QR(V i L i V ⊤ i Q t +V j L j V ⊤ j Q t ) 7:
end for 8:

V i ← Q max t ; L i ← diag(R max t ) 9: end loop the projection. Instead, AGPCA directly provides a compressed representation Y ≡ U ⋆⊤ Z that account for all dimensions of all observations.

Using the same ideas that allowed us to improve the efficiency in the DS case, we rely on the duality between the covariance and the Gram eigendecompositions. Recall that

C = 1 n ZZ ⊤ = U ⋆ L ⋆ U ⋆⊤ (17)
Introducing the SVD of the centered data Z = UL

1 2 V ⊤ , where ZZ ⊤ = ULU ⊤ and Z⊤ Z = VLV ⊤ , we have U ⋆ = U and L ⋆ = L/n. Consequently, Y = U ⋆⊤ Z = nL ⋆ 1 2 V ⊤ (18)
Define D = Z⊤ Z/D.

13

Clearly its eigendecomposition is V(L/D)V ⊤ = V(nL ⋆ /D)V ⊤ , and

D = 1 D (Z - 1 n Z11 ⊤ ) ⊤ (Z - 1 n Z11 ⊤ ) = I - 11 ⊤ n Z ⊤ Z D I - 11 ⊤ n ( 19
)
Observing that Equation 19only involves Z ⊤ Z/D, we remark that

Z ⊤ Z/D = ( i Z ⊤ i Z i )/( i D i
) is a weighted average of the locally-computable uncentered Gram matrices Z ⊤ i Z i . This average can be computed using the same protocol as for the DS case, defined by Equations [START_REF] Tipping | Mixtures of probabilistic principal component analyzers[END_REF][START_REF] Nikseresht | Gossip-based computation of a gaussian mixture model for distributed multimedia indexing[END_REF]. Using the factors V and L/D of D in Equation ( 18), we can obtain the compressed representations scaled by √ D. For most applications, this scaling factor has no impact. In applications for which the output basis needs to be orthonormal (and not only orthogonal), the scale can be recovered by gossiping D i with an initial weight on 1, provided the size of the network is known at each node.

Early-PCA aggregation in DC scenarios

Clearly, the Early-PCA aggregation scheme described in Section 3 also applies: We only need to compute and exchange the eigendecomposed forms

Z ⊤ i Z i = V i L i V ⊤ i .
Since local Gram matrices are n × n and are much larger than the local correlation matrices (of size D i × D i ≪ n × n), we use the same dual property as in Section 3 to obtain the Gram eigendecomposition from the local correlation matrices

Z i Z ⊤ i = U i L i U ⊤ i : Z i = U i L 1 2 i V ⊤ i ⇒ V i = Z ⊤ i U i L -1 2 i
The initial V i L i V ⊤ i can then be obtained at a low storage and computation cost. The corresponding emission and reception procedures are shown in Algorithms 3 and4.

Once the averaging protocol has provided all nodes with the eigendecomposition of Z ⊤ Z/D, we can locally apply Orthogonal Iterations to Equation ( 19) to obtain the eigendecomposition V(L ⋆ /D)V ⊤ of D at all nodes. Each node finally computes the full compressed data using Y = n(L ⋆ /D) 1 2 V ⊤ . This shows that using our AGPCA strategy, we can efficiently solve distributed PCA for both DS and DC scenarios by just swapping the roles of covariance and Gram matrices. A nice consequence is that theoretical results in one case hold in the other one.

Projection of subsequent observations

In DC scenarios, projecting future observations is not as straightforward as in the DS case, because new observations generate new entries in all nodes. Thus, all nodes have to participate in the projection of new observations. Remark the compressed representation y of a new observation z is computed by:

y = n(L ⋆ /D) -1 2 V ⊤ Z ⊤ z D ( 20
)
Observing that Z ⊤ z/D = i Z ⊤ i z i /D is also a weighted average, y can easily be recovered by simply Gossiping Z ⊤ i z i with initial weights D i and combining the result with the previously obtained factors V and L/D.

Theoretical analysis

In this section, we theoretically prove that AGPCA yields the same solution as a centralized PCA, provided assumptions are made on the rank of input data. We limit our analysis to the DS case, as the same arguments can easily be translated to the DC case by substituting the covariance matrices by the Gram matrices. Noting X ∈ R D the input data, this result is expressed in the following theorem: Theorem 1. If rank(X) ≤ q, then AGPCA yields the exact solution of a PCA run over X:

∀i, U ⊤ i X i = U ⋆⊤ X i (21)
PROOF. The centralized PCA uses the following eigendecomposition:

C = 1 n XX ⊤ -µµ ⊤ = U ⋆ L ⋆ U ⋆⊤ (22)
The proof sketch is as follows: We first show that the local decompositions

X i X ⊤ i = U i L i U ⊤
i perfectly reconstruct the local data X i . Then, we show our Gossip protocol aggregating the U i allows for a perfect reconstruction of i X i X ⊤ i at every node. The proof is completed by remarking AGPCA provides the eigendecomposition of the covariance matrices obtained from the reconstruction of i X i X ⊤ i . To prove local factorizations induce no loss of information, we use the assumption that rank(X) ≤ q. We thus have:

rank(X) ≤ q ⇒ rank(XX T ) ≤ q (23)
Let us denote H q = span(X) ⊂ R q . Since at any node i, X i is a subset of X, we have span(X i ) ⊆ H q , and consequently rank(X i X ⊤ i ) ≤ q. It follows that the initial local decompositions keeping only the q leading eigenvectors (the others being associated with null eigenvalues) perfectly reconstruct the local correlation matrices:

∀i, U i L i U ⊤ i = X i X ⊤ i = B i (24)
Considering the Gossip averaging, we first prove our update rules in Equations [START_REF] Tipping | Mixtures of probabilistic principal component analyzers[END_REF][START_REF] Nikseresht | Gossip-based computation of a gaussian mixture model for distributed multimedia indexing[END_REF] are equivalent to the Gossip updates in Equations [START_REF] Kempe | A decentralized algorithm for spectral analysis[END_REF][START_REF] Wiesel | Decomposable principal component analysis[END_REF]. Remark that a s and w s are explicitly updated using Equation (6-7). Checking validity of the sender update of B s (line 6 in Algorithm 1) is trivial:

B s (t + 1) = U s (t) 1 2 L s (t)U s (t) T = 1 2 B s (t).
On the receiver side, provided that span(U s (t)) ⊆ H q and span(U r (t)) ⊆ H q , using Othogonal Iterations to obtain the q leading eigenvectors of the weighted sum has two interesting properties: The reconstruction is perfect, that is

U r (t+1)L r (t+1)U r (t+1) = U r (t)L r (t)U r (t) ⊤ + 1 2 U s (t)L s (t)U s (t) ⊤ (25)
and the obtained basis is in the same subspace, i.e., U r (t+1) ⊆ H q . These properties come from the fact that span(U r ∪U s ) ⊆ H q , and thus rank(U r (t)L r (t)U r (t) ⊤ + U s (t)L s (t)U s (t) ⊤ /2) ≤ q Remarking that initially ∀i, span(U i (0)) ⊂ H q , it follows by induction that at any time t:

B r (t + 1) = U r (t+1)L r (t+1)U r (t+1) = U r (t)L r (t)U r (t) ⊤ + 1 2 U s (t)L s (t)U s (t) ⊤ = B r (t) + 1 2 B s (t)
That is, Equation ( 13) is equivalent to Equation [START_REF] Wiesel | Decomposable principal component analysis[END_REF]. Convergence of Equations (6-7) to the network average at all nodes has already been proved [START_REF] Kempe | Gossip-based computation of aggregate information[END_REF][START_REF] Iutzeler | Analysis of sum-weight-like algorithms for averaging in wireless sensor networks[END_REF]. For completeness though, we show that Equations (6-7) drive all C i to C. Let us introduce the Euclidean errors E i between covariances locally reconstructed by our protocol and the exact covariance matrix:

∀i, E i (t) = C i (t) -C 2 F = 1 w i (t) B i (t) - 1 w i (t) 2 a i (t)a i (t) ⊤ - 1 n XX ⊤ + µµ ⊤ 2 F ≤ 1 w i (t) B i (t) - 1 n XX ⊤ 2 F + 1 w i (t) 2 a i (t)a i (t) ⊤ -µµ ⊤ 2 F
From Equations ( 9), and (3-4), we have n = i w i (0), XX ⊤ = i B i (0) and µ = ( i a i (0))/( i w i (0)). Then, we obtain

E i (t) ≤ B i (t) w i (t) - j B j (0) j w j (0) 2 F + a i (t)a i (t) ⊤ w i (t) 2 - ( j a j (0))( j a j (0)) ⊤ ( j w j (0)) 2 2 F (26)
According to Equations (6-7), each single entry of (B i , a i , w i ) is identically updated. The following Lemma then shows that every quotient of B i or a i and w i tends to its network average (the proof is deferred to Appendix A):

Lemma 2. Under update rules in Equations [START_REF] Kempe | A decentralized algorithm for spectral analysis[END_REF][START_REF] Wiesel | Decomposable principal component analysis[END_REF]:

∀i, lim t→∞ B i (t) w i (t) = j B j (0) j w j (0)
and lim

t→∞ a i (t) w i (t) = j a j (0) j w j (0) (27)
In turn, it implies that both terms in Equation (26) tend to zero, and therefore entails the convergence of local estimates to the covariance matrix:

lim t→∞ C i (t) = C (28)

Experiments

In this section, we experimentally evaluate AGPCA, both on synthetic and natural data. Experiments for Distributed Samples and Distributed Coordinates scenarios are reported separately. Synthetic data were obtained by sampling n observations from a D-dimensional Gaussian distribution, which covariance matrix was arbitrarily generated with rank p ≪ D and such that all dimensions are correlated in R D . For natural data, we used the MNIST handwritten digits dataset, which contains 60000 grayscale images of 28 × 28 pixels, that is, D = 784 [START_REF] Lecun | Gradient-based learning applied to document recognition[END_REF]. Unless specified, the network is assumed fully connected, sender-receiver pairs being uniformly selected at random.

Distributed Samples scenario

In DS scenario, the objective of AGPCA is two-fold : (i) providing all nodes with the same projection matrix U, (ii) making U as close as possible to the solution U ⋆ of a centralized PCA solution computed over the aggregated network data X. We thus measure the average Euclidean reconstruction error between the locally reconstructed covariances C i and the exact C, defined by

E = C i (∞) - C 2 F / C 2 F .
We considered three cases: q = D, q = p, and q < p: • When q = D, the centralized PCA solution U ⋆ reconstructs the covariance with zero error. In this case, experiments on synthetic data show that AGPCA asymptotically provides all nodes with the exact U ⋆ . Figure 1 illustrates the convergence to this optimum versus time and show that both Early and Late PCA strategies allow a perfect reconstruction.

• When q = p ≪ D, Figure 2 confirms the theoretical result of Theorem 1, as convergence to the optimum still holds if the intrinsic dimension p of X is lower than q.

• By contrast, when q < p, the ideal projection itself accounts for less than 100% of the variance in the data. In such case, AGPCA still guarantees that all nodes converge to the same projection matrix U. This ensures that any observation is projected into the same output space whichever the node we consider. However, the error induced by the Orthogonal Iterations steps leads to a slight deviation from the optimal global PCA solution U ⋆ . Experiments reported in Figure 3 for various q reveal that even for high compression rates, the deviation from the optimal basis is still very low.

Similar experiments were conducted on the MNIST dataset, whose intrinsic dimensionality was shown in previous works to be much lower than D = 784. Figure 4 highlights that AGPCA achieves the same reconstruction error as a centralized PCA up to the numerical precision for q ≥ 75. For lower values, it is slightly suboptimal, but the deviation is kept under 1% until q = 3, and always under 2%. Figure 4 also shows that the convergence rate of the aggregation protocol is not impacted by switching from a Late PCA to an Early PCA scheme. Interestingly, lower values for q bring faster convergence.

Concerning communication costs, the efficiency of AGPCA is illustrated in Figure 5. This figure displays the number of messages each node has to emit in order to reach convergence (convergence is assumed when E improves by less than 0.01% between two message events). This number of messages per node appears to logarithmically scale with the number of nodes in the network. This ensures easy scaling to large networks.

Distributed Coordinates scenario

In a DC scenario, our accuracy criterion is still the Euclidean reconstruction error, but it is now computed as E = Y T i Y i -X T X . To evaluate AGPCA in such scenarios, we spread a synthetic dataset similar to Figure 1 by assigning a uniformly drawn number of dimensions D i to each node i, such that node 1 gets the D 1 first dimensions, node 2 gets the next D 2 and so forth, and such that i D i = D = 10000. In this scenario, Figure 6 shows that AGPCA behaves identically to the DS case, which is a logical consequence of using the same strategy.

Influence of network connectivity

Network topology usually has a great impact on speed, accuracy or even applicability of distributed learning algorithms. Gossip protocols have long been acclaimed for their fine behavior in a wide range of networking situations. It is worth highlighting that the Sum-Weight protocol used in AGPCA can be used to build asynchronous generalizations of broadcast, spanning tree, workers-master and scale-free connectivity:

• In an asynchronous broadcast scheme, senders still waken independently, but the same message is sent to all their neighbors. To obtain a broadcast protocol, we just need to replace the 1/2 coefficient in (6-7) by 1/N to respect mass conservation.

• In an asynchronous spanning tree communication scheme, the network has only N -1 links, which is the minimal number of links such that the network stays connected. In our experiments, we used the worst conditioned nondegenerate tree topology where nodes have at most 3 neighbors.

• In a workers-master scheme, N -1 nodes called workers hold local datasets and can only send messages to a single master node. This master is allowed to send messages to all workers, and may hold a local sample or not (if not, its estimates are initialized to 0). In our experiments, we considered two settings. In the first one, the workers and the master are selected for emission with the same probability, i.e., the master works at the same frequency as the workers. In the second one, the master is selected N times more often than workers, e.g., assuming it is a high-throughput server. Here, q = p = 30 (same data as Figure 6).

• An example of a scale-free network is the Barabási-Albert (BA) random graph model, which has a power-law distribution of node degrees. Nodes are connected to random neighbors such that the number of nodes having k neighbors scales in k -3 .

We run AGPCA for each of these connectivity models, varying the number of involved nodes (synthetic data similar to Figure 1). Figure 7 displays the number of messages per node to reach convergence (in the same sense as Figure 5) and compares it to the fully-connected setup we considered so far. Broadcast, workers-master and spanning tree schemes exhibit a poor scaling which appears linear with the number of nodes. By contrast, full and BA connectivity show a logarithmic dependence of the communication costs on the network size. In the workers-master case, we can recover a logarithmic scaling by allowing the master to communicate N -times faster than the workers (in practice, such performances at the master can be unrealistic).

The poor behavior of the broadcast scheme can be explained by remarking that senders always emit N -1 messages containing identical estimates. In the point-to-point case each message integrates the contribution of all preceding updates, thus making estimates mixing much faster. Concerning the spanning tree topology, its bad mixing properties have long been studied in the Gossip protocols and rumor spreading litterature. Due to its high average path length, information must flow, on average, through a much larger number of intermediary nodes to transit between any two peers.

To summarize, AGPCA, like asynchronous Gossip protocols in general, are best-suited for networks with maximally randomized communications and lowest average path lengths.

Conclusion

We presented an asynchronous and decentralized algorithm to solve PCA when data is spread over a network. Based on the integration of a dimensionality reduction operator into a Sum-Weight gossip averaging protocol, it is best suited for large setups with high extrinsic dimension, massive samples and large networks. Unlike other algorithms, it is applicable both in Distributed Samples and Distributed Coordinates scenarios, thanks to the duality between covariance and Gram matrices decompositions. Our theoretical and experimental studies show that it is formally equivalent to running a traditional PCA when the complete data has an intrinsic dimension lower than the output dimension, otherwise providing a low-error approximation of the optimum.

Perspectives include application to large-scale dimension reduction problems, where traditional methods are unapplicable, such as signatures compression in web-scale multimedia retrieval. Besides, while we considered a static scenario where data is provided all at once, Gossip protocols also enjoy nice dynamics when data and/or connectivity are time-evolving. Our approach could then be extended to deal with time-related phenomena, such as concept drift and dynamic networks. product of t realizations tends to have identical rows as t grows. This can be proved by introducing the coefficient of ergodicity τ (•) defined by ∀A ∈ R N ×N , τ (A) = 1 2 max i,j

A T (e ie j ) 1

In words, τ (A) corresponds to the maximum L 1 distance between any two rows of A. A has equal rows if and only if τ (A) = 0. When P(t) is the forward product of t realizations of a scrambling random matrix K, [START_REF] Seneta | Non-negative matrices and Markov chains[END_REF] states τ (P(t + 1)) = τ (P(t)K(t + 1)) ≤ τ (P(t))λ 2 (t),

where λ 2 (t) is the second largest eigenvalue of K(t + 1). Then, τ (P(t)) ≤ λ t 2 , where λ 2 is the maximal second largest eigenvalue over all realizations of K. K(t) being row-stochastic, we have 0 < λ 2 < 1. Hence, lim t→∞ τ (P(t)) = 0

Recalling that x ⊤ (t) = x ⊤ (0)P(t), if P(t) has equal rows we get ∀i,

x i (t) w i (t) = j x j (0)P ji (t) j w j (0)P ji (t) = P 1,i (t) j x j (0) P 1,i (t) j w j (0)

Simplifying P 1i (t) yields our final result: ∀i, lim t→∞ x i (t) w i (t) = j x j (0) j w j (0)

Figure 1 :

 1 Figure 1: Convergence of AGPCA compared to a Late PCA strategy and the centralized PCA solution on a synthetic dataset drawn from a Gaussian distribution (D = 200, p = 30, n = 10000), spread on N = 100 nodes. Here q is set to D.

Figure 2 :

 2 Figure 2: Convergence of AGPCA on the same data and network as Figure 1, but setting q = 30 = p ≪ D.

Figure 3 :Figure 4 :

 34 Figure 3: Convergence of AGPCA on the same data and network as Figure 1, but setting q = 10 < p.

Figure 5 :

 5 Figure5: Average number of message emissions per node required to converge against the number of nodes in the network, evaluated on MNIST with q = 50.

Figure 6 :

 6 Figure 6: Convergence of AGPCA in a Distributed Coordinates scenario on synthetic data similar to Figure 1 (D = 10000, n = 200, p = 30). Each node holds a number of dimensions of the complete data drawn uniformly in {1, . . . , 2DN -1}. Here, q = p = 30.

Figure 7 :

 7 Figure 7: Number of messages per node to reach convergence for various connectivity and networking models.Here, q = p = 30 (same data as Figure6).

When implementing AGPCA, care should be taken that iterations of Algorithm 1 and Algo-

rithm 2 be mutually exclusive, to ensure integrity of the concurrently updated variables.

Acknowledgements

This work is funded by the Culture 3D Cloud project as part of the french Funds for Digital Societies.

Appendix A. Proof of Lemma 2 PROOF. Let's independently consider any entry of B i or a i , denoted using the free variable x i . Let x = (x i ) i gather all the corresponding entries from each node i. Using the definition of x, we can rewrite Equations (6-7) into a matrix form:

x(t + 1) ⊤ = x(t) ⊤ K(t) where K(t) ≡ I + 1 2 e s (e re s ) ⊤ (A.1) K(t) is a random N × N matrix where (s, r) is the selected sender-receiver pair at time t. We can then write x(t) ⊤ = x(0) ⊤ P(t), where P(t) ≡ u≤t K(u).

Observe that ∀t, K(t) ⊤ 1 = 1 and P(t) ⊤ 1 = 1, meaning that K(t) and P(t) are row-stochastic. A row-stochastic matrix that is irreducible is said to be scrambling (see [START_REF] Rhodius | On the maximum of ergodicity coefficients, the dobrushin ergodicity coefficient, and products of stochastic matrices[END_REF]). Establishing irreducibility of K amounts to check if the network graph is strongly connected, i.e., for any two nodes (i, j) there exists a route from i to j through available connections. Assuming bidirectional communication is allowed between any pair of connected nodes, the only constraint to ensure K is irreducible is that the network graph is connected, which is a logical assumption in the distributed computing setup we consider. Consequently, K is a scrambling matrix. Scrambling matrices are shown to be weakly ergodic [START_REF] Seneta | Non-negative matrices and Markov chains[END_REF], that is, the


Jerome Fellus

David Picard

Philippe-Henri Gosselin

Asynchronous decentralized convex optimization through short-term gradient averaging

published or not. The documents may come

Asynchronous decentralized convex optimization through short-term gradient averaging

Jerome Fellus, David Picard,

Introduction

Recent successes of stochastic optimization methods based on averaged gradients have renewed the interest for gradient-based schemes to solve large convex problems [START_REF] Nesterov | Primal-dual subgradient methods for convex problems[END_REF][START_REF] Roux | A stochastic gradient method with an exponential convergence rate for finite training sets[END_REF]. However, very large problems such as training a SVM [START_REF] Bottou | Large-scale machine learning with stochastic gradient descent[END_REF] with Terabytes of input data are still computationally challenging. Meanwhile, the burgeoning of data-collecting devices (e.g., mobile phones, sensor networks, etc) leads to a logical spread of potential input data into many remote collections that can rarely be gathered into a single site. These computational and massive data spread challenges raise the need for distributed optimization algorithms. In spite of at least three decades of intensive research in this field [START_REF] Tsitsiklis | Problems in decentralized decision making and computation[END_REF], massively distributed optimization algorithms have proved inescapable in applicative contexts only a few years ago.

Stochastic Gradient Descent schemes (SGD) [START_REF] Bottou | Large-scale machine learning with stochastic gradient descent[END_REF], where only one training example is randomly drawn at each optimization step, have unlocked the problematic complexity of the original full-gradient descent method. Unfortunately, computing gradients on single examples entails a high instability that must be compensated by a decreasing step size to ensure convergence. Nesterov [START_REF] Nesterov | Primal-dual subgradient methods for convex problems[END_REF] noticed that SGD integrated new gradients with smaller weights than old ones, thus explaining its poor convergence rate in O( ln(t) √ t ). To circumvent this drawback, he proposed the Dual Averaging method (DA) where all past gradients are averaged and taken into account with the same weight when performing the update step. Although this leads to an improved convergence rate in O( 1 √ t ), this scheme still exhibits sub-linear convergence. In 2012, Le Roux et al [START_REF] Roux | A stochastic gradient method with an exponential convergence rate for finite training sets[END_REF] introduced the Stochastic Averaged Gradient method (SAG), achieving a nice exponential convergence rate in O(ρ t ), ρ < 1 (a rate comparable to the deterministic full gradient scheme) for strongly-convex objectives with a constant step size. They built on the intuition that, in case of finite datasets, only the most recent gradient with respect to each training vector should be taken into account in the averaged gradient (instead of all past ones in DA). Exponential convergence is then neither achieved when averaging one nor all previous gradients, but rather when taking into account some recent gradients and discarding older ones. Unfortunately, SAG's convergence proof [START_REF] Schmidt | Minimizing finite sums with the stochastic average gradient[END_REF] does not easily extends to the case where the number of averaged gradients is different from the number of training vectors. Whether exponential or quasi-exponential convergence still holds when averaging an arbitrary number of gradients remains an interesting open question. In this paper, we tackle it from an empirical perspective.

The main contribution of this paper is two-fold : We first introduce a Short-Term Averaged Gradient scheme to optimize a convex objective function, where the number of averaged gradients explicitly appears as a parameter. Secondly, we extend this scheme to distributed environments, and in particular decentralized networking contexts, using an asynchronous gossip averaging protocol. After formally stating our problem, we introduce our Short-Term Average gradient scheme and Gossip averaging. We then present our proposed Asynchronous Gossip Short-Term Averaged Gradient algorithm and report experimental results in centralized and distributed setups, before concluding.

Problem statement

Given a strongly-connected network G of N nodes, each node i holding a sample

X i = [x (i) 1 . . . x (i) ni ], x (i) j ∈ R D ,
our goal is to compute the unique minimizer y ⋆ ∈ R Q of the following strongly-convex objective function f :

∀y ∈ R Q , f (y) = N i=1 ni j=1 f (i) j (y), (1)
where each function f

(i)
j is a convex function associated with x (i) j . For this, each node i starts with a random local estimate y (i) and updates it using node-local computations on its own sample X (i) and communication with neighboring nodes. We aim at making all y (i) converge to y ⋆ . In addition, three constraints must be respected:

C1. No sample exchange -Only estimates and/or gradients can be exchanged.

C2. Asynchrony -Nodes must never wait for each other.

C3. Decentralization -All nodes and links must play the same role (i.e., all nodes run the same algorithm and links are randomly selected).

Short-Term Stochastic Averaged Gradient

First, consider the centralized case (N = 1, thus ( 1) is a traditional convex optimization problem). This problem can be solved using a Stochastic Averaged Gradient Descent scheme, that is, with the following general update rule:

∀i, g(t) = ∇f j (y(t)), j randomly drawn in {1 . . . n} (2)
y(t + 1) = y(t) -η t t τ =1 α t-τ +1 g(τ ), (3)
where (η t ) is a non-increasing sequence of step sizes and (α t-τ +1 ) t τ =1 a sequence of t weights associated with the t gradients seen so far (the most recent gradient takes weight α 1 and the oldest one takes α t ). (2-3) encompasses several existing algorithms:

Clearly, SGD [START_REF] Bottou | Large-scale machine learning with stochastic gradient descent[END_REF] corresponds to setting α 1 = 1 and ∀k > 1, α k = 0 since only the most recent gradient is kept. Variants integrating a momentum term [START_REF] Tseng | An incremental gradient (-projection) method with momentum term and adaptive stepsize rule[END_REF] with fixed weight β also fit to (2-3), with α k = β k .

In its initial statement, SAG [START_REF] Roux | A stochastic gradient method with an exponential convergence rate for finite training sets[END_REF] does not fit to (2-3), as it keeps the most recent gradient for each training example. If training vectors are uniformly drawn with replacement, these are not necessarily the same as the last n gradients. However, as highlighted in [START_REF] Shalev-Shwartz | Stochastic dual coordinate ascent methods for regularized loss[END_REF], the most efficient strategy to draw training examples is generally to select them without replacement. This can be implemented by performing deterministic passes over the full data and shuffling the dataset at the begin of each epoch. With this selection strategy, SAG is equivalent to (2-3) with α k = 1 n if k ≤ n and 0 otherwise. Dual Averaging [START_REF] Nesterov | Primal-dual subgradient methods for convex problems[END_REF] is a proximal method, i.e., y(t + 1) is not computed from y(t) but rather by directly projecting the averaged gradient using a proximal mapping. However, for simple (e.g., quadratic) proximal functions, one can derive from it an averaged gradient descent (2-3) where all gradients are given an equal weight α k = 1, ∀k ∈ {1 . . . t}.

SAG enjoys the fastest convergence rate of the above-cited algorithms. On the other hand, SAG memorizes as many gradients as training examples, resulting in a storage cost in O(nD), which can hardly be afforded in large scale setups. A natural idea is then to design an intermediary scheme, where only the M most recent gradients are kept in memory and averaged, where M is a user-defined parameter. We call this scheme Short-Term Averaged Gradient (STAG). Intuitively, for M close to 1, we get a behavior similar to SGD, and for M close to n, we get a SAG-like convergence.

Asynchronous gossip averaging

Now, consider a distributed case, that is, N > 1. To extend STAG to this setup, we would ideally provide all nodes with the average of the M most recent gradients seen over the network for each iteration. Obviously, C2 and C3 prevent us from computing this average by gathering local estimates in a master node. Instead, we rely on an asynchronous gossip averaging protocol, which is a special kind of Consensus protocols [START_REF] Bénézit | Weighted gossip: Distributed averaging using non-doubly stochastic matrices[END_REF]. Consensus protocols solve distributed averaging over a network by iteratively averaging estimates between neighboring nodes. Assuming each node i holds an arbitrary initial estimate s i (0) and a weight w i (0), and setting s(t) = (s 1 (t), . . . , s N (t)) and w(t) = (w 1 (t), . . . , w N (t)), consensus protocols compute the network average µ = N i si(0) N i wi(0) using the following update:

s(t + 1) ⊤ = s(t) ⊤ K(t + 1), w(t + 1) ⊤ = w(t) ⊤ K(t + 1), (4)
where K(t) is a constant primitive doubly stochastic matrix respecting G (i.e., (K(t)) ij = 0 only if i and j are connected in G). For constant K(t) = K, we have s(t) ⊤ = s(0) T K t and w(t) ⊤ = w(0) T K t . By Perron-Frobenius theorem and stochasticity of K (i.e., K1 = 1), the rows of K t tend to equality at rate O(λ t 2 ) where λ 2 is the sub-dominant eigenvalue of K. Consequently,

lim t→∞ s i (t) w i (t) = lim t→∞ N j s j (0)(K t ) ji N j w j (0)(K t ) ji = lim t→∞ (K t ) 1,i N j s j (0) (K t ) 1,i N j w j (0) = µ (5)
Consensus protocols have been successfully used in the distributed extension of Dual Averaging (DDA [START_REF] Duchi | Dual averaging for distributed optimization: convergence analysis and network scaling[END_REF]). However, constant and doubly-stochastic consensus matrix K are not compatible with C2 and C3. Hopefully, as shown in [START_REF] Bénézit | Weighted gossip: Distributed averaging using non-doubly stochastic matrices[END_REF], both assumptions can be relaxed : any random row-stochastic K(t) still satisfy (5) as long as t K(t) is primitive. Such randomized protocols, called asynchronous Gossip protocols, were recently applied to distributed SGD (for SVM training [START_REF] Hensel | Gadget svm: a gossip-based sub-gradient svm solver[END_REF]) and Dual Averaging (PS-DDA [START_REF] Tsianos | Push-sum distributed dual averaging for convex optimization[END_REF]). We define our asynchronous Gossip protocol by K(t) = 1 2 e i (e ⊤ ie ⊤ j ) where (i, j) is uniformly drawn among the edge set of G. This update rule appears very simple: at any time t, a sender node i is uniformly drawn, and sends half of its estimates to a uniformly chosen neighbor j which adds it to is own estimate:

s i (t + 1) = 1 2 s i (t) w i (t + 1) = 1 2 w i (t) , s j (t + 1) = s j (t) + 1 2 s i (t) w j (t + 1) = w j (t) + 1 2 w i (t) (6)
Remark that, contrarily to (4) that can require simultaneous messages between multiple nodes, update ( 6) only involves a single random link at a time, thus allowing a fully asynchronous (C2) and decentralized (C3) functioning.

The AGSTAG algorithm

In this section, we describe our proposed distributed optimization algorithm called AGSTAG (Asynchronous Gossip Short-Term Averaged Gradient). It is built on the STAG scheme described in section 3, associated with the asynchronous Gossip averaging protocol defined by [START_REF] Lecun | Gradient-based learning applied to document recognition[END_REF]. In AGSTAG, every nodes run the same local procedure (Algorithm 1). Each node i repeatedly perform the following steps without any synchronization with other nodes. First, it selects a training example (without replacement) to evaluate a new gradient. This gradient is added to a node-local gradients buffer in replacement of the oldest stored one, following a FIFO rule. This buffer thus always stores the M most recently evaluated local gradients. The difference between the oldest gradient and the newest one, respectively leaving and entering the buffer, is added to a local estimate z, while w counts the number of gradients currently stored in the buffer.

(z, w) is then iteratively exchanged with random neighboring nodes according to [START_REF] Lecun | Gradient-based learning applied to document recognition[END_REF].

Remark that this algorithm is asynchronous because there is no message back from the receiver to the sender (senders emit independently of received messages). Notice that the updates of z(t) imply that the global sum of the dual variables N i z (i) (t) is always equal to the global sum of the gradients stored in all queues. By [START_REF] Hensel | Gadget svm: a gossip-based sub-gradient svm solver[END_REF], during the course of AGSTAG, the local dual variables z (i) (t) behave as estimators of the network-wide average of all currently stored gradients, thus mimicking a centralized STAG, with the difference that there are now N M gradients globally averaged in the network instead of M . Consequently, M can be set N -times smaller to get the same behavior as the centralized scheme, resulting in a much lower memory cost.

Algorithm 1 AGSTAG (independently run at each networking node i) Parameters : (η t ) : sequence of step sizes. M : capacity of the gradients buffer.

S : number of messages sent at each iteration Init : y ∈ R Q at random ; z ∈ R Q ← 0 ; w ← 0 ; FIFO ← empty with M slots loop g ← ∇f

Experiments

We now present experimental results for both the STAG scheme and its distributed AGSTAG extension on the MNIST dataset [START_REF] Lecun | Gradient-based learning applied to document recognition[END_REF] (D = 784, n = 60000). Both algorithms are evaluated in a binary L 2 linear SVM training task (similar setup as [START_REF] Bottou | Large-scale machine learning with stochastic gradient descent[END_REF] with λ = 0.1). We only report results for the best-performing constant step-size and a complete network connectivity. The number of messages per iteration S was fixed to 10. Fig. 1 displays the convergence time of STAG in terms of the number of gradient evaluations needed to reach an ε-accurate minimizer to the objective f defined in [START_REF] Bénézit | Weighted gossip: Distributed averaging using non-doubly stochastic matrices[END_REF]. Observe that STAG still enjoys a low convergence time when M (the number of averaged gradient) is only half the size of the dataset. It also behaves well for small values of M and ε = 10 -2 because instability is contained within this threshold, but to obtain convergence under ε = 10 -3 we must use very low step sizes, implying much higher convergence times. This is not the case for large values of M , where STAG fastly converges up to any arbitrarily-low ε.

Fig. 2 shows the convergence slope of AGSTAG for various network sizes N , where training examples are uniformly spread over the nodes. At each node, we chose M = n i /10, that is, local gradients buffers are 90% smaller than local datasets. We can see that using a large network can drastically reduce the number of gradient evaluations per node. Besides, we observed that larger networks allowed greater constant step-sizes.

Conclusion

We introduced an asynchronous algorithm for distributed convex optimization called AGSTAG, that is made of two original parts: a Short-Term averaged gradient (STAG) optimization scheme and an asynchronous and decentralized Gossip averaging protocol. Unlike other methods, STAG turns the number of averaged gradients into a parameter, improving flexibility. The proposed asynchronous Gossip averaging protocol offers a simple and efficient communication mechanism to extend STAG to large scale distributed environments, as shown in our experimental evaluation.

jFig. 1 :Fig. 2 :

 12 Fig. 1: Convergence time of STAG versus gradient buffer size M

* This work is funded by the Culture 3D Cloud project.