Distributed heavy-ball: A generalization and acceleration of first-order methods with gradient tracking

I Introduction

We consider distributed optimization, where $n$ agents collaboratively solve the following problem:

Early work on this topic builds on the seminal work by Tsitsiklis in and includes Distributed Gradient Descent (DGD) and distributed dual averaging over undirected graphs. Leveraging push-sum consensus , Refs. extend the DGD framework to directed graphs. Based on a similar concept, Refs. propose Directed-Distributed Gradient Descent (D-DGD) for directed graphs that is based on surplus consensus . In general, the DGD-based methods achieve sublinear convergence at $\mathcal{O}\left(\frac{\log k}{\sqrt{k}}\right)$, where $k$ is the number of iterations, because of the diminishing step-size used in the iterations. The convergence rate of DGD can be improved with the help of a constant step-size but at the expense of an inexact solution . Follow-up work also includes augmented Lagrangians , which shows exact linear convergence for smooth and strongly-convex functions, albeit requiring higher computation at each iteration.

To improve convergence and retain computational simplicity, fast first-order methods that do not (explicitly) use a dual update have been proposed. Reference describes a distributed Nesterov-type method based on multiple consensus inner loops, at $\mathcal{O}\left(\frac{\log k}{k^{2}}\right)$ for smooth and convex functions, with bounded gradients. EXTRA uses the difference of two consecutive DGD iterates to achieve an $\mathcal{O}\left(\frac{1}{k}\right)$ rate for arbitrary convex functions and a $Q$-linear rate for strongly-convex functions. DEXTRA combines push-sum and EXTRA to achieve an $R$-linear rate over directed graphs given that a constant step-size is carefully chosen in some interval. Refs. apply an adapt-then-combine structure to EXTRA and generalize the symmetric weights in EXTRA to row-stochastic, over undirected graphs.

Noting that DGD-type methods are faster with a constant step-size, recent work uses a constant step-size and replaces the local gradient, at each agent in DGD, with an estimate of the global gradient. A method based on gradient tracking was first shown in over undirected graphs, which proposes Aug-DGM (that uses nonidentical step-sizes at the agents) with the help of dynamic consensus and shows convergence for smooth convex functions. When the step-sizes are identical, the convergence rate of Aug-DGM was derived to be $\mathcal{O}\left(\frac{1}{k}\right)$ for arbitrary convex functions and $R$-linear for strongly-convex functions in . ADD-OPT extends to directed graphs by combining push-sum with gradient tracking and derives a contraction in an arbitrary norm to establish an $R$-linear convergence rate when the global objective is smooth and strongly-convex. Ref. extends the analysis in to time-varying graphs and establishes an $R$-linear convergence using the small gain theorem . In contrast to the aforementioned methods , where the weights are doubly-stochastic for undirected graphs and column-stochastic for directed graphs, FROST uses row-stochastic weights, which have certain advantages over column-stochastic weights. Ref. unifies EXTRA and gradient tracking methods in a primal-dual framework over static undirected graphs. More recently, Ref. proposes distributed Nesterov over undirected graphs that also uses gradient tracking and shows a convergence rate of $\mathcal{O}((1-{c}{\mathcal{Q}^{-\frac{5}{7}}})^{k})$ for smooth, strongly-convex functions, where $\mathcal{Q}$ is the condition number of the global objective. Refs. , on the other hand, consider gradient tracking in distributed non-convex problems, while Ref. uses second-order information to accelerate the convergence.

Of significant relevance here is the $\mathcal{AB}$ algorithm , also appeared later in , which can be viewed as a generalization of distributed first-order methods with gradient tracking. In particular, the algorithms over undirected graphs in Refs. are a special case of $\mathcal{AB}$ because the doubly-stochastic weights therein are replaced by row- and column- stochastic weights. $\mathcal{AB}$ thus is naturally applicable to arbitrary directed graphs. Moreover, the use of both row- and column-stochastic weights removes the need for eigenvector estimationSimultaneous application of both row- and column-stochastic weights was first employed for average-consensus in and towards distributed optimization in , albeit without gradient tracking., required earlier in . Ref. derives an $R$-linear rate for $\mathcal{AB}$ when the objective functions are smooth and strongly-convex. In this paper, we provide an improved understanding of $\mathcal{AB}$ and extend it to the $\mathcal{AB}m$ algorithm, a distributed heavy-ball method, applicable to both undirected and directed graphs. We now summarize the main contributions:

We show that many of the existing accelerated first-order methods are either a special case of $\mathcal{AB}$ , or can be adapted from its equivalent forms .

We propose a distributed heavy-ball method, termed as $\mathcal{AB}m$, that combines $\mathcal{AB}$ with a heavy-ball (type) momentum term. To the best of our knowledge, this paper is the first to use a momentum term based on the heavy-ball method in distributed optimization.

$\mathcal{AB}m$ employs nonidentical step-sizes at the agents and thus its analysis naturally carries to nonidentical step-sizes in $\mathcal{AB}$ and to the related algorithms in .

We cast a unifying framework for consensus over arbitrary graphs that results from $\mathcal{AB}m$ and subsumes several well-known algorithms .

