Newton-Raphson Consensus for Distributed Convex Optimization

Introduction

Optimization is a pervasive concept underlying many aspects of modern life , and it also includes the management of distributed systems, i.e., artifacts composed by a multitude of interacting entities often referred to as “agents”. Examples are transportation systems, where the agents are both the vehicles and the traffic management devices (traffic lights), and smart electrical grids, where the agents are the energy producers-consumers and the power transformers-transporters.

Here we consider the problem of distributed optimization, i.e., the class of algorithms suitable for networked systems and characterized by the absence of a centralized coordination unit . Distributed optimization tools have received an increasing attention over the last years, concurrently with the research on networked control systems. Motivations comprise the fact that the former methods let the networks self-organize and adapt to surrounding and changing environments, and that they are necessary to manage extremely complex systems in an autonomous way with only limited human intervention. In particular we focus on unconstrained convex optimization, although there is a rich literature also on distributed constrained optimization such as Linear Programming .

The literature on distributed unconstrained convex optimization is extremely vast and a first taxonomy can be based whether the strategy uses or not the Lagrangian framework, see, e.g., [5, Chap. 5].

Among the distributed methods exploiting Lagrangian formalism, the most widely known algorithm is Alternating Direction Method of Multipliers (ADMM) , whose roots can be traced back to . Its efficacy in several practical scenarios is undoubted, see, e.g., and references therein. A notable size of the dedicated literature focuses on the analysis of its convergence performance and on the tuning of its parameters for optimal convergence speed, see, e.g., for Least Squares (LS) estimation scenarios, for linearly constrained convex programs, and for more general ADMM algorithms. Even if proved to be an effective algorithm, ADMM suffers from requiring synchronous communication protocols, although some recent attempts for asynchronous and distributed implementations have appeared .

On the other hand, among the distributed methods not exploiting Lagrangian formalisms, the most popular ones are the Distributed Subgradient Methods . Here the optimization of non-smooth cost functions is performed by means of subgradient based descent/ascent directions. These methods arise in both primal and dual formulations, since sometimes it is better to perform dual optimization. Subgradient methods have been exploited for several practical purposes, e.g., to optimally allocate resources in Wireless Sensor Networks , to maximize the convergence speeds of gossip algorithms , to manage optimality criteria defined in terms of ergodic limits . Several works focus on the analysis of the convergence properties of the DSM basic algorithm (see for a unified view of many convergence results). We can also find analyses for several extensions of the original idea, e.g., directions that are computed combining information from other agents and stochastic errors in the evaluation of the subgradients . Explicit characterizations can also show trade-offs between desired accuracy and number of iterations .

These methods have the advantage of being easily distributed, to have limited computational requirements and to be inherently asynchronous as shown in . However they suffer from low convergence rate since they require the update steps to decrease to zero as $1/t$ (being $t$ the time) therefore as a consequence the rate of convergence is sub-exponential. In fact, one of the current trends is to design strategies that improve the convergence rate of DSMs. For example, a way is to accelerate the convergence of subgradient methods by means of multi-step approaches, exploiting the history of the past iterations to compute the future ones . Another is to use Newton-like methods, when additional smoothness assumptions can be used. These techniques are based on estimating the Newton direction starting from the Laplacian of the communication graph. More specifically, distributed Newton techniques have been proposed in dual ascent scenarios . Since the Laplacian cannot be computed exactly, the convergence rates of these schemes rely on the analysis of inexact Newton methods . These Newton methods are shown to have super-linear convergence under specific assumptions, but can be applied only to specific optimization problems such as network flow problems.

Recently, several alternative approaches to ADMM and DSM have appeared. For example, in the authors construct contraction mappings by means of cyclic projections of the estimate of the optimum onto the constraints. A similar idea based on contraction maps is used in F-Lipschitz methods but it requires additional assumptions on the cost functions. Other methods are the control-based approach which exploits distributed consensus, the distributed randomized Kaczmarz method for quadratic cost functions, and distributed dual sub-gradient methods .

Statement of contributions

Here we propose a distributed Newton-Raphson optimization procedure, named Newton-Raphson Consensus (NRC), for the exact minimization of smooth multidimensional convex separable problems, where the global function is a sum of private local costs. With respect to the classification proposed before, the strategy exploits neither Lagrangian formalisms nor Laplacian estimation steps. More specifically, it is based on average consensus techniques and on the principle of separation of time-scales [46, Chap. 11]. The main idea is that agents compute and keep updated, by means of average consensus protocols, an approximated Newton-Raphson direction that is built from suitable Taylor expansions of the local costs. Simultaneously, agents move their local guesses towards the Newton-Raphson direction. It is proved that, if the costs satisfy some smoothness assumptions and the rate of change of the local update steps is sufficiently slow to allow the consensus algorithm to converge, then the NRC algorithm exponentially converges to the global minimizer.

