Distributed Learning in Non-Convex Environments -- Part I: Agreement at a Linear Rate

TL;DR

This paper analyzes the convergence behavior of distributed stochastic gradient algorithms in non-convex environments, showing that agents' estimates cluster around a network centroid at a linear rate.

Contribution

It extends existing non-convex distributed optimization analysis to stochastic settings, demonstrating clustering and convergence properties under gradient noise.

Findings

Agents' estimates cluster around the network centroid.

The diffusion strategy converges at a linear rate.

The approach handles persistent gradient noise in non-convex optimization.

Abstract

Driven by the need to solve increasingly complex optimization problems in signal processing and machine learning, there has been increasing interest in understanding the behavior of gradient-descent algorithms in non-convex environments. Most available works on distributed non-convex optimization problems focus on the deterministic setting where exact gradients are available at each agent. In this work and its Part II, we consider stochastic cost functions, where exact gradients are replaced by stochastic approximations and the resulting gradient noise persistently seeps into the dynamics of the algorithm. We establish that the diffusion learning strategy continues to yield meaningful estimates non-convex scenarios in the sense that the iterates by the individual agents will cluster in a small region around the network centroid. We use this insight to motivate a short-term model for network evolution over a finite-horizon. In Part II [2] of this work, we leverage this model to establish descent of the diffusion strategy through saddle points in O(1/$\mu$) steps and the return of approximately second-order stationary points in a polynomial number of iterations.

Tables

Table 1. TABLE I: Comparison of modeling assumptions and results for gradient-based methods. Statements marked with ⋆ are not explicitly stated but are implied by other conditions. The works marked with † establish global (asymptotic) convergence, which of course implies escape from saddle-points.

	Modeling conditions						Results
	Gradient	Hessian	Initialization	Perturbations	Step-size		Stationary	Saddle
Centralized
[20]	Lipschitz	—	—	SGD + Annealing	diminishing		$✓$	asymptotic†
[21]	Lipschitz & bounded⋆	Lipschitz	—	i.i.d. and bounded w.p. 1	constant		$✓$	polynomial
[22]	Lipschitz	—	Random	—	constant		$✓$	asymptotic
[23]	Lipschitz	Lipschitz	—	Selective & bounded w.p. 1	constant		$✓$	polynomial
[24]	Lipschitz	Lipschitz	—	SGD, bounded w.p. 1	alternating		$✓$	polynomial
[25]	Lipschitz	Lipschitz	—	Bounded variance, Lipschitz w.p. 1	constant		$✓$	polynomial
[26]	Lipschitz	Lipschitz	—	Bounded variance, Lipschitz w.p. 1	constant		$✓$	polynomial
[27]	Lipschitz	Lipschitz	—	Bounded variance, Lipschitz w.p. 1	constant		$✓$	polynomial
[28]	Lipschitz	Lipschitz	—	SGD, bounded w.p. 1	constant		$✓$	polynomial
[29]	Lipschitz	Lipschitz	—	SGD + Gaussian	constant		$✓$	polynomial
Decentralized
[30]	Lipschitz & bounded	—	—	—	constant		$✓$	—
[31]	Lipschitz	—	—	—	constant		$✓$	—
[32]	Lipschitz & bounded	—	—	i.i.d.	diminishing		$✓$	—
[33]	Lipschitz	Exists	Random	—	constant		$✓$	asymptotic
[34]	Bounded disagreement	—	—	SGD + Annealing	diminishing		$✓$	asymptotic†
This work	Bounded disagreement	Lipschitz	—	Bounded moments	constant		$✓$	polynomial

Equations

$$\underset{w}{\min }J(w),\mathrm{where}J(w)\triangleq \sum _{k=1}^N{p}_k{J}_k(w)$$

$${\mathbf{\varphi }}_{k,i}$$

$${\text{w}}_{k,i}$$

$$a_{\ell k}\ge 0,\sum _{\ell \in {\mathcal{N}}_k}{a}_{\ell k}=1,a_{\ell k}=0\mathrm{if}\ell \notin {\mathcal{N}}_k$$

$$Ap=p,\text{mathds}{1}^{\text{T}}p=1,p_k>0$$

$$\Vert \nabla J_k(x)-\nabla J_k(y)\Vert \le \delta \Vert x-y\Vert$$

$$\Vert \nabla J(x)-\nabla J(y)\Vert \le \delta \Vert x-y\Vert$$

$$J(y)\le J(x)+{\nabla J(x)}^{\text{T}}\left(y-x\right)+\frac{\delta }{2}{\Vert x-y\Vert }^2$$

$$-\delta I\le {\nabla }^{2}J(x)\le \delta I$$

$$\Vert \nabla J_k(x)-\nabla J_{\ell }(x)\Vert \le G$$

$$\Vert \nabla J_k(x)-\nabla J(x)\Vert =\Vert \sum _{\ell =1}^N{p}_{\ell }\left(\nabla J_k(x)-\nabla J_{\ell }(x)\right)\Vert$$

$$\le$$

$${\mathcal{F}}_i\triangleq \left\{{\text{bcw}}_0,{\text{bcw}}_1,\dots ,{\text{bcw}}_i\right\}$$

$$\text{E}\left\{\text{w}|\mathcal{F}\right\}=\text{w}$$

$$\text{E}\left\{\text{E}\left\{\text{x}|\mathcal{F}\right\}|\text{w}\in \mathcal{S}\right\}=\text{E}\left\{\text{x}|\text{w}\in \mathcal{S}\right\}$$

