Revisiting Randomized Gossip Algorithms: General Framework, Convergence Rates and Novel Block and Accelerated Protocols

TL;DR

This paper introduces a comprehensive framework for analyzing and designing randomized gossip algorithms, unifying existing methods and proposing new accelerated and block protocols with proven convergence rates.

Contribution

It develops a general framework that encompasses classical gossip algorithms, enabling the creation of faster, block, and accelerated protocols with theoretical guarantees.

Findings

Proposed novel block and accelerated gossip protocols.

Unified analysis of existing gossip algorithms within a new framework.

Extensive experimental validation on wireless network topologies.

Abstract

In this work we present a new framework for the analysis and design of randomized gossip algorithms for solving the average consensus problem. We show how classical randomized iterative methods for solving linear systems can be interpreted as gossip algorithms when applied to special systems encoding the underlying network and explain in detail their decentralized nature. Our general framework recovers a comprehensive array of well-known gossip algorithms as special cases, including the pairwise randomized gossip algorithm and path averaging gossip, and allows for the development of provably faster variants. The flexibility of the new approach enables the design of a number of new specific gossip methods. For instance, we propose and analyze novel block and the first provably accelerated randomized gossip protocols, and dual randomized gossip algorithms. From a numerical analysis viewpoint, our work is the first that explores in depth the decentralized nature of randomized iterative methods for linear systems and proposes them as methods for solving the average consensus problem. We evaluate the performance of the proposed gossip protocols by performing extensive experimental testing on typical wireless network topologies.

Tables

Table 1. Table 1 : Frequently used notation.

The Basics
$𝐀,b$	$m\times n$ matrix and $m\times 1$ vector defining the system $𝐀x=b$
$ℒ$	$\{x:𝐀x=b\}$ (solution set of the linear system)
$𝐁$	$n\times n$ symmetric positive definite matrix
${⟨x,y⟩}_{𝐁}$	$x^{\top }𝐁y$ ($𝐁$-inner product)
${\Vert x\Vert }_{𝐁}$	$\sqrt{{⟨x,x⟩}_{𝐁}}$ ($𝐁$-norm)
${𝐌}^{†}$	Moore-Penrose pseudoinverse of matrix $𝐌$
$𝐒$	a random real matrix with $m$ rows
$𝒟$	distribution from which matrix $𝐒$ is drawn ($𝐒\sim 𝒟$)
$𝐇$	$𝐒{\left({𝐒}^{\top }{\mathrm{𝐀𝐁}}^{-1}{𝐀}^{\top }𝐒\right)}^{†}{𝐒}^{\top }$
$𝐙$	${𝐀}^{\top }\mathrm{𝐇𝐀}$
$\mathrm{Range}\left(𝐌\right)$	range space of matrix $𝐌$
$\mathrm{Null}\left(𝐌\right)$	null space of matrix $𝐌$
$\mathbb{P}(\cdot )$	probability of an event
$𝔼[\cdot ]$	expectation
Projections
${\Pi }_{ℒ,𝐁}\left(x\right)$	projection of $x$ onto $ℒ$ in the $𝐁$-norm
${𝐁}^{-1}𝐙$	projection matrix, in the $𝐁$-norm, onto $\mathrm{Range}\left({𝐁}^{-1}{𝐀}^{\top }𝐒\right)$
Graphs
$𝒢=(𝒱,ℰ)$	an undirected graph with vertices $𝒱$ and edges $ℰ$
$n$	$=\left|𝒱\right|$ (number of vertices)
$m$	$=\left|ℰ\right|$ (number of edges)
$e=(i,j)\in ℰ$	edge of $𝒢$ connecting nodes $i,j\in 𝒱$
$d_i$	degree of node $i$
$c\in {ℝ}^n$	$=(c_1,\dots ,c_n)$; a vector of private values stored at the nodes of $𝒢$
$\overline{c}$	$\overline{c}=\frac{{\sum }_i^n{𝐁}_{ii}{c}_i}{{\sum }_i^n{𝐁}_{ii}}$ (the weighted average of the private values)
$𝐐\in {ℝ}^{m\times m}$	Incidence matrix of $𝒢$
$𝐋\in {ℝ}^{n\times n}$	$={𝐐}^{\top }𝐐$ (Laplacian matrix of $𝒢$)
$𝐃\in {ℝ}^{n\times n}$	$=\text{Diag}(d_1,d_2,\dots ,d_n)$ (Degree matrix of $𝒢$)
${𝐋}^{rw}\in {ℝ}^{n\times n}$	$={𝐃}^{-1}𝐋$ (random walk normalized Laplacian matrix of $𝒢$)
${𝐋}^{sym}\in {ℝ}^{n\times n}$	$={𝐃}^{-1/2}{\mathrm{𝐋𝐃}}^{-1/2}$ (symmetric normalized Laplacian matrix of $𝒢$)
$\alpha \left(𝒢\right)$	$={\lambda }_{\min }^{+}\left(𝐋\right)$ (algebraic connectivity of $𝒢$)
Eigenvalues
$𝐖$	${𝐁}^{-1/2}𝔼\left[𝐙\right]{𝐁}^{-1/2}$ (psd matrix)
${\lambda }_1,\dots ,{\lambda }_n$	eigenvalues of $𝐖$
${\lambda }_{\max },{\lambda }_{\min }^{+}$	largest and smallest nonzero eigenvalues of $𝐖$
Algorithms
$\omega$	relaxation parameter / stepsize
$\beta$	heavy ball momentum parameter
$\rho$	$1-\omega \left(2-\omega \right){\lambda }_{\min }^{+}$

Equations

$$\mathbf{H}:=\mathbf{S}({\mathbf{S}}^{\top }\mathbf{A}{\mathbf{B}}^{-1}{\mathbf{A}}^{\top }\mathbf{S}{)}^{†}{\mathbf{S}}^{\top }\text{and}\mathbf{Z}:={\mathbf{A}}^{\top }\mathbf{H}\mathbf{A},$$

$$\mathbf{W}={\mathbf{B}}^{-1/2}{\mathbf{A}}^{\top }\mathbb{E}[\mathbf{H}]\mathbf{A}{\mathbf{B}}^{-1/2}\stackrel{\text{eqref}ZETA}{=}{\mathbf{B}}^{-1/2}\mathbb{E}[\mathbf{Z}]{\mathbf{B}}^{-1/2},$$

$${\Pi }_{\mathcal{L},\mathbf{B}}(x):=\arg \underset{x^{\prime }\in \mathcal{L}}{\min }\Vert x^{\prime }-x{\Vert }_{\mathbf{B}}=x-{\mathbf{B}}^{-1}{\mathbf{A}}^{\top }(\mathbf{A}{\mathbf{B}}^{-1}{\mathbf{A}}^{\top }{)}^{†}(\mathbf{A}x-b).$$

$$x^{k+1}$$

$${\mathbf{S}}_k^{\top }\mathbf{A}x={\mathbf{S}}_k^{\top }b.$$

$$x^{k+1}=x^k-\omega \frac{{\mathbf{A}}_{i:}{x}^k-b_i}{\Vert {\mathbf{A}}_{i:}{\Vert }^2}{\mathbf{A}}_{i:}^{\top }.$$

$$x^{k+1}=x^k-\omega {\mathbf{A}}_{C:}^{\top }({\mathbf{A}}_{C:}{\mathbf{A}}_{C:}^{\top }{)}^{†}({\mathbf{A}}_{C:}{x}^k-b_C).$$

$$\underset{x=(x_1,\dots ,x_n)\in {\mathbb{R}}^n}{\min }\frac{1}{2}\Vert x-x^0{\Vert }_{\mathbf{B}}^2\text{subject to}\mathbf{A}x=b.$$

$$\mathrm{Null}({\mathbb{E}}_{\mathbf{S}\sim \mathcal{D}}[\mathbf{Z}])=\mathrm{Null}(\mathbf{A}).$$

$$\mathbb{E}[\Vert x^k-x^{\ast }{\Vert }_{\mathbf{B}}^2]\le {\rho }^k\Vert x^0-x^{\ast }{\Vert }_{\mathbf{B}}^2,$$

$$\rho :=1-\omega (2-\omega ){\lambda }_{\min }^{+}\in [0,1].$$

$$k\ge \frac{1}{1-\rho }\log \left(\frac{1}{ϵ}\right)\Rightarrow \mathbb{E}[\Vert x^k-x^{\ast }{\Vert }_{\mathbf{B}}^2]\le ϵ\Vert x^0-x^{\ast }{\Vert }_{\mathbf{B}}^2.$$

$$O\left(\frac{1}{1-\rho }\log \left(\frac{1}{ϵ}\right)\right).$$

$$\underset{x\in {\mathbb{R}}^d}{\min }\frac{1}{n}\sum _{i=1}^n{f}_i(x)$$

$$\bar{c}:=\frac{{\sum }_{i=1}^n{w}_i{c}_i}{{\sum }_{i=1}^n{w}_i},$$

$$\underset{x=(x_1,\dots ,x_n)\in {\mathbb{R}}^n}{\min }\frac{1}{2}\Vert x-c{\Vert }_{\mathbf{B}}^2\text{subject to}{x}_1=x_2=\cdots =x_n$$

$$x^{k+1}=\left[\mathbf{I}-\omega {\mathbf{A}}^{\top }{\mathbf{H}}_k\mathbf{A}\right]x^k=\left[\mathbf{I}-\omega {\mathbf{Z}}_k\right]x^k.$$

$$\frac{1}{n}{1}_n^{\top }{x}^{k+1}=\frac{1}{n}{1}_n^{\top }(\mathbf{I}-\omega {\mathbf{A}}^{\top }{\mathbf{H}}_k\mathbf{A})x^k=\frac{1}{n}{1}_n^{\top }\mathbf{I}{x}^k-\frac{1}{n}{1}_n^{\top }\omega {\mathbf{A}}^{\top }{\mathbf{H}}_k\mathbf{A}{x}^k\stackrel{\mathbf{A}{1}_n=0}{=}\frac{1}{n}{1}_n^{\top }{x}^k.$$

$$T_{ave}(\epsilon ):=\underset{c\in {\mathbb{R}}^n}{\sup }\inf \left\{k:\mathbb{P}\left(z^k>\epsilon z^0\right)\le \epsilon \right\}.$$

$$T_{ave}(ϵ)\le 3\frac{\log (1/\epsilon )}{\log (1/\rho )}\le 3\frac{\log (1/ϵ)}{1-\rho },$$

$$x^{k+1}=x^k-\omega \frac{e_i^{\top }(\mathbf{A}{x}^k-b)}{e_i^{\top }\mathbf{A}{\mathbf{B}}^{-1}{\mathbf{A}}^{\top }{e}_i}{\mathbf{B}}^{-1}{\mathbf{A}}^{\top }{e}_i=x^k-\omega \frac{{\mathbf{A}}_{i:}{x}^k-b_i}{\Vert {\mathbf{B}}^{-1/2}{\mathbf{A}}_{i:}^{\top }{\Vert }_2^2}{\mathbf{B}}^{-1}{\mathbf{A}}_{i:}^{\top }.$$

$$\Vert {\mathbf{B}}^{-1/2}{\mathbf{Q}}_{e:}^{\top }{\Vert }^2=\Vert {\mathbf{B}}^{-1/2}(f_i-f_j){\Vert }^2=\frac{1}{w_1}+\frac{1}{w_j}.$$

$$x^{k+1}$$

$$x_i^{k+1}=\left(1-\omega \frac{w_j}{w_j+w_i}\right)x_i^k+\omega \frac{w_j}{w_j+w_i}{x}_j^k\text{and}{x}_j^{k+1}=\omega \frac{w_i}{w_j+w_i}{x}_i^k+\left(1-\omega \frac{w_i}{w_j+w_i}\right)x_j^k.$$

$$x_i^{k+1}=\left(1-\frac{\omega }{2}\right)x_i^k+\frac{\omega }{2}{x}_j^k\text{and}{x}_j^{k+1}=\frac{\omega }{2}{x}_i^k+\left(1-\frac{\omega }{2}\right)x_j^k.$$

$$x_i^{k+1}=x_j^{k+1}=\frac{x_i^k+x_j^k}{2},$$

$$\Vert {\mathbf{B}}^{-1/2}{\mathbf{L}}_{i:}^{\top }{\Vert }^2={\Vert {\mathbf{B}}^{-1/2}(d_i{f}_i-\sum _{j\in {\mathcal{N}}_i}{f}_j)\Vert }^2=\frac{d_i^2}{w_i}+\sum _{j\in {\mathcal{N}}_i}\frac{1}{w_j}.$$

$$x^{k+1}$$

$$x_i^{k+1}=x_i^k-\frac{\omega \left(d_i{x}_i^k-{\sum }_{j\in {\mathcal{N}}_i}{x}_j^k\right)}{\frac{d_i^2}{w_i}+{\sum }_{j\in {\mathcal{N}}_i}\frac{1}{w_j}}\frac{d_i}{w_i},$$

$$x_j^{k+1}=x_j^k+\frac{\omega \left(d_i{x}_i^k-{\sum }_{\ell \in {\mathcal{N}}_i}{x}_{\ell }^k\right)}{\frac{d_i^2}{w_i}+{\sum }_{\ell \in {\mathcal{N}}_i}\frac{1}{w_{\ell }}}\frac{1}{w_j}.$$

Full text

From a numerical analysis viewpoint, our work is the first that explores in depth the decentralized nature of randomized iterative methods for linear systems and proposes them as methods for solving the average consensus problem.

We evaluate the performance of the proposed gossip protocols by performing extensive experimental testing on typical wireless network topologies.

Keywords randomized gossip algorithms $\cdot$ average consensus $\cdot$ weighted average consensus $\cdot$ stochastic methods $\cdot$ linear systems $\cdot$ randomized Kaczmarz $\cdot$ randomized block Kaczmarz $\cdot$ randomized coordinate descent $\cdot$ heavy ball momentum $\cdot$ Nesterov’s acceleration $\cdot$ duality $\cdot$ convex optimization $\cdot$ wireless sensor networks

**Mathematical Subject Classifications ** 93A14 $\cdot$ 68W15 $\cdot$ 68Q25 $\cdot$ 68W20 $\cdot$ 68W40 $\cdot$ 65Y20 $\cdot$ 90C15 $\cdot$ 90C20 $\cdot$ 90C25 $\cdot$ 15A06 $\cdot$ 15B52 $\cdot$ 65F10

1. 1 Introduction

2. 1.1 Main contributions

3. 1.2 Structure of the paper

4. 1.3 Notation

5. 2 Background - Technical Preliminaries

6. 2.1 Randomized iterative methods for linear systems

7. 2.2 Other related work

8. 3 Sketch and Project Methods as Gossip Algorithms

9. 3.1 Weighted average consensus

10. 3.2 Gossip algorithms through sketch and project framework

11. 3.2.1 Standard form and mass preservation

12. 3.2.2 $\varepsilon$-Averaging time

13. 3.3 Randomized Kaczmarz method as gossip algorithm

14. 3.3.1 AC system with incidence matrix $\mathbf{Q}$

15. 3.3.2 AC system with Laplacian matrix $\mathbf{L}$

16. 3.3.3 Details on complexity results

17. 3.4 Block gossip algorithms