On the analysis front, we show that $\mathcal{AB}$ (without momentum) converges faster as compared to the algorithms over directed graphs in , where separate iterations for eigenvector estimation are applied nonlinearly to the underlying algorithm. Towards $\mathcal{AB}m$, we establish a global $R$-linear convergence rate for smooth and strongly-convex objective functions when the largest step-size at the agents is positive and sufficiently small. This is in contrast to the earlier work on non-identical step-sizes within the framework of gradient tracking , which requires the heterogeneity among the step-sizes to be sufficiently small, i.e., the step-sizes are close to each other. We also acknowledge that similar to the centralized heavy-ball method , dating back to more than 50 years, and the recent work , a global acceleration can only be shown via numerical simulations. Following the standard practice, we provide simulations to verify that $\mathcal{AB}m$ has accelerated convergence, the effect of which is more pronounced when the global objective function is ill-conditioned.

We now describe the rest of the paper. Section II provides preliminaries, problem formulation, and introduces distributed heavy-ball, i.e., the $\mathcal{AB}m$ algorithm. Section III establishes the connection between $\mathcal{AB}$ and related algorithms. Section IV includes the main results on the convergence analysis, whereas Section V provides a family of average-consensus algorithms that result naturally from $\mathcal{AB}m$. Finally, Section VI provides numerical experiments and Section VII concludes the paper.

Basic Notation: We use lowercase bold letters to denote vectors and uppercase letters for matrices. The matrix, $I_{n}$, is the $n\times n$ identity, whereas $\mathbf{1}_{n}$ ($\mathbf{0}_{n}$) is the $n$-dimensional column vector of all ones (zeros). For an arbitrary vector, $\mathbf{x}$, we denote its $i$th element by $[\mathbf{x}]_{i}$ and its largest and smallest element by $[\mathbf{x}]_{\max}$ and $[\mathbf{x}]_{\min}$, respectively. We use $\mbox{diag}(\mathbf{x})$ to denote a diagonal matrix that has $\mathbf{x}$ on its main diagonal. For two matrices, $X$ and $Y$, $\mbox{diag}\left(X,Y\right)$ is a block-diagonal matrix with $X$ and $Y$ on its main diagonal, and $X\otimes Y$ denotes their Kronecker product. The spectral radius of a matrix, $X$, is represented by $\rho(X)$. For a primitive, row-stochastic matrix, $A$, we denote its left and right eigenvectors corresponding to the eigenvalue of $1$ by $\boldsymbol{\pi}_{r}$ and $\mathbf{1}_{n}$, respectively, such that $\boldsymbol{\pi}_{r}^{\top}\mathbf{1}_{n}=1$; similarly, for a primitive, column-stochastic matrix, $B$, we denote its left and right eigenvectors corresponding to the eigenvalue of $1$ by $\mathbf{1}_{n}$ and $\boldsymbol{\pi}_{c}$, respectively, such that $\mathbf{1}_{n}^{\top}\boldsymbol{\pi}_{c}=1$. For a matrix $X$, we denote $X_{\infty}$ as its infinite power (if it exists), i.e., $X_{\infty}=\lim_{k\rightarrow\infty}X^{k}.$ From the Perron-Frobenius theorem , we have $A_{\infty}=\mathbf{1}_{n}\boldsymbol{\pi}_{r}^{\top}$ and $B_{\infty}=\boldsymbol{\pi}_{c}\mathbf{1}_{n}^{\top}$. We denote $\left\|\cdot\right\|_{\mathcal{A}}$ and $\left\|\cdot\right\|_{\mathcal{B}}$ as some arbitrary vector norms, the choice of which will be clear in Lemma 1, while $\left\|\cdot\right\|$ denotes the Euclidean matrix and vector norms.

II Preliminaries and Problem Formulation

Consider $n$ agents connected over a directed graph, $\mathcal{G}=(\mathcal{V},\mathcal{E})$, where $\mathcal{V}=\{1,\cdots,n\}$ is the set of agents, and $\mathcal{E}$ is the collection of ordered pairs, $(i,j),i,j\in\mathcal{V}$, such that agent $j$ can send information to agent $i$, i.e., $j\rightarrow i$. We define $\mathcal{N}_{i}^{{\scriptsize\mbox{in}}}$ as the collection of in-neighbors of agent $i$, i.e., the set of agents that can send information to agent $i$. Similarly, $\mathcal{N}_{i}^{{\scriptsize\mbox{out}}}$ is the set of out-neighbors of agent $i$. Note that both $\mathcal{N}_{i}^{{\scriptsize\mbox{in}}}$ and $\mathcal{N}_{i}^{{\scriptsize\mbox{out}}}$ include agent $i$. The agents solve the following problem:

The graph, $\mathcal{G}$, is strongly-connected.

where $\mu_{i}\geq 0$ and $\sum_{i=1}^{n}\mu_{i}>0$.

It is well known that the best achievable convergence rate of the gradient descent algorithm,

is $\mathcal{O}((\tfrac{\mathcal{Q}-1}{\mathcal{Q}+1})^{k})$, where $\mathcal{Q}\triangleq\tfrac{l}{\mu}$ is the condition number of the objective function, $F$. Clearly, gradient descent is quite slow when $\mathcal{Q}$ is large, i.e., when the objective function is ill-conditioned. The seminal work by Polyak proposes the following heavy-ball method:

where $\beta\left(\mathbf{x}_{k}-\mathbf{x}_{k-1}\right)$ is interpreted as a “momentum” term, used to accelerate the convergence process. Polyak shows that with a specific choice of $\alpha$ and $\beta$, the heavy-ball method achieves a local accelerated rate of $\mathcal{O}((\tfrac{\mathcal{\sqrt{Q}}-1}{\mathcal{\sqrt{Q}}+1})^{k})$. By local, it is meant that the acceleration can only be analytically shown when $\|\mathbf{x}_{0}-\mathbf{x}^{*}\|$ is sufficiently small. Globally, i.e., for arbitrary initial conditions, only linear convergence is established, while an analytical characterization of the acceleration is still an open problem, see related work in . Numerical analysis and simulations are often employed to show global acceleration, i.e., it is possible to tune $\alpha$ and $\beta$ such that the heavy-ball method is faster than gradient descent .

II-B Distributed heavy-ball: The 𝒜​ℬ​m𝒜ℬ𝑚\mathcal{AB}m algorithm

where $\alpha_{i}\geq 0$ and $\beta_{i}\geq 0$ are respectively the local step-size and the momentum parameter adopted by agent $i$. The weights, $a_{ij}$’s and $b_{ij}$’s, are associated with the graph topology and satisfy the following conditions:

Note that the weight matrix, $A=\{a_{ij}\}$, in Eq. (2a) is RS (row-stochastic) and the weight matrix, $B=\{b_{ij}\}$ in Eq. (2b) is CS (column-stochastic), both of which can be implemented over undirected and directed graphs alike. Intuitively, Eq. (2b) tracks the average of local gradients, $\frac{1}{n}\sum_{i=1}^{n}\nabla f_{i}(\mathbf{x}^{i}_{k})$, see , and therefore Eq. (2a) asymptotically approaches the centralized heavy-ball, Eq. (1), as the descent direction $\mathbf{y}^{i}_{k}$ becomes the gradient of the global objective.

Vector form: For the sake of analysis, we now write $\mathcal{AB}m$ in vector form. We use the following notation:

We note here that when $\beta_{i}=0,\forall i$, $\mathcal{AB}m$ reduces to $\mathcal{AB}$ , albeit with two distinguishing features: (i) the algorithm in uses an identical step-size, $\alpha$, at each agent; and (ii) Eq. (2b) in is in an adapt-then-combine form.

III Connection with existing first-order methods

In this section, we provide a generalization of several existing methods that employ gradient tracking and show that $\mathcal{AB}$ lies at the heart of these approaches. To proceed, we rewrite the $\mathcal{AB}$ updates below (without momentum) .

Since $\mathcal{AB}$ uses both RS and CS weights simultaneously, it is natural to ask how are the optimization algorithms that require the weight matrices to be doubly-stochastic (DS) , or only CS , or only RS , are related to each other. We discuss this relationship next.

Optimization with DS weights: Refs. consider the following updates, termed as Aug-DGM in and DIGing in :

where $\mathcal{W}=W\otimes I_{p}$, and $W$ is a DS weight matrix. Clearly, to obtain DS weights, the underlying graph must be undirected (or balanced) and thus the algorithm in Eqs. (6) is not applicable to arbitrary directed graphs. That $\mathcal{AB}$ generalizes Eqs. (6) is straightforward as the DS weights naturally satisfy the RS requirement in the top update and the CS requirement in the bottom update, while the reverse is not true. Similarly, we note that a related algorithm, EXTRA , is given by

where the two weight matrices, $\mathcal{W}$ and $\widetilde{\mathcal{W}}$, must be symmetric and satisfy some other stringent requirements, see for details. Eliminating the $\mathbf{y}_{k}$-update in $\mathcal{AB}$, we note that $\mathcal{AB}$ can be written in the EXTRA format as follows:

It can be seen that the linear convergence of $\mathcal{AB}$ does not follow from the analysis in as $\mathcal{A+B}-I$ and $\mathcal{BA}$ are not necessarily symmetric. Analysis of the $\mathcal{AB}$ algorithm, therefore, generalizes that of EXTRA to non-doubly-stochastic and non-symmetric weight matrices.

Optimization with CS weights: We now relate $\mathcal{AB}$ to ADD-OPT/Push-DIGing that only require CS weights . Since $\mathcal{B}$ is already CS in $\mathcal{AB}$, it suffices to seek a state transformation that transforms $\mathcal{A}$ from RS to CS, while respecting the graph topology. To this aim, let us consider the following transformation on the $\mathbf{x}_{k}$-update in $\mathcal{AB}$: $\widetilde{\mathbf{x}}_{k}\triangleq\Pi_{r}\mathbf{x}_{k},$ where $\Pi_{r}\triangleq\mbox{diag}(n\boldsymbol{\pi}_{r})\otimes I_{p}$ and $\boldsymbol{\pi}_{r}$ is the left-eigenvector of the RS weight matrix, $A$, corresponding to the eigenvalue $1$. The resulting transformed $\mathcal{AB}$ is given by

where it is straightforward to show that $\mathcal{\widetilde{B}}=\Pi_{r}\mathcal{A}\Pi_{r}^{-1}$ is now CS and $\widetilde{\mathcal{B}}\left(\boldsymbol{\pi}_{r}\otimes I_{p}\right)=\boldsymbol{\pi}_{r}\otimes I_{p}$.