The main contribution of this work is to propose an algorithm that extends Newton-Raphson ideas in a distributed setting, thus being able to exploit second order information to speed up converge rate. By using singular perturbation theory we formally show that under suitable assumptions the convergence of the algorithm is exponential (linear in logspace). Differently, DSM algorithms have sublinear convergence rate even if the cost functions are smooth , although they are easy to implement and can be employed also for non-smooth cost functions and for constrained optimization. We also show by means of numerical simulations on real-world database benchmarks that the proposed algorithm exhibits faster convergence rates (in number of communications) than standard implementations of distributed ADMM algorithms , probably due to the second-order information embedded into the Newton-Raphson consensus. Although we have no theoretical guarantee of the superiority of the proposed algorithmic in terms of convergence rate, these simulations suggest that it is at least a potentially competitive algorithm. Moreover, one of the promising features of the NRC is that it is essentially based on average consensus algorithms, for which there exist robust implementations that encompass asynchronous communications, time-varying network topologies , directed graphs , and packet-losses effects.

Structure of the paper

7 collects the notation used through the whole paper, while

8 formulates the considered problem and provides some ancillary results that are then used to study the convergence properties of the main algorithm.

14 proposes the main optimization algorithm, provides convergence results and describes some strategies to trade-off communication and computational complexities with convergence speed.

16 compares, via numerical simulations, the performance of the proposed algorithm with several distributed optimization strategies available in the literature. Finally,

23 collects some final observations and suggests future research directions. We collect all the proofs in the Appendix.

Notation

As in [46, p. 116], we say that a function $V$ is a Lyapunov function for a specific dynamics if $V$ is continuously differentiable and satisfies $V(0)=0$, $V(x)>0$ for $x\neq 0$, and $\dot{V}(x)\leq 0$.

Problem formulation and preliminary results

Our main contribution is to characterize the convergence properties of the distributed Newton-Raphson (NR) scheme proposed in

14. In doing so we both exploit standard singular perturbation analysis tools [46, Chap. 11] and a set of ancillary results, collected for readability in this section.

The logical flow of these ancillary results is the following:

10.2 claims that, under suitable assumptions, forward-Euler discretizations of stable continuous dynamics lead to stable discrete dynamics. This basic result enables reasoning on continuous-time systems. Then, Sections 10.3 and 11.1 respectively claim that single- and multi-agent continuous-time NR dynamics satisfy these discretization assumptions. Sections 12.1 and 13.1 then generalize these dynamics by introducing perturbation terms that mimic the behavior of the proposed main optimization algorithm, and characterize their stability properties. Summarizing, the ancillary results characterize the stability properties of systems that are progressive approximations of the dynamics under investigation.

is a well-defined global cost. We assume that the aim of the agents is to cooperate and distributedly compute the minimizer of $\overline{f}$, namely

We now enforce the following simplifying assumptions, valid throughout the rest of the paper:

The scalar $c$ is assumed to be known by all the agents a-priori. Assumption 1 ensures that $x^{\ast}$ in (2) exists and is unique. The strictly positive definite Hessian is moreover a mild sufficient condition to guarantee that the minimum $x^{*}$ defined in (2) will be globally exponentially stable under the continuous and discrete Newton-Raphson dynamics described in the following Theorem 3. We also notice that, for the subsequent Theorems 2 and 3, in principle just the average function $\overline{f}$ needs to have specific properties, and thus no conditions for the single $f_{i}$’s are required (that for example might be even non convex). For the convergence of the distributed NR scheme we will nonetheless enforce the more restrictive Assumptions 5 and 9, not presented now for readability issues. In the rest of this section, in order to simplify notation, we will considerer, without loss of generality, the following translated cost functions:

so that the origin becomes the minimizer of the averaged cost function $\overline{f^{\prime}}(x)$, i.e. $\overline{f^{\prime}}(0)=0$.

2 Stability of discretized dynamics

This subsection aims to show that, under suitable assumptions, forward-Euler discretization of suitable exponentially stable continuous-time dynamics maintains the same global exponential stability properties.

for system (4) the origin is globally exponentially stable;

for the following forward-Euler discretization of system (4),

there exists a positive scalar $\overline{\varepsilon}$ such that for every $\varepsilon\in(0,\overline{\varepsilon})$ the origin is globally exponentially stable.

3 Stability of single-agent NR dynamics

This subsection shows that the results of

10.2 apply to continuous NR dynamics, i.e., that forward-Euler discretizations maintain global exponential stability propertiesWe notice that other asymptotic properties of continuous time NR methods are available in the literature, e.g., ..

i.e., Theorem 2 applies to dynamics (8) and (9).

For suitable choices of $\overline{\overline{h^{\prime}}}(x)$ the dynamics (8) corresponds to continuous versions of well known descent dynamics. Indeed, the correspondences are

Under Assumption 1, the origin is a globally exponentially stable point for dynamics (8). Moreover there exists $\overline{\varepsilon}>0$ such that the origin is a globally exponentially stable point also for dynamics (9) for all $\varepsilon<\overline{\varepsilon}$.

The previous lemma and theorems do not require $\overline{\overline{h^{\prime}}}(x)$ to be differentiable. However, differentiability may be used to linearize the system dynamics and obtain explicit rates of convergence. In fact, the linearized dynamics around the origin is given by