18. 4 Faster and Provably Accelerated Randomized Gossip Algorithms

19. 4.1 Gossip algorithms with heavy ball momentum

20. 4.1.1 Sketch and project with heavy ball momentum

21. 4.1.2 Randomized Kaczmarz gossip with heavy ball momentum

22. 4.1.3 Connections with existing fast randomized gossip algorithms

23. 4.1.4 Randomized block Kaczmarz gossip with heavy ball momentum

24. 4.1.5 Mass preservation

25. 4.2 Provably accelerated randomized gossip algorithms

26. 4.2.1 Accelerated Kaczmarz methods using Nesterov’s momentum

27. 4.2.2 Theoretical guarantees of AccRK

28. 4.2.3 Accelerated randomized gossip algorithms

29. 5 Dual Randomized Gossip Algorithms

30. 5.1 Dual problem and SDSA

31. 5.2 Randomized Newton method as a dual gossip algorithm

32. 6 Further Connections Between Methods for Solving Linear Systems and Gossip Algorithms

33. 7 Numerical Evaluation

34. 7.1 Convergence on weighted average consensus

35. 7.2 Benefit of block variants

36. 7.3 Accelerated gossip algorithms

37. 7.3.1 Impact of momentum parameter on mRK

38. 7.3.2 Comparison of mRK and shift-Register algorithm [40]

39. 7.3.3 Impact of momentum parameter on mRBK

40. 7.3.4 Performance of AccGossip

41. 7.4 Relaxed randomized gossip without momentum

42. 8 Conclusion and Future Direction of Research

43. 9 Acknowledgments

44. A Missing Proofs

45. A.1 Proof of Theorem 4

46. A.2 Proof of Theorem 5

47. B Notation Glossary

1 Introduction

Average consensus is a fundamental problem in distributed computing and multi-agent systems. It comes up in many real world applications such as coordination of autonomous agents, estimation, rumour spreading in social networks, PageRank and distributed data fusion on ad-hoc networks and decentralized optimization. Due to its great importance there is much classical [81, 13] and recent [88, 87, 6] work on the design of efficient algorithms/protocols for solving it.

In the average consensus (AC) problem we are given an undirected connected network ${\cal G}=({\cal V},{\cal E})$ with node set ${\cal V}=\{1,2,\dots,n\}$ and edges ${\cal E}$. Each node $i\in{\cal V}$ “knows” a private value $c_{i}\in\mathbb{R}$. The goal of AC is for every node to compute the average of these private values, $\bar{c}:=\frac{1}{n}\sum_{i}c_{i}$, in a distributed fashion. That is, the exchange of information can only occur between connected nodes (neighbors).

One of the most attractive classes of protocols for solving the average consensus problem are gossip algorithms. The development and design of gossip algorithms was studied extensively in the last decade. The seminal 2006 paper of Boyd et al. [6] motivated a fury of subsequent research and gossip algorithms now appear in many applications, including distributed data fusion in sensor networks [88], load balancing [10] and clock synchronization [19]. For a survey of selected relevant work prior to 2010, we refer the reader to the work of Dimakis et al. [14]. For more recent results on randomized gossip algorithms we suggest [94, 40, 63, 43, 55, 3]. See also [15, 4, 64, 29, 30].

1.1 Main contributions

In this work, we connect two areas of research which until now have remained remarkably disjoint in the literature: randomized iterative (projection) methods for solving linear systems and randomized gossip protocols for solving the average consensus. This connection enables us to make contributions by borrowing from each body of literature to the other and using it we propose a new framework for the design and analysis of novel efficient randomized gossip protocols.

The main contributions of our work include:

• •

RandNLA. We show how classical randomized iterative methods for solving linear systems can be interpreted as gossip algorithms when applied to special systems encoding the underlying network and explain in detail their decentralized nature. Through our general framework we recover a comprehensive array of well-known gossip protocols as special cases. In addition our approach allows for the development of novel block and dual variants of all of these methods. From a randomized numerical linear algebra (RandNLA) viewpoint our work is the first that explores in depth, the decentralized nature of randomized iterative methods for solving linear systems and proposes them as efficient methods for solving the average consensus problem (and its weighted variant).

• •

Weighted AC. The methods presented in this work solve the more general weighted average consensus (Weighted AC) problem (Section 3.1) popular in the area of distributed cooperative spectrum sensing networks. The proposed protocols are the first randomized gossip algorithms that directly solve this problem with finite-time convergence rate analysis. In particular, we prove linear convergence of the proposed protocols and explain how we can obtain further acceleration using momentum. To the best of our knowledge, the existing decentralized protocols that solve the weighted average consensus problem show convergence but without convergence analysis.

• •

Acceleration. We present novel and provably accelerated randomized gossip protocols. In each step, of the proposed algorithms, all nodes of the network update their values using their own information but only a subset of them exchange messages. The protocols are inspired by the recently proposed accelerated variants of randomized Kaczmarz-type methods and use momentum terms on top of the sketch and project update rule (gossip communication) to obtain better theoretical and practical performance. To the best of our knowledge, our accelerated protocols are the first randomized gossip algorithms that converge to a consensus with a provably accelerated linear rate without making any further assumptions on the structure of the network. Achieving an accelerated linear rate in this setting using randomized gossip protocols was an open problem.

• •

Duality. We reveal a hidden duality of randomized gossip algorithms, with the dual iterative process maintaining variables attached to the edges of the network. We show how the randomized coordinate descent and randomized Newton methods work as edge-based dual randomized gossip algorithms.

• •

Experiments. We corroborate our theoretical results with extensive experimental testing on typical wireless network topologies. We numerically verify the linear convergence of the our protocols for solving the weighted AC problem. We explain the benefit of using block variants in the gossip protocols where more than two nodes update their values in each iteration. We explore the performance of the proposed provably accelerated gossip protocols and show that they significantly outperform the standard pairwise gossip algorithm and existing fast pairwise gossip protocols with momentum. An experiment showing the importance of over-relaxation in the gossip setting is also presented.

This paper contains a synthesis and a unified presentation of the randomized gossip protocols proposed in Loizou and Richtárik [46], Loizou and Richtárik [49] and Loizou et al. [45]. In [46], building upon the results from [26], a connection between the area of randomized iterative methods for linear systems and gossip algorithms was established and block gossip algorithm were developed. Then in [49] and [45] faster and provably accelerated gossip algorithms were proposed using the heavy ball momentum and Nesterov’s acceleration technique, respectively. This paper expands upon these results and presents proofs for theorems that are referenced in the above papers. We also conduct several new experiments.

We believe that this work could potentially open up new avenues of research in the area of decentralized gossip protocols.

1.2 Structure of the paper

This work is organized as follows. Section 2 introduces the necessary background on basic randomized iterative methods for linear systems that will be used for the development of randomized gossip protocols. Related work on the literature of linear system solvers, randomized gossip algorithms for averaging and gossip algorithms for consensus optimization is presented. In Section 3 the more general weighted average consensus problem is described and the connections between the two areas of research (randomized projection methods for linear systems and gossip algorithms) is established. In particular we explain how methods for solving linear systems can be interpreted as gossip algorithms when applied to special systems encoding the underlying network and elaborate in detail their distributed nature. Novel block gossip variants are also presented. In Section 4 we describe and analyze fast and provably accelerated randomized gossip algorithms. In each step of these protocols all nodes of the network update their values but only a subset of them exchange their private values. Section 5 describes dual randomized gossip algorithms that operate with values that are associated to the edges of the network and Section 6 highlights further connections between methods for solving linear systems and gossip algorithms. Numerical evaluation of the new gossip protocols is presented in Section 7. Finally, concluding remarks are given in Section 8.

1.3 Notation

For convenience, a table of the most frequently used notation is included in the Appendix B. In particular, with boldface upper-case letters denote matrices; $\mathbf{I}$ is the identity matrix. By $\|\cdot\|$ and $\|\cdot\|_{F}$ we denote the Euclidean norm and the Frobenius norm, respectively. For a positive integer number $n$, we write $[n]:=\{1,2,\dots,n\}$. By ${\cal L}$ we denote the solution set of the linear system $\mathbf{A}x=b$, where $\mathbf{A}\in\mathbb{R}^{m\times n}$ and $b\in\mathbb{R}^{m}$.

We shall often refer to specific matrix expressions involving several matrices. In order to keep these expressions brief throughout the paper it will be useful to define the following two matrices:

[TABLE]

depending on a random matrix ${\bf S}\in\mathbb{R}^{m\times q}$ drawn from a given distribution ${\cal D}$ and on an $n\times n$ positive definite matrix ${\bf B}$ which defines the geometry of the space. In particular we define $\mathbf{B}-$inner product in $\mathbb{R}^{n}$ via $\langle x,z\rangle_{\bf B}:=\langle{\bf B}x,z\rangle=x^{\top}\mathbf{B}z$ and an induced $\mathbf{B}-$norm, $\|x\|_{\bf B}:=(x^{\top}{\bf B}x)^{1/2}$. By $\mathbf{A}_{i:}$ and $\mathbf{A}_{:j}$ we denote the $i^{th}$ row and the $j^{th}$ column of matrix $\mathbf{A}$, respectively. By $\dagger$ we denote the Moore-Penrose pseudoinverse.

The complexity of all gossip protocols presented in this paper is described by the spectrum of matrix

[TABLE]

where the expectation is taken over $\mathbf{S}\sim{\cal D}$. With $\lambda_{\min}^{+}$ and $\lambda_{\max}$ we indicate the smallest nonzero and the largest eigenvalue of matrix $\mathbf{W}$, respectively.

Vector $x^{k}=(x^{k}_{1},\dots,x^{k}_{n})\in\mathbb{R}^{n}$ represents the vector with the private values of the $n$ nodes of the network at the $k^{th}$ iteration while with $x_{i}^{k}$ we denote the value of node $i\in[n]$ at the $k^{th}$ iteration. ${\cal N}_{i}\subseteq{\cal V}$ denotes the set of nodes that are neighbors of node $i\in{\cal V}$. By $\alpha({\cal G})$ we denote the algebraic connectivity of graph ${\cal G}$.

Throughout the paper, $x^{*}$ is the projection of $x^{0}$ onto ${\cal L}$ in the $\mathbf{B}$-norm. We write $x^{*}=\Pi_{{\cal L},\mathbf{B}}(x^{0})$. An explicit formula for the projection of $x$ onto set ${\cal L}$ is given by

[TABLE]

Finally, with $\mathbf{Q}\in\mathbb{R}^{|{\cal E}|\times n}$ we define the incidence matrix and with $\mathbf{L}\in\mathbb{R}^{n\times n}$ the Laplacian matrix of the network. Note that it holds that $\mathbf{L}=\mathbf{Q}^{\top}\mathbf{Q}$. Further, with $\mathbf{D}$ we denote the degree matrix of the graph. That is, $\mathbf{D}=\text{{Diag}}(d_{1},d_{2},\dots,d_{n})\in\mathbb{R}^{n\times n}$ where $d_{i}$ is the degree of node $i\in{\cal V}$.

2 Background - Technical Preliminaries

Solving linear systems is a central problem in numerical linear algebra and plays an important role in computer science, control theory, scientific computing, optimization, computer vision, machine learning, and many other fields. With the advent of the age of big data, practitioners are looking for ways to solve linear systems of unprecedented sizes. In this large scale setting, randomized iterative methods are preferred mainly because of their cheap per iteration cost and because they can easily scale to extreme dimensions.

2.1 Randomized iterative methods for linear systems

Kaczmarz-type methods are very popular for solving linear systems $\mathbf{A}x=b$ with many equations. The (deterministic) Kaczmarz method for solving consistent linear systems was originally introduced by Kaczmarz in 1937 [34]. Despite the fact that a large volume of papers was written on the topic, the first provably linearly convergent variant of the Kaczmarz method—the randomized Kaczmarz Method (RK)—was developed more than 70 years later, by Strohmer and Vershynin [77]. This result sparked renewed interest in design of randomized methods for solving linear systems [56, 57, 18, 51, 93, 58, 74, 41]. More recently, Gower and Richtárik [25] provide a unified analysis for several randomized iterative methods for solving linear systems using a sketch-and-project framework. We adopt this framework in this paper.

In particular, the analysis in [25] was done under the assumption that matrix $\mathbf{A}$ has full column rank. This assumption was lifted in [26], and a duality theory for the method developed. Later, in [73], it was shown that the sketch and project method of [25] can be interpreted as stochastic gradient descent applied to a suitable stochastic optimization problem and relaxed variants of the proposed methods have been presented.

The sketch-and-project algorithm [25, 73] for solving a consistent linear system $\mathbf{A}x=b$ has the form

[TABLE]

where in each iteration, matrix ${\bf S}_{k}\in\mathbb{R}^{m\times q}$ is sampled afresh from an arbitrary distribution ${\cal D}$.222We stress that there is no restriction on the number of columns of matrix ${\bf S}_{k}$ ($q$ can be varied) In [25] it was shown that many popular randomized algorithms for solving linear systems, including RK method and randomized coordinate descent method (a.k.a Gauss-Seidel method) can be cast as special cases of the above update by choosing an appropriate combination of the distribution ${\cal D}$ and the positive definite matrix $\mathbf{B}$.

In the special case that $\omega=1$ (no relaxation), the update rule of equation (3) can be equivalently written as follows:

[TABLE]

This equivalent presentation of the method justifies the name Sketch and Project. In particular, the method is a two step procedure: (i) Draw random matrix ${\bf S}_{k}\in\mathbb{R}^{m\times q}$ from distribution ${\cal D}$ and formulate a sketched system $\mathbf{S}_{k}^{\top}\mathbf{A}x=\mathbf{S}_{k}^{\top}b$, (ii) Project the last iterate $x^{k}$ into the solution set of the sketched system.

A formal presentation of the Sketch and Project method is shown in Algorithm 1.

In this work, we are mostly interested in two special cases of the sketch and project framework— the randomized Kaczmarz (RK) method and its block variant, the randomized block Kaczmarz (RBK) method. In addition, in the following sections we present novel scaled and accelerated variants of these two selected cases and interpret their gossip nature. In particular, we focus on explaining how these methods can solve the average consensus problem and its more general version, the weighted average consensus (subsection 3.1).

Let $e_{i}\in\mathbb{R}^{m}$ be the $i^{\text{th}}$ unit coordinate vector in $\mathbb{R}^{m}$ and let $\mathbf{I}_{:C}$ be column submatrix of the $m\times m$ identity matrix with columns indexed by $C\subseteq[m]$. Then RK and RBK methods can be obtained as special cases of the update rule (3) as follows:

• •

RK: Let $\mathbf{B}=\mathbf{I}$ and $\mathbf{S}_{k}=e_{i}$, where $i\in[m]$ is chosen independently at each iteration, with probability $p_{i}>0$. In this setup the update rule (3) simplifies to

[TABLE]

• •

RBK: Let $\mathbf{B}=\mathbf{I}$ and $\mathbf{S}=\mathbf{I}_{:C}$, where set $C\subseteq[m]$ is chosen independently at each iteration, with probability $p_{C}\geq 0$. In this setup the update rule (3) simplifies to

[TABLE]