In order to implement the above equations, two different CS matrices ($\widetilde{\mathcal{B}}$ and $\mathcal{B}$) suffice, as long as they are primitive and respect the graph topology. The second update requires the right-eigenvector of the CS matrix used in the first update, i.e., $\widetilde{\mathcal{B}}$. Since this eigenvector is not known locally to any agent, ADD-OPT/Push-DIGing propose learning this eigenvector with the following iterations: $\mathbf{w}_{k+1}=\widetilde{\mathcal{B}}\mathbf{w}_{k},\mathbf{w}_{0}=\mathbf{1}_{np}$. The algorithms provided in essentially implement Eqs. (8), albeit with two differences: (i) the same CS weight matrix is used in all updates; and, (ii) the division in Eq. (8b) is replaced by the estimated component, $\mathbf{w}_{k+1}^{i}$, of the left-eigenvector at each agent. This nonlinearity causes stability issues in ADD-OPT/Push-DIGing, whereas their convergence compared to $\mathcal{AB}$ is slower because such an eigenvector estimation is not needed in the latter on the account of using the RS weights. Furthermore, the local step-sizes are now given by $n\alpha[\boldsymbol{\pi}_{r}]_{i}$ that shows that ADD-OPT/Push-DIGing should work with nonidentical step-sizes.

Optimization with RS weights: The state transformation technique discussed above also leads to an algorithm from $\mathcal{AB}$ that only requires RS weights. Since $\mathcal{A}$ in $\mathcal{AB}$ is RS, a transformation now is imposed on the $\mathbf{y}_{k}$-update and is given by $\widetilde{\mathbf{y}}_{k}\triangleq\Pi_{c}^{-1}\mathbf{y}_{k}$, where $\Pi_{c}\triangleq\mbox{diag}(\boldsymbol{\pi}_{c})\otimes I_{p},$ and $\boldsymbol{\pi}_{c}$ is the right-eigenvector of the CS weight matrix, $B$, corresponding to the eigenvalue $1$. Equivalently, $\mathcal{AB}$ is given by

where $\widetilde{\mathcal{A}}=\Pi_{c}^{-1}\mathcal{B}\Pi_{c}$ is now RS and $\left(\boldsymbol{\pi}_{c}^{\top}\otimes I_{p}\right)\widetilde{\mathcal{A}}=\boldsymbol{\pi}_{c}^{\top}\otimes I_{p}$. Since the above form of $\mathcal{AB}$ cannot be implemented because $\boldsymbol{\pi}_{c}$ is not locally known, an eigenvector estimation is used in FROST and the division in Eq. (9b) is replaced with the appropriate estimated component of $\boldsymbol{\pi}_{c}$. The observations on different weight matrices in the two updates, nonidentical step-sizes, stability, and convergence made earlier for ADD-OPT/Push-DIGing are also applicable here.

In conclusion, the $\mathcal{AB}$ algorithm has various equivalent representations and several already-known protocols can in fact be derived from these representations. In a similar way, $\mathcal{AB}m$ leads to protocols that add momentum to Aug-DGM, ADD-OPT/Push-DIGing, and FROST. We will revisit the relationship and equivalence cast here in Sections V and VI. In Section V, we will show that both $\mathcal{AB}$ and $\mathcal{AB}m$ naturally provide a non-trivial class of average-consensus algorithms, a special case of which are and surplus consensus . In Section VI, we will compare these algorithms numerically.

IV Convergence Analysis

We now start the convergence analysis of the proposed distributed heavy-ball method, $\mathcal{AB}m$. In the following, we first provide some auxiliary results borrowed from the literature.

The following lemma establishes contractions with RS and CS matrices under arbitrary norms ; note thacontraction in the Euclidean norm is not applicable unless the weight matrix is DS as in . A similar result was first presented in for CS matrices, and later in for RS matrices.

where $0<\sigma_{\mathcal{A}}<1$ and $0<\sigma_{\mathcal{B}}<1$ are some constants.

The next lemma from states that the sum of $\mathbf{y}_{k}^{i}$’s preserves the sum of local gradients. This is a direct consequence of the dynamic consensus employed with CS weights in the $\mathbf{y}_{k}$-update of $\mathcal{AB}m$.

$(\mathbf{1}_{n}^{\top}\otimes I_{p})\mathbf{y}_{k}=(\mathbf{1}_{n}^{\top}\otimes I_{p})\nabla\mathbf{f}(\mathbf{x}_{k}),\forall k$.

The next lemma is standard in the convex optimization theory . It states that the distance to the optimizer contracts at each step in the standard gradient descent method.

Let $F$ be $\mu$-strongly-convex and $l$-smooth. For $0<\alpha<\frac{2}{l}$, we have

where $\sigma_{F}=\max\left(\left|1-\mu\alpha\right|,\left|1-l\alpha\right|\right)$.

Finally, we provide a result from nonnegative matrix theory.

IV-B Main results