In particular, for the NR descent it holds that $\overline{\overline{h^{\prime}}}(x)=\nabla^{2}\overline{f^{\prime}}(x)$. Thus in this case $F(0)=-I$, since $\nabla\overline{f^{\prime}}(0)=0$, and this says that the linearized continuous time NR dynamics is $\dot{x}=-x$, independent of the cost $\overline{f^{\prime}}(x)$ and whose rate of convergence is unitary and uniform along any direction.

We now generalize (8) by considering $N$ coupled dynamical systems that, when starting at the very same initial condition, behave like $N$ decoupled systems (8). This novel dynamics is the core of the slow-dynamics embedded in the main algorithm presented in

14. In this section we also include additional assumptions to show that the generalization of (8) presented here preserves global exponential stability and some other additional properties.

so that $h^{\prime}_{i}(x)=h^{\prime}_{i}(x)^{T}$ for all $x$. Moreover let

and $g^{\prime}(\bm{x}),\overline{g^{\prime}}(\bm{x}),\overline{\overline{g^{\prime}}}(\overline{x})$ be defined accordingly as for $h^{\prime}_{i}$.

The definitions of $h^{\prime}_{i}$ and $g^{\prime}_{i}$ are instrumental to generalize the NR dynamics (8) to the distributed case. Indeed, let

(with the existence of $\overline{h^{\prime}}(\bm{x})^{-1}$ guaranteed by the following Assumption 5). It is easy to verify that the previous functions satisfy the following properties:

i.e., as the combination of $N$ independent dynamical systems that are driven by the same forcing term $\psi(\bm{x})$.

i.e., we obtain dynamics (7), that is, thanks to Theorem 3 and the assumption that $\overline{h^{\prime}}(\bm{x})$ is invertible, globally exponentially stable.

The question is then whether dynamics (17) is exponentially stable also in the general case where the $x_{i}(0)$’s may not be identical. To characterize this case we assume some additional global properties:

Using the previous assumptions we can now prove global stability of dynamics (17):

Under Assumptions 1 and 5, and for a suitable positive scalar $\eta$,

As in Lemma 4, combining Theorem 6 with Theorem 2 it is possible to claim that (17) and its discrete-time counterpart are globally exponentially stable.

Consider then the perturbed version of the multi-agent NR dynamics (17),

where the division is a Hadamard division, as recalled in

7. Direct inspection of dynamics (24) then shows that

The next lemma provides perturbations interconnection bounds that will be used in Theorem 12.

Under Assumptions 1 and 5 there exist positive scalars $a_{x}$, $a_{\Delta}$ s.t., for all $\bm{x}$ and $\bm{\chi}$,

Moreover, $\bm{x}^{eq}(\xi)=\mathds{1}_{N}\otimes x^{eq}(\xi)$, with

which corresponds to the translated version of the original perturbed system $\phi_{x}(\bm{x},\bm{\xi})$, which has now the property that the origin is an equilibrium point, i.e., $\phi^{\prime}_{x}(\bm{0},\xi)=0,\forall\|\xi\|\leq r$.

To prove the global exponential stability of (30) we need the flow $\phi^{\prime}_{x}$ to satisfy a global Lipschitz condition:

With these assumptions we can prove that the origin is a globally exponentially stable equilibrium for dynamics (30):

$V_{\textrm{PNR}}(\bm{x})$ defined in (21) is a Lyapunov function for (30);

Again, as in Lemma 4, combining Theorem 10 with Theorem 2 it is possible to claim that (30) and its discrete-time counterpart are globally exponentially stable.

2 Quadratic Functions

Before presenting the main algorithm, we show that quadratic costs satisfy all the previous assumptions. In fact, let us consider then

Based on this definition we have the following result:

satisfy Assumptions 1, 5 and 9 for $h^{\prime}_{i}(x)=\nabla^{2}f^{\prime}_{i}(x)$.

Newton-Raphson Consensus

In this section we provide an algorithm to distributively compute the minimizer of the function $x^{*}$ defined in (2). The algorithm will be shown to converge to $x^{*}$ even if $x^{*}\neq 0$. The proof of convergence will be based on the results derived in the previous sections via a suitable translation of the argument of the cost functions, which basically reduces the problem to the special case $x^{*}=0$.

Consider then Algorithm 1, where g\big{(}\bm{x}(-1)\big{)}=\bm{0} and h\big{(}\bm{x}(-1)\big{)}=\bm{0} in the initialization step should be intended as initialization of suitable registers and not as operations involving the quantity $\bm{x}(-1)$.

Intuitively, the algorithm functions as follows: if the dynamics of the $x_{i}(k)$s is sufficiently slow w.r.t. the dynamics of the $y_{i}(k)$s and $z_{i}(k)$s, then the two latter quantities tend to reach consensus. Then, the more these quantities reach consensus, the more the products $\left[z_{i}(k)\right]_{c}^{-1}y_{i}(k)$ exhibit these two specific characteristics: i) being the same among the various agent; ii) representing Newton descent directions. Thus, the more the $y_{i}(k)$s and $z_{i}(k)$s in Algorithm 1 are sufficiently close, the more the various $x_{i}(k)$s are driven by the same forcing term, that makes them converge to the same value, equal to the optimum $x^{\ast}$.