In several papers [26, 73, 48, 50], it was shown that even in the case of consistent linear systems with multiple solutions, Kaczmarz-type methods converge linearly to one particular solution: the projection (on $\mathbf{B}$-norm) of the initial iterate $x^{0}$ onto the solution set of the linear system. This naturally leads to the formulation of the best approximation problem:

[TABLE]

where $\mathbf{A}\in\mathbb{R}^{m\times n}$. In the rest of this manuscript, $x^{*}$ denotes the solution of (7) and we write $x^{*}=\Pi_{{\cal L},\mathbf{B}}(x^{0})$.

Exactness.

An important assumption that is required for the convergence analysis of the randomized iterative methods under study is exactness. That is:

[TABLE]

The exactness property is of key importance in the sketch and project framework, and should be seen as an assumption on the distribution ${\cal D}$ and not on matrix $\mathbf{A}$.

Clearly, an assumption on the distribution ${\cal D}$ of the random matrices $\mathbf{S}$ should be required for the convergence of Algorithm 1. For an instance, if ${\cal D}$ is such that, $\mathbf{S}=e_{1}$ with probability 1, where $e_{1}\in\mathbb{R}^{m}$ be the $1^{\text{st}}$ unit coordinate vector in $\mathbb{R}^{m}$, then the algorithm will select the same row of matrix $\mathbf{A}$ in each step. For this choice of distribution it is clear that the algorithm will not converge to a solution of the linear system. The exactness assumption guarantees that this will not happen.

For necessary and sufficient conditions for exactness, we refer the reader to [73]. Here it suffices to remark that the exactness condition is very weak, allowing ${\cal D}$ to be virtually any reasonable distribution of random matrices. For instance, a sufficient condition for exactness is for the matrix ${\mathbb{E}\left[{\bf H}\right]}$ to be positive definite [26]. This is indeed a weak condition since it is easy to see that this matrix is symmetric and positive semidefinite without the need to invoke any assumptions; simply by design.

A much stronger condition than exactness is ${\mathbb{E}\left[{\bf Z}\right]}\succ 0$ which has been used for the analysis of the sketch and project method in [25]. In this case, the matrix $\mathbf{A}\in\mathbb{R}^{m\times n}$ of the linear system requires to have full column rank and as a result the consistent linear system has a unique solution.

The convergence performance of the Sketch and Project method (Algorithm 1) under the exactness assumption for solving the best approximation problem is described by the following theorem.

Theorem 1** ([73]).**

Let assume exactness and let $\{x^{k}\}_{k=0}^{\infty}$ be the iterates produced by the sketch and project method (Algorithm 1) with step-size $\omega\in(0,2)$. Set, $x^{*}=\Pi_{{\cal L},\mathbf{B}}(x^{0})$. Then,

[TABLE]

where

[TABLE]

Recall that $\lambda_{\min}^{+}$ denotes the minimum nonzero eigenvalue of matrix $\mathbf{W}:=\mathbb{E}[\mathbf{B}^{-1/2}\mathbf{A}^{\top}\mathbf{H}\mathbf{A}\mathbf{B}^{-1/2}]$.

In other words, using standard arguments, from Theorem 1 we observe that for a given $\epsilon\in(0,1)$ we have that:

[TABLE]

We say that the iteration complexity of sketch and project method is,

[TABLE]

2.2 Other related work

On Sketch and Project Methods.

Variants of the sketch-and-project methods have been recently proposed for solving several other problems. [22] and [27] use similar ideas for the development of linearly convergent randomized iterative methods for computing/estimating the inverse and pseudoinverse of a large matrix, respectively. A limited memory variant of the stochastic block BFGS method for solving the empirical risk minimization problem arising in machine learning was proposed by  [23]. Tu et al. [82] utilize the sketch-and-project framework to show that breaking block locality can accelerate block Gauss-Seidel methods. In addition, they develop an accelerated variant of the method for a specific distribution ${\cal D}$. A sketch and project method with the heavy ball momentum was studied in [48, 47] and an accelerated (in the sense of Nesterov) variant of the method proposed in [24] for the more general Euclidean setting and applied to matrix inversion and quasi-Newton updates. Inexact variants of Algorithm 1 have been proposed in [50]. As we have already mentioned, in [73], through the development of stochastic reformulations, a stochastic gradient descent interpretation of the sketch and project method has been proposed. Recently, using a different stochastic reformulation, [21] performed a tight convergence analysis of stochastic gradient descent in a more general convex setting. The analysis proposed in [21] recovers the linear convergence rate of sketch and project method (Theorem 1) as special case.

Gossip algorithms for average consensus

The problem of average consensus has been extensively studied in the automatic control and signal processing literature for the past two decades [14], and was first introduced for decentralized processing in the seminal work [81]. A clear connection between the rate of convergence and spectral characteristics of the underlying network topology over which message passing occurs was first established in [6] for pairwise randomized gossip algorithms.

Motivated by network topologies with salient properties of wireless networks (e.g., nodes can communicate directly only with other nearby nodes), several methods were proposed to accelerate the convergence of gossip algorithms. For instance, [5] proposed averaging among a set of nodes forming a path in the network (this protocol can be seen as special case of our block variants in Section 3.4). Broadcast gossip algorithms have also been analyzed [4] where the nodes communicate with more than one of their neighbors by broadcasting their values.

While the gossip algorithms studied in [6, 5, 4] are all first-order (the update of $x^{k+1}$ only depends on $x^{k}$), a faster randomized pairwise gossip protocol was proposed in [7] which suggested to incorporate additional memory to accelerate convergence. The first analysis of this protocol was later proposed in [40] under strong conditions. It is worth to mention that in the setting of deterministic gossip algorithms theoretical guarantees for accelerated convergence were obtained in [65, 35]. In Section 4 we propose fast and provably accelerated randomized gossip algorithms with memory and compare them in more detail with the fast randomized algorithm proposed in [7, 40].

Gossip algorithms for multiagent consensus optimization.

In the past decade there has been substantial interest in consensus-based mulitiagent optimization methods that use gossip updates in their update rule [55, 90, 76]. In multiagent consensus optimization setting , $n$ agents or nodes, cooperate to solve an optimization problem. In particular, a local objective function $f_{i}:\mathbb{R}^{d}\rightarrow\mathbb{R}$ is associated with each node $i\in[n]$ and the goal is for all nodes to solve the optimization problem

[TABLE]

by communicate only with their neighbors. In this setting gossip algorithms works in two steps by first executing some local computation followed by communication over the network [55]. Note that the average consensus problem with $c_{i}$ as node $i$ initial value can be case as a special case of the optimization problem (11) when the function values are $f_{i}(x)=(x-c_{i})^{2}$.

Recently there has been an increasing interest in applying mulitagent optimization methods to solve convex and non-convex optimization problems arising in machine learning [80, 39, 1, 2, 9, 36, 32]. In this setting most consensus-based optimization methods make use of standard, first-order gossip, such as those described in [6], and incorporating momentum into their updates to improve their practical performance.

3 Sketch and Project Methods as Gossip Algorithms

In this section we show how by carefully choosing the linear system in the constraints of the best approximation problem (7) and the combination of the parameters of the Sketch and Project method (Algorithm 1) we can design efficient randomized gossip algorithms. We show that the proposed protocols can actually solve the weighted average consensus problem, a more general version of the average consensus problem described in Section 1. In particular we focus, on a scaled variant of the RK method (5) and on the RBK (6) and understand the convergence rates of these methods in the consensus setting, their distributed nature and how they are connected with existing gossip protocols.

3.1 Weighted average consensus

In the weighted average consensus (Weighted AC) problem we are given an undirected connected network ${\cal G}=({\cal V},{\cal E})$ with node set ${\cal V}=\{1,2,\dots,n\}$ and edges ${\cal E}$. Each node $i\in{\cal V}$ holds a private value $c_{i}\in\mathbb{R}$ and its weight $w_{i}$. The goal of this problem is for every node to compute the weighted average of the private values,

[TABLE]

in a distributed fashion. That is, the exchange of information can only occur between connected nodes (neighbors).

Note that in the special case when the weights of all nodes are the same ($w_{i}=r$ for all $i\in[n]$) the weighted average consensus is reduced to the standard average consensus problem. However, there are more special cases that could be interesting. For instance the weights can represent the degree of the nodes ($w_{i}=d_{i}$) or they can denote a probability vector and satisfy $\sum_{i}^{n}w_{i}=1$ with $w_{i}>0$.

It can be easily shown that the weighted average consensus problem can be expressed as optimization problem as follows:

[TABLE]

where matrix $\mathbf{B}=\textbf{Diag}(w_{1},w_{2},\dots,w_{n})$ is a diagonal positive definite matrix (that is $w_{i}>0$ for all $i\in[n]$) and $c=(c_{1},\dots,c_{n})^{\top}$ the vector with the initial values $c_{i}$ of all nodes $i\in{\cal V}$. The optimal solution of this problem is $x^{*}_{i}=\frac{\sum_{i=1}^{n}w_{i}c_{i}}{\sum_{i=1}^{n}w_{i}}$ for all $i\in[n]$ which is exactly the solution of the weighted average consensus.

As we have explained, the standard average consensus problem can be cast as a special case of weighted average consensus. However, in the situation when the nodes have access to global information related to the network, such as the size of the network (number of nodes $n=|{\cal V}|$) and the sum of the weights $\sum_{i=1}^{n}w_{i}$, then any algorithm that solves the standard average consensus can be used to solve the weighted average consensus problem with the initial private values of the nodes changed from $c_{i}$ to $\frac{nw_{i}c_{i}}{\sum_{i=1}^{n}w_{i}}$.

The weighted AC problem is popular in the area of distributed cooperative spectrum sensing networks [33, 66, 91, 92]. In this setting, one of the goals is to develop decentralized protocols for solving the cooperative sensing problem in cognitive radio systems. The weights in this case represent a ratio related to the channel conditions of each node/agent [33]. The development of methods for solving the weighted AC problem is an active area of research (check [33] for a recent comparison of existing algorithms). However, to the best of our knowledge, existing analysis for the proposed algorithms focuses on showing convergence and not on providing convergence rates. Our framework allows us to obtain novel randomized gossip algorithms for solving the weighted AC problem. In addition, we provide a tight analysis of their convergence rates. In particular, we show convergence with a linear rate. See Section 7.1 for an experiment confirming linear convergence of one of our proposed protocols on typical wireless network topologies.

3.2 Gossip algorithms through sketch and project framework

We propose that randomized gossip algorithms should be viewed as special case of the Sketch and Project update to a particular problem of the form (7). In particular, we let $c=(c_{1},\dots,c_{n})$ be the initial values stored at the nodes of ${\cal G}$, and choose $\mathbf{A}$ and $b$ so that the constraint $\mathbf{A}x=b$ is equivalent to the requirement that $x_{i}=x_{j}$ (the value stored at node $i$ is equal to the value stored at node $j$) for all $(i,j)\in{\cal E}$.

Definition 2**.**

We say that $\mathbf{A}x=b$ is an “average consensus (AC) system” when $\mathbf{A}x=b$ iff $x_{i}=x_{j}$ for all $(i,j)\in{\cal E}$.

It is easy to see that $\mathbf{A}x=b$ is an AC system precisely when $b=0$ and the nullspace of $\mathbf{A}$ is $\{t1_{n}:t\in\mathbb{R}\}$, where $1_{n}$ is the vector of all ones in $\mathbb{R}^{n}$. Hence, $\mathbf{A}$ has rank $n-1$. Moreover in the case that $x^{0}=c$, it is easy to see that for any AC system, the solution of (7) necessarily is $x^{*}=\bar{c}\cdot 1_{n}$ — this is why we singled out AC systems. In this sense, any algorithm for solving (7) will “find” the (weighted) average $\bar{c}$. However, in order to obtain a distributed algorithm we need to make sure that only “local” (with respect to ${\cal G}$) exchange of information is allowed.

It can be shown that many linear systems satisfy the above definition.

For example, we can choose $b=0$ and $\mathbf{A}=\mathbf{Q}\in\mathbb{R}^{|{\cal E}|\times n}$ to be the incidence matrix of ${\cal G}$. That is, $\mathbf{Q}\in\mathbb{R}^{|{\cal E}|\times n}$ such that $\mathbf{Q}x=0$ directly encodes the constraints $x_{i}=x_{j}$ for $(i,j)\in{\cal E}$. That is, row $e=(i,j)\in{\cal E}$ of matrix $\mathbf{Q}$ contains value $1$ in column $i$, value $-1$ in column $j$ (we use an arbitrary but fixed order of nodes defining each edge in order to fix $\mathbf{Q}$) and zeros elsewhere. A different choice is to pick $b=0$ and ${\bf A}={\bf L}=\mathbf{Q}^{\top}\mathbf{Q}$, where ${\bf L}$ is the Laplacian matrix of network ${\cal G}$. Depending on what AC system is used, the sketch and project methods can have different interpretations as gossip protocols.

In this work we mainly focus on the above two AC systems but we highlight that other choices are possible333Novel gossip algorithms can be proposed by using different AC systems to formulate the average consensus problem. For example one possibility is using the random walk normalized Laplacian $\mathbf{L}^{rw}=\mathbf{D}^{-1}\mathbf{L}$. For the case of degree-regular networks the symmetric normalized Laplacian matrix $\mathbf{L}^{sym}=\mathbf{D}^{-1/2}\mathbf{L}\mathbf{D}^{-1/2}$ can also being used.. In Section 4.2 for the provably accelerated gossip protocols we also use a normalized variant ($\|\mathbf{A}_{i:}\|^{2}=1$) of the Incidence matrix.

3.2.1 Standard form and mass preservation

Assume that $\mathbf{A}x=b$ is an AC system. Note that since $b=0$, the general sketch-and-project update rule (3) simplifies to:

[TABLE]

This is the standard form in which randomized gossip algorithms are written. What is new here is that the iteration matrix $\mathbf{I}-\omega\mathbf{Z}_{k}$ has a specific structure which guarantees convergence to $x^{*}$ under very weak assumptions (see Theorem 1). Note that if $x^{0}=c$, i.e., the starting primal iterate is the vector of private values (as should be expected from any gossip algorithm), then the iterates of (13) enjoy a mass preservation property (the proof follows the fact that $\mathbf{A}1_{n}=0$):

Theorem 3** (Mass preservation).**

If $\mathbf{A}x=b$ is an AC system, then the iterates produced by (13) satisfy: $\frac{1}{n}\sum_{i=1}^{n}x_{i}^{k}=\bar{c}$, for all $k\geq 0$.

Proof.

Let fix $k\geq 0$ then,

[TABLE]

∎

3.2.2 $\varepsilon$-Averaging time

Let $z^{k}:=\|x^{k}-x^{*}\|$. The typical measure of convergence speed employed in the randomized gossip literature, called $\varepsilon$-averaging time and here denoted by $T_{ave}(\varepsilon)$, represents the smallest time $k$ for which $x^{k}$ gets within $\varepsilon z^{0}$ from $x^{*}$, with probability greater than $1-\varepsilon$, uniformly over all starting values $x^{0}=c$. More formally, we define

[TABLE]

This definition differs slightly from the standard one in that we use $z^{0}$ instead of $\|c\|$.

Inequality (9), together with Markov inequality, can be used to give a bound on $K(\varepsilon)$, formalized next:

Theorem 4**.**

Assume $\mathbf{A}x=b$ is an AC system. Let $x^{0}=c$ and $\mathbf{B}$ be positive definite diagonal matrix. Assume exactness. Then for any $0<\varepsilon<1$ we have

[TABLE]

where $\rho$ is defined in (10).

Proof.

See Appendix A.1. ∎

Note that under the assumptions of the above theorem, $\mathbf{W}=\mathbf{B}^{-1/2}\mathbb{E}[\mathbf{Z}]\mathbf{B}^{-1/2}$ only has a single zero eigenvalue, and hence $\lambda_{\min}^{+}(\mathbf{W})$ is the second smallest eigenvalue of $\mathbf{W}$. Thus, $\rho$ is the second largest eigenvalue of $\mathbf{I}-\mathbf{W}$. The bound on $K(\varepsilon)$ appearing in Thm 4 is often written with $\rho$ replaced by $\lambda_{2}(\mathbf{I}-\mathbf{W})$ [6].

In the rest of this section we show how two special cases of the sketch and project framework, the randomized Kaczmarz (RK) and its block variant, randomized block Kaczmatz (RBK) work as gossip algorithms for the two AC systems described above.

3.3 Randomized Kaczmarz method as gossip algorithm

As we described before the sketch and project update rule (3) has several parameters that should be chosen in advance by the user. These are the stepsize $\omega$ (relaxation parameter), the positive definite matrix $\mathbf{B}$ and the distribution ${\cal D}$ of the random matrices $\mathbf{S}$.

In this section we focus on one particular special case of the sketch and project framework, a scaled/weighted variant of the randomized Kaczmarz method (RK) presented in (5), and we show how this method works as gossip algorithm when applied to special systems encoding the underlying network. In particular, the linear systems that we solve are the two AC systems described in the previous section where the matrix is either the incidence matrix $\mathbf{Q}$ or the Laplacian matrix $\mathbf{L}$ of the network.

As we described in (5) the standard RK method can be cast as special case of the sketch and project update (3) by choosing $\mathbf{B}=\mathbf{I}$ and $\mathbf{S}=e_{i}$. In this section, we focus on a small modification of this algorithm and we choose the positive definite matrix $\mathbf{B}$ to be $\mathbf{B}=\textbf{Diag}(w_{1},w_{2},\dots,w_{n})$, the diagonal matrix of the weights presented in the weighted average consensus problem.

Scaled RK:

Let us have a general consistent linear system $\mathbf{A}x=b$ with $\mathbf{A}\in\mathbb{R}^{m\times n}$. Let us also choose $\mathbf{B}=\textbf{Diag}(w_{1},w_{2},\dots,w_{n})$ and $\mathbf{S}_{k}=e_{i}$, where $i\in[m]$ is chosen in each iteration independently, with probability $p_{i}>0$. In this setup the update rule (3) simplifies to

[TABLE]

This small modification of RK allow us to solve the more general weighted average consensus presented in Section 3.1 (and at the same time the standard average consensus problem if $\mathbf{B}=r\mathbf{I}$ where $r\in\mathbb{R}$). To the best of our knowledge, even if this variant is special case of the general Sketch and project update, was never precisely presented before in any setting.

3.3.1 AC system with incidence matrix $\mathbf{Q}$

Let us represent the constraints of problem (12) as linear system with matrix $\mathbf{A}=\mathbf{Q}\in\mathbb{R}^{|{\cal E}|\times n}$ be the Incidence matrix of the graph and right had side $b=0$. Lets also assume that the random matrices $\mathbf{S}\sim{\cal D}$ are unit coordinate vectors in $\mathbb{R}^{m}=\mathbb{R}^{|{\cal E}|}$.

Let $e=(i,j)\in{\cal E}$ then from the definition of matrix $\mathbf{Q}$ we have that $\mathbf{Q}_{e:}^{\top}=f_{i}-f_{j}$ where $f_{i},f_{j}$ are unit coordinate vectors in $\mathbb{R}^{n}$. In addition, from the definition the diagonal positive definite matrix $\mathbf{B}$ we have that

[TABLE]

Thus in this case the update rule (14) simplifies:

[TABLE]

From (16) it can be easily seen that only the values of coordinates $i$ and $j$ update their values. These coordinates correspond to the private values $x^{k}_{i}$ and $x^{k}_{j}$ of the nodes of the selected edge $e=(i,j)$. In particular the values of $x^{k}_{i}$ and $x^{k}_{j}$ are updated as follows:

[TABLE]

Remark 1**.**

In the special case that $\mathbf{B}=r\mathbf{I}$ where $r\in\mathbb{R}$ (we solve the standard average consensus problem) the update of the two nodes is simplified to

[TABLE]

If we further select $\omega=1$ then this becomes:

[TABLE]

which is the update of the standard pairwise randomized gossip algorithm first presented and analyzed in [6].

3.3.2 AC system with Laplacian matrix $\mathbf{L}$

The AC system takes the form $\mathbf{L}x=0$, where matrix $\mathbf{L}\in\mathbb{R}^{n\times n}$ is the Laplacian matrix of the network. In this case, each row of the matrix corresponds to a node. Using the definition of the Laplacian, we have that $\mathbf{L}_{i:}^{\top}=d_{i}f_{i}-\sum_{j\in{\cal N}_{i}}f_{j}$, where $f_{i},f_{j}$ are unit coordinate vectors in $\mathbb{R}^{n}$ and $d_{i}$ is the degree of node $i\in{\cal V}$.

Thus, by letting $\mathbf{B}=\textbf{Diag}(w_{1},w_{2},\dots,w_{n})$ to be the diagonal matrix of the weights we obtain:

[TABLE]

In this case, the update rule (14) simplifies to:

[TABLE]

From (20), it is clear that only coordinates $\{i\}\cup{\cal N}_{i}$ update their values. All the other coordinates remain unchanged. In particular, the value of the selected node $i$ (coordinate $i$) is updated as follows:

[TABLE]

while the values of its neighbors $j\in{\cal N}_{i}$ are updated as:

[TABLE]

Remark 2**.**

Let $\omega=1$ and $\mathbf{B}=r\mathbf{I}$ where $r\in\mathbb{R}$ then the selected nodes update their values as follows:

[TABLE]

That is, the selected node $i$ updates its value to the average of its neighbors and itself, while all the nodes $j\in{\cal N}_{i}$ update their values using the current value of node $i$ and all nodes in ${\cal N}_{i}$.

In a wireless network, to implement such an update, node $i$ would first broadcast its current value to all of its neighbors. Then it would need to receive values from each neighbor to compute the sums over ${\cal N}_{i}$, after which node $i$ would broadcast the sum to all neighbors (since there may be two neighbors $j_{1},j_{2}\in{\cal N}_{i}$ for which $(j_{1},j_{2})\notin{\cal E}$). In a wired network, using standard concepts from the MPI library, such an update rule could be implemented efficiently by defining a process group consisting of $\{i\}\cup{\cal N}_{i}$, and performing one Broadcast in this group from $i$ (containing $x_{i}$) followed by an AllReduce to sum $x_{\ell}$ over $\ell\in{\cal N}_{i}$. Note that the terms involving diagonal entries of $\mathbf{B}$ and the degrees $d_{i}$ could be sent once, cached, and reused throughout the algorithm execution to reduce communication overhead.

3.3.3 Details on complexity results

Recall that the convergence rate of the sketch and project method (Algorithm 1) is equivalent to:

[TABLE]

where $\omega\in(0,2)$ and $\mathbf{W}=\mathbf{B}^{-1/2}\mathbf{A}^{\top}\mathbb{E}[\mathbf{H}]\mathbf{A}\mathbf{B}^{-1/2}$ (from Theorem 1). In this subsection we explain how the convergence rate of the scaled RK method (14) is modified for different choices of the main parameters of the method.

Let us choose $\omega=1$ (no over-relaxation). In this case, the rate is simplified to $\rho=1-\lambda_{\min}^{+}$.

Note that the different ways of modeling the problem (AC system) and the selection of the main parameters (weight matrix $\mathbf{B}$ and distribution ${\cal D}$) determine the convergence rate of the method through the spectrum of matrix $\mathbf{W}$.

Recall that in the $k^{th}$ iterate of the scaled RK method (14) a random vector $\mathbf{S}_{k}=e_{i}$ is chosen with probability $p_{i}>0$. For convenience, let us choose444Similar probabilities have been chosen in [25] for the convergence of the standard RK method ($\mathbf{B}=\mathbf{I}$). The distribution ${\cal D}$ of the matrices $\mathbf{S}$ used in equation (24) is common in the area of randomized iterative methods for linear systems and is used to simplify the analysis and the expressions of the convergence rates. For more choices of distributions we refer the interested reader to [25]. It is worth to mention that the probability distribution that optimizes the convergence rate of the RK and other projection methods can be expressed as the solution to a convex semidefinite program [25, 11].:

[TABLE]

Then we have that:

[TABLE]

and

[TABLE]

Incidence Matrix:

Let us choose the AC system to be the one with the incidence matrix $\mathbf{A}=\mathbf{Q}$. Then $\|\mathbf{B}^{-1/2}\mathbf{Q}^{\top}\|^{2}_{F}=\sum_{i=1}^{n}\frac{d_{i}}{b_{i}}$ and we obtain

[TABLE]

If we further have $\mathbf{B}=\mathbf{D}$, then $\mathbf{W}=\frac{\mathbf{D}^{-1/2}\mathbf{L}\mathbf{D}^{-1/2}}{n}$ and the convergence rate simplifies to:

[TABLE]

If $\mathbf{B}=r\mathbf{I}$ where $r\in\mathbb{R}$ (solve the standard average consensus problem), then $\mathbf{W}=\frac{\mathbf{L}}{\|\mathbf{Q}\|^{2}_{F}}=\frac{\mathbf{L}}{\sum_{i=1}^{n}d_{i}}=\frac{\mathbf{L}}{2m}$ and the convergence rate simplifies to

[TABLE]

The convergence rate (27) is identical to the rate proposed for the convergence of the standard pairwise gossip algorithm in [6]. Recall that in this special case the proposed gossip protocol has exactly the same update rule with the algorithm presented in [6], see equation (18).

Laplacian Matrix:

If we choose to formulate the AC system using the Laplacian matrix $\mathbf{L}$, that is $\mathbf{A}=\mathbf{L}$, then $\|\mathbf{B}^{-1/2}\mathbf{L}^{\top}\|^{2}_{F}=\sum_{i=1}^{n}\frac{d_{i}(d_{i}+1)}{\mathbf{B}_{ii}}$ and we have:

[TABLE]

If $\mathbf{B}=\mathbf{D}$, then the convergence rate simplifies to:

[TABLE]

If $\mathbf{B}=r\mathbf{I}$, where $r\in\mathbb{R}$, then $\mathbf{W}=\frac{\mathbf{L}^{2}}{\|\mathbf{L}\|^{2}_{F}}=\frac{\mathbf{L}^{2}}{\sum_{i=1}^{n}d_{i}(d_{i}+1)}$ and the convergence rate simplifies to

[TABLE]

3.4 Block gossip algorithms

Up to this point we focused on the basic connections between the convergence analysis of the sketch and project methods and the literature of randomized gossip algorithms. We show how specific variants of the randomized Kaczmarz method (RK) can be interpreted as gossip algorithms for solving the weighted and standard average consensus problems.

In this part we extend the previously described methods to their block variants related to randomized block Kaczmarz (RBK) method (6). In particular, in each step of the sketch and project method (3), the random matrix $\mathbf{S}$ is selected to be a random column submatrix of the $m\times m$ identity matrix corresponding to columns indexed by a random subset $C\subseteq[m]$. That is, $\mathbf{S}=\mathbf{I}_{:C}$, where a set $C\subseteq[m]$ is chosen in each iteration independently, with probability $p_{C}\geq 0$ (see equation (6)). Note that in the special case that set $C$ is a singleton with probability 1 the algorithm is simply the randomized Kaczmarz method of the previous section.

To keep things simple, we assume that $\mathbf{B}=\mathbf{I}$ (standard average consensus, without weights) and choose the stepsize $\omega=1$. In the next section, we will describe gossip algorithms with heavy ball momentum and explain in detail how the gossip interpretation of RBK change in the more general case of $\omega\in(0,2)$.

Similar to the previous subsections, we formulate the consensus problem using either $\mathbf{A}=\mathbf{Q}$ or $\mathbf{A}=\mathbf{L}$ as the matrix in the AC system. In this setup, the iterative process (3) has the form:

[TABLE]

which, as explained in the introduction, can be equivalently written as:

[TABLE]

Essentially in each step of this method the next iterate is evaluated to be the projection of the current iterate $x^{k}$ onto the solution set of a row subsystem of $\mathbf{A}x=0$.

AC system with Incidence Matrix:

In the case that $\mathbf{A}=\mathbf{Q}$ the selected rows correspond to a random subset $C\subseteq{\cal E}$ of selected edges. While (28) may seem to be a complicated algebraic (resp. variational) characterization of the method, due to our choice of $\mathbf{A}=\mathbf{Q}$ we have the following result which gives a natural interpretation of RBK as a gossip algorithm (see also Figure 1).

Theorem 5** (RBK as Gossip algorithm: RBKG).**

Consider the AC system with the constraints being expressed using the Incidence matrix $\mathbf{Q}$. Then each iteration of RBK (Algorithm (28)) works as gossip algorithm as follows:

Select a random set of edges $C\subseteq{\cal E}$, 2. 2.

Form subgraph ${\cal G}_{k}$ of ${\cal G}$ from the selected edges 3. 3.

For each connected component of ${\cal G}_{k}$, replace node values with their average.

Proof.

See Appendix A.2. ∎

Using the convergence result of general Theorem 1 and the form of matrix $\mathbf{W}$ (recall that in this case we assume $\mathbf{B}=\mathbf{I}$, $\mathbf{S}=\mathbf{I}_{:C}\sim{\cal D}$ and $\omega=1$), we obtain the following complexity for the algorithm:

[TABLE]

For more details on the above convergence rate of randomized block Kaczmarz method with meaningfully bounds on the rate in a more general setting we suggest the papers [57, 58].

There is a very closed relationship between the gossip interpretation of RBK explained in Theorem 5 and several existing randomized gossip algorithms that in each step update the values of more than two nodes. For example the path averaging algorithm porposed in [5] is a special case of RBK, when set $C$ is restricted to correspond to a path of vertices. That is, in path averaging, in each iteration a path of nodes is selected and the nodes that belong to it update their values to their exact average. A different example is the recently proposed clique gossiping [44] where the network is already divided into cliques and through a random procedure a clique is activated and the nodes of it update their values to their exact average. In [6] a synchronous variant of gossip algorithm is presented where in each step multiple node pairs communicate exactly at the same time with the restriction that these simultaneously active node pairs are disjoint.

It is easy to see that all of the above algorithms can be cast as special cases of RBK if the distribution ${\cal D}$ of the random matrices is chosen carefully to be over random matrices $\mathbf{S}$ (column sub-matrices of Identity) that update specific set of edges in each iteration. As a result our general convergence analysis can recover the complexity results proposed in the above works.