The convergence analysis of $\mathcal{AB}m$ is based on deriving a contraction relationship between the following four quantities: (i) $\|\mathbf{x}_{k+1}-\mathcal{A}_{\infty}\mathbf{x}_{k+1}\|_{\mathcal{A}}$, the consensus error in the network; (ii) $\|\mathcal{A}_{\infty}\mathbf{x}_{k+1}-\mathbf{1}_{n}\otimes\mathbf{x}^{*}\|$, the optimality gap; (iii) $\left\|\mathbf{x}_{k+1}-\mathbf{x}_{k}\right\|$, the state difference; and (iv) $\|\mathbf{y}_{k+1}-\mathcal{B}_{\infty}\mathbf{y}_{k+1}\|_{\mathcal{B}}$, the (biased) gradient estimation error. We will establish an LTI-system inequality where the state vector is the collection of these four quantities and then develop the convergence properties of the corresponding system matrix. Before we proceed, note that since all vector norms on finite-dimensional vector spaces are equivalent , there exist positive constants $c_{\mathcal{A}\mathcal{B}},c_{\mathcal{B}\mathcal{A}},c_{2\mathcal{A}},c_{\mathcal{A}2},c_{2\mathcal{B}},c_{\mathcal{B}2}$ such that

We also define $\overline{\alpha}\triangleq\left[\boldsymbol{\alpha}\right]_{\max}$ and $\overline{\beta}\triangleq\left[\boldsymbol{\beta}\right]_{\max}$. In the following, we first provide an upper bound on the estimate, $\mathbf{y}_{k}$, of the gradient of the global objective that will be useful in deriving the aforementioned LTI system.

The following inequality holds, $\forall k$:

Recall that $\mathcal{B}_{\infty}=(\boldsymbol{\pi}_{c}\otimes I_{p})(\mathbf{1}_{n}^{\top}\otimes I_{p}).$ We have

We next bound $\left\|\mathcal{B}_{\infty}\mathbf{y}_{k}\right\|$:

where the first inequality uses Jensen’s inequality and the last inequality uses the fact that $\left\|\mathcal{B}_{\infty}\right\|=\sqrt{n}\|\boldsymbol{\pi}_{c}\|$. The lemma follows by plugging Eq. (IV-B) into Eq. (12). ∎

In the next Lemmas 6-9, we derive the relationships among the four quantities mentioned above. We start with a bound on $\|\mathbf{x}_{k+1}-\mathcal{A}_{\infty}\mathbf{x}_{k+1}\|_{\mathcal{A}}$, the consensus error in the network.

The following inequality holds, $\forall k$:

First, note that $\mathcal{A}_{\infty}\mathcal{A}=\mathcal{A}_{\infty}$. Following the $\mathbf{x}_{k}$-update of $\mathcal{AB}m$ in Eq. (4a) and using the one-step contraction property of $\mathcal{A}$ from Lemma 1, we have:

Next, we derive a bound for $\left\|\mathcal{A}_{\infty}\mathbf{x}_{k+1}-\mathbf{1}_{n}\otimes\mathbf{x}^{*}\right\|$, which can be interpreted as the optimality gap between the network accumulation state, $\mathcal{A}_{\infty}\mathbf{x}_{k}$, and the global minimizer, $\mathbf{1}_{n}\otimes\mathbf{x}^{*}$.

The following inequality holds, $\forall k$, when

$0<\boldsymbol{\pi}_{r}^{\top}\mbox{diag}(\boldsymbol{\alpha})\boldsymbol{\pi}_{c}<\frac{2}{nl}$:

where $\lambda=\max\left\{\left|1-\mu n\boldsymbol{\pi}_{r}^{\top}\mbox{diag}(\boldsymbol{\alpha})\boldsymbol{\pi}_{c}\right|,\left|1-ln\boldsymbol{\pi}_{r}^{\top}\mbox{diag}(\boldsymbol{\alpha})\boldsymbol{\pi}_{c}\right|\right\}.$

Recall the $\mathbf{x}_{k}$-update of $\mathcal{AB}m$ in Eq. (4a), we have that

where in the last inequality, we use $\mathcal{B}_{\infty}\mathbf{y}_{k}=\mathcal{B}_{\infty}\nabla\mathbf{f}\left(\mathbf{x}_{k}\right)$ adapted from Lemma 2. Since the last two terms in Eq. (IV-B) match the last two terms in Eq. (7), what is left is to bound the first term. Before we proceed, define

Now we bound the first term in Eq. (IV-B). We have

and we bound $s_{1}$ and $s_{2}$ next. Using Lemma 3, we have that if $0<\boldsymbol{\pi}_{r}^{\top}\mbox{diag}(\boldsymbol{\alpha})\boldsymbol{\pi}_{c}<\frac{2}{nl}$,

where $\lambda=\max\left\{\left|1-\mu n\boldsymbol{\pi}_{r}^{\top}\mbox{diag}(\boldsymbol{\alpha})\boldsymbol{\pi}_{c}\right|,\left|1-ln\boldsymbol{\pi}_{r}^{\top}\mbox{diag}(\boldsymbol{\alpha})\boldsymbol{\pi}_{c}\right|\right\}.$ We next bound $s_{2}$. Since $\nabla F(\widetilde{\mathbf{x}}_{k})=\frac{1}{n}(\mathbf{1}_{n}^{\top}\otimes I_{p})\nabla\mathbf{f}(\widetilde{\mathbf{x}}_{k})$,

and the lemma follows from Eqs. (IV-B), (IV-B), and (IV-B). ∎