We now characterize the convergence properties of Algorithm 1. Let us define

Consider the dynamics defined by Algorithm 1 with possibly nonzero initial conditions. If $\xi^{y}=0$ and $\xi^{z}=0$, then under Assumptions 1 and 5 there exists a positive scalar $\overline{\varepsilon}>0$ such that Theorem 2 holds, i.e., the algorithm can be considered a forward-Euler discretization of a globally exponentially stable continuous dynamics. Thus the local estimates $x_{i}(k)$ produced by the algorithm exponentially converge to the global minimizer, i.e.,

Consider now that, due to finite-precision issues, the quantities $\xi^{y}$ and $\xi^{z}$ may be non-null. Non-null initial $\xi^{y}$ and $\xi^{z}$ will make the proposed algorithm converge to a point that, in general does not coincide with the global optimum $x^{\ast}$. Nonetheless in this case the computed solution, as a function of the initial conditions, is a smooth function and thus small errors in the initial conditions do not produce dramatic errors in the computation of the optimum:

s.t. the local estimates exponentially converge to it, i.e.,

We notice that Theorem 13 ensures global convergence properties w.r.t. the initial conditions $x_{i}(0)$’s by requiring Assumptions 1, 5 and 9, while for the same convergence properties Theorem 12 requires only Assumptions 1 and 5. The difference is that Theorem 13 considers a non-null perturbation $\xi$ and Assumption 9 is needed to cope with this additional perturbation term.

The Assumptions 1, 5 and 9 are not needed if only local convergence is ought. In fact, local differentiability, and therefore local Lipschitzianity, of the cost functions $f_{i}(x)$ at the minimizer $x^{*}$ is sufficient to guarantee that Assumptions 5 and 9 are locally valid. As so, the proof that the equilibrium point is a locally exponentially stable point is exactly the same, with the difference that all bounds and inequalities are local. This observation is summarized in the following theorem.

for all $\varepsilon\in(0,\overline{\varepsilon})$ and initial conditions

Numerical simulations suggest that the algorithm is robust w.r.t. numerical errors and quantization noise. We also notice that Theorem 12 guarantees the existence of a critical value $\overline{\varepsilon}$ but does not provide indications on its value. This is a known issue in all the systems dealing with separation of time scales. A standard rule of thumb is then to let the rate of convergence of the fast dynamics be sufficiently faster than the one of the slow dynamics, typically 2-10 times faster. In our algorithm the fast dynamics inherits the rate of convergence of the consensus matrix $P$, given by its spectral gap $\sigma(P)$, i.e., its spectral radius $\rho(P)=1-\sigma(P)$. The rate of convergence of the slow dynamics is instead governed by (18), which is nonlinear and therefore possibly depending on the initial conditions. However, close to the equilibrium point the dynamic behavior is approximately given by \dot{\overline{x}}(t)\approx-\big{(}\overline{x}(t)-x^{*}\big{)}, thus, since $x_{i}(k)\approx\overline{x}(\varepsilon k)$, then the convergence rate of the algorithm approximately given by $1-\varepsilon$.

Thus we aim to let $1-\rho(P)\gg 1-(1-\varepsilon),$ which provides the rule of thumb

which is suitable for generic cost functions. We then notice that, although the spectral gap $\sigma(P)$ might not be known in advance, it is possible to distributedly estimate it, see, e.g., . However, such rule of thumb might be very conservative. In fact, if all the $f_{i}$’s are quadratic and are, w.l.o.g. s.t. $\nabla^{2}f_{i}\geq cI$, then one can set $\varepsilon=1$ and neglect the thresholding $\left[\cdot\right]_{c}$, so that the procedure reduces to

where $\bm{x}(k)\mathrel{\mathop{:}}=\left[x_{1}^{T}(k)\;,\ldots,\;x_{N}^{T}(k)\right]^{T}$, $\bm{y}(k)\mathrel{\mathop{:}}=\left[y_{1}^{T}(k)\;,\ldots,\;y_{N}^{T}(k)\right]^{T}$, $\bm{z}(k)\mathrel{\mathop{:}}=\left[z_{1}(k)\;,\ldots,\;z_{N}(k)\right]^{T}$. Thus:

Consider Algorithm 1 with arbitrary initial conditions $x_{i}(0)$, quadratic cost functions $f_{i}=\frac{1}{2}\left(x-d_{i}\right)^{T}A_{i}\left(x-d_{i}\right)$ with $A_{i}>0$ and $\varepsilon=1$. Then $\left\|x_{i}(k)-x^{\ast}\right\|\leq\alpha\left(\rho(P)\right)^{k}$ for all $k,i$ and for a suitable positive $\alpha$.