Finally, as we mentioned, in the special case in which set $C$ is always a singleton, Algorithm (28) reduces to the standard randomized Kaczmarz method. This means that only a random edge is selected in each iteration and the nodes incident with this edge replace their local values with their average. This is the pairwise gossip algorithm of Boyd er al. [6] presented in equation (18). Theorem 5 extends this interpretation to the case of the RBK method.

AC system with Laplacian Matrix:

For this choice of AC system the update is more complicated. To simplify the way that the block variant work as gossip we make an extra assumption. We assume that the selected rows of the constraint $\mathbf{I}_{:C}^{\top}\mathbf{L}x=0$ in update (29) have no-zero elements at different coordinates. This allows to have a direct extension of the serial variant presented in Remark 2. Thus, in this setup, the RBK update rule (28) works as gossip algorithm as follows:

$|C|$ nodes are activated (with restriction that the nodes are not neighbors and they do not share common neighbors) 2. 2.

For each node $i\in C$ we have the following update:

[TABLE]

The above update rule can be seen as a parallel variant of update (23). Similar to the convergence in the case of Incidence matrix, the RBK for solving the AC system with a Laplacian matrix converges to $x^{*}$ with the following rate (using result of Theorem 1):

[TABLE]

4 Faster and Provably Accelerated Randomized Gossip Algorithms

The main goal in the design of gossip protocols is for the computation and communication to be done as quickly and efficiently as possible. In this section, our focus is precisely this. We design randomized gossip protocols which converge to consensus fast with provable accelerated linear rates. To the best of our knowledge, the proposed protocols are the first randomized gossip algorithms that converge to consensus with an accelerated linear rate.

In particular, we present novel protocols for solving the average consensus problem where in each step all nodes of the network update their values but only a subset of them exchange their private values. The protocols are inspired from the recently developed accelerated variants of randomized Kaczmarz-type methods for solving consistent linear systems where the addition of momentum terms on top of the sketch and project update rule provides better theoretical and practical performance.

In the area of optimization algorithms, there are two popular ways to accelerate an algorithm using momentum. The first one is using the Polyak’s heavy ball momentum [67] and the second one is using the theoretically much better understood momentum introduced by Nesterov [59, 61]. Both momentum approaches have been recently proposed and analyzed to improve the performance of randomized iterative methods for solving linear systems.

To simplify the presentation, the accelerated algorithms and their convergence rates are presented for solving the standard average consensus problem ($\mathbf{B}=\mathbf{I}$). Using a similar approach as in the previous section, the update rules and the convergence rates can be easily modified to solve the more general weighted average consensus problem. For the protocols in this section we use the incidence matrix $\mathbf{A}=\mathbf{Q}$ or its normalized variant to formulate the AC system.

4.1 Gossip algorithms with heavy ball momentum

The recent works [47, 48] propose and analyze heavy ball momentum variants of several stochastic optimization algorithms for solving stochastic quadratic optimization problems and linear systems. One of the proposed algorithms is the sketch and project method (Algorithm 1) with heavy ball momentum. In particular, the authors focus on explaining how the method can be interpreted as SGD with momentum—also known as the stochastic heavy ball method (SHB). SHB is a well known algorithm in the optimization literature for solving stochastic optimization problems, and extremely popular in areas such as deep learning [78, 79, 37, 85]. However, even if SHB is used extensively in practice, its theoretical convergence behavior is not well understood. [47, 48] were the first works that prove linear convergence of SHB in any setting.

In this subsection we focus on the sketch and project method with heavy ball momentum. We present the main theorems showing its convergence performance as presented in [47, 48] and explain how special cases of the general method work as gossip algorithms when are applied to a special system encoding the underlying network.

4.1.1 Sketch and project with heavy ball momentum

The update rule of the sketch and project method with heavy ball momentum as proposed and analyzed in [47, 48] is formally presented in the following algorithm:

Using, $\mathbf{B}=\mathbf{I}$ and the same choice of distribution ${\cal D}$ as in equations (5) and (6) we can now obtain momentum variants of the RK and RBK as special case of the above algorithm as follows:

• •

RK with momentum (mRK):

[TABLE]

• •

RBK with momentum (mRBK):

[TABLE]

In [48], two main theoretical results describing the behavior of Algorithm 2 (and as a result also the special cases mRK and mRBK) have been presented555Note that in [48] the analysis have been made on the same framework of Theorem 1 with general positive definite matrix $\mathbf{B}$..

Theorem 6** (Theorem 1, [48]).**

Choose $x^{0}=x^{1}\in\mathbb{R}^{n}$. Let $\{x^{k}\}_{k=0}^{\infty}$ be the sequence of random iterates produced by Algorithm 2 and let assume exactness. Assume $0<\omega<2$ and $\beta\geq 0$ and that the expressions $a_{1}:=1+3\beta+2\beta^{2}-(\omega(2-\omega)+\omega\beta)\lambda_{\min}^{+}$ and $a_{2}:=\beta+2\beta^{2}+\omega\beta\lambda_{\max}$ satisfy $a_{1}+a_{2}<1$. Set $x^{*}=\Pi_{{\cal L},\mathbf{B}}(x^{0})$. Then

[TABLE]

where $q=\frac{1}{2}(a_{1}+\sqrt{a_{1}^{2}+4a_{2}})$ and $\delta=q-a_{1}$. Moreover, $a_{1}+a_{2}\leq q<1$.

Theorem 7** (Theorem 4, [48]).**

Assume exactness. Let $\{x^{k}\}_{k=0}^{\infty}$ be the sequence of random iterates produced by Algorithm 2, started with $x^{0}=x^{1}\in\mathbb{R}^{n}$, with relaxation parameter (stepsize) $0<\omega\leq 1/\lambda_{\max}$ and momentum parameter $(1-\sqrt{\omega\lambda_{\min}^{+}})^{2}<\beta<1$. Let $x^{*}=\Pi_{{\cal L},\mathbf{I}}(x^{0})$. Then there exists a constant $C>0$ such that for all $k\geq 0$ we have

[TABLE]

Using Theorem 7 and by a proper combination of the stepsize $\omega$ and the momentum parameter $\beta$, Algorithm 2 enjoys an accelerated linear convergence rate in mean [48].

Corollary 1**.**

1. (i)

If $\omega=1$ and $\beta=(1-\sqrt{0.99\lambda_{\min}^{+}})^{2}$, then the iteration complexity of Algorithm 2 becomes: $O\left(\sqrt{1/\lambda_{\min}^{+}}\log(1/\epsilon)\right)$. 2. (ii)

If $\omega=1/\lambda_{\max}$ and $\beta=(1-\sqrt{0.99\lambda_{\min}^{+}/\lambda_{\max}})^{2}$, then the iteration complexity of Algorithm 2 becomes: $O\left(\sqrt{\lambda_{\max}/\lambda_{\min}^{+}}\log(1/\epsilon)\right)$.

Having presented Algorithm 2 and its convergence analysis results, let us now describe its behavior as a randomized gossip protocol when applied on the AC system $\mathbf{A}x=0$ with $\mathbf{A}=\mathbf{Q}\in|{\cal E}|\times n$ (incidence matrix of the network).

Note that since $b=0$ (from the AC system definition), method (32) can be simplified to:

[TABLE]

In the rest of this section we focus on two special cases of (36): RK with heavy ball momentum (equation (33) with $b_{i}=0$) and RBK with heavy ball momentum (equation (34) with $b_{C}=0$).

4.1.2 Randomized Kaczmarz gossip with heavy ball momentum

As we have seen in previous section when the standard RK is applied to solve the AC system $\mathbf{Q}x=0$, one can recover the famous pairwise gossip algorithm [6]. Algorithm 3 describes how a relaxed variant of randomized Kaczmarz with heavy ball momentum ($0<\omega<2$ and $0\leq\beta<1$) behaves as a gossip algorithm. See also Figure (2) for a graphical illustration of the method.

Remark 3**.**

In the special case that $\beta=0$ (zero momentum) only the two nodes of edge $e=(i,j)$ update their values. In this case the two selected nodes do not update their values to their exact average but to a convex combination that depends on the stepsize $\omega\in(0,2)$. To obtain the pairwise gossip algorithm of [6], one should further choose $\omega=1$.

Distributed Nature of the Algorithm: Here we highlight a few ways to implement mRK in a distributed fashion.

• •

Pairwise broadcast gossip: In this protocol each node $i\in{\cal V}$ of the network ${\cal G}$ has a clock that ticks at the times of a rate 1 Poisson process. The inter-tick times are exponentially distributed, independent across nodes, and independent across time. This is equivalent to a global clock ticking at a rate $n$ Poisson process which wakes up an edge of the network at random. In particular, in this implementation mRK works as follows: In the $k^{th}$ iteration (time slot) the clock of node $i$ ticks and node $i$ randomly contact one of its neighbors and simultaneously broadcast a signal to inform the nodes of the whole network that is updating (this signal does not contain any private information of node $i$). The two nodes $(i,j)$ share their information and update their private values following the update rule of Algorithm 3 while all the other nodes update their values using their own information. In each iteration only one pair of nodes exchange their private values.

• •

Synchronous pairwise gossip: In this protocol a single global clock is available to all nodes. The time is assumed to be slotted commonly across nodes and in each time slot only a pair of nodes of the network is randomly activated and exchange their information following the update rule of Algorithm 3. The remaining not activated nodes update their values using their own last two private values. Note that this implementation of mRK comes with the disadvantage that it requires a central entity which in each step requires to choose the activated pair of nodes666We speculate that a completely distributed synchronous gossip algorithm that finds pair of nodes in a distributed manner without any additional computational burden can be design following the same procedure proposed in Section III.C of [6]..

• •

Asynchronous pairwise gossip with common counter: Note that the update rule of the selected pair of nodes $(i,j)$ in Algorithm 3 can be rewritten as follows:

[TABLE]

[TABLE]

In particular observe that the first part of the above expressions $x_{i}^{k}+\beta(x_{i}^{k}-x_{i}^{k-1})$ (for the case of node $i$) is exactly the same with the update rule of the non activate nodes at $k^{th}$ iterate (check step 5 of Algorithm 3) . Thus, if we assume that all nodes share a common counter that keeps track of the current iteration count and that each node $i\in{\cal V}$ remembers the iteration counter $k_{i}$ of when it was last activated, then step 5 of Algorithm 3 takes the form:

• –

$x_{i}^{k+1}=i_{k}\left[x_{i}^{k}+\beta(x_{i}^{k}-x_{i}^{k-1})\right]+\frac{\omega}{2}(x_{j}^{k}-x_{i}^{k}),$

• –

$x_{j}^{k+1}=j_{k}\left[x_{j}^{k}+\beta(x_{j}^{k}-x_{j}^{k-1})\right]+\frac{\omega}{2}(x_{i}^{k}-x_{j}^{k}),$

• –

$k_{i}=k_{j}=k+1,$

• –

Any other node $\ell$: $x_{\ell}^{k+1}=x_{\ell}^{k},$

where $i_{k}=k-k_{i}$ ($j_{k}=k-k_{j}$) denotes the number of iterations between the current iterate and the last time that the $i^{th}$ ($j^{th}$) node is activated. In this implementation only a pair of nodes communicate and update their values in each iteration (thus the justification of asynchronous), however it requires the nodes to share a common counter that keeps track the current iteration count in order to be able to compute the value of $i_{k}=k-k_{i}$.

4.1.3 Connections with existing fast randomized gossip algorithms

In the randomized gossip literature there is one particular method closely related to our approach. It was first proposed in [7] and its analysis under strong conditions was presented in [40]. In this work local memory is exploited by installing shift registers at each agent. In particular we are interested in the case of two registers where the first stores the agent’s current value and the second the agent’s value before the latest update. The algorithm can be described as follows. Suppose that edge $e=(i,j)$ is chosen at time $k$. Then,

• •

Node $i$: $x_{i}^{k+1}=\omega(\frac{x_{i}^{k}+x_{j}^{k}}{2})+(1-\omega)x_{i}^{k-1},$

• •

Node $j$: $x_{i}^{k+1}=\omega(\frac{x_{i}^{k}+x_{j}^{k}}{2})+(1-\omega)x_{j}^{k-1},$

• •

Any other node $\ell$: $x_{\ell}^{k+1}=x_{\ell}^{k},$

where $\omega\in[1,2)$. The method was analyzed in [40] under a strong assumption on the probabilities of choosing the pair of nodes, that as the authors mentioned, is unrealistic in practical scenarios, and for networks like the random geometric graphs. At this point we should highlight that the results presented in [48] hold for essentially any distribution ${\cal D}$ 777The only restriction is the exactness condition to be satisfied. See Theorem 6. and as a result in the proposed gossip variants with heavy ball momentum such problem cannot occur.

Note that, in the special case that we choose $\beta=\omega-1$ in the update rule of Algorithm 3 is simplified to:

• •

Node $i$: $x_{i}^{k+1}=\omega(\frac{x_{i}^{k}+x_{j}^{k}}{2})+(1-\omega)x_{i}^{k-1},$

• •

Node $j$: $x_{i}^{k+1}=\omega(\frac{x_{i}^{k}+x_{j}^{k}}{2})+(1-\omega)x_{j}^{k-1},$

• •

Any other node $\ell$: $x_{\ell}^{k+1}=\omega x_{\ell}^{k}+(1-\omega)x_{\ell}^{k-1}.$

In order to apply Theorem 6, we need to assume that $0<\omega<2$ and $\beta=\omega-1\geq 0$ which also means that $\omega\in[1,2)$. Thus for $\omega\in[1,2)$ and momentum parameter $\beta=\omega-1$ it is easy to see that our approach is very similar to the shift-register algorithm. Both methods update the selected pair of nodes in the same way. However, in Algorithm 3 the not selected nodes of the network do not remain idle but instead update their values using their own previous information.

By defining the momentum matrix $\mathbf{M}=\text{{Diag}}(\beta_{1},\beta_{2},\dots,\beta_{n})$, the above closely related algorithms can be expressed, in vector form, as:

[TABLE]

In particular, in mRK every diagonal element of matrix $\mathbf{M}$ is equal to $\omega-1$, while in the algorithm of [7, 40] all the diagonal elements are zeros except the two values that correspond to nodes $i$ and $j$ that are equal to $\beta_{i}=\beta_{j}=\omega-1$.

Remark 4**.**

The shift register algorithm of [40] and Algorithm 3 of this work can be seen as the two limit cases of the update rule (37). As we mentioned, the shift register method [40] uses only two non-zero diagonal elements in $\mathbf{M}$, while our method has a full diagonal. We believe that further methods can be developed in the future by exploring the cases where more than two but not all elements of the diagonal matrix $\mathbf{M}$ are non-zero. It might be possible to obtain better convergence if one carefully chooses these values based on the network topology. We leave this as an open problem for future research.

4.1.4 Randomized block Kaczmarz gossip with heavy ball momentum

Recall that Theorem 5 explains how RBK (with no momentum and no relaxation) can be interpreted as a gossip algorithm. In this subsection by using this result we explain how relaxed RBK with momentum works. Note that the update rule of RBK with momentum can be rewritten as follows:

[TABLE]

and recall that $x^{k+1}=\left(\mathbf{I}-\mathbf{A}_{C:}^{\top}(\mathbf{A}_{C:}\mathbf{A}_{C:}^{\top})^{\dagger}\mathbf{A}_{C:}\right)x^{k}$ is the update rule of the standard RBK (28).