The next step is to bound the state difference, $\left\|\mathbf{x}_{k+1}-\mathbf{x}_{k}\right\|$.

The following inequality holds, $\forall k$:

Note that $\mathcal{A}\mathcal{A}_{\infty}=\mathcal{A}_{\infty}$ and hence $\mathcal{A}\mathcal{A}_{\infty}-\mathcal{A}_{\infty}$ is a zero matrix. Following the $\mathbf{x}_{k}$-update of $\mathcal{AB}m$, we have:

The final step in formulating the LTI system is to write $\left\|\mathbf{y}_{k+1}-\mathcal{B}_{\infty}\mathbf{y}_{k+1}\right\|$, the biased gradient estimation error, in terms of the other three quantities. We call this biased to make a distinction with the unbiased gradient estimation error: $\left\|\mathbf{y}_{k+1}-\mathcal{W}_{\infty}\mathbf{y}_{k+1}\right\|$, where $\mathcal{W}$ is doubly-stochastic.

The following inequality holds, $\forall k$:

Note that $\mathcal{B}_{\infty}\mathcal{B}=\mathcal{B}_{\infty}$. From Eq. (4b), we have:

where in the inequality above we use the contraction property of $\mathcal{B}$ from Lemma 1. The proof follows by applying the result of Lemma 8 to the inequality above. ∎

With the help of the Lemmas 6-9, we now present the main result of this paper, i.e., the $\mathcal{AB}m$ algorithm converges to the global minimizer at a global $R$-linear rate.

Let $0<\boldsymbol{\pi}_{r}^{\top}\mbox{diag}(\boldsymbol{\alpha})\boldsymbol{\pi}_{c}<\frac{2}{nl}$, then the following LTI inequality holds entry-wise:

and the constants $a_{i}$’s in the above expression are

When the largest step-size, $\overline{\alpha}$, satisfies

and when the largest momentum parameter, $\overline{\beta}$, satisfies

where $\delta_{1},\delta_{2},\delta_{3},\delta_{4}$ are arbitrary constants such that

then $\rho(J_{\boldsymbol{\alpha},\overline{\beta}})<1$ and thus $\|\mathbf{x}_{k}-\mathbf{1}_{n}\otimes\mathbf{x}^{*}\|$ converges to zero linearly at the rate of $\mathcal{O}(\rho(J_{\boldsymbol{\alpha},\overline{\beta}}))^{k}$.

It is straightforward to verify Eq. (18) by combining Lemmas 6-9. The next step is to find the range of $\overline{\alpha}$ and $\overline{\beta}$ such that $\rho(J_{\boldsymbol{\alpha},\overline{\beta}})<1$. In the light of Lemma 4, we solve for a positive vector $\boldsymbol{\delta}=[\delta_{1},\delta_{2},\delta_{3},\delta_{4}]^{\top}$ and the range of $\overline{\alpha}$ and $\overline{\beta}$ such that the following inequality holds:

which is equivalent to the following four conditions:

Recall $\lambda$ in Lemma 7, when $\overline{\alpha}<\frac{1}{nl\boldsymbol{\pi}_{r}^{\top}\boldsymbol{\pi}_{c}}$, we have

Therefore, the third condition in Eq. (37) is satisfied when

For the right hand side of the Eq. (36), (40), (38) and (39) to be positive, each one of these equations needs to satisfy the conditions we give below.

We first choose arbitrary positive constants, $\delta_{3}$ and $\delta_{4}$, then pick $\delta_{1}$ satisfying Eqs. (45) and (48), and finally choose $\delta_{2}$ according to Eq. (42). Note that $\delta_{1},\delta_{2},\delta_{3},$ and $\delta_{4}$ are chosen to ensure that the upper bounds on $\overline{\alpha}$ are all positive. Subsequently, from Eqs. (41), (45), and (48), together with the requirement that $\overline{\alpha}<\frac{1}{nl\boldsymbol{\pi}_{r}^{\top}\boldsymbol{\pi}_{c}}$, we obtain the upper bound on the largest step-size, $\overline{\alpha}$. Finally, the original four conditions in Eqs. (36), (40), (38) and (39) lead to an upper bound on $\overline{\beta}$, and the theorem follows. ∎

Remark 1: In Theorem 1, we have established the $R$-linear rate of $\mathcal{AB}m$ when the largest step-size, $\overline{\alpha}$, and the largest momentum parameter, $\overline{\beta}$, respectively follow the upper bounds described in Eq. (29) and Eq. (30). Note that $\delta_{1},\delta_{2},\delta_{3},\delta_{4}$ therein are tunable parameters and only depend on the network topology and the objective functions. The upper bounds for $\overline{\alpha}$ and $\overline{\beta}$ may not be computable for arbitrary directed graphs as the contraction coefficients, $\sigma_{\mathcal{A}}$, $\sigma_{\mathcal{B}}$, and the norm equivalence constants may be unknown. However, when the graph is undirected, we can obtain computable bounds for $\overline{\alpha}$ and $\overline{\beta}$, as developed in for example. The upper bound on $\overline{\beta}$ also implies that if the step-sizes are relatively large, only small momentum parameters can be picked to ensure stability.