Thus, if the cost functions are close to be quadratic then the overall rate of convergence is limited by the rate of convergence of the embedded consensus algorithm. Moreover, the values of $\varepsilon$ that still guarantee convergence can be much larger than those dictated by the rule of thumb (32).

10.3, by selecting different structures for $h_{i}(x)$ one can obtain different procedures with different convergence properties and different computational/communication requirements. Plausible choices for $h_{i}$ are the ones in (13c), and the correspondences are the following:

$\bullet$ $h_{i}(x)=\nabla^{2}f_{i}(x)$ $\rightarrow$ Newton-Raphson Consensus (NRC): in this case it is possible to rewrite the main algorithm and show that, for sufficiently small $\varepsilon$, $x_{i}(k)\approx\overline{x}(\varepsilon k)$, where $\overline{x}(t)$ evolves according to the continuous-time Newton-Raphson dynamics

which can be shown to converge to the global optimum $x^{\ast}$ with a convergence rate that in general is slower than the Newton-Raphson when the global cost function is skewed.

$\bullet$ $h_{i}(x)=I$ $\rightarrow$ Gradient Descent Consensus (GDC): this choice is motivated in frameworks where the computation of the local second derivatives $\displaystyle\left.\frac{\partial^{2}f_{i}}{\partial x_{m}^{2}}\right|_{x}$ is expensive (with $x_{m}$ indicating here the $m$-th component of $x$), or where the second derivatives simply might not be continuous. With this choice the main algorithm reduces to a distributed gradient-descent procedure. In fact, for sufficiently small $\varepsilon$, $x_{i}(k)\approx\overline{x}(\varepsilon k)$ with $\overline{x}(t)$ evolving according to the continuous-time dynamics

which one again is guaranteed to converge to the global optimum $x^{\ast}$.

The following Table 1 summarizes the various costs of the previously proposed strategies.

We remark that $\overline{\varepsilon}$ in Theorem 12 depends also on the particular choice for $h_{i}$. The list of choices for $h_{i}$ given above is not exhaustive. For example, future directions are to implement distributed quasi-Newton procedures. To this regard, we recall that approximations of the Hessians that do not maintain symmetry and positive definiteness or are bad conditioned require additional modification steps, e.g., through Cholesky factorizations .

Finally, we notice that in scalar scenarios JC and NRC are equivalent, while GDC corresponds to algorithms requiring just the knowledge of first derivatives.

Numerical Examples

18.1 we analyze the effects of different choices of $\varepsilon$ on the NRC on regular graphs and exponential cost functions. We then propose two machine learning problems in

20.1, used in Sections 20.2 and 20.3, and numerically compare the convergence performance of the NRC, JC, GDC algorithms and other distributed convex optimization algorithms on random geometric graphs.

Notice that we will use cost functions that may not satisfy Assumptions 1, 5 and 9 to highlight the fact that the algorithm seems to have favorable numerical properties and large basins of stability even if the assumptions needed for global stability are not satisfied.

Consider a ring network of $S=30$ agents that communicate only to their left and right neighbors through the consensus matrix

so that the spectral radius $\rho(P)\approx 0.99$, implying a spectral gap $\sigma(P)\approx 0.01$. Consider also scalar costs of the form $f_{i}(x)=c_{i}e^{a_{i}x}+d_{i}e^{-b_{i}x},$ $i=1,\ldots,N,$ with $a_{i},b_{i}\sim\mathcal{U}\left[0,0.2\right]$, $c_{i},d_{i}\sim\mathcal{U}\left[0,1\right]$ and where $\mathcal{U}$ indicates the uniform distribution.

Figure 1 compares the evolution of the local states $x_{i}$ of the continuous system (43) for different values of $\varepsilon$. When $\varepsilon$ is not sufficiently small, then the trajectories of $x_{i}(t)$ are different even if they all start from the same initial condition $x_{i}(0)=0$. As $\varepsilon$ decreases, the difference between the two time scales becomes more evident and all the trajectories $x_{i}(k)$ become closer to the trajectory given by the slow NR dynamics $\overline{x}(\varepsilon k)$ given in (18) and guaranteed to converge to the global optimum $x^{\ast}$.

In Figure 2 we address the robustness of the proposed algorithm w.r.t. the choice of the initial conditions. In particular, Figure 2(a) shows that if $\alpha=\beta=0$ then the local states $x_{i}(t)$ converge to the optimum $x^{\ast}$ for arbitrary initial conditions $x_{i}(0)$. Figure 2(b) considers, besides different initial conditions $x_{i}(0)$, also perturbed initial conditions $v(0)$, $w(0)$, $y(0)$, $z(0)$ leading to non null $\alpha$’s and $\beta$’s. More precisely we apply Algorithm 1 to different random initial conditions s.t. $\alpha,\beta\sim\mathcal{U}\left[-\sigma,\sigma\right]$. Figure 2(b) shows the boxplots of the errors $x_{i}(+\infty)-x^{\ast}$ for different $\sigma$’s based on 300 Monte Carlo runs with $\varepsilon=0.01$ and $N=30$.