Thus, in analogy to the standard RBK, in the $k^{th}$ step, a random set of edges is selected and $q\leq n$ connected components are formed as a result. This includes the connected components that belong to both sub-graph ${\cal G}_{k}$ and also the singleton connected components (nodes outside the ${\cal G}_{k}$). Let us define the set of the nodes that belong in the $r\in[q]$ connected component at the $k^{th}$ step ${\cal V}_{r}^{k}$, such that ${\cal V}=\cup_{r\in[q]}{\cal V}_{r}^{k}$ and $|{\cal V}|=\sum_{r=1}^{q}|{\cal V}_{r}^{k}|$ for any $k>0$.

Using the update rule (38), Algorithm 4 shows how mRBK is updating the private values of the nodes of the network (see also Figure 3 for the graphical interpretation).

Note that in the update rule of mRBK the nodes that are not attached to a selected edge (do not belong in the sub-graph ${\cal G}_{k}$) update their values via $x_{\ell}^{k+1}=x_{\ell}^{k}+\beta(x_{\ell}^{k}-x_{\ell}^{k-1})$. By considering these nodes as singleton connected components their update rule is exactly the same with the nodes of sub-graph ${\cal G}_{k}$. This is easy to see as follows:

[TABLE]

Remark 5**.**

In the special case that only one edge is selected in each iteration ($\mathbf{S}_{k}\in\mathbb{R}^{m\times 1}$) the update rule of mRBK is simplified to the update rule of mRK. In this case the sub-graph ${\cal G}_{k}$ is the pair of the two selected edges.

Remark 6**.**

In previous section we explained how several existing gossip protocols for solving the average consensus problem are special cases of the RBK (Theorem 5). For example two gossip algorithms that can be cast as special cases of the standard RBK are the path averaging proposed in [5] and the clique gossiping [44]. In path averaging, in each iteration a path of nodes is selected and its nodes update their values to their exact average ($\omega=1$). In clique gossiping, the network is already divided into cliques and through a random procedure a clique is activated and the nodes of it update their values to their exact average ($\omega=1$). Since mRBK contains the standard RBK as a special case (when $\beta=0$), we expect that these special protocols can also be accelerated with the addition of momentum parameter $\beta\in(0,1)$.

4.1.5 Mass preservation

One of the key properties of some of the most efficient randomized gossip algorithms is mass preservation. That is, the sum (and as a result the average) of the private values of the nodes remains fixed during the iterative procedure ($\textstyle\sum_{i=1}^{n}x_{i}^{k}=\sum_{i=1}^{n}x_{i}^{0},\quad\forall k\geq 1$). The original pairwise gossip algorithm proposed in [6] satisfied the mass preservation property, while exisiting fast gossip algorithms [7, 40] preserving a scaled sum. In this subsection we show that mRK and mRBK gossip protocols presented above satisfy the mass preservation property. In particular, we prove mass preservation for the case of the block randomized gossip protocol (Algorithm 4) with momentum. This is sufficient since the randomized Kaczmarz gossip with momentum (mRK), Algorithm 3 can be cast as special case.

Theorem 8**.**

Assume that $x^{0}=x^{1}=c$. That is, the two registers of each node have the same initial value. Then for the Algorithms 3 and 4 we have $\sum_{i=1}^{n}x_{i}^{k}=\sum_{i=1}^{n}c_{i}$ for any $k\geq 0$ and as a result, $\frac{1}{n}\sum_{i=1}^{n}x_{i}^{k}=\bar{c}$.

Proof.

We prove the result for the more general Algorithm 4. Assume that in the $k^{th}$ step of the method $q$ connected components are formed. Let the set of the nodes of each connected component be ${\cal V}_{r}^{k}$ so that ${\cal V}=\cup_{r=\{1,2,...q\}}{\cal V}_{r}^{k}$ and $|{\cal V}|=\sum_{{r}=1}^{q}|{\cal V}_{r}^{k}|$ for any $k>0$. Thus:

[TABLE]

Let us first focus, without loss of generality, on connected component $r\in[q]$ and simplify the expression for the sum of its nodes:

[TABLE]

By substituting this for all $r\in[q]$ into the right hand side of (41) and from the fact that ${\cal V}=\cup_{r\in[q]}{\cal V}_{r}^{k}$, we obtain:

[TABLE]

Since $x^{0}=x^{1}$, we have $\sum_{i=1}^{n}x_{i}^{0}=\sum_{i=1}^{n}x_{i}^{1}$, and as a result $\sum_{i=1}^{n}x_{i}^{k}=\sum_{i=1}^{n}x_{i}^{0}$ for all $k\geq 0$. ∎

4.2 Provably accelerated randomized gossip algorithms

In the results of this subsection we focus on one specific case of the Sketch and Project framework, the RK method (5). We present two accelerated variants of the randomized Kaczmarz (RK) where the Nesterov’s momentum is used, for solving consistent linear systems and we describe their theoretical convergence results. Based on these methods we propose two provably accelerated gossip protocols, along with some remarks on their implementation.

4.2.1 Accelerated Kaczmarz methods using Nesterov’s momentum

There are two different but very similar ways to provably accelerate the randomized Kaczmarz method using Nesterov’s acceleration. The first paper that proves *asymptotic convergence *with an accelerated linear rate is [41]. The proof technique is similar to the framework developed by Nesterov in [60] for the acceleration of coordinate descent methods. In [82, 24] a modified version for the selection of the parameters was proposed and a non-asymptotic accelerated linear rate was established. In Algorithm 5, pseudocode of the Accelerated Kaczmarz method (AccRK) is presented where both variants can be cast as special cases, by choosing the parameters with the correct way.

There are two options for selecting the parameters of the AccRK for solving consistent linear systems with normalized matrices, which we describe next.

From [41]: Choose $\lambda\in[0,\lambda_{\min}^{+}(\mathbf{A}^{\top}\mathbf{A})]$ and set $\gamma_{-1}=0$. Generate the sequence $\{\gamma_{k}:k=0,1,\dots,K+1\}$ by choosing $\gamma_{k}$ to be the largest root of

[TABLE]

and generate the sequences $\{\alpha_{k}:k=0,1,\dots,K+1\}$ and $\{\beta_{k}:k=0,1,\dots,K+1\}$ by setting

[TABLE] 2. 2.

From [24]: Let

[TABLE]

Choose the three sequences to be fixed constants as follows: $\beta_{k}=\beta=1-\sqrt{\frac{\lambda_{\min}^{+}(\mathbf{W})}{\nu}}$,  $\gamma_{k}=\gamma=\sqrt{\frac{1}{\lambda_{\min}^{+}(\mathbf{W})\nu}}$,   $\alpha_{k}=\alpha=\frac{1}{1+\gamma\nu}\in(0,1)$ where $\mathbf{W}=\frac{\mathbf{A}^{\top}\mathbf{A}}{m}$.

4.2.2 Theoretical guarantees of AccRK

The two variants (Option 1 and Option 2) of AccRK are closely related, however their convergence analyses are different. Below we present the theoretical guarantees of the two options as presented in [41] and [24].

Theorem 9** ([41]).**

Let $\{x^{k}\}_{k=0}^{\infty}$ be the sequence of random iterates produced by Algorithm 5 with the Option 1 for the parameters. Let $\mathbf{A}$ be normalized matrix and let $\lambda\in[0,\lambda_{\min}^{+}(\mathbf{A}^{\top}\mathbf{A})]$. Set $\sigma_{1}=1+\frac{\sqrt{\lambda}}{2m}$ and $\sigma_{2}=1-\frac{\sqrt{\lambda}}{2m}$. Then for any $k\geq 0$ we have that:

[TABLE]

Corollary 2** ([41]).**

Note that as $k\rightarrow\infty$, we have that $\sigma^{k}_{2}\rightarrow 0$. This means that the decrease of the right hand side is governed mainly by the behavior of the term $\sigma_{1}$ in the denominator and as a result the method converge asymptotically with a decrease factor per iteration: $\sigma_{1}^{-2}=(1+\frac{\sqrt{\lambda}}{2m})^{-2}\approx 1-\frac{\sqrt{\lambda}}{m}.$ That is, as $k\rightarrow\infty$:

[TABLE]

Thus, by choosing $\lambda=\lambda_{\min}^{+}$ and for the case that $\lambda_{\min}^{+}$ is small, Algorithm 5 will have significantly faster convergence rate than RK. Note that the above convergence results hold only for normalized matrices $\mathbf{A}\in\mathbb{R}^{m\times n}$, that is matrices that have $\|\mathbf{A}_{i:}\|=1$ for any $i\in m$.

Using Corollary 2, Algorithm 5 with the first choice of the parameters converges linearly with rate $\left(1-\sqrt{\lambda}/m\right)$. That is, it requires $O\left(m/\sqrt{\lambda}\log(1/\epsilon)\right)$ iterations to obtain accuracy $\mathbb{E}[\|x^{k}-x^{*}\|^{2}]\leq\epsilon 4\lambda\|x^{0}-x^{*}\|^{2}_{(\mathbf{A}^{\top}\mathbf{A})^{\dagger}}$.

Theorem 10** ([24]).**

Let $\mathbf{W}=\frac{\mathbf{A}^{\top}\mathbf{A}}{m}$ and let assume exactness888Note that in this setting $\mathbf{B}=\mathbf{I}$, which means that $\mathbf{W}=\mathbb{E}[\mathbf{Z}]$, and the exactness assumption takes the form ${\rm Null}({\bf W})={\rm Null}({\bf A})$.. Let $\{x^{k},y^{k},v^{k}\}$ be the iterates of Algorithm 5 with the Option 2 for the parameters. Then

[TABLE]

where $\Psi^{k}=\mathbb{E}\left[\|v^{k}-x^{*}\|_{{\bf W}^{\dagger}}^{2}+\frac{1}{\mu}\|x^{k}-x^{*}\|^{2}\right]$.

The above result implies that Algorithm 5 converges linearly with rate $1-\sqrt{\lambda_{\min}^{+}(\mathbf{W})/\nu}$, which translates to a total of $O\left(\sqrt{\nu/\lambda_{\min}^{+}(\mathbf{W})}\log(1/\epsilon)\right)$ iterations to bring the quantity $\Psi^{k}$ below $\epsilon>0$. It can be shown that $1\leq\nu\leq 1/\lambda_{\min}^{+}(\mathbf{W})$, (Lemma 2 in [24]) where $\nu$ is as defined in (43). Thus, $\sqrt{\frac{1}{\lambda_{\min}^{+}(\mathbf{W})}}\leq\sqrt{\frac{\nu}{\lambda_{\min}^{+}(\mathbf{W})}}\leq\frac{1}{\lambda_{\min}^{+}(\mathbf{W})},$ which means that the rate of AccRK (Option 2) is always better than that of the RK with unit stepsize which is equal to $O\left(\frac{1}{\lambda_{\min}^{+}(\mathbf{W})}\log(1/\epsilon)\right)$ (see Theorem 1).

In [24], Theorem 10 has been proposed for solving more general consistent linear systems (the matrix $\mathbf{A}$ of the system is not assumed to be normalized). In this case $\mathbf{W}=\mathbb{E}[\mathbf{Z}]$ and the parameter $\nu$ is slightly more complicated than the one of equation (43). We refer the interested reader to [24] for more details.

Comparison of the convergence rates:

Before describe the distributed nature of the AccRK and explain how it can be interpreted as a gossip algorithm, let us compare the convergence rates of the two options of the parameters for the case of general normalized consistent linear systems ($\|\mathbf{A}_{i:}\|=1$ for any $i\in[m]$).

Using Theorems 9 and 10, it is clear that the iteration complexity of AccRK is

[TABLE]

and

[TABLE]

for the Option 1 and Option 2 for the parameters, respectively.

In the following derivation we compare the iteration complexity of the two methods.

Lemma 11**.**

Let matrices $\mathbf{C}\in\mathbb{R}^{n\times n}$ and $\mathbf{C}_{i}\in\mathbb{R}^{n\times n}$ where $i\in[m]$ be positive semidefinite, and satisfying $\sum_{i=1}^{m}\mathbf{C}_{i}=\mathbf{C}$. Then

[TABLE]

Proof.

From the definition of the matrices it holds that $\mathbf{C}_{i}\preceq\mathbf{C}$ for any $i\in[m]$. Using the properties of Moore-Penrose pseudoinverse, this implies that

[TABLE]

Therefore

[TABLE]

From the definition of the matrices by taking the sum over all $i\in[m]$ we obtain:

[TABLE]

which completes the proof. ∎

Let us now choose $\mathbf{C}_{i}={\bf A}_{i:}^{\top}{\bf A}_{i:}$ and $\mathbf{C}={\bf A}^{\top}{\bf A}$. Note that from their definition the matrices are positive semidefinite and satisfy $\sum_{i=1}^{m}{\bf A}_{i:}^{\top}{\bf A}_{i:}={\bf A}^{\top}{\bf A}$. Using Lemma 11 it is clear that:

[TABLE]

or in other words, for any vector $v\notin{\rm Null}{(\mathbf{A})}$ we set the inequality

[TABLE]

Multiplying both sides by $m$, we set:

[TABLE]

Using the above derivation, it is clear from the definition of the parameter $\nu$ (43), that $\nu\leq m.$ By combining our finding with the bounds already obtained in [24] for the parameter $\nu$, we have that:

[TABLE]

Thus, by comparing the two iteration complexities of equations (44) and (45) it is clear that Option 2 for the parameters [24] is always faster in theory than Option 1 [41]. To the best of our knowledge, such comparison of the two choices of the parameters for the AccRK was never presented before.

4.2.3 Accelerated randomized gossip algorithms

Having presented the complexity analysis guarantees of AccRK for solving consistent linear systems with normalized matrices, let us now explain how the two options of AccRK behave as gossip algorithms when they are used to solve the linear system $\mathbf{A}x=0$ where $\mathbf{A}\in\mathbb{R}^{|{\cal E}|\times n}$ is the normalized incidence matrix of the network. That is, each row $e=(i,j)$ of $\mathbf{A}$ can be represented as $(\mathbf{A}_{e:})^{\top}=\frac{1}{\sqrt{2}}(e_{i}-e_{j})$ where $e_{i}$ (resp.$e_{j}$) is the $i^{th}$ (resp. $j^{th}$) unit coordinate vector in $\mathbb{R}^{n}$.

By using this particular linear system, the expression $\frac{\mathbf{A}_{i:}y^{k}-b_{i}}{\|\mathbf{A}_{i:}\|^{2}}\mathbf{A}_{i:}^{\top}$ that appears in steps 7 and 8 of AccRK takes the following form when the row $e=(i,j)\in{\cal E}$ is sampled:

[TABLE]

Recall that with $\mathbf{L}$ we denote the Laplacian matrix of the network. For solving the above AC system (see Definition 2), the standard RK requires $O\left(\left(\frac{2m}{\lambda_{\min}^{+}(\mathbf{L})}\right)\log(1/\epsilon)\right)$ iterations to achieve expected accuracy $\epsilon>0$. To understand the acceleration in the gossip framework this should be compared to the

[TABLE]

of AccRK (Option 1) and the

[TABLE]

of AccRK (Option 2).

