---
abstract: 'This work concerns the analysis and design of distributed first-order optimization algorithms over time-varying graphs. The goal of such algorithms is to optimize a global function that is the average of local functions using only local computations and communications. Several different algorithms have been proposed that achieve linear convergence to the global optimum when the local functions are strongly convex. We provide a unified analysis that yields the worst-case linear convergence rate as a function of the condition number of the local functions, the spectral gap of the graph, and the parameters of the algorithm. The framework requires solving a small semidefinite program whose size is fixed; it does not depend on the number of local functions or the dimension of their domain. The result is a computationally efficient method for distributed algorithm analysis that enables the rapid comparison, selection, and tuning of algorithms. Finally, we propose a new algorithm, which we call SVL, that is easily implementable and achieves a faster worst-case convergence rate than all other known algorithms.'
author:
- 'Akhil Sundararajan$^{1,2}$'
- 'Bryan Van Scoy$^{1}$'
- 'Laurent Lessard$^{1,2}$'
bibliography:
- 'distralg\_tcns.bib'
title: 'Analysis and Design of First-Order Distributed Optimization Algorithms over Time-Varying Graphs'
---

Introduction {#sec:intro}
============

In distributed optimization, a network of agents, such as computing nodes, robots, or mobile sensors, work collaboratively to optimize a global objective. Specifically, each agent $i\in\{1,\dots,n\}$ has access to a local function $f_i$ and must minimize the average of all agents’ local functions $$\label{eq:distrop-problem}
\min_{x\in {\mathbb{R}}^d} f(x),
\quad\text{where }f(x) {\colonequals}\frac{1}{n}\sum_{i=1}^n f_i(x),$$

by querying its local gradient ${\nabla\! f}_i$, exchanging information with neighboring agents, and performing local computations.

This work aims to study the reliability of distributed optimization algorithms in the presence of a *time-varying* communication graph. Such a scenario could occur if communication links fail due to interference, mobile agents move out of range, or an adversary is jamming communications.

Distributed optimization is relevant in many application areas. For example, in large-scale machine learning [@forero2010consensus; @johansson2008distributed], $n$ could represent the number of computing units available for training a large data set. Each $f_i$ then denotes the loss function corresponding to the training examples assigned to unit $i$. Another example is sensor networks [@rabbat04], where each sensor may have a limited power budget, communication bandwidth, or sensing capability. [The goal is to aggregate all]{} local data without having a single point of failure. Other applications include distributed spectrum sensing [@dist_spectrum_sensing] and resource allocation across geographic regions [@dist_power_control].

Distributed optimization generalizes both average consensus and centralized optimization, as we now explain.

#### Consensus

If each agent uses the initial value $x_i^0$ and local objective $f_i(x) = {\lVert{x-x_i^0}\rVert}^2$, distributed optimization reduces to *average consensus* [@tsitsiklis; @xiaoboyd04]. The unique optimizer of  is then the average of all initial states: $x^\star=\frac{1}{n}\sum_{i=1}^n x_i^0$. Using a *gossip* update of the form $x_i^{k+1} = \sum_{i=1}^n {W}_{ij} x_j^k$ where ${W}$ is carefully chosen[,]{} such methods converge exponentially: ${\lVert{x_i^k -x^\star}\rVert} \le \rho^k$ with $\rho \in (0,1)$ that depends on ${W}$ [@FDLA]. This is called a *linear rate* in the optimization community.

#### Optimization

If $n=1$ or if all $f_i$ are identical, we recover the standard centralized optimization setup. Linear convergence can be guaranteed in certain cases. For example, the gradient descent method $x_i^{k+1} = x_i^k - \alpha {\nabla\!}f_i(x_i^k)$ achieves linear convergence if $f_i$ is continuously differentiable, smooth, and strongly convex (formally stated in Assumption \[assumption:local\_functions\]) [@NesterovBook].

A linear convergence rate for the general case was first achieved by the *exact first-order algorithm* (EXTRA) [@EXTRA]. This algorithm requires storing the previous state in memory:

\[eq:EXTRA\] $$\begin{aligned}
    x_i^1 &= \sum_{j=1}^n {W}_{ij}\,x_j^0 - \alpha\, {\nabla\! f}_i(x_i^0),\qquad x_i^0\text{ arbitrary}, \\
    \label{eq:EXTRAb}
    x_i^{k+2} &= x_i^{k+1}+\sum_{j=1}^n {W}_{ij}\,x_j^{k+1} - \sum_{j=1}^n \widetilde{{W}}_{ij}\,x_j^k -\alpha\,\bigl({\nabla\! f}_i(x_i^{k+1})-{\nabla\! f}_i(x_i^k)\bigr)
    \end{aligned}$$

where ${W}$ and $\widetilde{{W}}$ are gossip matrices that satisfy certain technical conditions and $\alpha$ is sufficiently small. [Several additional linear-rate algorithms have since been proposed, including: AugDGM [@AugDGM], DIGing [@DIGing; @QuLi], Exact Diffusion [@ExactDiffusion1; @ExactDiffusion2], NIDS [@NIDS], and a unified method [@Unification]. Each of these methods have updates similar to  in that they require agents to store previous iterates or gradients.]{}

Although linear convergence rates were obtained for the algorithms above, each algorithm differs in the nature and strength of its convergence analysis guarantees. For example, some works show (non-constructively) the existence of a linear rate [@AB] whereas others provide specific tuning recommendations with associated analytic rate bounds (which may be conservative) [@EXTRA; @NIDS]. Numerical simulations are also frequently used [@ABN], but can be misleading because algorithm performance depends on the graph topology, choice of functions, algorithm initialization, and algorithm tuning.

The present work makes an effort to systematize the analysis and design of distributed optimization algorithms. We now summarize our main contributions.

**Analysis framework.** We present a universal analysis framework that [provides an upper bound on the]{} worst-case linear convergence rate $\rho$ of a wide range of distributed algorithms as a function of the parameters $\kappa$ (local function conditioning) and $\sigma$ (network connectedness). Our main result, Theorem \[thm:main-result\], is a semidefinite program (SDP) parameterized by [ $(\kappa,\sigma)$]{} whose solution yields an upper bound on $\rho$. The SDP has a small fixed size that does not depend on the number of agents $n$ or the dimension of the function domains [and is efficiently solvable]{}. Our SDP yields robust performance guarantees when the graph is allowed to vary [(even adversarially)]{} at each iteration. [Fig. \[fig:money\_plot\] compares the worst-case linear rate $\rho$ for 8 different algorithms.]{}

**Algorithm design.** [We present a]{} new distributed algorithm, which we name SVL (the authors’ initials). SVL is derived by optimizing the SDP from our analysis framework and provides the fastest known convergence rate to date [for this time-varying graph setting]{}. [The rate]{} depends explicitly on $\kappa$ and $\sigma$, so no tuning is required if these parameters are known or estimated in advance. [When the graph is well-connected, SVL recovers the performance of gradient descent, which is optimal in this time-varying graph setting.]{}

**Worst-case examples.** Although our analysis technique only provides *upper bounds* on the worst-case convergence rate [for]{} distributed algorithms, we outline a computationally tractable optimization procedure that finds numerically matching lower bounds [by constructing worst-case trajectories, suggesting the bounds found via our analysis technique are tight.]{}

[Distributed algorithms that achieve *accelerated*[@QuLi_accelerated; @ABN; @DHB] or *optimal* [@Scaman17] linear rates have also been proposed. It turns out such methods are not guaranteed to achieve acceleration when the graph is time-varying. We discuss this phenomenon in Section \[sec:lower\_bounds\], where we derive lower bounds for the time-varying setting.]{}

The paper is organized as follows. We [describe notation and assumptions]{} in Section \[prelim\]. We state and prove our main result for certifying worst-case rate bounds in Section \[main\_result\]. We present our SVL algorithm and discuss interpretations in Section \[sec:algo\_design\]. Finally, we demonstrate the tightness of our bounds by generating worst-case trajectories in Section \[numerics\].

Preliminaries {#prelim}
=============

Notation
--------

Let $I_n$ be the identity matrix in ${\mathbb{R}}^{n\times n}$. The symbol $1_n$ denotes the column vector of all ones in ${\mathbb{R}}^n$. $\Pi {\colonequals}\frac{1}{n}1_n 1_n^{\mathsf{T}}$ is the projection matrix onto $1_n$. We will sometimes omit subscripts when dimensions are clear from context. Unless otherwise indicated, Greek letters denote scalar parameters, lower-case letters denote column vectors, and upper-case letters denote matrices. Exceptions include the scalars $m$ and $L$, which we use in Assumption \[assumption:local\_functions\] to conform with convention. The symbol $\otimes$ denotes the Kronecker matrix product. ${\lVert{x}\rVert}$ denotes the standard Euclidean norm of a vector $x$, and ${\lVert{A}\rVert}{\colonequals}\sup_{x\ne 0} {\lVert{Ax}\rVert}/{\lVert{x}\rVert}$ is the spectral norm of a matrix $A$. Unless otherwise indicated, subscripts refer to individual agents while superscripts refer to iteration count. For brevity, we write the symmetric quadratic form $x^{\mathsf{T}}Q x$ as ${\begin{bmatrix}\star\end{bmatrix}}^{\mathsf{T}}Q x$.

Define the graph $\mathcal{G}{\colonequals}({\mathcal{V}},{\mathcal{E}})$ where ${\mathcal{V}}{\colonequals}\{1,\dots,n\}$ is the set of agents and ${\mathcal{E}}$ is the set of pairs of agents $(i,j)$ that are connected. ${\mathcal{L}}\in {\mathbb{R}}^{n\times n}$ is a *Laplacian matrix* associated with $\mathcal{G}$ if ${\mathcal{L}}1_n = 0$ and ${\mathcal{L}}_{ij} = 0$ if $(i,j)\notin {\mathcal{E}}$. The *spectral gap* of ${\mathcal{L}}$ is defined as the second-smallest eigenvalue magnitude of ${\mathcal{L}}$. Since we consider time-varying graphs, we let ${\mathcal{L}}^k$ denote a Laplacian [matrix associated with]{} $\mathcal{G}^k$. [We denote a symbol on agent $i$ at iteration $k$ by $x_i^k$ along with its associated fixed point $x_i^\star$. For all such symbols, we denote their aggregation over all agents as $$x^k {\colonequals}{\begin{bmatrix}x_1^k \\ \vdots \\ x_n^k\end{bmatrix}} \quad\text{and}\quad
    x^\star = {\begin{bmatrix}x_1^\star \\ \vdots \\ x_n^\star\end{bmatrix}}.$$ We denote the associated local and global error coordinates as $\tilde x_i^k {\colonequals}x_i^k - x_i^\star$ and $\tilde x^k {\colonequals}x^k - x^\star$, respectively. ]{}

Function and Graph Assumptions
------------------------------

[We assume that the local function gradients satisfy the following *sector bound*.]{}

\[assumption:local\_functions\] Given $0<m\le L$, the the local objective functions $f_i$ are continuously differentiable and each satisfy [ $$\bigl( {\nabla\! f}_i(y)-{\nabla\! f}_i(y_\text{opt}) - m\,(y-y_\text{opt})\bigr)^{\mathsf{T}}\bigl( {\nabla\! f}_i(y)-{\nabla\! f}_i(y_\text{opt}) - L\,(y-y_\text{opt})\bigr) \le 0$$ for all $y\in{\mathbb{R}}^d$, where $y_\text{opt}$ satisfies $\sum_{i=1}^n {\nabla\! f}_i(y_\text{opt}) = 0$.]{}

One way to satisfy Assumption \[assumption:local\_functions\] is if the local functions $f_i$ are $L$-Lipschitz continuous and $m$-strongly convex, though in general, Assumption \[assumption:local\_functions\] is much weaker.

We define the condition ratio as $\kappa{\colonequals}L/m$. This quantity captures how much the curvature of the objective function varies. If $f$ is twice differentiable, $\kappa$ is an upper bound on the condition number of the Hessian $\nabla^2 f$. [In general, as $\kappa\to\infty$, the functions become poorly conditioned and more difficult to optimize using first-order methods.]{}

The graph associated with the network of agents can change at each step of the algorithm, so we assume the following about the sequence of graph Laplacian matrices $\{{\mathcal{L}}^k\}$.

\[assumption:Laplacian\] The following properties hold at each step of the algorithm.

1.  The graph is connected: there always exists a path between any two nodes in $\mathcal{G}^k$. This implies that the zero eigenvalue of ${\mathcal{L}}^k$ has a multiplicity of one for all $k$.\[it:a1\]

2.  The graph is balanced: every node has equal in-degree and out-degree. This means that $1_n^{\mathsf{T}}{\mathcal{L}}^k = 0$ for all $k$.\[it:a2\]

3.  The spectral gap of the time-varying graph is uniformly bounded. In particular, we assume there exists $\sigma \in [0,1)$ such that ${\lVert{ I-\Pi-{\mathcal{L}}^k }\rVert} \leq \sigma$ for all $k$. Since the spectral radius of a matrix is always upper-bounded by its spectral norm, this implies that $\sigma$ is a uniform bound on the spectral gap of each Laplacian matrix in $\{{\mathcal{L}}^k\}$.\[it:a3\]

\[rem:B-connected\] The assumption that $\mathcal{G}^k$ must be connected for all $k$ is a strong assumption. Works that consider directed or time varying graphs typically make weaker assumptions, such as a *joint spectrum property* or *$B$-connectedness* [@DIGing]. Nevertheless, our setting (which is equivalent to $B$-connectedness with $B=1$) is still weaker than assuming a constant graph. Indeed, NIDS [@NIDS] converges for any $\sigma$ when the graph is constant, but in Section \[sec:worst-case\], we construct a sequence of graphs that drives NIDS to instability.

Algorithm Form
--------------

In this paper, we consider the broad class of distributed optimization algorithms that satisfy the algebraic equations

\[alg\] $$\begin{gathered}
    {\begin{bmatrix} x^{k+1}_i \\ y^k_i \\ z^k_i \end{bmatrix}}
    = {\begin{bmatrix} A & B_u & B_v \\ C_y & D_{yu} & D_{yv} \\ C_z & D_{zu} & D_{zv} \end{bmatrix}}
    {\begin{bmatrix} x^k_i \\ u^k_i \\ v^k_i\end{bmatrix}}, \label{alg1} \\
    u^k_i = {\nabla\! f}_i (y^k_i), \qquad
    v^k_i = \sum_{j=1}^n{\mathcal{L}}^k_{ij} z^k_j, \label{alg2}\\
    \sum_{j=1}^n \left( F_x x^k_j + F_u u^k_j \right) = 0. \label{alg3}\end{gathered}$$

Equation  describes how agent $i$’s state $x_i^k$ evolves with iteration $k$. The local gradient ${\nabla\!}f_i$ is evaluated at $y_i^k$ and the quantity $z_i^k$ is transmitted to neighboring agents in . Finally, we allow for linear state-input invariants to be enforced in . Such invariants typically arise from requiring a particular initialization for the algorithm.

The matrices $A$, $D_{yu}$, and $D_{zv}$ are square, and the other matrices have compatible dimensions. The dimension of $A$ is the number of local states on each agent, the dimension of $D_{yu}$ is one, and the dimension of $D_{zv}$ is the number of variables that each agent transmits with neighbors at each iteration.

To simplify notation, we assume the objective function is one-dimensional ($d=1$). We can recover the general $d$ case by replacing each scalar symbol with a $1\times d$ row vector (e.g., $u_i^k\in{\mathbb{R}}^{1\times d}$) and interpreting each local gradient ${\nabla\! f}_i$ as a map from ${\mathbb{R}}^{1\times d}$ to ${\mathbb{R}}^{1\times d}$.

Not all instances of are efficiently implementable. For example, if ${D_{yu}\ne 0}$, then $y_i^k$ depends on $u_i^k$, which then depends on $y_i^k$. Such circular dependencies arise naturally in proximal algorithms, where an inner optimization problem must be solved at each iteration. For instance, given a convex differentiable $f$ and parameter $\lambda > 0$, the proximal algorithm $$x^{k+1} = \mathbf{prox}_{\lambda f}(x^k) {\colonequals}\operatorname*{\arg\min}_{x} \bigl( \lambda f(x) + \tfrac{1}{2}{\lVert{x-x^k}\rVert}^2\bigr)$$ satisfies the optimality condition $\lambda {\nabla\!}f(x^{k+1}) + x^{k+1} - x^k = 0$ and can therefore be expressed in the form of  as follows: $$\begin{aligned}
        x^{k+1} &= x^k - \lambda u^k, &
        y^k     &= x^k - \lambda u^k, &
        u^k     &= {\nabla\! f}(y^k).
        \end{aligned}$$ In the forthcoming analysis, we treat implementability and analysis separately. That is, we derive convergence rate bounds for general algorithms of the form , regardless of whether they can be efficiently implemented. However, we note that a sufficient condition for avoiding circular dependencies is if the feedthrough term satisfies $$\begin{aligned}
\label{implementable}
        {\begin{bmatrix}D_{yu} & D_{yv} \\ D_{zu} & D_{zv}\end{bmatrix}} = {\begin{bmatrix}0 & D_{yv} \\ 0 & 0\end{bmatrix}} \quad\text{or}\quad
        {\begin{bmatrix}0 & 0 \\ D_{zu} & 0\end{bmatrix}}.
    \end{aligned}$$

Putting a distributed optimization algorithm into the form of  is a straightforward algebraic exercise, which we now demonstrate for two recently proposed algorithms. These algorithms are parameterized by a stepsize $\alpha$ and a gossip matrix $W$. To relate the gossip matrix to the Laplacian matrix, we set $W = I-\mu \mathcal{L}$ for some scalar $\mu\ne 0$. This provides an additional tuning parameter, and is akin to the method of successive overrelaxation used in the numerical solutions of linear systems of equations [@sor_book].

#### EXTRA.

The EXTRA algorithm  has a state that depends on two previous timesteps. Using the authors’ recommendation of $\widetilde{W} = \tfrac{1}{2}(I+W)$ together with $W = I -\mu {\mathcal{L}}^k$, the equations become $$\begin{aligned}
    x^1 &= x^0 - \alpha {\nabla\!}f(x^0) - \mu {\mathcal{L}}^k x^0, \\
    x^{k+2} &= 2x^{k+1}-x^k -\alpha \left( {\nabla\!}f(x^{k+1}) - {\nabla\!}f(x^k) \right) - \mu {\mathcal{L}}^k \left( x^{k+1} - \tfrac{1}{2} x^k \right).\end{aligned}$$ Define the state $(x^{k+1},x^k,{\nabla\!}f(x^k))$. The outputs are now functions of the state: $y^k {\colonequals}x^{k+1}$ and $z^k {\colonequals}x^{k+1}-\tfrac{1}{2}x^k$. Finally, summing across agents (left-multiplying by $1^{\mathsf{T}}$) and using $1^{\mathsf{T}}{\mathcal{L}}^k = 0$, we find that $1^{\mathsf{T}}\left( x^{k+1} - x^k + \alpha {\nabla\!}f(x^k) \right)$ is independent of $k$, and identically zero thanks to how $x^1$ is initialized. The parameters that characterize EXTRA are shown below and in Table \[table::algo\_params\]. $$\begin{aligned}
\left[\begin{array}{c:c:c}
A & B_u & B_v \\ \hdashline
C_y & D_{yu} & D_{yv} \\ \hdashline
C_z & D_{zu} & D_{zv} \\ \hdashline
F_x & F_u & 
\end{array}\right]
&=
\left[\begin{array}{ccc:c:c}
    2 & -1 & \alpha & -\alpha & -\mu \\
    1 & 0 & 0 & 0 & 0 \\
    0 & 0 & 0 & 1 & 0 \\ \hdashline
    1 & 0 & 0 & 0 & 0 \\ \hdashline
    1 & -\tfrac{1}{2} & 0 & 0 & 0 \\ \hdashline
    1 & -1 & \alpha & 0 &
\end{array}\right].\end{aligned}$$

#### DIGing.

The DIGing algorithm [@DIGing; @QuLi], is an example of a *gradient tracking* algorithm. It begins with an arbitrary $x^0$ and has two update equations: $$\begin{aligned}
s^0 &= {\nabla\!}f(x^0), \\
x^{k+1} &= W x^k -\alpha s^k, \\
s^{k+1} &= \widetilde{W}s^k  + {\nabla\!}f (x^{k+1}) - {\nabla\!}f(x^k).\end{aligned}$$ Using the authors’ recommendation of $\widetilde{W} = W$, defining $W = I -\mu {\mathcal{L}}^k$ as before, and defining the state as $(x^{k},s^k,{\nabla\!}f(x^k))$, we find that the output is $y^k {\colonequals}x^{k+1}$, two quantities must be communicated between agents, $z^k {\colonequals}(x^k, s^k)$, and the invariant is $1^{\mathsf{T}}(s^k - {\nabla\! f}(x^k)) = 0$. The parameters that characterize DIGing are shown in Table \[table::algo\_params\].

A similar derivation can be applied to a variety of algorithms. Table \[table::algo\_params\] summarizes the parameterizations for 8 recently proposed algorithms.

\[table::algo\_params\]

  -- -- ----------------------------------------------------------------- ------------------------- --
         ${\addtolength{\arraycolsep}{-2pt}}\left[\begin{array}{cc:c:c}   EXTRA [@EXTRA]            
                            1 & \beta & -\alpha & -\gamma \\                                        
                              0 & 1 & 0 & -1 \\ \hdashline                                          
                            1 & 0 & 0 & -\delta \\ \hdashline                                       
                               1 & 0 & 0 & 0 \\ \hdashline                                          
                                      0 & 1 & 0 &                                                   
                                 \end{array}\right]$                                                

         ${\addtolength{\arraycolsep}{-2pt}}\left[\begin{array}{cc:c:c}   NIDS [@NIDS]              
                               2 & -1 & -\alpha & -\mu \\                                           
                      1 & 0 & -\alpha & -\tfrac12\mu \\ \hdashline                                  
                         1 & 0 & -\tfrac12 \mu & 0 \\ \hdashline                                    
                               1 & 0 & 0 & 0 \\ \hdashline                                          
                                      1 & -1 & 0 &                                                  
                                 \end{array}\right]$                                                

         ${\addtolength{\arraycolsep}{-2pt}}\left[\begin{array}{cc:c:cc}  DIGing [@DIGing; @QuLi]   
                           1 & -\alpha & -\alpha & -\mu & 0 \\                                      
                           0 & 1 & 0 & 0 & -\mu \\ \hdashline                                       
                             1 & 0 & 0 & 0 & 0 \\ \hdashline                                        
                                  1 & 0 & 0 & 0 & 0 \\                                              
                      -\tfrac{L+m}{2} & 1 & 1 & 0 & 0\\ \hdashline                                  
                                     0 & 1 & 0 & &                                                  
                                 \end{array}\right]$                                                

         ${\addtolength{\arraycolsep}{-2pt}}\left[\begin{array}{cc:c:cc}  AugDGM [@AugDGM]          
                           1 & -\alpha & -\alpha & -\mu & 0 \\                                      
                           0 & 1 & 0 & 0 & -\mu \\ \hdashline                                       
                             1 & 0 & 0 & 0 & 0 \\ \hdashline                                        
                                  1 & 0 & 0 & 0 & 0 \\                                              
                           -L & 1 & 1 & L\mu & 0\\ \hdashline                                       
                                     0 & 1 & 0 & &                                                  
                                 \end{array}\right]$                                                
  -- -- ----------------------------------------------------------------- ------------------------- --

Existence of a Fixed Point
--------------------------

Not all instances of algorithm  solve the distributed optimization problem . For an algorithm to be valid, (i) there must exist a fixed point corresponding to the optimal solution, and (ii) the iterates must converge to the fixed point. We address convergence to a fixed point in our main result of Section \[main\_result\]. In this section, however, we provide simple conditions for verifying the existence of such a fixed point.

A distributed algorithm of the form  has a fixed point $(x^\star, y^\star, z^\star, u^\star, v^\star)$ corresponding to the optimal solution of  for all functions satisfying Assumption \[assumption:local\_functions\] and all graphs satisfying Assumption \[assumption:Laplacian\] if the following conditions hold.

\[eq:fixedpoint\_conditions\]

-   **Consensus and Optimality:** All agents must achieve consensus on the point at which the gradient is evaluated, and the point must be a stationary (first-order optimal) point of $f$. This means that the fixed point must satisfy $y_1^\star = \ldots = y_n^\star$ and $u_1^\star + \dots + u_n^\star = 0$, or in vector form, $$\label{eq:fixedpoint_consensus_optimality}
            (I-\Pi)\, y^\star = 0 \quad\text{and}\quad
            1^{\mathsf{T}}u^\star = 0.$$

-   **Robustness to Graph:** The fixed point must not depend on the sequence of graphs $\{{\mathcal{L}}^k\}$, so $z_1^\star = \ldots = z_n^\star$ and $v_1^\star = \dots = v_n^\star = 0$, or in vector form, $$(I-\Pi)\, z^\star = 0 \quad\text{and}\quad
             v^\star = 0.$$

-   **Robustness to Functions:** The fixed point must satisfy $y_1^\star = \ldots = y_n^\star = y_\text{opt}$ and $u_i^\star = {\nabla\! f}_i(y_\text{opt})$, where $y_\text{opt}$ is the optimizer of . For these to hold for any objective function $f$, we need $$1^{\mathsf{T}}y^\star \text{ and } (I-\Pi)\, u^\star \text{ unconstrained}.$$

The following proposition characterizes algorithms with such a fixed point, which we prove in Appendix \[app:fixedpoint\].

\[prop:fixedpoint\] An algorithm of the form  has a fixed point $(x^\star,y^\star,z^\star,u^\star,v^\star)$ that satisfies the conditions in  if and only if

\[eq:fixedpoint\] $$\begin{gathered}
        \label{eq:fixedpoint1}
        \mathrm{null}(A-I) \cap \mathrm{row}(C_y) \cap \mathrm{null}(F_x) \neq \{0\}\\
        \label{eq:fixedpoint2}
        \text{and}\quad
        {\begin{bmatrix} B_u \\ D_{yu} \\ D_{zu} \end{bmatrix}} \in \mathrm{col}\!\left({\begin{bmatrix} A-I \\ C_y  \\ C_z \end{bmatrix}}\right).
    \end{gathered}$$

Here, “null”, “col”, and “row” denote the nullspace, column space, and row space, respectively. Both EXTRA and DIGing as derived above satisfy the conditions in  and therefore have a fixed point corresponding to the optimal solution of .

\[rem:convergence\] Proposition \[prop:fixedpoint\] guarantees that any instance of algorithm  [satisfying ]{} has a desirable fixed point in the presence of a time-varying graph; all agents agree on a common stationary point of . However, Proposition \[prop:fixedpoint\] does not ensure that the algorithm necessarily converges to this fixed point, nor does it characterize the rate of convergence. These questions will be explored in Section \[main\_result\].

Lower Bounds on Worst-Case Convergence Rates {#sec:lower_bounds}
--------------------------------------------

We now construct simple lower bounds on the worst-case asymptotic convergence rate of the iterates for any valid algorithm of the form . We do so by separately considering the two specific instances discussed in Section \[sec:intro\]

#### Consensus

Consider the scalar local quadratic functions $f_i(y) = \tfrac{L}{2}\,(y-r_i)^2$. Then Assumption \[assumption:local\_functions\] holds with $m=L$ and $y_\text{opt}=\tfrac{1}{n} \sum_{i=1}^n r_i$.

#### Optimization

Consider the case $n=1$. For the graph to satisfy Assumption \[assumption:Laplacian\], the Laplacian matrix must be ${\mathcal{L}}^k=0$, which has spectral gap $\sigma=0$.

In both cases above, the algorithm reduces to a linear system in feedback with sector-bounded nonlinearity: in the sector $(1-\sigma,1+\sigma)$ for consensus and $(m,L)$ for optimization. Further, the linear part of the system is strictly proper (since the algorithm is implementable) and must contain an integrator (due to the fixed-point conditions). Then using the lower bound for such systems in [@lowerbound], we obtain the following.

There does *not* exist an algorithm of the form  that satisfies the implementability conditions  and fixed-point conditions  and such that, for all objective functions and Laplacian matrices satisfying Assumptions \[assumption:local\_functions\] and \[assumption:Laplacian\], there exists a constant $c>0$ such that the bound $\|x_i^k-y_\text{opt}\|\le c\,\rho_\text{lb}^k$ holds for all agents $i\in\{1,\ldots,n\}$ and all iterations $k\ge 0$, where $\rho_\text{lb} = \max\bigl\{\tfrac{\kappa-1}{\kappa+1}, \, \sigma\bigr\}$.

\[rem:no\_free\_lunch\] These lower bounds, which are achieved by ordinary gradient descent, imply that accelerated algorithms such as the recently proposed SSDA [@Scaman17] or distributed versions of heavy-ball [@DHB] or Nesterov acceleration [@QuLi_accelerated; @ABN] do not in fact achieve accelerated rates in the worst case in our time-varying setting.

Main Result {#main_result}
===========

Our main theorem, Theorem \[thm:main-result\], consists of a small convex semidefinite program (SDP) whose feasibility guarantees the linear convergence of a distributed algorithm in the form of . The algorithm parameters, problem data $(\kappa, \sigma)$, and candidate linear rate $\rho$ all appear as parameters in the SDP. Furthermore, the SDP has a fixed size that does not depend on $n$ (the number of agents) or $d$ (the dimension of the domain of $f$) and can thus be efficiently solved using a variety of established solvers.

\[thm:main-result\] Consider the distributed optimization problem  solved using algorithm . Suppose Assumptions \[assumption:local\_functions\] and \[assumption:Laplacian\] hold and further assume the algorithm satisfies the fixed point conditions . Define the matrices $$\begin{gathered}
        M_0 {\colonequals}{\begin{bmatrix}-2mL & L+m \\ L+m & -2\end{bmatrix}} \quad\text{and}\quad
        M_1 {\colonequals}{\begin{bmatrix}\sigma^2-1 & 1 \\ 1 & -1\end{bmatrix}}.
    \end{gathered}$$ Let $\Psi$ be a matrix whose columns form a basis for the nullspace of ${\begin{bmatrix} F_x & F_u\end{bmatrix}}$. If there exist $P \succ 0$, $Q\succ 0$, and $R\succeq 0$ of appropriate sizes such that

\[eq:sdp-main\] $$\begin{aligned}
        \Psi^{\mathsf{T}}\left[\begin{array}{cc}
            A & B_u\\ I & 0 \\ \hdashline
            C_y & D_{yu}\\ 0 & I
        \end{array}\right]^{\mathsf{T}}\left[\begin{array}{cc:c}
            P & 0 & 0 \\
            0 & -\rho^2 P & 0 \\ \hdashline
            0 & 0 & M_0
        \end{array}\right]
        \left[\begin{array}{cc}
            A & B_u\\ I & 0 \\ \hdashline
            C_y & D_{yu}\\ 0 & I
        \end{array}\right]
        \Psi &\preceq 0 \label{eq:lb-rho} \\
        \left[\begin{array}{ccc}
            A & B_u & B_v \\
            I & 0 & 0 \\ \hdashline
            C_y & D_{yu} & D_{yv} \\
            0 & I & 0 \\ \hdashline
            C_z & D_{zu} & D_{zv} \\
            0 & 0 & I
        \end{array}\right]^{\mathsf{T}}\left[\begin{array}{cc:c:c}
            Q & 0 & 0 & 0 \\
            0 & -\rho^2 Q & 0 & 0 \\ \hdashline
            0 & 0 & M_0 & 0 \\ \hdashline
            0 & 0 & 0 & M_1\otimes R
        \end{array}\right]
        \left[\begin{array}{ccc}
            A & B_u & B_v \\
            I & 0 & 0 \\ \hdashline
            C_y & D_{yu} & D_{yv} \\
            0 & I & 0 \\ \hdashline
            C_z & D_{zu} & D_{zv} \\
            0 & 0 & I
        \end{array}\right]
            &\preceq 0
        \label{eq:lmi-dis}
    \end{aligned}$$

then there exists a constant $c>0$ independent of $i$ and $k$ such that for all agents $i\in\{1,\dots,n\}$ and all iterations $k\ge 0$, $$\label{eq:bound}
        {\lVert{ x_i^k - x_i^\star }\rVert} \le c\,\rho^k$$ for some fixed point $(x_i^\star, y_i^\star, z_i^\star, u_i^\star, v_i^\star)$ that satisfies .

[For fixed algorithm parameters $A,B_u,B_v,C_y,C_z,D_{yu}$, $D_{yv},D_{zu},D_{zv},F_x,F_u$, function parameters $m$ and $L$, graph parameter $\sigma$, and candidate rate $\rho$, the SDP  is a linear matrix inequality (LMI) in the variables $(P,Q,R)$, and therefore convex. Indeed, and are decoupled and their feasibility may be checked separately.]{} To find the best (smallest) upper bound, we observe that feasibility of  for some $\rho_0$ implies feasibility for all $\rho \ge \rho_0$. A bisection search on $\rho$ is then guaranteed to find the minimal $\rho$, even though  is not jointly convex in $(P,Q,R,\rho)$. [While our result is only a sufficient condition for convergence, we provide empirical evidence in Section \[sec:worst-case\] that suggests that it is in fact tight.]{}

\[rem:bounds\] Our main theorem provides conditions under which the *state* converges to a fixed point linearly with rate $\rho$. However, when the algorithm also satisfies the conditions in  for being efficiently implementable, then under the conditions of Theorem \[thm:main-result\], there exist constants $c_u$, $c_v$, $c_y$, and $c_z$ such that for all agents $i$ and all iterations $k$, $$\begin{aligned}
        \|u_i^k - u_i^\star\| &\le c_u\,\rho^k, &
        \|y_i^k - y_i^\star\| &\le c_y\,\rho^k, &
        \|v_i^k - v_i^\star\| &\le c_v\,\rho^k, &
        \|z_i^k - z_i^\star\| &\le c_z\,\rho^k,
    \end{aligned}$$ for some fixed point $(x_i^\star, y_i^\star, z_i^\star, u_i^\star, v_i^\star)$ that satisfies . In particular, the output sequence $y_i^k$ of each agent converges to the optimizer $y_\text{opt}$ of  linearly with rate $\rho$.

The core idea behind Theorem \[thm:main-result\] is to posit a quadratic Lyapunov candidate of the form [$$\begin{aligned}
\label{eq:V}
V^k{\colonequals}(x^k - x^\star)^{\mathsf{T}}\bigl( \Pi\otimes P + (I-\Pi)\otimes Q\bigr) (x^k - x^\star)\end{aligned}$$ for some appropriate choice of $P,Q \succ 0$.]{} Feasibility of  can be shown to imply $V^{k+1}\leq \rho^2 V^k$, which ensures linear convergence of the distributed optimization algorithm when $\rho<1$. A preliminary (and less concise) version of Theorem \[thm:main-result\] appeared in [@sundararajan_allerton]. [The proof of Theorem \[thm:main-result\] is given in Appendix \[proof\_main\].]{}

Algorithm Design {#sec:algo_design}
================

We now use Theorem \[thm:main-result\] to design a distributed optimization algorithm, which we name SVL. Our guiding principle is to seek the fastest possible rate bound guarantee while keeping the algorithm as simple as possible. Therefore, we seek an algorithm with two states that only requires one state to be communicated at every timestep. Inspired by our previous work in which we developed a canonical form for distributed algorithms over time-invariant graphs [@canform], we restrict our search to algorithms of the form  with $$\begin{aligned}
\label{eq:SVL-template}
    \left[\begin{array}{c:c:c}
    A & B_u & B_v \\ \hdashline
    C_y & D_{yu} & D_{yv} \\ \hdashline
    C_z & D_{zu} & D_{zv} \\ \hdashline
    F_x & F_u &
    \end{array}\right]
    &=
    \left[\begin{array}{cc:c:c}
        1 & \beta & -\alpha & -\gamma \\
        0 & 1 & 0 & -1 \\ \hdashline
        1 & 0 & 0 & -\delta \\ \hdashline
        1 & 0 & 0 & 0 \\ \hdashline
        0 & 1 & 0 & 
    \end{array}\right].\end{aligned}$$ As long as $\beta \ne 0$, this algorithm satisfies the fixed point conditions of Proposition \[prop:fixedpoint\]. Moreover, the update equations satisfy  and therefore do not contain circular dependencies, so we can implement the algorithm in a straightforward fashion as in Algorithm \[alg:canonical\]. To motivate the structure of our algorithm, we show how it corresponds to an inexact version of the alternating direction method of multipliers (ADMM), as well as how it reduces to well-known consensus and optimization algorithms in special cases. But first, we show how to use the SDP  to choose the algorithm parameters.

[ **Initialization:** Let ${\mathcal{L}}^k \in{\mathbb{R}}^{n\times n}$ be a Laplacian matrix. [Agents $i\in\{1,\ldots,n\}$ choose initial local state $x_i^0 \in {\mathbb{R}}^d$ arbitrarily and $w_i^0\in{\mathbb{R}}^d$ such that $\sum_{i=1}^{n}w_i^0=0$ (e.g. $w_i^0=0$).]{} **Local communication** $~~\hphantom{w_i^{k+1}}\mathllap{v_{i}^k} = \sum_{j=1}^n {\mathcal{L}}^k_{ij}\,x_j^k \hfill \text{(C.1)}$ **Local gradient computation** $~~\hphantom{w_i^{k+1}}\mathllap{y_i^k} = x_i^k - \delta\,v_{i}^k \hfill \text{(C.2)}$ $~~\hphantom{w_i^{k+1}}\mathllap{u_i^k} = {\nabla\! f}_i(y_i^k) \hfill \text{(C.3)}$ **Local state update** $~~\hphantom{w_i^{k+1}}\mathllap{x_i^{k+1}} = x_i^k + \beta\,w_i^k - \alpha\,u_i^k - \gamma\, v_{i}^k \hfill \text{(C.4)}$ $~~\hphantom{w_i^{k+1}}\mathllap{w_i^{k+1}} = w_i^k - v_{i}^k \hfill \text{(C.5)}$ ]{}

Choosing the Algorithm Parameters
---------------------------------

The problem of minimizing the worst-case convergence rate $\rho$ over the algorithm parameters $(\alpha,\beta,\gamma,\delta)$ and SDP solution [$(P,Q,R)$]{} subject to the SDP being feasible is difficult due to the nonlinear matrix [inequalities ]{}. Instead, we show that for a particular choice of $(\alpha,\gamma,\delta)$, the remaining parameters $(\beta,\rho)$ can be chosen such that the [SDP is feasible, where the]{} matrix in  is rank one. We have performed extensive numerical optimizations of the SDP, suggesting that the optimal parameters do in fact have this structure. We now state our main design result, which describes [the convergence rate of]{} the SVL algorithm. We prove the result in Appendix \[sec:SVL-proof\].

\[thm:SVL\] Consider applying Algorithm \[alg:canonical\] to the distributed optimization problem , and suppose Assumptions \[assumption:local\_functions\] and \[assumption:Laplacian\] hold [with $0<m<L$ and $0\le\sigma < 1$]{}. Define $\eta{\colonequals}1+\rho-\kappa\,(1-\rho)$ and choose the parameters $$\begin{aligned}
\label{eq:paramopt}
  \alpha &= \frac{1-\rho}{m}, &
  \gamma &= 1+\beta, &
  \delta &= 1,
  \end{aligned}$$ where $\beta$ and ${\rho\in\bigl[\tfrac{L-m}{L+m},1\bigr)}$ satisfy the constraints

\[eq:sol-opt\] $$\begin{aligned}
  \bigl(2\beta-(1-\rho)(\kappa+1)\bigr)(\beta-1+\rho^2) &< 0, \label{eq:beta-constraint}\\
  \rho^2\,\bbbl(\frac{\beta-1+\rho^2}{\beta-1+\rho}\bbbr)
    \bbbl(\frac{2-\eta-2\beta}{2\rho^2\beta - (1-\rho^2)\eta}\bbbr) 
    \bbbl(\frac{(2\rho^2+\eta)\beta - (1-\rho^2)\eta}{(1+\rho)(\eta-2\eta\rho+2\rho^2)-(2\rho^2+\eta)\beta}\bbbr) &= \sigma^2. \label{eq:sigopt}
  \end{aligned}$$

Then there exists a constant $c>0$ independent of $i$ and $k$ such that for all agents $i\in\{1,\dots,n\}$ and all iterations $k\ge 0$, ${\lVert{ y_i^k - y_\text{opt} }\rVert} \le c\,\rho^k $ where $y_\text{opt}\in{\mathbb{R}}^d$ is the optimizer of .

Theorem \[thm:SVL\] provides conditions on parameters $(\alpha,\beta,\gamma,\delta)$ of Algorithm \[alg:canonical\] such that the algorithm converges with rate at least $\rho$. The theorem, however, does not address the problem of optimizing the convergence rate since $\beta$ and $\rho$ must only be chosen to satisfy the constraints . This is because the optimal parameters do not admit a closed-form solution for the convergence rate $\rho$ as a function of the spectral gap $\sigma$ and function parameters $m$ and $L$. However, we now provide a systematic method for computing the optimal parameters.

The parameters must satisfy , but this equation does not have a closed-form solution for $\rho$. Instead, we consider fixing the [rate $\rho$]{} and maximizing the corresponding spectral gap. We can then choose $\beta$ to maximize $\sigma^2$ in . Setting the derivative equal to zero, we find that the value of $\beta$ which maximizes $\sigma^2$ for a fixed convergence rate $\rho$ satisfies $$\frac{\textrm{d}\sigma^2}{\textrm{d}\beta}=0 \quad\implies\quad 
\bigl(\beta\bigl(1-\kappa+2\rho(1+\rho)\bigr)-\eta(1-\rho^2)\bigr)
\bigl(s_0 + s_1\beta + s_2\beta^2 + s_3\beta^3\bigr)=0,$$ where the coefficients $s_i$ are given by $$\begin{aligned}
s_0 &{\colonequals}\eta\,\bigl(1-\rho^2\bigr)^2 \bigl(\eta-(3-\eta)\eta\rho+2(1-\eta)\rho^2+2\rho^3\bigr), \\
s_1 &{\colonequals}-\bigl(1-\rho^2\bigr) \Bigl(\eta^3\rho+4\rho^5-2\eta\rho^2(2\rho^2+\rho-3)+\eta^2\,\bigl(4 \rho^3-4 \rho^2-6 \rho+3\bigr)\Bigr), \\
s_2 &{\colonequals}3\eta(1-\rho)^2 (1+\rho) (2\rho^2+\eta), \\
s_3 &{\colonequals}(2\rho^2+\eta) (2\rho^3-\eta).\end{aligned}$$ Solving the first factor for $\beta$, we find that it does not satisfy the inequality  and is therefore not a valid solution. The optimal $\beta$ must then make the second factor zero. Therefore, we can do a bisection search over $\rho$, where at each iteration of the bisection search we solve the cubic equation $$\begin{aligned}
\label{eq:cubic}
s_0 + s_1\beta + s_2\beta^2 + s_3\beta^3 = 0\end{aligned}$$ to find the unique real solution $\beta$ that satisfies . Substituting this value for $\beta$ into  we can solve for $\sigma$. If this value is less than $\sigma$, we increase $\rho$; otherwise, we decrease $\rho$. We then repeat this procedure until $\sigma$ is sufficiently close to the spectral gap. We summarize this procedure for finding the parameters $\beta$ and $\rho$ that optimize the worst-case convergence rate in Algorithm \[alg:paramopt\]; we refer to Algorithm \[alg:canonical\] using these parameters along with those in  as SVL.

**Initialization:** Let $0<m<L$, $0\le\sigma<1$, and ${\varepsilon}>0$. Define $\kappa{\colonequals}L/m$. Set $\rho_1=0$ and $\rho_2=1$. $\rho = (\rho_1+\rho_2)/2$ Let $\beta$ be the unique real solution to  that satisfies . Using this value of $\beta$, let $\hat{\sigma}$ denote the solution to . $\rho_1 = \rho$ $\rho_2 = \rho$\
$\rho,\beta$

Using this procedure for computing the worst-case convergence rate of SVL, Fig. \[fig:optalg\] displays $\rho$ as a function of the spectral gap $\sigma$ and the centralized gradient rate $\tfrac{\kappa-1}{\kappa+1}$. One of the remarkable aspects of the SVL algorithm is that it actually achieves the same worst-case convergence rate as *centralized* gradient descent if the spectral gap is sufficiently small. In this case, there is sufficient mixing among the agents so that the convergence rate is limited by the difficulty of the optimization problem and not the problem of having agents agree on the solution (i.e., consensus). This corresponds to the horizontal lines for small values of $\sigma$ in the [top]{} panel of Fig. \[fig:optalg\]. Viewed another way, the convergence rate is limited by the difficulty of the optimization problem when the problem is ill-conditioned (i.e., $\kappa$ is large), which corresponds to the curves approaching the straight line at $\rho=\tfrac{\kappa-1}{\kappa+1}$ in the [bottom]{} panel of Fig. \[fig:optalg\].

We conjecture that the SVL parameters $(\alpha,\beta,\gamma,\delta)$ produce the fastest worst-case convergence rate over all algorithms in the form of Algorithm \[alg:canonical\] that is certifiable using Theorem \[thm:main-result\]. However, we make no formal claims of optimality of the SVL algorithm in this paper.

![Worst-case [linear rate $\rho$]{} of SVL in Theorem \[thm:SVL\] [as a function of $\kappa$ and $\sigma$. Top plot: as $\kappa\to 1$ (quadratic objective), we obtain $\rho = \sigma$ (optimal linear consensus rate). Bottom plot: as $\sigma\to 0$ (fully connected graph), we obtain $\rho = \tfrac{\kappa-1}{\kappa+1}$ (optimal centralized gradient rate).]{} []{data-label="fig:optalg"}](SVL_plot)

Interpretation of SVL as Inexact ADMM {#sec:admm}
-------------------------------------

To motivate the structure of SVL, we show how SVL can be interpreted as an inexact version of the alternating direction method of multipliers (ADMM). Using the formulation in [@ADMM Section 7.1], the problem  can be solved using ADMM:

$$\begin{aligned}
x_i^{k+1} &= \operatorname*{\arg\min}_{x} f_{i}(x)+(x-y_i^k)^{\mathsf{T}}z_i^k+\tfrac{\beta}{2} {\lVert{x-y_i^k}\rVert}^2 \label{eq:admm1} \\
y_i^{k+1} &= \frac{1}{n} \sum_{j=1}^n x_j^{k+1} \label{eq:admm2} \\
z_i^{k+1} &= z_i^k + \beta\,(x_i^{k+1}-y_i^{k+1}) \label{eq:admm3}\end{aligned}$$

where $(x_i^k,y_i^k,z_i^k)$ are the variables associated with agent $i$ at time $k$, and $\beta$ is the ADMM parameter. To implement this algorithm, however, each agent must solve the local optimization problem  *exactly* as well as compute the *exact* average  at each iteration. Instead, we consider a variant where the computations and communications are *inexact*. Specifically, we replace the exact minimization  with a single gradient step with initial condition $y_i^k$ and stepsize $\alpha>0$, and we replace the exact averaging step  with a single gossip step using the Laplacian matrix ${\mathcal{L}}^k$. This gives the following inexact version of ADMM: $$\begin{aligned}
x_i^{k+1} &= y_i^k - \alpha\,\bigl({\nabla\! f}(y_i^k)+z_i^k\bigr) \\
y_i^{k+1} &= x_i^{k+1} - \sum_{j=1}^n {\mathcal{L}}_{ij}^{k+1}\,x_j^{k+1} \\
z_i^{k+1} &= z_i^k + \beta\,(x_i^{k+1}-y_i^{k+1})\end{aligned}$$ Defining the state $w_i^k {\colonequals}-\tfrac{\alpha}{\beta} z_i^{k-1}$, this algorithm is equivalent to Algorithm \[alg:canonical\] with $\gamma=1+\beta$ and $\delta=1$. In other words, SVL corresponds to an inexact version of ADMM, where $\alpha$ is the stepsize of the gradient step and $\beta$ is the ADMM parameter. See [@DADMM1; @DADMM2] for other distributed ADMM variants.

![image](algo_comparison)

Special Cases {#sec:special-cases}
-------------

We now show how the SVL algorithm reduces to well-known consensus and optimization algorithms in special cases.

#### $n=1$:

[With]{} only [one]{} agent, the distributed optimization problem  is equivalent to centralized optimization. In this case, the Laplacian matrix is simply the scalar ${\mathcal{L}}^k=0$, so $v_1^k=0$ for all $k\ge 0$. Algorithm \[alg:canonical\] then simplifies to $$x_1^{k+1} = x_1^k - \alpha\,{\nabla\! f}(x_1^k), \qquad x_1^0\text{ arbitrary,}$$ which is ordinary gradient descent with stepsize $\alpha$. The fastest possible gradient rate of $\rho=\frac{\kappa-1}{\kappa+1}$ is achieved when $\alpha=\frac{2}{L+m}$.

#### $\kappa=1$:

When the condition ratio is unity (i.e., $m=L$), the distributed optimization problem  is equivalent to average consensus. In this case, the parameters of SVL are simply $\alpha=\tfrac{1}{L}$, $\beta=1$, $\gamma=2$, and $\delta=1$. Also, the objective functions are quadratic, so we may assume without loss of generality that they have the form ${f_i^k}(x) = \tfrac{L}{2} \|x-{r_i^k}\|^2$, where ${r_i^k}\in{\mathbb{R}}^d$ is a parameter on agent $i\in\{1,\ldots,n\}$ [at iteration $k$]{}. The SVL algorithm then simplifies to $$x_i^{k+1} = x_i^k - \sum_{j=1}^n {\mathcal{L}}_{ij}^k\,x_j^k + \bigl(r_i^k-r_i^{k-1}\bigr), \qquad x_i^0 = r_i^0,$$ which is a dynamic average consensus algorithm since the reference signals are continually injected into the dynamics [@kia2019tutorial]. [When the objective functions are constant, the $r_i$ terms cancel from the iterations and only affect the initial conditions.]{} [This case]{} is [referred to as *static*]{} average consensus [@tsitsiklis], [and]{} the worst-case rate of convergence is $\rho=\sigma$ [@xiaoboyd04].

Numerical Results {#numerics}
=================

In this section, we compare the worst-case performance of SVL with [that of]{} other first-order distributed algorithms.

Algorithm Comparison (Upper Bounds) {#sec:upper_bounds}
-----------------------------------

[Theorem \[thm:main-result\] provides an upper bound on the worst-case convergence rate. We used this result to compare all algorithms]{} in Table \[table::algo\_params\], [including SVL]{}. The results are shown in Fig. \[fig:money\_plot\]. For each algorithm, we used a bisection search to find the smallest rate $\rho$ that yielded a feasible solution to the SDP . We implemented the SDP in Julia [@julia] with the JuMP [@jump] modeling package and the Mosek [interior point]{} solver [@mosek]. [In an outer loop,]{} we performed a parameter search for each algorithm [to find the step size $\alpha$ and overrelaxation parameter $\mu$]{} that yielded the smallest possible $\rho$. Specifically, we used Brent’s method and the Nelder–Mead method, respectively, as implemented in the Optim package [@optim] as $\sigma$ ranged from 0 to 1.

As shown in Fig. \[fig:money\_plot\], [optimizing over $\mu$ further improves worst-case performance]{}. Our proposed SVL algorithm outperforms all methods we tested. [Also shown in Fig. \[fig:money\_plot\] is the lower bound described in Section \[sec:lower\_bounds\], namely $\rho \ge \max\{\tfrac{\kappa-1}{\kappa+1},\sigma\}$, which holds for any distributed algorithm.]{}

![image](trajectories)

Approximate Worst-Case Examples (Lower Bounds) {#sec:worst-case}
----------------------------------------------

[In an effort to show that the upper bounds for each algorithm in Fig 2 were likely tight, we searched for signals $\{x^k,u^k,v^k,y^k,z^k\}$ that satisfied  for some choice of $f_i$ and ${\mathcal{L}}^k$ satsifying Assumptions \[assumption:local\_functions\] and \[assumption:Laplacian\], respectively.]{}

[We first solved a relaxed version of the problem, where we replaced Assumptions \[assumption:local\_functions\] and \[assumption:Laplacian\] by the weaker conditions  and , respectively. We used the following greedy heuristic. For a given algorithm and rate $\rho$, we solved  to obtain $(P,Q,R)$. At each time step $k$, we then maximized the Lyapunov increment $V^{k+1}-\rho^2 V^k$, where $V^k$ is defined in . We solved the following optimization problem for $k\ge 0$.]{}

[$$\label{opt:low_bound}
    \begin{aligned}
    \operatorname*{maximize}_{u^k_i,v^k_i\in{\mathbb{R}}^d} \qquad &V^{k+1}-\rho^2 V^k \\
    \text{such that}\qquad & \text{\eqref{alg1}, \eqref{alg3}, and \eqref{quad:graph} hold,} \\
    & \text{\eqref{quad:func} holds for $i=1,\dots,n$,} \\
    & 1^{\mathsf{T}}v^k = 0.
\end{aligned}$$ For $k=0$, we also included $x^0$ as an optimization variable and the normalization $V^0=1$. For $k\ge 1$, we solved  using the $x^k$ found at the previous iteration and warm-starting $u^k,v^k$. We used the Ipopt [@Ipopt] local solver with default settings since  is a nonconvex quadratically constrained quadratic program. Note that we must choose parameters $n$ and $d$.]{}

[Our relaxed heuristic using $n=d=2$ was successful in constructing trajectories that matched the worst-case bounds from . To illustrate, we simulated EXTRA, NIDS, DIGing, and SVL with $\kappa=10$ and a few values of $\sigma$ in Fig. \[fig:wc-ex\]. For each trajectory, we plotted $\|y^k-y^\star\|$ together with the corresponding upper bound $\rho$ found from Theorem \[thm:main-result\]. We obtained similar results for the other algorithms from Table \[table::algo\_params\].]{}

[Since we used the relaxation  to construct $z^k$ and $v^k$, there is no guarantee that there will exist a *linear* Laplacian ${\mathcal{L}}^k$ such that $v^k = {\mathcal{L}}^k z^k$. However, finding whether such an ${\mathcal{L}}^k$ exists amounts to solving a convex optimization problem: $$\label{eq:constructL}
\begin{aligned}
\operatorname*{minimize}_{{\mathcal{L}}^k \in {\mathbb{R}}^{n\times n}} \qquad & {\lVert{I-\Pi-{\mathcal{L}}^k}\rVert} \\
\text{such that}\qquad &({\mathcal{L}}^k \otimes I) z^k = v^k, \\
&{\mathcal{L}}^k 1 = 0,\quad 1^{\mathsf{T}}{\mathcal{L}}^k = 0.
\end{aligned}$$ If  is feasible and its optimal value is less than or equal to $\sigma$, then the associated ${\mathcal{L}}^k$ is a valid Laplacian matrix at timestep $k$. While there is no guarantee that will even be feasible, we reasoned that since there are $n^2$ variables and $2n+ndc$ linear constraints, where $d$ and $c$ are the number of rows of $C_y$ and $C_z$, respectively, we could increase our chances of finding feasible ${\mathcal{L}}^k$ with $n$ large and $d$ and $c$ small.]{}

[In Figure \[fig:NIDS-unstable\], we show a successful construction for the NIDS algorithm, which has $c=1$. We solved  with $n=15$ and $d=1$, and solved  at each timestep. An optimal cost for  of $\sigma$ was always achieved.]{} [This result indicates that the upper bound for NIDS in Fig. \[fig:money\_plot\] is likely tight, and that NIDS is not robustly stable in the time-varying setting. In other words, the *network-independent* rate bound enjoyed by NIDS in the constant-graph setting [@NIDS Thm. 2] does not carry over to the time-varying setting.]{}

![[Worst-case trajectories for NIDS found by solving  and successfully solving  to construct a sequence of Laplacians $\{{\mathcal{L}}^k\}$. Simulations were performed using optimized $\alpha$, $\mu=1$, $\kappa=10$, $n=15$, and $d=1$ for $\sigma\in\{0.3,0.6,0.9\}$. Trajectories were plotted with their accompanying rate bounds (dashed lines) from Theorem \[thm:main-result\] and translated to improve clarity.]{}[]{data-label="fig:NIDS-unstable"}](NIDS-unstable)

[There may be other approaches to finding a worst-case ${\mathcal{L}}^k$ that perform better. For example, one might try alternating convex optimizations or including ${\mathcal{L}}^k$ directly as an optimization variable in a nonlinear program.]{}

Conclusion
==========

We presented a universal analysis framework for a broad class of first-order distributed optimization algorithms over time-varying graphs. The framework provides worst-case certificates of linear convergence via semidefinite programming, and we show empirically that our rate bounds are likely tight. Optimizing the SDP from our analysis framework, we designed a novel distributed algorithm, SVL, which outperforms [all known algorithms in this time-varying setting]{}.

Appendix
========

Proof of Proposition \[prop:fixedpoint\] {#app:fixedpoint}
----------------------------------------

Suppose  holds, and denote the optimizer of  by $y_\text{opt}$. Then there exist vectors $p$ and $q$ such that $$\begin{aligned}
    \left\{\begin{aligned}
        0 &= (A-I)\,p \\
        y_\text{opt} &= C_y p \\
        0 &= F_x p
    \end{aligned}\right.
    \qquad\text{and}\qquad
    \left\{\begin{aligned}
        B_u &= (A-I)\,q \\
        D_{yu} &= C_y q \\
        D_{zu} &= C_z q.
    \end{aligned}\right.\end{aligned}$$ For all $i\in\{1,\ldots,n\}$, use these vectors to define the points $$\begin{aligned}
    x_i^\star &= p - q\,{\nabla\! f}_i(y_\text{opt}), &
    y_i^\star &= y_\text{opt}, &
    z_i^\star &= C_z p, \\
    u_i^\star &= {\nabla\! f}_i(y_\text{opt}), &
    v_i^\star &= 0.\end{aligned}$$ This is a fixed point of algorithm , and the fixed point satisfies the conditions in  since $y_\text{opt}$ is the optimizer of .

Now suppose $(x^\star,y^\star,z^\star,u^\star,v^\star)$ is a fixed point of  satisfying . Let $p = (1/n) \sum_{i=1}^n x_i^\star$. Since $1^{\mathsf{T}}u^\star = 0$, $v^\star = 0$, and $1^{\mathsf{T}}y^\star$ is unconstrained, we have from  and  that $p\ne 0$ is in the set . Now let $v$ be any nonzero vector such that $v^{\mathsf{T}}1=0$. Then from , we have that $$0 = {\begin{bmatrix}A-I \\ C_y \\ C_z\end{bmatrix}} (v^{\mathsf{T}}x^\star) + {\begin{bmatrix}B_u \\ D_{yu} \\ D_{zu}\end{bmatrix}} (v^{\mathsf{T}}u^\star).$$ Since this must hold for arbitrary $v^{\mathsf{T}}u^\star$, this implies . 

Proof of Theorem \[thm:main-result\] {#proof_main}
------------------------------------

Assumptions \[assumption:local\_functions\] and \[assumption:Laplacian\] lead to quadratic inequalities that will be useful in proving our main result. These are stated in the following propositions.

\[prop:functions\] Suppose Assumption \[assumption:local\_functions\] holds for the local objective functions $f_i$. [Let $(y_i^k,u_i^k)$ satisfy , and let $(y_i^\star,u_i^\star)$ be a fixed point that satisfies . Then $$\label{quad:func}
        {\begin{bmatrix}\tilde y^k \\ \tilde u^k\end{bmatrix}}^{\mathsf{T}}(M_0\otimes I) {\begin{bmatrix}\tilde y^k \\ \tilde u^k\end{bmatrix}} \ge 0.$$ ]{}

[[**Proof.**]{} Using the definition of $M_0$, the quadratic form is $${\begin{bmatrix}\tilde y^k \\ \tilde u^k\end{bmatrix}}^{\mathsf{T}}(M_0\otimes I) {\begin{bmatrix}\tilde y^k \\ \tilde u^k\end{bmatrix}}
        = -2\sum_{i=1}^n (\tilde u_i^k - m\tilde y_i^k)^{\mathsf{T}}(\tilde u_i^k - L\tilde y_i^k).$$ Since the fixed point satisfies , Assumption \[assumption:local\_functions\] implies that this is nonnegative with $y_\text{opt} = y_1^\star = \ldots = y_n^\star$. [ ]{}]{}

\[prop:graph\] Suppose Assumption \[assumption:Laplacian\] holds for the graph $\mathcal{G}^k$ at each iteration. [Let $(z_i^k,v_i^k)$ satisfy , and let $(z_i^\star,v_i^\star)$ be a fixed point that satisfies . Then for all $R\succeq 0$, $$\label{quad:graph}
        {\begin{bmatrix}\tilde z^k \\ \tilde v^k\end{bmatrix}}^{\mathsf{T}}\bigl( M_1\otimes (I-\Pi)\otimes R\bigr) {\begin{bmatrix}\tilde z^k \\ \tilde v^k\end{bmatrix}} \ge 0.$$ ]{}

From the definition of the matrix norm and Assumption \[assumption:Laplacian\], we have that $$\begin{aligned}
    \sigma \ge {\bl\lVert{I-\Pi-{\mathcal{L}}^k}\br\rVert}
    &= {\bl\lVert{(I-\Pi-{\mathcal{L}}^k)(I-\Pi)}\br\rVert} \\
    &= \max_{y\in {\mathbb{R}}^{n}, y\ne 0} \frac{{\bl\lVert{(I-\Pi-{\mathcal{L}}^k)(I-\Pi)y}\br\rVert}}{{\lVert{y}\rVert}}.\end{aligned}$$ Without loss of generality, $y = \Pi \eta + (I-\Pi)\,\phi$, where $\eta$ and $\phi$ are arbitrary. By orthogonality, ${\lVert{y}\rVert}^2 = {\lVert{\Pi \eta}\rVert}^2 + {\lVert{(I-\Pi)\,\phi}\rVert}^2$. Substituting the decomposition of $y$ into the above [inequality]{}, $$\begin{aligned}
    \sigma &\ge
    \max_{\phi,\eta\in{\mathbb{R}}^n, y\ne 0} \frac{{\bl\lVert{(I-\Pi-{\mathcal{L}}^k)(I-\Pi)\,\phi}\br\rVert}}{\sqrt{{\lVert{\Pi \eta}\rVert}^2 + {\lVert{(I-\Pi)\,\phi}\rVert}^2 }} \\
    &= \max_{\phi\in{\mathbb{R}}^n, y\ne 0} \frac{{\bl\lVert{(I-\Pi-{\mathcal{L}}^k)(I-\Pi)\,\phi}\br\rVert}}{{\lVert{(I-\Pi)\,\phi}\rVert}} \\
    &= \max_{\phi\in{\mathbb{R}}^n, y\ne 0} \frac{{\bl\lVert{(I-\Pi)(\phi-{\mathcal{L}}^k \phi)}\br\rVert}}{{\lVert{(I-\Pi)\,\phi}\rVert}},\end{aligned}$$ where the last two steps follow because the maximum is attained with $\eta=0$, and ${\mathcal{L}}^k \Pi = \Pi {\mathcal{L}}^k = {\mathbf{0}}$. Squaring both sides and rewriting as a quadratic form yields $$\begin{aligned}
\label{ineq_simple}
    {\begin{bmatrix}\phi \\ {\mathcal{L}}^k\phi\end{bmatrix}}^{\mathsf{T}}\bigl(M_1\otimes (I-\Pi)\bigr) {\begin{bmatrix}\phi \\ {\mathcal{L}}^k\phi\end{bmatrix}} \ge 0\end{aligned}$$ [for all $\phi\in{\mathbb{R}}^n$. Now let $p$ denote the dimension of $z_i^k$. Then since $R\succeq 0$, it has the decomposition $$R = \sum_{\ell=1}^p \mu_\ell \, w_\ell w_\ell^{\mathsf{T}},$$ where $w_\ell\in{\mathbb{R}}^p$ and $\mu_\ell\ge 0$. Then using that $\tilde v^k = ({\mathcal{L}}^k\otimes I_p)\,\tilde z^k$, the quadratic form is $$\begin{aligned}
    {\begin{bmatrix}\tilde z^k \\ \tilde v^k\end{bmatrix}}^{\mathsf{T}}\bigl(M_1\otimes (I-\Pi)\otimes R\bigr) {\begin{bmatrix}\tilde z^k \\ \tilde v^k\end{bmatrix}}
    &= \sum_\ell \mu_\ell {\begin{bmatrix}\star\end{bmatrix}}^{\mathsf{T}}\bigl(M_1\otimes (I-\Pi)\bigr) {\begin{bmatrix}(I\otimes w_\ell^{\mathsf{T}})\,\tilde z^k \\[1pt] (I\otimes w_\ell^{\mathsf{T}})\,\tilde v^k\end{bmatrix}} \\
    &= \sum_\ell \mu_\ell {\begin{bmatrix}\star\end{bmatrix}}^{\mathsf{T}}\bigl(M_1\otimes (I-\Pi)\bigr) {\begin{bmatrix}(I\otimes w_\ell^{\mathsf{T}})\,\tilde z^k \\[1pt] {\mathcal{L}}^k\,(I\otimes w_\ell^{\mathsf{T}})\,\tilde z^k\end{bmatrix}},\end{aligned}$$ which is nonnegative from  with $\phi \leftarrow (I\otimes w_\ell^{\mathsf{T}})\,\tilde z^k$. ]{}

Let $(x^k,y^k,z^k,u^k,v^k)$ denote a trajectory of algorithm . Since the algorithm satisfies the fixed point conditions  (by assumption), we have from Proposition \[prop:fixedpoint\] that there exists a fixed point $(x^\star,y^\star,z^\star,u^\star,v^\star)$ satisfying . The global optimizer is unique from Assumption \[assumption:local\_functions\], so the fixed point conditions  imply that $y_1^\star = \ldots = y_n^\star = y_\text{opt}$ with $y_\text{opt}$ the optimizer of .

Since the trajectory satisfies the invariant  and the columns of $\Psi$ form a basis for the nullspace of ${\begin{bmatrix}F_x & F_u\end{bmatrix}}$, there exists a vector $\tilde s^k$ such that $$\Psi\, \tilde s^k = \frac{1}{\sqrt{n}}\sum_{i=1}^n {\begin{bmatrix}\tilde x_i^k \\ \tilde u_i^k\end{bmatrix}}.$$ Multiplying the matrix in  on the right and left by $\tilde s^k$ and its transpose, respectively, we obtain the consensus inequality

\[eq:inequalities\] $$\label{eq:cons_ineq}
        (\tilde x^{k+1})^{\mathsf{T}}(\Pi\otimes P)\, \tilde x^{k+1}
        - \rho^2\, (\tilde x^k)^{\mathsf{T}}(\Pi\otimes P)\, \tilde x^k \\
        + {\begin{bmatrix}\tilde y^k \\ \tilde u^k\end{bmatrix}}^{\mathsf{T}}(M_0\otimes\Pi) {\begin{bmatrix}\tilde y^k \\ \tilde u^k\end{bmatrix}} \leq 0.$$ Now choose the vectors $w_2,\ldots,w_n\in{\mathbb{R}}^n$ such that the matrix ${\begin{bmatrix}1_n/\sqrt{n} & w_2 & \ldots & w_n\end{bmatrix}}$ is orthonormal. Then we can multiply the matrix in  on the right and left by the weighted sum $$\sum_{i=1}^n (w_\ell)_i {\begin{bmatrix}\tilde x_i^k \\ \tilde u_i^k \\ \tilde v_i^k\end{bmatrix}}$$ and its transpose, respectively, and sum over $\ell\in\{2,\ldots,n\}$ to obtain the disagreement inequality $$\begin{gathered}
\label{eq:dis_ineq}
    (\tilde x^{k+1})^{\mathsf{T}}\bigl((I-\Pi)\otimes Q\bigr) \tilde x^{k+1} - \rho^2 (\tilde x^k)^{\mathsf{T}}\bigl((I-\Pi)\otimes Q\bigr)\tilde x^k \\
    + {\begin{bmatrix}\tilde y^k \\ \tilde u^k\end{bmatrix}}^{\mathsf{T}}\bigl(M_0\otimes (I-\Pi)\bigr) {\begin{bmatrix}\tilde y^k \\ \tilde u^k\end{bmatrix}}
    + {\begin{bmatrix}\tilde z^k \\ \tilde v^k\end{bmatrix}}^{\mathsf{T}}\bigl( M_1\otimes (I-\Pi)\otimes R\bigr) {\begin{bmatrix}\tilde z^k \\ \tilde v^k\end{bmatrix}} \le 0,
    \end{gathered}$$

where we used that $\{w_i\}_{i=1}^n$ form an orthonormal basis for ${\mathbb{R}}^n$. Summing the inequalities in , we obtain $$V^{k+1} - \rho^2 V^k
    + {\begin{bmatrix}\tilde y^k \\ \tilde u^k\end{bmatrix}}^{\mathsf{T}}(M_0\otimes I) {\begin{bmatrix}\tilde y^k \\ \tilde u^k\end{bmatrix}}
    + {\begin{bmatrix}\tilde z^k \\ \tilde v^k\end{bmatrix}}^{\mathsf{T}}\bigl( M_1\otimes (I-\Pi)\otimes R\bigr) {\begin{bmatrix}\tilde z^k \\ \tilde v^k\end{bmatrix}} \le 0,$$ where $V^k$ is defined in . The quadratic forms in the last two terms are nonnegative from Propositions \[prop:functions\] and \[prop:graph\], which implies $V^{k+1}\le\rho^2\,V^k$. We then apply this inequality iteratively to obtain $V^k \le \rho^{2k}\,V_0$ for all $k\ge 0$. Now define $$T {\colonequals}\Pi\otimes P + (I-\Pi)\otimes Q,$$ and note that $T\succ 0$ since $P$ and $Q$ are positive definite. Then letting $\operatorname{cond}(T) = \lambda_\text{max}(T)/\lambda_\text{min}(T)$ denote the condition number of $T$, we have the bound $$\begin{aligned}
    \|x_i^k\!-\!x_i^\star\|^2 &\le \|x^k\!-\!x^\star\|^2
    \le \operatorname{cond}(T)\,V^k
    \le \rho^{2k}\operatorname{cond}(T)\,V^0.\end{aligned}$$ Therefore, the bound  holds with $c = \sqrt{\operatorname{cond}(T)\,V^0}$. 

Proof of Theorem \[thm:SVL\] {#sec:SVL-proof}
----------------------------

[Substituting the template  into the LMI  reduces to $$P_{11} {\begin{bmatrix}1-\rho^2 & -\alpha \\ -\alpha & \alpha^2\end{bmatrix}} + M_0 \preceq 0,$$ which is satisfied with $\alpha = (1-\rho)/m$ and $P_{11}=\tfrac{m\,(L-m)}{\rho\,(1-\rho)}$. Note that this LMI is known to describe the convergence rate of centralized gradient descent; see [@lessard16 Section 4.4].]{}

Now consider the potential solution [to  given by]{} $$\begin{aligned}
  {Q} &= \frac{t_3}{\alpha^2\rho^2} {\begin{bmatrix}1+\rho^2\frac{t_1}{t_4} & -1 \\ -1 & 1\end{bmatrix}}
  \quad\text{and}\quad
  {R} = \frac{t_5}{\alpha^2 t_2},\quad\text{where} \hspace{-5cm}\\
    t_1 &{\colonequals}2\,(1-\beta)-\eta, &
    t_2 &{\colonequals}\beta-1+\rho^2, \\
    t_3 &{\colonequals}\beta\,(\eta+2\rho^2)-\eta\,(1-\rho^2), &
    t_4 &{\colonequals}2\beta\rho^2 - \eta\,(1-\rho^2), \\
    t_5 &{\colonequals}(1-\beta-\rho)\bigl(\beta\,(\eta+2\rho^2)-(1-\rho^2) (1-\kappa+2\kappa\rho)\bigr),\hspace{-5cm}\\
    t_6 &{\colonequals}\bigl(2-\alpha(L+m)\bigr)(1\!-\!\rho^2)^2 \!-\! \bigl(2(1\!-\!\rho^4)-\alpha(L+m)\bigr)\beta.\hspace{-5cm}\end{aligned}$$ Using these values along with the value for $\sigma^2$ in , the matrix in  is equal to the rank-one matrix $-\tfrac{1}{t_2 t_4} z z^{\mathsf{T}}$, where $$\begin{aligned}
    z {\colonequals}\frac{1}{\alpha\rho} {\begin{bmatrix}
        t_6 \\
        -t_2 t_3 \\
        \alpha\,t_2\,\bigl(2-\alpha\,(L+m)\bigr) \\
        \beta\,\bigl(t_3-\alpha\rho^2(L+m)\bigr)\end{bmatrix}}.\end{aligned}$$ In order for this to be a valid solution, we must have $t_3 > 0$ and $t_1/t_4 > 0$ (so that ${Q}\succ 0$), $t_5/t_2\ge 0$ (so that ${R\succeq 0}$), and $t_2 t_4>0$ (so that  holds). All of these inequalities hold if and only if  holds. Therefore, the SDP has a rank-one solution using the parameters in  if $\beta$ and $\rho$ satisfy . The convergence bound then follows from Theorem \[thm:main-result\] [and Remark \[rem:bounds\].]{}