1 Optimization problems

where $E_{i}$ is the set of emails available to agent $i$, $E=\cup_{i=1}^{N}E_{i}$, and $\gamma$ is a global regularization parameter. In the following numerical experiments we consider $|E|=5000$ emails from the spam-nonspam UCI repository, available at http://archive.ics.uci.edu/ml/datasets/Spambase, randomly assigned to 30 different users communicating as in graph of Figure 4. For each email we consider $3$ features (the frequency of words “make”, “address”, “all”) so that the corresponding optimization problem is 4-dimensional.

In both the previous problems the optimum, in the following indicated for simplicity with $x^{\ast}$, has been computed with a centralized NR with the termination rule “stop when in the last 5 steps the norm of the guessed $x^{\ast}$ changed less than $10^{-9}\%$”.

2 Comparison of the NRC, JC and GDC algorithms

In Figure 20 we analyze the performance of the three proposed NRC, JC and GDC algorithms defined by the various choices for $h_{i}(x)$ in Algorithm 1 in terms of the relative MSE

for the classification and regression optimization problem described above. The consensus matrix $P$ has been by selecting the Metropolis-Hastings weights which are consistent with the communication graph . Panels 3(a) and 3(c) report the MSE obtained at a specific iteration ($k=40$) by the various algorithms, as a function of $\varepsilon$. These plots thus inspect the sensitivity w.r.t. the choice of the tuning parameters. Consistently with the theorems in the previous section, the GDC and JC algorithms are stable only for $\varepsilon$ sufficiently small, while NRC exhibit much larger robustness and best performance for $\varepsilon=1$. Panels 3(b) and 3(d) instead report the evolutions of the relative MSE as a function of the number of iterations $k$ for the optimally tuned algorithms.

We notice that the differences between NRC and JC are evident but not resounding, due to the fact that the Jacobi approximations are in this case a good approximation of the analytical Hessians. Conversely, GDC presents a slower convergence rate which is a known drawback of gradient descent algorithms.

3 Comparisons with other distributed convex optimization algorithms

We now compare Algorithm 1 and its accelerated version, referred as Fast Newton-Raphson Consensus (FNRC) and described in detail below in Algorithm 2), with three popular distributed convex optimization methods, namely the DSM, the Distributed Control Method (DCM) and the ADMM, described respectively in Algorithm 3, 4 and 5. The following discussion provides some details about these strategies.

$\bullet$ FNRC is an accelerated version of Algorithm 1 that inherits the structure of the so called second order diffusive schedules, see, e.g., , and exploits an additional level of memory to speed up the convergence properties of the consensus strategy. Here the weights multiplying the $g_{i}$’s and $h_{i}$’s are necessary to guarantee exact tracking of the current average, i.e., \sum_{i}y_{i}(k)=\sum_{i}g_{i}\big{(}x(k-1)\big{)} for all $k$. As suggested in , we set the $\varphi$ that weights the gradient and the memory to $\displaystyle\varphi=\frac{2}{1+\sqrt{1-\rho(P)^{2}}}$. This guarantees second order diffusive schedules to be faster than first order ones (even if this does not automatically imply the FNRC to be faster than the NRC). This setting can be considered a valid heuristic to be used when $\rho(P)$ is known. For the graph in Figure 4, $\varphi\approx 1.4730$.

$\bullet$ DSM, as proposed in , alternates consensus steps on the current estimated global minimum $x_{i}(k)$ with subgradient updates of each $x_{i}(k)$ towards the local minimum. To guarantee the convergence, the amplitude of the local subgradient steps should appropriately decrease. Algorithm 3 presents a synchronous DSM implementation, where $\varrho$ is a tuning parameter and $P$ is the matrix of Metropolis-Hastings weights.

$\bullet$ DCM, as proposed in , differentiates from the gradient searching because it forces the states to the global optimum by controlling the subgradient of the global cost. This approach views the subgradient as an input/output map and uses small gain theorems to guarantee the convergence property of the system. Again, each agents $i$ locally computes and exchanges information with its neighbors, collected in the set $\mathcal{N}_{i}\mathrel{\mathop{:}}=\{j\,|\,(i,j)\in\mathcal{E}\}$. DCM is summarized in Algorithm 4, where $\mu,\nu>0$ are parameters to be designed to ensure the stability property of the system. Specifically, $\mu$ is chosen in the interval $\displaystyle 0<\mu<\frac{2}{2\max_{i=\{1,\ldots,N\}}{|\mathcal{N}_{i}|}+1}$ to bound the induced gain of the subgradients. Also here the parameters have been manually tuned for best convergence rates.

$\bullet$ ADMM, instead, requires the augmentation of the system through additional constraints that do not change the optimal solution but allow the Lagrangian formalism. There exist different implementations of ADMM in distributed contexts, see, e.g., [7, 60, 12, pp. 253-261]. For simplicity we consider the following formulation,

where the auxiliary variables $z_{(i,j)}$ correspond to the different links in the network, and where the local Augmented Lagrangian is given by