Algorithm 6 describes in a single framework how the two variants of AccRK of Section 4.2.1 behave as gossip algorithms when are used to solve the above linear system. Note that each node $\ell\in{\cal V}$ of the network has two local registers to save the quantities $v^{k}_{\ell}$ and $x^{k}_{\ell}$. In each step using these two values every node $\ell\in{\cal V}$ of the network (activated or not) computes the quantity $y^{k}_{\ell}=\alpha_{k}v^{k}_{\ell}+(1-\alpha_{k})x^{k}_{\ell}$. Then in the $k^{th}$ iteration the activated nodes $i$ and $j$ of the randomly selected edge $e=(i,j)$ exchange their values $y^{k}_{i}$ and $y^{k}_{j}$ and update the values of $x^{k}_{i}$, $x^{k}_{j}$ and $v^{k}_{i}$, $v^{k}_{j}$ as shown in Algorithm 6. The rest of the nodes use only their own $y^{k}_{\ell}$ to update the values of $v^{k}_{i}$ and $x^{k}_{i}$ without communicate with any other node.

The parameter $\lambda^{+}_{\min}(\mathbf{L})$ can be estimated by all nodes in a decentralized manner using the method described in [8]. In order to implement this algorithm, we assume that all nodes have synchronized clocks and that they know the rate at which gossip updates are performed, so that inactive nodes also update their local values. This may not be feasible in all applications, but when it is possible (e.g., if nodes are equipped with inexpensive GPS receivers, or have reliable clocks) then they can benefit from the significant speedup achieved.

5 Dual Randomized Gossip Algorithms

An important tool in optimization literature is duality. In our setting, instead of solving the original minimization problem (primal problem) one may try to develop dual in nature methods that have as a goal to directly solve the dual maximization problem. Then the primal solution can be recovered through the use of optimality conditions and the development of an affine mapping between the two spaces (primal and dual).

Like in most optimization problems (see [54, 75, 70, 17]), dual methods have been also developed in the literature of randomized algorithms for solving linear systems [86, 26, 48]. In this section, using existing dual methods and the connection already established between the two areas of research (methods for linear systems and gossip algorithms), we present a different viewpoint that allows the development of novel dual randomized gossip algorithms.

Without loss of generality we focus on the case of $\mathbf{B}=\mathbf{I}$ (no weighted average consensus). For simplicity, we formulate the AC system as the one with the incidence matrix of the network ($\mathbf{A}=\mathbf{Q}$) and focus on presenting the distributed nature of dual randomized gossip algorithms with no momentum. While we focus only on no-momentum protocols, we note that accelerated variants of the dual methods could be easily obtained using tools from Section 4.

5.1 Dual problem and SDSA

The Lagrangian dual of the best approximation problem (7) is the (bounded) unconstrained concave quadratic maximization problem [26]:

[TABLE]

A direct method for solving the dual problem – Stochastic Dual Subspace Accent (SDSA)– was first proposed in [26]. SDSA is a randomized iterative algorithm for solving (49), which updates the dual vectors $y^{k}$ as follows:

[TABLE]

where $\mathbf{S}_{k}$ is a matrix chosen in an i.i.d. fashion throughout the iterative process from an arbitrary but fixed distribution (which is a parameter of the method) and $\lambda^{k}$ is a vector chosen afterwards so that $D(y^{k}+\mathbf{S}_{k}\lambda)$ is maximized in $\lambda$. In particular, SDSA proceeds by moving in random subspaces spanned by the random columns of ${\bf S}_{k}$.

In general, the maximizer in $\lambda$ is not unique. In SDSA, $\lambda^{k}$ is chosen to be the least-norm maximizer, which leads to the iteration

[TABLE]

It can be shown that the iterates $\{x^{k}\}_{k\geq 0}$ of the sketch and project method (Algorithm 1) can be arised as affine images of the iterates $\{y^{k}\}_{k\geq 0}$ of the dual method (50) as follows [26, 48]:

[TABLE]

In [26] the dual method was analyzed for the case of unit stepsize ($\omega=1$). Later in [48] the analysis extended to capture the cases of $\omega\in(0,2)$.

It can be shown,[26, 48], that if $\mathbf{S}_{k}$ is chosen randomly from the set of unit coordinate/basis vectors in $\mathbb{R}^{m}$, then the dual method (50) is randomized coordinate descent [38, 72], and the corresponding primal method is RK (5). More generally, if $\mathbf{S}_{k}$ is a random column submatrix of the $m\times m$ identity matrix, the dual method is the randomized Newton method [69], and the corresponding primal method is RBK (6). In Section 5.2 we shall describe the more general block case in more detail.

The basic convergence result for the dual iterative process (50) is presented in the following theorem. We set $y^{0}=0$ so that $x^{0}=c$ and the affine mapping of equation (51) is satisfied. Note that, SDSA for solving the dual problem (49) and sketch and project method for solving the best approximation problem (7) (primal problem) converge exactly with the same convergence rate.

Theorem 12** (Complexity of SDSA [26, 48]).**

Assume exactness and let $y^{0}=0$. Let $\{y^{k}\}$ be the iterates produced by SDSA (equation (50)) with step-size $\omega\in(0,2)$ . Then,

[TABLE]

where $\rho:=1-\omega(2-\omega)\lambda_{\min}^{+}\left(\mathbf{W}\right)\in(0,1)$.

5.2 Randomized Newton method as a dual gossip algorithm

In this subsection we bring a new insight into the randomized gossip framework by presenting how the dual iterative process that is associated to RBK method solves the AC problem with $\mathbf{A}=\mathbf{Q}$ (incidence matrix). Recall that the right hand side of the linear system is $b=0$. For simplicity, we focus on the case of $\mathbf{B}=\mathbf{I}$ and $\omega=1$.

Under this setting ($\mathbf{A}=\mathbf{Q}$, $\mathbf{B}=\mathbf{I}$ and $\omega=1$) the dual iterative process (50) takes the form:

[TABLE]

and from Theorem 12 converges to a solution of the dual problem as follows:

[TABLE]

Note that the convergence rate is exactly the same with the rate of the RBK under the same assumptions (see (30)).

This algorithm is a randomized variant of the Newton method applied to the problem of maximizing the quadratic function $D(y)$ defined in (49). Indeed, in each iteration we perform the update $y^{k+1}=y^{k}+\mathbf{I}_{C:}\lambda^{k}$, where $\lambda^{k}$ is chosen greedily so that $D(y^{k+1})$ is maximized. In doing so, we invert a random principal submatrix of the Hessian of $D$, whence the name.

Randomized Newton Method (RNM) was first proposed by Qu et al. [69]. RNM was first analyzed as an algorithm for minimizing smooth strongly convex functions. In [26] it was also extended to the case of a smooth but weakly convex quadratics. This method was not previously associated with any gossip algorithm.

The most important distinction of RNM compared to existing gossip algorithms is that it operates with values that are associated to the edges of the network. To the best of our knowledge, it is the first randomized dual gossip method999Since our first paper [46] appeared online, papers [29] and [32] propose dual randomized gossip algorithms as well. In [29] the authors focus on presenting privacy preserving dual randomized gossip algorithms and [32] solve the dual problem (49) using accelerated coordinate descent.. In particular, instead of iterating over values stored at the nodes, RNM uses these values to update “dual weights” $y^{k}\in\mathbb{R}^{m}$ that correspond to the edges ${\cal E}$ of the network. However, deterministic dual distributed averaging algorithms were proposed before [71, 20]. Edge-based methods have also been proposed before; in particular in [84] an asynchronous distributed ADMM algorithm presented for solving the more general consensus optimization problem with convex functions.

Natural Interpretation. In iteration $k$, RNM (Algorithm (53)) executes the following steps:

1. Select a random set of edges ${\cal S}_{k}\subseteq{\cal E}$, 2) Form a subgraph ${\cal G}_{k}$ of ${\cal G}$ from the selected edges, 3) The values of the edges in each connected component of ${\cal G}_{k}$ are updated: their new values are a linear combination of the private values of the nodes belonging to the connected component and of the adjacent edges of their connected components. (see also example of Figure 4).

Dual Variables as Advice. The weights $y^{k}$ of the edges have a natural interpretation as advice that each selected node receives from the network in order to update its value (to one that will eventually converge to the desired average).

Consider RNM performing the $k^{th}$ iteration and let ${\cal V}_{r}$ denote the set of nodes of the selected connected component that node $i$ belongs to. Then, from Theorem 5 we know that $x_{i}^{k+1}=\sum_{i\in{\cal V}_{r}}x_{i}^{k}/|{\cal V}_{r}|$. Hence, by using (51), we obtain the following identity:

[TABLE]

Thus in each step $(\mathbf{A}^{\top}y^{k+1})_{i}$ represents the term (advice) that must be added to the initial value $c_{i}$ of node $i$ in order to update its value to the average of the values of the nodes of the connected component $i$ belongs to.

Importance of the dual perspective: It was shown in [69] that when RNM (and as a result, RBK, through the affine mapping (51)) is viewed as a family of methods indexed by the size $\tau=|{\cal S}|$ (we choose ${\cal S}$ of fixed size in the experiments), then $\tau\to 1/(1-\rho)$, where $\rho$ is defined in (10), decreases superlinearly fast in $\tau$. That is, as $\tau$ increases by some factor, the iteration complexity drops by a factor that is at least as large. Through preliminary numerical experiments in Section 7.2 we experimentally show that this is true for the case of AC systems as well.

6 Further Connections Between Methods for Solving Linear Systems and Gossip Algorithms

In this section we highlight some further interesting connections between linear systems solvers and gossip protocols for average consensus:

• •

Eavesdrop gossip as special case of Kaczmarz-Motzkin method. In [83] greedy gossip with eavesdropping (GGE), a novel randomized gossip algorithm for distributed computation of the average consensus problem was proposed and analyzed. In particular it was shown that that greedy updates of GGE lead to rapid convergence. In this protocol, the greedy updates are made possible by exploiting the broadcast nature of wireless communications. During the operation of GGE, when a node decides to gossip, instead of choosing one of its neighbors at random, it makes a greedy selection, choosing the node which has the value most different from its own. In particular the method behaves as follows:

At the $k^{th}$ iteration of GGE, a node $i_{k}$ is chosen uniformly at random from $[n]$. Then, $i_{k}$ identifies a neighboring node $j_{k}\in{\cal N}_{i}$ satisfying:

[TABLE]

which means that the selected node $i_{k}$ identifies a neighbor that currently has the most different value from its own. This choice is possible because each node $i\in{\cal V}$ maintains not only its own local variable $x_{i}^{k}$, but also a copy of the current values at its neighbors $x_{j}^{k}$ for $j\in{\cal N}_{i}$. In the case that node $i_{k}$ has multiple neighbors whose values are all equally (and maximally) different from its current value, it chooses one of these neighbors at random. Then node $i_{k}$ and $j_{k}$ update their values to:

[TABLE]

In the area of randomized methods for solving large linear system there is one particular method, the Kaczmarz-Motzkin algorithm [12, 28] that can work as gossip algorithm with the same update as the GGE when is use to solve the homogeneous linear system with matrix the Incidence matrix of the network.

Update rule of Kaczmarz-Motzkin algorithm (KMA) [12, 28]:

Choose sample of $d_{k}$ constraints, $P_{k}$, uniformly at random from among the rows of matrix $\mathbf{A}$. 2. 2.

From among these $d_{k}$ constraints, choose $t_{k}=\text{argmax}_{i\in P_{k}}\mathbf{A}_{i:}x^{k}-b_{i}.$ 3. 3.

Update the value: $x^{k+1}=x^{k}-\frac{\mathbf{A}_{t_{k}:}x^{k}-b_{i}}{\|\mathbf{A}_{t_{k}:}\|^{2}}\mathbf{A}_{t_{k}:}^{\top}.$

It is easy to verify that when the Kaczmarz-Motzkin algorithm is used for solving the AC system with $\mathbf{A}=\mathbf{Q}$ (incidence matrix) and in each step of the method the chosen constraints $d_{k}$ of the linear system correspond to edges attached to one node it behaves exactly like the GGE. From numerical analysis viewpoint an easy way to choose the constraints $d_{k}$ that are compatible to the desired edges is in each iteration to find the indexes of the non-zeros of a uniformly at random selected column (node) and then select the rows corresponding to these indexes.

Therefore, since GGE [83] is a special case of the KMA (when the later applied to special AC system with Incidence matrix) it means that we can obtain the convergence rate of GGE by simply use the tight conergence analysis presented in [12, 28] 101010Note that the convergence theorems of [12, 28] use $d_{k}=d$. However, with a small modification in the original proof the theorem can capture the case of different $d_{k}$.. In [83] it was mentioned that analyzing the convergence behavior of GGE is non-trivial and not an easy task. By establishing the above connection the convergence rates of GGE can be easily obtained as special case of the theorems presented in [12].

In Section 4 we presented provably accelerated variants of the pairwise gossip algorithm and of its block variant. Following the same approach one can easily develop accelerated variants of the GGE using the recently proposed analysis for the accelerated Kaczmarz-Motzkin algorithm presented in [52].

• •

Inexact Sketch and Project Methods:

Recently in [50], several inexact variants of the sketch and project method (1) have been proposed. As we have already mentioned the sketch and project method is a two step procedure algorithm where first the sketched system is formulated and then the last iterate $x^{k}$ is exactly projected into the solution set of the sketched system. In [50] the authors replace the exact projection with an inexact variant and suggested to run a different algorithm (this can be the sketch and project method itself) in the sketched system to obtain an approximate solution. It was shown that in terms of time the inexact updates can be faster than their exact variants.

In the setting of randomized gossip algorithms for the AC system with Incidence matrix ($\mathbf{A}=\mathbf{Q}$) , $\mathbf{B}=\mathbf{I}$ and $\omega=1$ a variant of the inexact sketch and project method will work as follows (similar to the update proved in Theorem 5):

Select a random set of edges $C\subseteq{\cal E}$. 2. 2.

Form subgraph ${\cal G}_{k}$ of ${\cal G}$ from the selected edges. 3. 3.

Run the pairwise gossip algorithm of [6] (or any variant of the sketch and project method) on the subgraph ${\cal G}_{k}$ until an accuracy $\epsilon$ is achieved (reach a neighborhood of the exact average).

• •

Non-randomized gossip algorithms as special cases of Kaczmarz methods:

In the gossip algorithms literature there are efficient protocols that are not randomized[53, 31, 42, 89]. Typically, in these algorithms the pairwise exchanges between nodes it happens in a deterministic, such as predefined cyclic, order. For example, $T$-periodic gossiping is a protocol which stipulates that each node must interact with each of its neighbours exactly once every $T$ time units. It was shown that under suitable connectivity assumptions of the network $\mathcal{G}$, the $T$-periodic gossip sequence will converge at a rate determined by the magnitude of the second largest eigenvalue of the stochastic matrix determined by the sequence of pairwise exchanges which occurs over a period. It has been shown that if the underlying graph is a tree, the mentioned eigenvalue is constant for all possible $T$-periodic gossip protocols.