Remark 2: The nonidentical step-sizes in gradient tracking methods have previously been studied in . These works rely on some notion of heterogeneity among the step-sizes, defined respectively as the relative deviation of the step-sizes from their average, $\frac{\|(I-W)\boldsymbol{\alpha}\|}{\|W\boldsymbol{\alpha}\|}$, in , and as the ratio of the largest to the smallest step-size, ${[\boldsymbol{\alpha}]_{\max}}/{[\boldsymbol{\alpha}]_{\min}}$, in . The authors then show that when the heterogeneity is sufficiently small and when the largest step-size follows a bound that is a function of the heterogeneity, the proposed algorithms converge to the global minimizer. It is worth noting that sufficiently small step-sizes do not guarantee sufficiently small heterogeneity in both of the above definitions. In contrast, the upper bound on the largest step-size in this paper, Eq. (29), is independent of any notion of heterogeneity and only depends on the objective functions and the network topology. Each agent therefore locally picks a sufficiently small step-size without any coordination. Based on the discussion in Section III, our approach thus improves the analysis in . Besides, Eq. (29) allows the existence of zero step-sizes among the agents as long as the largest step-size is positive and is sufficiently small.

Remark 3: To show that $\mathcal{AB}m$ has an $R$-linear rate for sufficiently small $\overline{\alpha}$ and $\overline{\beta}$, one can alternatively use matrix perturbation analysis as in (Theorem 1). However, it does not provide explicit upper bounds on $\overline{\alpha}$ and $\overline{\beta}$ in closed form.

V Average-Consensus from 𝒜​ℬ​m𝒜ℬ𝑚\mathcal{AB}m

In this section, we show that $\mathcal{AB}m$ subsumes a novel average-consensus algorithm over strongly-connected directed graphs. To show this, we choose the objective functions as

Clearly, the minimization of $\widetilde{F}=\sum_{i=1}^{n}\widetilde{f}_{i}$ is now achieved at $\mathbf{x}^{*}=\tfrac{1}{n}\sum_{i=1}^{n}\boldsymbol{\upsilon}_{i}$. The $\mathcal{AB}m$ algorithm, Eq. (4), thus naturally leads to the following average-consensus algorithm, termed as $\mathcal{AB}m$-$\mathcal{C}$, with $\nabla\mathbf{f}(\mathbf{x}_{k+1})-\nabla\mathbf{f}(\mathbf{x}_{k})=\mathbf{x}_{k+1}-\mathbf{x}_{k}$; for the sake of simplicity, we choose $\alpha_{i}=\alpha,\beta_{i}=\beta,\forall i$:

Its local implementation at each agent $i$ is given by:

where $\mathbf{x}^{i}_{0}=\boldsymbol{\upsilon}_{i}$ and $\mathbf{y}_{i}^{0}=0,~{}\forall i$.

From the analysis of $\mathcal{AB}m$, an $R$-linear convergence of $\mathcal{AB}m$-$\mathcal{C}$ to the average of $\boldsymbol{\upsilon}_{i}$’s is clear from Theorem 1. It may be possible to make concrete rate statements by studying the spectral radius of the following system matrix:

However, this analysis is beyond the scope of this paper. We note that when $\beta=0$, the above equations still converge to the average of $\boldsymbol{\upsilon}_{i}$’s according to Theorem 1. What is surprising is that, with $\beta=0$, $\mathcal{AB}m$-$\mathcal{C}$ reduces to

which is surplus consensus , after a state transformation with $\mbox{diag}\left(I,-I\right)$; in fact, any state transformation of the form $\mbox{diag}(I,\widetilde{I})$ applies here as long as $\widetilde{I}$ is diagonal (to respect the graph topology) and invertible. More importantly, compared with surplus consensus , $\mathcal{AB}m$-$\mathcal{C}$ uses information from the past iterations. This history information is in fact the momentum from a distributed optimization perspective, which may lead to accelerated convergence as we will numerically show in Section VI.

Following this discussion, choosing the local functions as $\widetilde{f}_{i}$’s in , or in ADD-OPT , or in FROST , we get average-consensus with only DS, CS, or RS weights. The protocol that results directly from $\mathcal{AB}$ is surplus consensus, while the one resulting directly from FROST was presented in . With the analysis provided in Section III, we see that the algorithm in is in fact related to surplus consensus after a state transformation. Clearly, accelerated average-consensus based exclusively on either row- or column-stochastic weights can be abstracted from the discussion herein, after adding a momentum term.

VI Numerical Experiments

We now provide numerical experiments to illustrate the theoretical findings described in this paper. To this aim, we use two different graphs: an undirected graph, $\mathcal{G}_{1}$, and a directed graph, $\mathcal{G}_{2}$. Both graphs have $n=500$ agents and are generated using nearest neighbor rules and then we add less than $0.05\%$ random links. The number of edges in all cases is less $4\%$ of the total possible edges. Since the graphs are randomly generated across experiments, two sample graphs are shown in Fig. 1, without the self-edges and random links for visual clarity. We generate DS weights using the Laplacian method: $W=I-\tfrac{1}{\max_{i}\deg_{i}+1}L$, where $L$ is the graph Laplacian and $\deg_{i}$ is the degree of node $i$. Additionally, we generate RS and CS weights with the uniform weighting strategy: $a_{ij}=\tfrac{1}{|\mathcal{N}_{j}^{{\scriptsize\mbox{in}}}|}$ and $b_{ij}=\tfrac{1}{|\mathcal{N}_{j}^{{\scriptsize\mbox{out}}}|},\forall i,j$. We note that both weighting strategies are applicable to undirected graphs, while only the uniform strategy can be used over directed graphs.