with $\delta$ a tuning parameter (see for a discussion on how to tune it) and the $y_{(i,j)}$’s Lagrange multipliers.

The computational, communication and memory costs of these algorithms is reported in Table 2. Notice that the computational and memory costs of ADMM algorithms depends on how nodes minimize the local augmented Lagrangian $L_{i}(x_{i},k)$. E.g., in our simulations the step has been performed through a dedicated Newton-Raphson procedure with associated $O\left(n^{3}\right)$ computational costs and $O\left(n^{2}\right)$ memory costs.

Figure 22 then compares the previously cited algorithms as did in Figure 20. The first panel thus reports the relative MSE of the various algorithms at a given number of iterations ($k=40$) as a function of the parameters. The second panel instead reports the temporal evolution of the relative MSE for the case of optimal tuning.

We notice that the DCM and the DSM are both much slower, in terms of communications iterations, than the NRC, FNRC and ADMM. Moreover, both the NRC and its accelerated version converge faster than the ADMM, even if not tuned at their best. These numerical examples seem to indicate that the proposed NRC might be a viable alternative to the ADMM, although further comparisons are needed to strengthen this claim. Moreover, a substantial potential advantage of NRC as compared to ADMM is that the former can be readily adapted to asynchronous and time-varying graphs, as preliminary made in . Moreover, as in the case of the FNRC, the strategy can implement any improved linear consensus algorithm.

Conclusion

We proposed a novel distributed optimization strategy suitable for convex, unconstrained, multidimensional, smooth and separable cost functions. The algorithm does not rely on Lagrangian formalisms and acts as a distributed Newton-Raphson optimization strategy by repeating the following steps: agents first locally compute and update second order Taylor expansions around the current local guesses and then they suitably combine them by means of average consensus algorithms to obtain a sort of approximated Taylor expansion of the global cost. This allows each agent to infer a local Newton direction, used to locally update the guess of the global minimum.

Importantly, the average consensus protocols and the local updates steps have different time-scales, and the whole algorithm is proved to be convergent only if the step-size is sufficiently slow. Numerical simulations based on real-world databases show that, if suitably tuned, the proposed algorithm is faster then ADMMs in terms of number of communication iterations, although no theoretical proof is provided.

The set of open research paths is extremely vast. We envisage three main avenues. The first one is to study how the agents can dynamically and locally tune the speed of the local updates w.r.t. the consensus process, namely how to tune their local step-size $\varepsilon_{i}$. In fact large values of $\varepsilon$ gives faster convergence but might lead to instability. A second one is to let the communication protocol be asynchronous: in this regard we notice that some preliminary attempts can be found in . A final branch is about the analytical characterization of the rate of convergence of the proposed strategies, a theoretical comparison with ADMMs, and the extensions to non-smooth convex functions.

References

A

$\bullet$ Discrete to continuous dynamics) The dynamics of Algorithm 1 can be written in state space as

with suitable initial conditions. (42) can then be interpreted as the forward-Euler discretization of

with $\bm{\chi}\mathrel{\mathop{:}}=\left(\bm{v}^{\prime},\bm{w}^{\prime},\bm{y}^{\prime},\bm{z}^{\prime}\right)$, so that (43) becomes

where $g^{\perp}\left(\bm{x}\right)\mathrel{\mathop{:}}=g(\bm{x})-\mathds{1}_{N}\otimes\overline{g}(\bm{x})$ (equivalent definition for $h^{\perp}$). Notice that (44) has the origin as an equilibrium point. Moreover this dynamics exploits the function $\phi_{x}$ defined in (24), with $\bm{\chi}^{y}=\bm{y}^{\prime}+\bm{v}^{\prime}$, and $\bm{\chi}^{z}=\bm{z}^{\prime}+\bm{w}^{\prime}$.

The next step is to exploit the structure of $K$ (more precisely, the fact that it contains an average consensus matrix) to reduce the dynamics, i.e., to eliminate the dynamics of the average since the latter does not change in time. To this aim, we analyze the behavior of the average of the $y_{i}^{\prime}$s, i.e., the behavior of $\left(\mathds{1}_{N}^{T}\otimes I_{n}\right)\dot{\bm{y}}^{\prime}$. To this point, consider the third equation in (44). Recalling that $(A\otimes B)(C\otimes D)=AB\otimes CD$, and exploiting the fact that $\mathds{1}_{N}^{T}P^{\prime}=0$, we notice that $\left(\mathds{1}_{N}^{T}\otimes I_{n}\right)K=0$. Moreover, from the definitions of $g$ and $\overline{g}$,

Since $N=\mathds{1}_{N}^{T}\mathds{1}_{N}$, it follows also that

for all $t\geq 0$, i.e., $\mathds{1}^{T}\bm{y}^{\prime}(t)=\mathds{1}^{T}\bm{y}^{\prime}(0)\equiv 0$. Similarly it is possible to show that $\bm{z}^{\prime}(t)\equiv 0$. This eventually implies that