In this work we focus only on randomized gossip protocols. However we speculate that the above non-randomized gossip algorithms would be able to express as special cases of popular non-randomized projection methods for solving linear systems [68, 62, 16]. Establishing connections like that is an interesting future direction of research and can possibly lead to the development of novel block and accelerated variants of many non-randomized gossip algorithms, similar to the protocols we present in Sections 3 and 4.

7 Numerical Evaluation

In this section, we empirically validate our theoretical results and evaluate the performance of the proposed randomized gossip algorithms. The section is divided into four main parts, in each of which we highlight a different aspect of our contributions.

In the first experiment, we numerically verify the linear convergence of the Scaled RK algorithm (see equation (14)) for solving the weighted average consensus problem presented in Section 3.1. In the second part, we explain the benefit of using block variants in the gossip protocols where more than two nodes update their values in each iteration (protocols presented in Section 3.4). In the third part, we explore the performance of the faster and provably accelerated gossip algorithms proposed in Section 4. In the last experiment, we numerically show that relaxed variants of the pairwise randomized gossip algorithm converge faster than the standard randomized pairwise gossip with unit stepsize (no relaxation). This gives a specific setting where the phenomenon of over-relaxation of iterative methods for solving linear systems is beneficial.

In the comparison of all gossip algorithms we use the relative error measure $\|x^{k}-x^{*}\|_{\mathbf{B}}^{2}/\|x^{0}-x^{*}\|_{\mathbf{B}}^{2}$ where $x^{0}=c\in\mathbb{R}^{n}$ is the starting vector of the values of the nodes and matrix $\mathbf{B}$ is the positive definite diagonal matrix with weights in its diagonal (recall that in the case of standard average consensus this can be simply $\mathbf{B}=\mathbf{I}$). Depending on the experiment, we choose the values of the starting vector $c\in\mathbb{R}^{n}$ to follow either a Gaussian distribution or uniform distribution or to be integer values such that $c_{i}=i\in\mathbb{R}$. In the plots, the horizontal axis represents the number of iterations except in the figures of subsection 7.2, where the horizontal axis represents the block size.

In our implementations we use three popular graph topologies from the area of wireless sensor networks. These are the cycle (ring graph), the 2-dimension grid and the random geometric graph (RGG) with radius $r=\sqrt{\log(n)/n}$. In all experiments we formulate the average consensus problem (or its weighted variant) using the incidence matrix. That is, $\mathbf{A}=\mathbf{Q}$ is used as the AC system. Code was written in Julia 0.6.3.

7.1 Convergence on weighted average consensus

As we explained in Section 3, the sketch and project method (Algorithm 1) can solve the more general weighted AC problem. In this first experiment we numerically verify the linear convergence of the Scaled RK algorithm (14) for solving this problem in the case of $\mathbf{B}=\mathbf{D}$. That is, the matrix $\mathbf{B}$ of the weights is the degree matrix $\mathbf{D}$ of the graph ($\mathbf{B}_{ii}=d_{i}$, $\forall i\in[n]$). In this setting the exact update rule of the method is given in equation (17), where in order to have convergence to the weighted average the chosen nodes are required to share not only their private values but also their weight $\mathbf{B}_{ii}$ (in our experiment this is equal to the degree of the node $d_{i}$). In this experiment the starting vector of values $x^{0}=c\in\mathbb{R}^{n}$ is a Gaussian vector. The linear convergence of the algorithm is clear in Figure 5.

7.2 Benefit of block variants

We devote this experiment to evaluate the performance of the randomized block gossip algorithms presented in Sections 3.4 and 5. In particular, we would like to highlight the benefit of using larger block size in the update rule of randomized Kaczmarz method and as a result through our established connection of the randomized pairwise gossip algorithm [6] (see equation (18)).

Recall that in Section 5 we show that both RBK and RNM converge to the solution of the primal and dual problems respectively with the same rate and that their iterates are related via a simple affine transform (51). In addition note that an interesting feature of the RNM [69], is that when the method viewed as algorithm indexed by the size $\tau=|C|$, it enjoys superlinear speedup in $\tau$. That is, as $\tau$ (block size) increases by some factor, the iteration complexity drops by a factor that is at least as large (see Section 5.2). Since RBK and RNM share the same rates this property naturally holds for RBK as well.

We show that for a connected network ${\cal G}$, the complexity improves superlinearly in $\tau=|C|$, where $C$ is chosen as a subset of ${\cal E}$ of size $\tau$, uniformly at random (recall the in the update rule of RBK the random matrix is $\mathbf{S}=\mathbf{I}_{:C}$). Similar to the rest of this section in comparing the number of iterations for different values of $\tau$, we use the relative error $\varepsilon=\|x^{k}-x^{*}\|^{2}/\|x^{0}-x^{*}\|^{2}$. We let $x^{0}_{i}=c_{i}=i$ for each node $i\in{\cal V}$ (vector of integers). We run RBK until the relative error becomes smaller than $0.01$. The blue solid line in the figures denotes the actual number of iterations (after running the code) needed in order to achieve $\varepsilon\leq 10^{-2}$ for different values of $\tau$. The green dotted line represents the function $f(\tau):=\frac{\ell}{\tau}$, where $\ell$ is the number of iterations of RBK with $\tau=1$ (i.e., the pairwise gossip algorithm). The green line depicts linear speedup; the fact that the blue line (obtained through experiments) is below the green line points to superlinear speedup. In this experiment we use the Cycle graph with $n=30$ and $n=100$ nodes (Figure 6) and the $4\times 4$ two dimension grid graph (Figure 7). Note that, when $|C|=m$ the convergence rate of the method becomes $\rho=0$ and as a result it converges in one step.

7.3 Accelerated gossip algorithms

We devote this subsection to experimentally evaluate the performance of the proposed accelerated gossip algorithms: mRK (Algorithm 3), mRBK (Algorithm 4) and AccGossip with the two options of the parameters (Algorithm 6). In particular we perform four experiments. In the first two we focus on the performance of the mRK and how the choice of stepsize (relaxation parameter) $\omega$ and heavy ball momentum parameter $\beta$ affect the performance of the method. In the next experiment we show that the addition of heavy ball momentum can be also beneficial for the performance of the block variant mRBK. In the last experiment we compare the standard pairwise gossip algorithm (baseline method) from [6], the mRK and the AccGossip and show that the probably accelerated gossip algorithm, AccGossip outperforms the other algorithms and converge as predicted from the theory with an accelerated linear rate.

7.3.1 Impact of momentum parameter on mRK

As we have already presented in the standard pairwise gossip algorithm (equation (18)) the two selected nodes that exchange information update their values to their exact average while all the other nodes remain idle. In our framework this update can be cast as special case of mRK when $\beta=0$ and $\omega=1$.

In this experiment we keep the stepsize fixed and equal to $\omega=1$ which means that the pair of the chosen nodes update their values to their exact average and we show that by choosing a suitable momentum parameter $\beta\in(0,1)$ we can obtain faster convergence to the consensus for all networks under study. The momentum parameter $\beta$ is chosen following the suggestions made in [48] for solving general consistent linear systems. See Figure 8 for more details. It is worth to point out that for all networks under study the addition of a heavy ball momentum term is beneficial in the performance of the method.

7.3.2 Comparison of mRK and shift-Register algorithm [40]

In this experiment we compare mRK with the shift register gossip algorithm (pairwise momentum method, abbreviation: Pmom) analyzed in [40]. We choose the parameters $\omega$ and $\beta$ of mRK in such a way in order to satisfy the connection established in Section 4.1.3. That is, we choose $\beta=\omega-1$ for any choice of $\omega\in(1,2)$. Observe that in all plots of Figure 9 mRK outperforms the corresponding shift-register algorithm.

7.3.3 Impact of momentum parameter on mRBK

In this experiment our goal is to show that the addition of heavy ball momentum accelerates the RBK gossip algorithm presented in Section 3.4. Without loss of generality we choose the block size to be equal to $\tau=5$. That is, the random matrix $\mathbf{S}_{k}\sim{\cal D}$ in the update rule of mRBK is a $m\times 5$ column submatrix of the indetity $m\times m$ matrix. Thus, in each iteration $5$ edges of the network are chosen to form the subgraph ${\cal G}_{k}$ and the values of the nodes are updated according to Algorithm 4. Note that similar plots can be obtained for any choice of block size. We run all algorithms with fixed stepsize $\omega=1$. From Figure 10, it is obvious that for all networks under study, choosing a suitable momentum parameter $\beta\in(0,1)$ gives faster convergence than having no momentum, $\beta=0$.

7.3.4 Performance of AccGossip

In the last experiment on faster gossip algorithms we evaluate the performance of the proposed provably accelerated gossip protocols of Section 4.2. In particular we compare the standard RK (pairwise gossip algorithm of [6]) the mRK (Algorithm 3) and the AccGossip (Algorithm 6) with the two options for the selection of the parameters presented in Section 4.2.1.

The starting vector of values $x^{0}=c$ is taken to be a Gaussian vector. For the implementation of mRK we use the same parameters with the ones suggested in the stochastic heavy ball (SGB) setting in [48]. For the AccRK (Option 1) we use $\lambda=\lambda_{\min}^{+}(\mathbf{A}^{\top}\mathbf{A})$ and for AccRK (Option 2) we select $\nu=m$111111For the networks under study we have $m<\frac{1}{\lambda_{\min}^{+}(\mathbf{W})}$. Thus, by choosing $\nu=m$ we select the pessimistic upper bound of the parameter (48) and not its exact value (43). As we can see from the experiments, the performance is still accelerated and almost identical to the performance of AccRK (Option 1) for this choice of $\nu$.. From Figure 11 it is clear that for all networks under study the two randomized gossip protocols with Nesterov momentum are faster than both the pairwise gossip algorithm of [6] and the mRK/SHB (Algorithm 3). To the best of our knowledge Algorithm 6 (Option 1 and Option 2) is the first randomized gossip protocol that converges with provably accelerated linear rate and as we can see from our experiment its faster convergence is also obvious in practice.

7.4 Relaxed randomized gossip without momentum

In the area of randomized iterative methods for linear systems it is know that over-relaxation (using of larger step-sizes) can be particularly helpful in practical scenarios. However, to the best of our knowledge there is not theoretical justification of why this is happening.

In our last experiment we explore the performance of relaxed randomized gossip algorithms ($\omega\neq 1$) without momentum and show that in this setting having larger stepsize can be particularly beneficial.

As we mentioned before (see Theorem 1) the sketch and project method (Algorithm 1) converges with linear rate when the step-size (relaxation parameter) of the method is $\omega\in(0,2)$ and the best theoretical rate is achieved when $\omega=1$. In this experiment we explore the performance of the standard pairwise gossip algorithm when the step-size of the update rule is chosen in $(1,2)$. Since there is no theoretical proof of why over-relaxation can be helpful we perform the experiments using different starting values of the nodes. In particular we choose the values of vector $c\in R^{n}$ to follow (i) Gaussian distribution, (ii) Uniform Distribution and (iii) to be integers values such that $c_{i}=i\in R$. Our findings are presented in Figure 12. Note that for all networks under study and for all choices of starting values having larger stepsize, $\omega\in(1,2)$ can lead to better performance. Interesting observation from Figure 12 is that the stepsizes $\omega=1.8$ and $\omega=1.9$ give the best performance (among the selected choices of stepsizes) for all networks and for all choices of starting vector $x^{0}=c$.

8 Conclusion and Future Direction of Research

In this work we present a general framework for the analysis and design of randomized gossip algorithms. Using tools from numerical linear algebra and the area of randomized projection methods for solving linear systems we propose novel serial, block and accelerated gossip protocols for solving the average consensus and weighted average consensus problems.

We believe that this work could open up several future avenues for research. Using similar approach with the one presented in this manuscript, many popular projection methods can be interpreted as gossip algorithms when used to solve linear systems encoding the underlying network. This can lead to the development of novel distributed protocols for average consensus.

In addition, we speculate that the gossip protocols presented in this work can be extended to the more general setting of multi-agent consensus optimization where the goal is to minimize the average of convex or non-convex functions $\frac{1}{n}\sum_{i=1}^{n}f_{i}(x)$ in a decentralized way [55].

9 Acknowledgments

The authors would like to thank Mike Rabbat for useful discussions related to the literature of gossip algorithms and for his comments during the writing of this paper.

Appendix A Missing Proofs

A.1 Proof of Theorem 4

Proof.

Let $z^{k}:=\|x^{k}-x^{*}\|$ , $x^{0}=c$ is the starting point and $\rho$ is as defined in (10). From Theorem 1 we know that sketch and project method converges with

[TABLE]

where $x^{*}$ is the solution of (7). Inequality (55), together with Markov inequality can be used to give the following bound

[TABLE]

Therefore, as long as $k$ is large enough so that $\rho^{k}\leq\varepsilon^{3}$, we have $\mathbb{P}\left(z^{k}/z^{0}\geq\varepsilon^{2}\right)\leq\varepsilon$. That is, if

[TABLE]

then:

[TABLE]

Hence, an upper bound for value $T_{ave}(\epsilon)$ can be obtained as follows,

[TABLE]

where in last inequality we use $1/\log(1/\rho)\leq 1/1-\rho$ which is true because $\rho\in(0,1)$. ∎

A.2 Proof of Theorem 5

Proof.

The following notation conventions are used in this proof. With $q_{k}$ we indicate the number of connected components of subgraph ${\cal G}_{k}$, while with ${\cal V}_{r}$ we denote the set of nodes of each connected component $q_{k}$ $(r\in\{1,2,\dots,q_{k}\})$. Finally, $|{\cal V}_{r}|$ shows the cardinality of set ${\cal V}_{r}$. Notice that, if ${\cal V}$ is the set of all nodes of the graph then ${\cal V}=\underset{r=\{1,2,...q\}}{\cup}{\cal V}_{r}$ and $|{\cal V}|=\sum\limits_{{r}=1}^{q}|{\cal V}_{r}|$.

Note that from equation (29), the update of RBK for $\mathbf{A}=\mathbf{Q}$(Incidence matrix) can be expressed as follows:

[TABLE]

Notice that $\mathbf{I}_{:C}^{\top}\mathbf{Q}$ is a row submatrix of matrix $\mathbf{Q}$ with rows those that correspond to the random set $C\subseteq{\cal E}$ of the edges. From the expression of matrix $\mathbf{Q}$ we have that

[TABLE]

Now, using this, it can be seen that the constraint $\mathbf{I}_{:C}^{\top}\mathbf{Q}x=0$ of problem (7) is equivalent to $q$ equations (number of connected components) where each one of them forces the values $x_{i}^{k+1}$ of the nodes $i\in{\cal V}_{r}$ to be equal. That is, if we use $z_{r}$ to represent the value of all nodes that belong in the connected component $r$ then:

[TABLE]

and the constrained optimization problem (7) can expressed as unconstrained as follows:

[TABLE]

where $z=(z_{1},z_{2},\dots,z_{q})\in\mathbb{R}^{q}$ is the vector of all values $z_{r}$ when $r\in\{1,2,\dots,q\}$. Since our problem is unconstrained the minimum of equation (60) is obtained when $\nabla\phi^{k}(z)=0$.

By evaluating partial derivatives of (60) we obtain:

[TABLE]

As a result,

[TABLE]

Thus from (59), the value of each node $i\in{\cal V}_{r}$ will be updated to

[TABLE]

. ∎

Appendix B Notation Glossary