$$\text{mathds}{I}_{\mathcal{S}}(\text{w})=\{\begin{array}{l}1,\mathrm{if}\text{w}\in \mathcal{S}\\ 0,\mathrm{otherwise}.\end{array}$$

$$\text{E}\left\{\text{x}\text{mathds}{I}_{\mathcal{S}}(\text{w})\right\}=$$

$$+\text{E}\left\{\text{x}\text{mathds}{I}_{\mathcal{S}}(\text{w})|\text{w}\notin \mathcal{S}\right\}\cdot \mathrm{Pr}\left\{\text{w}\notin \mathcal{S}\right\}$$

$$=$$

$$\text{E}\left\{\text{x}|\text{w}\in \mathcal{S}\right\}=\frac{\text{E}\left\{\text{x}\text{mathds}{I}_{\mathcal{S}}(\text{w})\right\}}{\mathrm{Pr}\left\{\text{w}\in \mathcal{S}\right\}}$$

$$\text{E}\left\{\text{E}\left\{\text{x}|\mathcal{F}\right\}\text{mathds}{I}_{\mathcal{S}}(\text{w})\right\}$$

$$=$$

$$+\text{E}\left\{\text{E}\left\{\text{x}|\mathcal{F}\right\}\text{mathds}{I}_{\mathcal{S}}(\text{w})|\text{w}\notin \mathcal{S}\right\}\cdot \mathrm{Pr}\left\{\text{w}\notin \mathcal{S}\right\}$$

$$=$$

$$\text{E}\left\{\text{E}\left\{\text{x}|\mathcal{F}\right\}|\text{w}\in \mathcal{S}\right\}\stackrel{\text{eqref}eq:equi{v}_{e}xp2}{=}\frac{\text{E}\left\{\text{E}\left\{\text{x}|\mathcal{F}\right\}\text{mathds}{I}_{\mathcal{S}}(\text{w})\right\}}{\mathrm{Pr}\left\{\text{w}\in \mathcal{S}\right\}}$$

$$\stackrel{(a)}{=}$$

$$\stackrel{\text{eqref}eq:equi{v}_{e}xp}{=}$$

$${\text{s}}_{k,i}({\text{w}}_{k,i-1})={\widehat{\nabla J}}_k({\text{w}}_{k,i-1})-\nabla J_k({\text{w}}_{k,i-1})$$

$$\text{E}\left\{{\text{s}}_{k,i}({\text{w}}_{k,i-1})|{\mathcal{F}}_{i-1}\right\}$$

$$\text{E}\left\{\Vert {\text{s}}_{k,i}({\text{w}}_{k,i-1}){\Vert }^4|{\mathcal{F}}_{i-1}\right\}$$

$$\text{E}\left\{{\text{s}}_{k,i}({\text{w}}_{k,i-1}){\text{s}}_{\ell ,i}({\text{w}}_{\ell ,i-1}{)}^{\text{T}}|{\mathcal{F}}_{i-1}\right\}=0$$

$$\text{E}\left\{\Vert {\text{s}}_{k,i}({\text{w}}_{k,i-1}){\Vert }^2|{\mathcal{F}}_{i-1}\right\}$$

$$=$$

$$\le$$

$$J_k(w)=\frac{1}{S}\sum _{s=1}^{S}Q(w,x_{k,s})$$

$$w^o\triangleq \arg \underset{w}{\min }\frac{1}{N}\sum _{k=1}^N\left(\frac{1}{S}\sum _{s=1}^{S}Q(w,x_{k,s})\right)$$

Full text

I Introduction

The broad objective of distributed adaptation and learning is the solution of global, stochastic optimization problems by networked agents through localized interactions and in the absence of information about the statistical properties of the data. When constant, rather than diminishing, step-sizes are employed, the resulting algorithms are adaptive in nature and are able to adapt to drifts in the data statistics. In this work, we consider a collection of $N$ agents, where each agent $k$ is equipped with a stochastic risk of the form $J_{k}(w)=\E_{x}Q_{k}(w;\x_{k})$ with $Q_{k}(w;\x_{k})$ referring to the loss function, $w\in\mathds{R}^{M}$ denoting a parameter vector, and $\x_{k}$ referring to the stochastic data. The expectation is over the probability distribution of the data. The objective of the network is to seek the Pareto solution:

[TABLE]

where the $p_{k}$ are positive weights that are normalized to add up to one and will be specified further below; in particular, in the special case when the $\{p_{k}\}$ are identical, they can be removed from (1). Algorithms for the solution of (1) have been studied extensively over recent years both with inexact [3, 4, 5, 6] and exact [7, 8, 9] gradients. Here, we focus on the following diffusion strategy, which has been shown in previous works to provide enhanced performance and stability guarantees under constant step-size learning and adaptive scenarios [10, 4]:

[TABLE]

where $\widehat{\nabla J}_{k}(\cdot)$ denotes a stochastic approximation for the true local gradient ${\nabla J}_{k}(\cdot)$. The intermediate estimate $\boldsymbol{\phi}_{k,i}$ is obtained at agent $k$ by taking a stochastic gradient update relative to the local cost $J_{k}(\cdot)$. The intermediate estimates are then fused across local neighborhoods where $a_{\ell k}$ are convex combination weights satisfying:

[TABLE]

The symbol ${\cal N}_{k}$ denotes the set of neighbors of agent $k$. {assumption}[Strongly-connected graph] We shall assume that the graph described by the weighted combination matrix $A=[a_{\ell k}]$ is strongly-connected [4]. This means that there exists a path with nonzero weights between any two agents in the network and, moreover, at least one agent has a nontrivial self-loop, $a_{kk}>0$. ∎

It then follows from the Perron-Frobenius theorem [11, 12, 4] that $A$ has a single eigenvalue at one while all other eigenvalues are strictly inside the unit circle, so that $\rho(A)=1$. Moreover, if we let $p$ denote the right-eigenvector of $A$ that is associated with the eigenvalue at one, and if we normalize the entries of $p$ to add up to one, then it also holds that all entries of $p$ are strictly positive, i.e.,

[TABLE]

where the $\{p_{k}\}$ denote the individual entries of the Perron vector, $p$.

I-A Related Works

The performance of the diffusion algorithm (2a)–(2b) has been studied extensively in differentiable settings [10, 5], with extensions to multi-task [13], constrained [14], and non-differentiable [15] environments. A common assumption in these works, along with others studying the behavior of distributed optimization algorithms in general, is that of convexity (or strong-convexity) of the aggregate risk $J(w)$. While many problems of interest such as least-squares estimation [4], logistic regression [4], and support vector machines [16] are convex, there has been increased interest in the optimization of non-convex cost functions. Such problems appear frequently in the design of robust estimators [17] and the training of more complex machine learning architectures such as those involving dictionary learning [18] and artificial neural networks [19].

Motivated by these applications, recent works have pursued the study of optimization algorithms for non-convex problems, both in the centralized and distributed settings [20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37]. While some works focus on establishing convergence to a stationary point [30, 31], there has been growing interest in examining the ability of gradient descent implementations to escape from saddle points, since such points represent bottlenecks to the underlying learning problem [19]. We defer a detailed discussion on the plethora of related works on second-order guarantees [38, 20, 21, 22, 23, 33, 24, 25, 27, 26, 28, 29, 34] to Part II [2], where we will be establishing the ability of the diffusion strategy (2a)–(2b) to escape strict-saddle points efficiently. For ease of reference, the modeling conditions and results from this and related works are summarized in Table I.

The key contributions of Parts I and II this work are three-fold. To the best of our knowledge, we present the first analysis establishing efficient (i.e., polynomial) escape from strict-staddle points in the distributed setting. Second, we establish that the gradient noise process is sufficient to ensure efficient escape without the need to alter it by adding artificial forms of perturbations, interlacing steps with small and large step-sizes, or imposing a dispersive noise assumption. Third, relative to the existing literature on centralized non-convex optimization, where the focus is mostly on deterministic or finite-sum optimization, our modeling conditions are specifically tailored to the scenario of learning from stochastic streaming data. In particular, we only impose bounds on the gradient noise variance in expectation, rather than assume a bound with probability one [24, 28] or a sub-Gaussian distribution [29]. Furthermore, we assume that any Lipschitz conditions only hold on the expected stochastic gradient approximation, rather than for every realization, with probability one [25, 27, 26].

I-B Preview of Results

We first establish that in non-convex environments, as was already shown earlier in [5] for convex environments, the evolution of the individual iterates $\w_{k,i}$ at the agents continues to be well-described by the evolution of the weighted centroid vector $\sum_{k=1}^{N}p_{k}\w_{k,i}$ in the sense that the iterates from across the network will cluster around this centroid after sufficient iterations. We subsequently consider two cases separately and establish descent in both of them. The first case corresponds to the region where the gradient at the network centroid is large and establish that descent can occur in one iteration. The second and more challenging case occurs when the gradient norm is small, but there is a sufficiently negative eigenvalue in the Hessian matrix. We establish Part II [2] that the recursion will continue to descend along the aggregate cost at a rate of $O(\mu)$ per $O(1/\mu)$ iterations. Combined with the first result, this descent relation allows us to provide guarantees about the second-order optimality of the returned iterates.

The flow of the argument is summarized in Fig. 1. We decompose $\mathds{R}^{M}$ into the set of approximate first-order stationary points, i.e., those with $\|\nabla J(w)\|^{2}\leq O(\mu)$ and the complement, i.e., the large-gradient regime. For the large-gradient regime, descent is established in Theorem II-C. We proceed to further decompose the set of approximate first-order stationary points into those that are $\tau$-strict-saddle, i.e., those that have a Hessian with significant negative eigenvalue $\lambda_{\min}\left(\nabla^{2}J(w)\right)\leq-\tau$, and the complement, which are approximately second-order stationary points. For $\tau$-strict-saddle points we establish descent in Part II [2, Theorem 1]. Finally, in Part II [2, Theorem 2], we conclude that the centroid will reach an approximately second-order stationary point in a polynomial number of iterations.

II Evolution Analysis

We shall perform the analysis under the following common assumptions on the gradients and their approximations. {assumption}[Lipschitz gradients] For each $k$, the gradient $\nabla J_{k}(\cdot)$ is Lipschitz, namely, for any $x,y\in\mathds{R}^{M}$:

[TABLE]

In light of (1) and Jensen’s inequality, this implies for the aggregate cost:

[TABLE]

∎

The Lipschitz gradient conditions (5) and (6) imply bounds on the both the function value and the Hessian matrix (when it exists), which will be used regularly throughout the derivations. In particular, we have for the function values:

[TABLE]

For the Hessian matrix we have [4]:

[TABLE]

{assumption}

[Bounded gradient disagreement] For each pair of agents $k$ and $\ell$, the gradient disagreement is bounded, namely, for any $x\in\mathds{R}^{M}$:

[TABLE]

∎

This assumption is similar to the one used in [34] under a diminishing step-size with annealing. Note that condition (9) is weaker than the more common assumption of bounded gradients. Condition (9) is automatically satisfied in cases where the expected risks $J_{k}(\cdot)$ are common (though agents still may see different realizations of data), or in the case of centralized stochastic gradient descent where the number of agents is one. This condition is also satisfied whenever agent-specific risks with bounded gradients are regularized by common regularizers with potentially unbounded gradients, as is common in many machine learning applications. Observe that (9) implies a similar condition on the deviation from the centralized gradient via Jensen’s inequality:

[TABLE]

{definition}

[Filtration] We denote by $\boldsymbol{\mathcal{F}}_{i}$ the filtration generated by the random processes $\w_{k,j}$ for all $k$ and $j\leq i$:

[TABLE]

where $\bcw_{j}\triangleq\mathrm{col}\left\{\w_{1,j},\ldots,\w_{k,j}\right\}$ contains the iterates across the network at time $j$. Informally, $\boldsymbol{\mathcal{F}}_{i}$ captures all information that is available about the stochastic processes $\w_{k,j}$ across the network up to time $i$. ∎

Throughout the following derivations, we will frequently rely on appropriate conditionings to make the analysis tractable. A frequent theme will be the exchange of conditioning on filtrations by conditioning on events. To this end, the following lemma will be used repeatedly. {lemma}[Conditioning] Suppose $\w\in\mathds{R}^{M}$ is a random variable measurable by $\boldsymbol{\mathcal{F}}$. In other words, $\w$ is deterministic conditioned on $\boldsymbol{\mathcal{F}}$ and

[TABLE]

Then,

[TABLE]

for any deterministic set $\mathcal{S}\subseteq\mathds{R}^{M}$ and random $\x\in\mathds{R}^{M}$.

Proof:

Denote by $\mathds{I}_{\mathcal{S}}(\w)$ the random indicator function:

[TABLE]

Since $\w$ is measurable by $\boldsymbol{\mathcal{F}}$, then $\mathds{I}_{\mathcal{S}}(\w)$ is measurable by $\boldsymbol{\mathcal{F}}$ as well. In other words, the event $\w\in\mathcal{S}$ is deterministic conditioned on $\boldsymbol{\mathcal{F}}$. Furthermore, for the random variable $\x\mathds{I}_{\mathcal{S}}(\w)$, we have:

[TABLE]

Rearranging yields:

[TABLE]

Similarly, for the random variable $\E\left\{\x|\boldsymbol{\mathcal{F}}\right\}\mathds{I}_{\mathcal{S}}(\w)$, we have:

[TABLE]

It then follows that:

[TABLE]

where in step $(a)$ we pulled $\mathds{I}_{\mathcal{S}}(\w)$ into the inner expectation, since it is deterministic conditioned on $\boldsymbol{\mathcal{F}}$ and $(b)$ follows from the law of total expectation. ∎

{assumption}

[Gradient noise process] For each $k$, the gradient noise process is defined as

[TABLE]

and satisfies

[TABLE]

for some non-negative constant $\sigma^{4}$. We also assume that the gradient noise pocesses are pairwise uncorrelated over the space conditioned on $\boldsymbol{\mathcal{F}}_{i-1}$, i.e.:

[TABLE]

∎ Property (20a) means that the gradient noise construction is unbiased on average. Property (20b) means that the fourth-moment of the gradient noise is bounded. These properties are automatically satisfied for several costs of interest [4, 10]. Note, that the bound on the fourth-order moment, in light of Jensen’s intequality, immediately implies:

[TABLE]

While our primary interest is in the development of algorithms that allow for learning from streaming data, we remark briefly that the results obtained in this work are equally applicable to empirical risk minimization via stochastic gradient descent, by assuming that the streaming data is selected according to a particular distribution.

{example}

[Empirical Risk Minimization] Suppose the costs $J_{k}(\cdot)$ are empirical based on locally collected data ${\left\{x_{k,s}\right\}}_{s=1}^{S}$ and take the form:

[TABLE]

In empirical risk minimization (ERM) problems, we are interested in finding a vector $w^{o}$ that minimizes the following empirical risk over the data across the entire network:

[TABLE]

If we introduce the uniformly-distributed random variable $\x_{k}=x_{k,s}$ with probability $\frac{1}{S}$ for all $s$, then the cost (24) is equivalent to solving:

[TABLE]

which is of the same form as (1) with $p_{k}=\frac{1}{N}$. The resulting gradient noise process satisfies the assumptions imposed in this work under appropriate conditions on the risk $Q(\cdot,\cdot)$. This observation has been leveraged to accurately quantify the performance of stochastic gradient descent, as well as mini-batch and importance sampling generalizations, for emprical minimization of convex risks in [7]. ∎

II-A Network basis transformation

In analyzing the dynamics of the distributed algorithm (2a)–(2b), it is useful to introduce the following extended quantities by collecting variables from across the network:

[TABLE]

where $\otimes$ denotes the Kronecker product operation. We can then write the diffusion recursion (2a)–(2b) compactly as

[TABLE]

By construction, the combination matrix $A$ is left-stochastic and primitive and hence admits a Jordan decomposition of the form $A=V_{\epsilon}JV_{\epsilon}^{-1}$ with [4, 5]:

[TABLE]

where $J_{\epsilon}$ is a block Jordan matrix with the eigenvalues $\lambda_{2}(A)$ through $\lambda_{N}(A)$ on the diagonal and $\epsilon$ on the first lower sub-diagonal. The extended matrix $\mathcal{A}$ then satisfies $\mathcal{A}=\mathcal{V}_{\epsilon}\mathcal{J}\mathcal{V}_{\epsilon}^{-1}$ with $\mathcal{V}_{\epsilon}={V}_{\epsilon}\otimes I_{N}$, $\mathcal{J}=J\otimes I_{N}$, $\mathcal{V}_{\epsilon}^{-1}={V}_{\epsilon}^{-1}\otimes I_{N}$. The spectral properties of $A$ and its corresponding eigendecomposition have been exploited extensively in the study of the diffusion learning strategy in the convex setting [4, 5], and will continue to be useful in non-convex scenarios.

Multiplying both sides of (29) by $\left(p^{\T}\otimes I\right)$ from the left, we obtain in light of (4):

[TABLE]

Letting $\w_{c,i}\triangleq\sum_{k=1}^{K}p_{k}\w_{k,i}=\left(p^{\T}\otimes I\right)\bcw_{i}$ and exploiting the block-structure of the gradient term, we find:

[TABLE]

Note that $\w_{c,i}$ is a convex combination of iterates across the network and can be viewed as a weighted centroid. The recursion for $\w_{c,i}$ is reminiscent of a stochastic gradient step associated with the aggregate cost $\sum_{k=1}^{N}p_{k}J_{k}(w)$ with the exact gradients $\nabla J_{k}(\cdot)$ replaced by stochastic approximations $\widehat{\nabla J}_{k}(\cdot)$ and with the stochastic gradients evaluated at $\w_{k,i-1}$, rather than $\w_{c,i-1}$. In fact, we can write:

[TABLE]

where we defined the perturbation terms:

[TABLE]

We use the subscript $i-1$ for $\boldsymbol{d}_{i-1}$ to emphasize that it depends on data up to time $i-1$, in contrast to $\s_{i}$ which is also dependent on the most recent data from time $i$. Observe that $\boldsymbol{d}_{i-1}$ arises from the disagreement within the network, and in particular that if each $\w_{k,i-1}$ remains close to the network centroid $\w_{c,i-1}$, this perturbation will be small in light of the Lipschitz condition (5) on the gradients. The second perturbation term $\s_{i}$ arises from the noise introduced by stochastic gradient approximations at each agent. We now establish that recursion (38) will continue to exhibit some of the desired properties of (centralized) gradient descent, despite the presence of persistent and coupled perturbation terms.

II-B Network disagreement

To begin with, we study more closely the evolution of the individual estimates $\w_{k,i}$ relative to the network centroid $\w_{c,i}$. Multiplying (29) by $\mathcal{V}_{R}^{\T}\triangleq\left({V}_{R}^{\T}\otimes I\right)$ from the left yields in light of (35):

[TABLE]

Then, for the deviation from the network centroid:

[TABLE]

so that the deviation from the centroid can be easily recovered from $\mathcal{V}_{R}^{\T}\bcw_{i}$ in (II-B). Proceeding with (II-B), we find:

[TABLE]

where $(a)$ follows from the sub-multiplicative property of norms, and $(b)$ follows from Jensen’s inequality ${\|a+b\|}^{4}\leq\frac{1}{\alpha^{3}}\|a\|^{4}+\frac{1}{\left(1-\alpha\right)^{3}}\|b\|^{4}$ with

[TABLE]

for sufficiently small $\epsilon$ due to Assumption 3, where $\lambda_{2}\triangleq\rho\left(A-\mathds{1}p^{\T}\right)$. We observe that the term ${\left\|\mathcal{V}_{R}^{\T}\bcw_{i}\right\|}^{4}$ contracts at an exponential rate given by $\left\|J_{\epsilon}^{\T}\right\|\approx\lambda_{2}$ for small $\epsilon$, also known as the mixing rate of the graph. Iterating this relation and applying Assumptions 3–II, we obtain the following result.

{theorem}

[Network disagreement (4th order)] Under assumptions 3–II, the network disagreement is bounded after sufficient iterations $i\geq i_{o}$ by:

[TABLE]

where

[TABLE]

and $o(\mu^{4})$ denotes a term that is higher in order than $\mu^{4}$.

Proof:

Appendix A. ∎

Note again, that Jensen’s inequality immediately implies for the second-order moment:

[TABLE]

where $(a)$ follows from (45) and sub-additivity of the square root, i.e. $\sqrt{x+y}\leq\sqrt{x}+\sqrt{y}$. This result establishes that, for every agent $k$, we have after sufficient iterations $i\geq i_{o}$:

[TABLE]

or, by Markov’s inequality [39]:

[TABLE]

and hence $\w_{k,i}$ will be arbitrarily close to $\w_{c,i}$ with arbitrarily high probability for all agents. This result has two implications. First, it allows us to use the network centroid $\w_{c,i}$ as a proxy for all iterates $\w_{k,i}$ in the network, since all agents will cluster around the network centroid after sufficient iterations. Second, it allows us to bound the perturbation terms encountered in (38). {lemma}[Perturbation bounds (2nd and 4th order)] Under assumptions 3–II and for sufficiently small step-sizes $\mu$, the perturbation terms are bounded as:

[TABLE]

after sufficient iterations $i\geq i_{0}$.

Proof:

Appendix B. ∎

{definition}

[Sets] To simplify the notation in the sequel, we introduce following sets:

[TABLE]

where $\tau$ is a small positive parameterm, $c_{1}$ and $c_{2}$ are constants:

[TABLE]

and $0<\pi<1$ is a parameter to be chosen. Note that $\mathcal{G}^{C}=\mathcal{H}\cup\mathcal{M}$. We also define the probabilities $\pi^{\mathcal{G}}_{i}\triangleq\mathrm{Pr}\left\{\w_{c,i}\in\mathcal{G}\right\}$, $\pi^{\mathcal{H}}_{i}\triangleq\mathrm{Pr}\left\{\w_{c,i}\in\mathcal{H}\right\}$ and $\pi^{\mathcal{M}}_{i}\triangleq\mathrm{Pr}\left\{\w_{c,i}\in\mathcal{M}\right\}$. Then for all $i$, we have $\pi^{\mathcal{G}}_{i}+\pi^{\mathcal{H}}_{i}+\pi^{\mathcal{M}}_{i}=1$. ∎ The definitions (52)–(55) decompose the parameter-space $\mathds{R}^{M}$ into two disjoint sets $\mathcal{G}$ and $\mathcal{G}^{C}$, and further sub-divides $\mathcal{G}^{C}$ into $\mathcal{H}$ and $\mathcal{M}$. The set $\mathcal{G}$ denotes the set all points $w$ where the norm of the gradient is large, while $\mathcal{G}^{C}=\mathcal{H}\cup\mathcal{M}$ denotes the set of all points where the norm of the gradient is small, i.e., approximately first-order stationary points. In a manner similar to related works on the escape from strict-saddle points, we further decompose the set $\mathcal{G}^{C}$ of approximate first-order stationary points into those points $w\in\mathcal{H}$ that do have a significant negative eigenvalue, and those in $\mathcal{M}$ that do not [21, 23]. Points in the parameter space that have a small gradient norm and no significant negative eigenvalue are referred to as second-order stationary points, while points in $\mathcal{H}$ are known as strict saddle-points due to the presence of a strictly negative eigenvalue in the Hessian matrix. In the sequel, we will establish descent for centroids in $\mathcal{G}$ in Theorem II-C and centroids in $\mathcal{H}$ in Part II [2, Theorem 1], and hence the approach of a point in $\mathcal{M}$ with high probability after a polynomial number of iterations in Part II [2, Theorem 2]. Second-order stationary points are generally more likely to be “good” minimizers than first-order stationary points, which could even correspond to local maxima. Furthermore, for a certain class of cost functions, known as “strict-saddle” functions, second-order stationary points always correspond to local minimia for sufficiently small $\tau$ [21].

II-C Evolution of the network centroid

Having established in (47), that after sufficient iterations, all agents in the network will have contracted around the centroid in a small cluster for small step-sizes, we can now leverage $\w_{c,i}$ as a proxy for all $\w_{k,i}$. From Assumption II and (7), we have the following bound:

[TABLE]

From (38), we then obtain:

[TABLE]

This relation, along with (38) and the results from Lemma 49, allow us to establish the following theorem. {theorem}[Descent relation] Beginning at $\w_{c,i-1}$ in the large gradient regime $\mathcal{G}$, we can bound:

[TABLE]

as long as $\pi_{i-1}^{\mathcal{G}}=\mathrm{Pr}\left\{\w_{c,i-1}\in\mathcal{G}\right\}\neq 0$ where the relevant constants are listed in definition II-B. On the other hand, beginning at $\w_{c,i-1}\in\mathcal{M}$, we can bound:

[TABLE]

as long as $\pi_{i-1}^{\mathcal{M}}=\mathrm{Pr}\left\{\w_{c,i-1}\in\mathcal{M}\right\}\neq 0$.

Proof:

Appendix D. ∎

Relation (II-C) guarantees a lower bound on the expected improvement when the gradient norm at the current iterate is suffiently large, i.e. $\w_{c,i-1}\in\mathcal{G}$ is not an approximately first-order stationary point. On the other hand, when $\w_{c,i-1}\in\mathcal{M}$, inequality (II-C) it establishes an upper bound on the expected ascent. The respective bounds can be balanced by appropriately choosing $\pi$, which will be leveraged in Part II [2]. We are left to treat the third possibility, namely $\w_{c,i-1}\in\mathcal{H}$. In this case, since the norm of the gradient is small, it is no longer possible to guarantee descent in a single iteration. We shall study the dynamics in more detail in the sequel.

II-D Behavior around stationary points

In the vicinity of saddle-points, the norm of the gradient is not sufficiently large to guarantee descent at every iteration as indicated by (II-C). Instead, we will study the cumulative effect of the gradient, as well as perturbations, over several iterations. For this purpose, we introduce the following second-order condition on the cost functions, which is common in the literature [4, 21, 23]. {assumption}[Lipschitz Hessians] Each $J_{k}(\cdot)$ is twice-differentiable with Hessian $\nabla^{2}J_{k}(\cdot)$ and, there exists $\rho\geq 0$ such that:

[TABLE]

By Jensen’s inequality, this implies that $J(\cdot)=\sum_{k=1}^{N}p_{k}J_{k}(\cdot)$ also satisfies:

[TABLE]

∎

Let $i^{\star}$ denote an arbitraty point in time. We use $i^{\star}$ in order to emphasize approximately first-order stationary points, where the norm of the gradient is small. Such first-order stationary points $\w_{c,i^{\star}}\in\mathcal{G}^{C}$ could either be in the set of second-order stationary points $\mathcal{M}$ or in the set of strict-saddle points $\mathcal{H}$. Our objective is to show that when $\w_{c,i^{\star}}\in\mathcal{H}$, we can guarantee descent after several iterations. To this end, starting at $i^{\star}$, we have for $i\geq 0$:

[TABLE]

Subsequent analysis will rely on an auxilliary model, referred to as a short-term model. It will be seen that this model is more tractable and evolves “close” to the true recursion under the second-order smoothness condition on the Hessian matrix (63) and as long as the iterates remain close to a stationary point. A similar approach has been introduced and used to great advantage in the form of a “long-term model” to derive accurate mean-square deviation performance expressions for strongly-convex costs in [4, 10, 40, 41]. The approach was also used to provide a “quadratic approximation” to establish the ability of stochastic gradient based algorithms to escape from strict saddle-points in the single-agent case under i.i.d. perturbations in [21].

For the driving gradient term in (64), we have from the mean-value theorem [4]:

[TABLE]

where

[TABLE]

Subtracting (64) from $\w_{c,i^{\star}}$, we obtain:

[TABLE]

We introduce short-hand notation for the deviation:

[TABLE]

Note that $\widetilde{\w}_{i}^{i^{\star}}$ denotes the deviation of the network centroid $\w_{c,i^{\star}+i}$ at time $i^{\star}+i$ from the initial, approximately first-order stationary point $\w_{c,i^{\star}}$. Establishing escape from saddle-points is equivalent to establishing the growth of $\widetilde{\w}_{i}^{i^{\star}}$ whenever $\w_{c,i^{\star}}\in\mathcal{H}$. We hence expect the deviation to grow over time, but would like to establish that $\w_{c,i^{\star}+i}$ moves away from $\w_{c,i^{\star}}$ in a direction of descent. We can then write more compactly:

[TABLE]

The time-varying nature of $\boldsymbol{H}_{i^{\star}+i}$ makes this recursion difficult to study. We hence introduce the following auxilliary recursion, initialized at $\w_{c,i^{\star}}^{\prime}=\w_{c,i^{\star}}$, where $\boldsymbol{H}_{i^{\star}+i}$ is replaced by $\nabla^{2}J(\w_{c,i^{\star}})$ and the perturbation term $\mu\boldsymbol{d}_{i^{\star}+i}$ is omitted:

[TABLE]

or, more compactly, with $\widetilde{\w}_{i}^{\prime}{}^{i^{\star}}\triangleq\w_{c,i^{\star}}-\w_{c,i^{\star}+i}^{\prime}$

[TABLE]

Of course, this second model is only useful in studying the behavior of the original recursion (64) if the iterates generated by both models remain close to each other, which we shall prove to be true. Specifically, if we write:

[TABLE]

then $\boldsymbol{u}_{i^{\star}+i+1}$ will be shown to be negligible in some sense. Results along this line have been established in the centralized and distributed contexts for strongly-convex costs [4, 10] and in the centralized setting for strict saddle points [21]. We show here that this conclusion holds more generally in the vicinity of $O(\mu)$-first-order stationary points. Before establishing deviation bounds, we establish a short lemma which will be used repeatedly. {lemma}[A limiting result] For $T,\mu,\delta>0$ and $k\in\mathds{Z}_{+}$ with $\mu<\frac{1}{\delta}$, we have:

[TABLE]

Proof:

Appendix C. ∎

{lemma}

[Deviation bounds] Suppose $\mathrm{Pr}\left\{\w_{c,i^{\star}}\in\mathcal{H}\right\}\neq 0$. Then, the following quantities are conditionally bounded:

[TABLE]

for $i\leq\frac{T}{\mu}$, where $T$ denotes an arbitrary constant that is independent of the step-size $\mu$.

Proof:

Appendix E. ∎

These deviation bounds establish that, beginning at a strict-saddle point $\w_{c,i^{\star}}$ at time $i^{\star}$ the iterates will remain close to $\w_{c,i^{\star}}$ for the next $O(1/\mu)$ iterations. Consequently, the stort-term model will be sufficiently accurate for the next $O(1/\mu)$ iterations. We will establish formally in Part II [2] that the small-deviation bounds in Lemma II-D ensure descent of the true recursion can be inferred by studying only the evolution of the short-term model, which is significantly more tractable.

III Application: Robust Regression

Consider a scenario where each agent $k$ in the network observes streaming realizations $\left\{\boldsymbol{\gamma}(k,i),\boldsymbol{h}_{k,i}\right\}$ from the linear model $\boldsymbol{\gamma}(k)=\boldsymbol{h}_{k}^{\mathsf{T}}w^{o}+\boldsymbol{v}(k)$ where $\boldsymbol{\gamma}(k)$ denotes scalar observations and $\boldsymbol{v}(k)$ denotes measurement noise. One common approach for estimating $w^{o}$ in a distributed setting is via least-mean-square error estimation, resulting in the local cost functions:

[TABLE]

The resulting problem is convex and has been studied extensively in the literature. While effective under the assumption of Gaussian noise, and similar well-behaved noise conditions, this approach is susceptible to outliers caused by heavy-tailed distributions for $\boldsymbol{v}(k)$ [17]. This is caused by the fact that the quadratic risk penalizes errors proportionally to their squared norm, and as such has a tendency to over-correct outliers, even if they are rare. Several alternative robust cost functions have been suggested in the literature. We consider two in particular in order to illustrate the advantages of allowing for non-convex costs in the context of robust estimation, namely the Huber loss $Q_{k}^{\mathrm{H}}(w;\x_{k})$ and Tukey’s biweight loss $Q_{k}^{\mathrm{B}}(w;\x_{k})$ [17]. For ease of notation, let $\boldsymbol{e}(w)\triangleq\boldsymbol{\gamma}(k)-\boldsymbol{h}_{k}^{\mathsf{T}}w$. Then:

[TABLE]

where $c_{H},c_{B}$ are tuning constants. The Huber cost is merely convex (and not strongly-convex), while the Tukey loss is non-convex. Both losses satisfy assumptions 3–II imposed in this work. In particular, since the Huber risk $J_{k}^{\mathrm{H}}(w)$ has a unique, local minimum, which also happens to be locally strongly-convex, we can conclude that despite the absence of strong-convexity, the algorithm will converge to within $O(\mu)$ of the global minimum. The Tukey loss on the other hand, is non-convex, and is therefore a more challenging problem. The setting for the simulation results is shown in Figure 2.

Performance is illustrated in Fig. 3. We first show the performance of each cost in the nominal scenario, where $\boldsymbol{v}(k)\sim\mathcal{N}(0,\sigma_{v}^{2})$. We observe that the distributed strategies outperform the non-cooperative ones, and that despite differences in the rate of convergence, there is negligible difference in the performance of the mean-square-error, Huber and Tukey variations. In the presence of outliers, modeled as a bimodal distribution with $\boldsymbol{v}(k)\sim(1-\epsilon)\mathcal{N}(0,\sigma_{v}^{2})+\epsilon\mathcal{N}(10,\sigma_{v}^{2})$ and $\epsilon=0.1$, the performance of the mean-square-error solution dramatically deteriorates, as is to be expected in the presence of deviations from the nominal model.

Appendix A Proof of Lemma II-B

Starting from (II-B), taking norms of both sides and computing the fourth power, we find:

[TABLE]

where step $(a)$ follows from convexity of $\|\cdot\|^{4}$ and Jensen’s inequality, i.e. ${\left\|a+b\right\|}^{4}=\frac{1}{\alpha^{3}}{\left\|a\right\|}^{4}+\frac{1}{{(1-\alpha)}^{3}}{\left\|b\right\|}^{4}$. To begin with, we study the stochastic gradient term in some greater detail. We have:

[TABLE]

For the first term we have:

[TABLE]

where $(a)$ follows from the fact that (35) implies $V_{R}^{\T}\mathds{1}=0$, $(b)$ follows from the sub-multiplicity of norms and $(c)$ expands ${\|\cdot\|}^{2}$. For the gradient noise term we find under expectation:

[TABLE]

where $(a)$ follows from Cauchy-Schwarz, which implies ${\left(\sum_{k=1}^{N}x_{k}\right)}^{2}\leq N\sum_{k=1}^{N}x_{k}^{2}$. Plugging these relations back into (82), we obtain:

[TABLE]

We can iterate, starting from $i=0$, to obtain:

[TABLE]

where $(a)$ follows from $\sum_{n=1}^{i}\|J_{\epsilon}^{\T}\|^{n-1}\leq\sum_{n=1}^{\infty}\|J_{\epsilon}^{\T}\|^{n-1}={\left(1-\|J_{\epsilon}^{\T}\|\right)}^{-1}$, and $(b)$ holds whenever:

[TABLE]

Finally, we have from (42) under (88):

[TABLE]

where $(a)$ follows from the sub-multiplicative property of norms.We conclude that all agents in the network will contract around the centroid vector $\left(\mathds{1}p^{\T}\otimes I\right)\bcw_{i}$ after sufficient iterations.

Appendix B Proof of Lemma 49

We begin by studying the perturbation term $\boldsymbol{s}_{i}$. We have:

[TABLE]

where $(a)$ follows from $\sum_{k=1}^{N}p_{k}=1$ and Jensen’s inequality and $(b)$ follows from the fourth-order moment condition in Assumption II. For the second perturbation term, we have

[TABLE]

where $(a)$ again follows from Jensen’s inequality, $(b)$ follows from the Lipschitz gradient condition in Assumption II, and we introduced $\bcw_{c,i-1}\triangleq\mathds{1}\otimes\w_{c,i-1}$. Result (50) follows by applying (89) to (91).

Appendix C Proof of Lemma 72

For the natural logarithm of the expression, we have:

[TABLE]

Since the logarithm is continuous over $\mathds{R}_{+}$, we have:

[TABLE]

We examine the fraction inside the limit more closely. Since both the numerator and denominator of the fraction approach zero as $\mu\to 0$, we apply L’Hôpital’s rule:

[TABLE]

Hence, we find:

[TABLE]

Appendix D Proof of Theorem II-C

We begin with (II-C) and take expectations conditioned on $\bcw_{i-1}$ to obtain:

[TABLE]

where cross-terms were removed in $(a)$ due to the conditional zero-mean condition (20a), $(b)$ follows from ${\|a+b\|}^{2}\leq 2{\|a\|}^{2}+2{\|b\|}^{2}$ and from $-2a^{\T}b\leq\|a\|^{2}+\|b\|^{2}$ and $(c)$ is a result of grouping terms and Lemma 49.

Note that (96) continues to be random due to the conditioning on $\bcw_{i-1}$, but that it holds for every choice of $\bcw_{i-1}$ with probability $1$. Furthermore, since $\w_{c,i-1}=\sum_{k=1}^{N}p_{k}\w_{k,i-1}$, the centroid $\w_{c,i-1}$ is deterministic conditioned on $\bcw_{i-1}$. As such, the event $\w_{c,i-1}\in\mathcal{G}$ is deterministing conditioned on $\bcw_{i-1}$, and (96) holds for every $\w_{c,i-1}\in\mathcal{G}$. We can hence take expectations over $\w_{c,i-1}\in\mathcal{G}$ and apply Lemma II to find:

[TABLE]

In step $(a)$ we applied definition II-B, and in particular, that from (52) ${\left\|\nabla J(\w_{c,i-1})\right\|}^{2}\geq\mu\frac{c_{2}}{c_{1}}\left(1+\frac{1}{\pi}\right)$ whenever $\w_{c,i-1}\in\mathcal{G}$, which implies:

[TABLE]

We also collected constants into $c_{1}$ and $c_{2}$ defined in (56)–(57) for brevity. Step $(b)$ is obtained by grouping terms. Note that from lemma 49, we have a bound on $\E{\|\boldsymbol{d}_{i-1}\|}^{2}$, but not on the partial expectation conditioned over $\w_{c,i-1}\in\mathcal{G}$. We can decompose the full expectation:

[TABLE]

which implies

[TABLE]

so that we obtain for (D):

[TABLE]

Similarly:

[TABLE]

where $(a)$ follows from the fact that ${\left\|{\nabla J(\w_{c,i-1})}\right\|}^{2}\geq 0$ with probability $1$ and $(b)$ made use of the same argument that led to (101).

Appendix E Proof of Lemma II-D

We refer to (69). Suppose $i\leq\frac{T}{\mu}$, where $T$ is an arbitrary constant independent of $\mu$. We then have for $i\geq 0$:

[TABLE]

where $(a)$ follows from the conditional zero-mean property of the gradient noise term in Assumption II, $(b)$ follows from Jensen’s inequality

[TABLE]

with $\alpha=\mu\delta<1$ and $(c)$ follows from the same inequality with $\alpha=\frac{1}{2}$. Step $(d)$ follows from the sub-multiplicative property of norms along with $-\delta I\leq\nabla^{2}J(\w_{c,i^{\star}})\leq\delta I$, which follows from the Lipschitz gradient condition in Assumption II. Since $\w_{c,i^{\star}}$ is deterministic conditioned on $\boldsymbol{\mathcal{F}}_{i^{\star}+i}$ we can now take expectations over $\w_{c,i^{\star}}\in\mathcal{H}$ to obtain:

[TABLE]

where $(a)$ follows from the perturbation bounds in Lemma 49 and the starting assumption that $\w_{c,i^{\star}}$ is an $O(\mu)$-square stationary point. Note that, at time $i=0$, we have:

[TABLE]

and hence the initial deviation is zero, by definition. Iterating, starting at $i=0$ yields:

[TABLE]

where the last line follows from Lemma 72 after noting that:

[TABLE]

This establishes (74). We proceed to establish a bound on the fourth-order moment. Using the inequality [4]:

[TABLE]

we have:

[TABLE]

where in step $(a)$ we dropped cross-terms due to the conditional zero-mean property of the gradient noise in Assumption II, step $(b)$ follows from the fourth-order conditions on the gradient noise in Assumption II along with the perturbation bounds in Lemma 49, and $(c)$ follows from Jensen’s inequality, i.e. ${\|a+b+c\|}^{2}\leq 3{\|a\|}^{2}+3{\|b\|}^{2}+3{\|c\|}^{2}$. Taking expectations over $\w_{c,i^{\star}}\in\mathcal{H}$ on both sides and collecting constant factors along with $\mu$ in appropriate $O(\cdot)$ terms:

[TABLE]

Finally, from Jensen’s inequality, we find for $0<\alpha<1$:

[TABLE]

and hence for $\alpha=1-\mu\delta$ and $0<\mu<\frac{1}{\delta}$:

[TABLE]

Hence,

[TABLE]

Recall again that $\widetilde{\w}_{0}^{i^{\star}}=0$ and therefore iterating yields:

[TABLE]

where in $(a)$ we expanded:

[TABLE]

and the last step follows from Lemma 72. This establishes (76). Eq. (75) then follows from Jensen’s inequality via:

[TABLE]

We now study the difference between the short-term model (71) and the true recursion (69). We have:

[TABLE]

Before proceeding, note that the difference between the Hessians in the driving term can be bounded as:

[TABLE]

where $(a)$ follows Jensen’s inequality and $(b)$ follows form the Lipschitz Hessian assumption II-D. Returning to (118) and taking norms yields:

[TABLE]

where $(a)$ again follows from Jensen’s inequality (104) with $\alpha=1-\mu\delta$ and $(b)$ follows from the same inequality with $\alpha=\frac{1}{2}$. Taking expecations over $\w_{c,i^{\star}}\in\mathcal{H}$ yields:

[TABLE]

where $(a)$ follows from the bound on the network disagreement in Lemma II-D.

Since both the true and the short-term model are initialized at $\w_{c,i^{\star}}$, we have $\w_{c,i^{\star}+0}-\w_{c,i^{\star}+0}^{\prime}=0$. Iterating and applying the same argument as above leads to:

[TABLE]

which is (77).