that means, recalling that $\bm{y}^{\prime}=\bm{y}^{\prime\parallel}+\bm{y}^{\prime\perp}$ and $\bm{z}^{\prime}=\bm{z}^{\prime\parallel}+\bm{z}^{\prime\perp}$, that (44) can be equivalently rewritten as

where now $\bm{\chi}^{\prime}\mathrel{\mathop{:}}=\left(\bm{v},\bm{w},\bm{y}^{\prime\perp},\bm{z}^{\prime\perp}\right)$ and where the novel initial conditions for the changed variables are

$\bullet$ c) the boundary layer of (45) is computed by setting $\bm{x}^{\prime}(t)=\bm{x}^{\prime}$. Since a constant $\bm{x}^{\prime}$ implies $\dot{\bm{x}}^{\prime}=\phi_{x}=0$, this boundary layer reduces to a linear system globally exponentially converging to the origin. Notice that this implies that, in the original coordinates system,

In the novel coordinates system we thus consider, as a Lyapunov function, $\frac{1}{2}\|\bm{\chi}^{\prime}\|^{2}$.

$\bullet$ d) the reduced system of (45) is computed by plugging $\bm{\chi}^{\prime}=\bm{0}$ into the equations (i.e., by setting $\bm{v}^{\prime}(t)=\bm{0}$, $\bm{w}^{\prime}(t)=\bm{0}$, $\bm{y}^{\prime\perp}(t)=\bm{0}$, $\bm{z}^{\prime\perp}(t)=\bm{0}$). Defining then

where $\psi$ and $\phi_{\textrm{PNR}}$ are the functions defined in (15) and (17), respectively. Thus the reduced system, thanks to Theorem 6, admits $\bm{x}^{\ast}$ as a global exponentially stable equilibrium, and admits $V_{\textrm{PNR}}$ in (21) as a Lyapunov function.

$\bullet$ e) we now notice that the interconnection of the boundary layer and reduced systems maintains the global stability, since their Lyapunov functions are quadratic type. Thus (see [46, pp. 453]) the global system is asymptotically globally stable. To check that forward-Euler discretizations of the system preserve these stability properties we then consider as a global Lyapunov function the function

that is clearly positive definite for every $d\in(0,1)$, and prove that inequalities (5c) of Theorem 2 are satisfied.

proof that (5a) holds: from (22a) and the structure of $V$ it follows immediately that

proof that (5c) holds: applying (20d) and (26a) to (45) it follows that (5c) holds with

proof that (5b) holds: the part relative to the slow dynamics is already characterized by (31a). For the part relative to the fast dynamics, since $\frac{\partial\frac{1}{2}\|\bm{\chi}\|^{2}}{\partial\bm{\chi}}=\bm{\chi}^{T}$ to check that (5b) corresponds to check the negativity of the terms

These terms can then be majorized using (20d) and (26a). E.g., the third term can be majorized with

where $\sigma(P)$ is the spectral gap of $P$. Applying similar concepts also to the other terms it follows that (5b) holds with

The proof is identical to the one of Theorem 12 with the exception that the substitution is now $\bm{x}^{\prime\prime}=\bm{x}-\bm{x}^{\ast}-\mathds{1}_{N}\otimes\Psi(\xi^{y},\xi^{z})$. Indeed one can prove the stability of the novel system using the same Lyapunov function of Theorem 12. Notice that we are ensured that there exists a sufficiently small neighborhood of the origin for which the function $\Psi$ exists due to the smoothness conditions in Assumption 1. $\diamondsuit$

The proof is the local version of the one in Theorem 13. Indeed the local versions of Assumptions 1, 5 and 9 always hold, i.e., they hold when considering $\bm{x}$ s.t. $\|\bm{x}\|\leq r^{\prime}$, and one can thus repeat that reasonings using local perspectives. $\diamondsuit$

Consider for simplicity the scalar case. Let $y^{*}\mathrel{\mathop{:}}=\frac{1}{N}\sum_{i}A_{i}d_{i}$ and $z^{*}\mathrel{\mathop{:}}=\frac{1}{N}\sum_{i}A_{i}$, so that $\displaystyle x^{*}=\frac{y^{*}}{z^{*}}$. Since $\bm{y}(k+1)=P\bm{y}(k)$ and $\bm{z}(k+1)=P\bm{z}(k)$, given the assumptions on $P$, there exist positive $\alpha_{y},\alpha_{z}$ independent of $\bm{x}(0)$ s.t. $|y_{i}(k)-y^{*}|\leq\alpha_{y}\left(\rho(P)\right)^{k}$ and $|z_{i}(k)-z^{*}|\leq\alpha_{z}\left(\rho(P)\right)^{k}$. The claim thus follows considering that $x_{i}(k)=\frac{y_{i}(k)}{\left[z_{i}(k)\right]_{c}}$ and that, since the elements of $P$ are non negative, all the $z_{i}(k)$ are non smaller than $c$ for all $k\geq 0$ (i.e., the operator $\left[\cdot\right]_{c}$ is always performing as the identity operator). $\diamondsuit$