where $\frac{\lambda}{2}\|\mathbf{b}\|_{2}^{2}$ is a regularization term used to prevent over-fitting of the data. The feature vectors, $\mathbf{c}_{ij}$’s, are randomly generated from a Gaussian distribution with zero mean and the binary labels are randomly generated from a Bernoulli distribution. We plot the average of residuals at each agent, $\frac{1}{n}\sum_{i=1}^{n}\|\mathbf{x}_{i}(k)-\mathbf{x}^{*}\|_{2}$, and first compare the performance of the following over undirected graphs in Fig. 2 (Left): (i) $\mathcal{AB}m$with RS and CS weights; (ii) $\mathcal{AB}m$with DS weights; (iii) distributed optimization based on gradient tracking from , with DS weights; (iv) EXTRA from ; and, (v) centralized gradient descent.

Next, we compare the performance similarly over directed graphs in Fig. 2 (Right). Here, the algorithms with doubly-stochastic weights are not applicable, and instead we compare $\mathcal{AB}m$ with $\mathcal{AB}$ , ADD-OPT/Push-DIGing , and centralized gradient descent. The weight matrices are chosen as we discussed before and the algorithm parameters are hand-tuned for best performance (except for gradient descent where the optimal step-size is given by $\alpha=\tfrac{2}{\mu+l}$). We note that momentum improves the convergence when compared to applicable algorithms without momentum, while ADD-OPT/Push-DIGing are much slower because of the eigenvector estimation, see Section III for details.

VI-B Distributed Quadratic Programming

For small condition numbers, we note that gradient descent is quite fast and the distributed algorithms suffer from a relatively slower fusion over the graphs. Recall that the optimal convergence rate of gradient decent is $\mathcal{O}((\tfrac{\mathcal{Q}-1}{\mathcal{Q}+1})^{k})$. When the condition number is large, gradient descent is quite conservative allowing fusion to catch up. Finally, we note that $\mathcal{AB}m$, with momentum, outperforms the centralized gradient descent when the condition number is large. This observation is consistent with the existing literature, see e.g., .

VI-C 𝒜​ℬ​m𝒜ℬ𝑚\mathcal{AB}m and Average-Consensus

We now provide numerical analysis and simulations to show that $\mathcal{AB}m$-$\mathcal{C}$, in Eq. (49), possibly achieves acceleration when compared with surplus-consensus, in Eq. (56). To explain our choice of $\alpha$ and $\beta$, we first note that the power limit of the system matrix in Eq. (56), denoted as $\mathcal{H}$, is :

where $\mathcal{W}_{\infty}=(\frac{1}{n}\mathbf{1}_{n}\mathbf{1}_{n}^{\top})\otimes I_{p}$. It is straightforward to show that $\mathcal{H}^{k}-\mathcal{H_{\infty}}=\left(\mathcal{H}-\mathcal{H}_{\infty}\right)^{k}.$ Similarly, for the augmented system matrix, $\widetilde{\mathcal{H}}$, in Eq. (49), we observe that

and it can be verified that $\mathcal{\widetilde{H}}^{k}-\mathcal{\widetilde{H}_{\infty}}=(\mathcal{\widetilde{H}}-\mathcal{\widetilde{H}}_{\infty})^{k}.$ We therefore use grid search to choose the optimal $\alpha^{*}$ in $\mathcal{H}$ and the optimal $\widetilde{\alpha}^{*}$ and $\widetilde{\beta}^{*}$ in $\mathcal{\widetilde{H}}$, which respectively minimize $\rho(\mathcal{H}-\mathcal{H}_{\infty})$ and $\rho(\mathcal{\widetilde{H}}-\mathcal{\widetilde{H}}_{\infty})$. Numerically, we observe that it may be possible for the minimum of $\rho(\mathcal{\widetilde{H}}-\mathcal{\widetilde{H}}_{\infty})$ to be smaller than that of $\rho\left(\mathcal{H}-\mathcal{H}_{\infty}\right)$. The convergence speed comparison between $\mathcal{AB}m$-$\mathcal{C}$ and surplus consensus is shown in Fig 5 over a directed graph, $\mathcal{G}_{2}$.

VII Conclusions

In this paper, we provide a framework for distributed optimization that removes the need for doubly-stochastic weights and thus is naturally applicable to both undirected and directed graphs. Using a state transformation based on the non-$\mathbf{1}_{n}$ eigenvector, we show that the underlying algorithm, $\mathcal{AB}$, based on a simultaneous application of both RS and CS weights, lies at the heart of several algorithms studied earlier that rely on eigenvector estimation when using only CS (or only RS) weights. We then propose the distributed heavy-ball method, termed as $\mathcal{AB}m$, that combines $\mathcal{AB}$ with a heavy-ball (type) momentum term. To the best of our knowledge, this paper is the first to use a momentum term based on the heavy-ball method in distributed optimization. We show that $\mathcal{AB}m$ subsumes a novel average-consensus algorithm as a special case that unifies earlier attempts over directed graphs, with potential acceleration due to the momentum term.

References