# Towards Parameter-free Distributed Optimization: a Port-Hamiltonian Approach

Rodrigo Aldana-López, Alessandro Macchelli, Giuseppe Notarstefano,
Rosario Aragüés, Carlos Sagüés 111This work was supported via projects
PID2021-124137OB-I00 and TED2021-130224B-I00 funded by
MCIN/AEI/10.13039/501100011033, by ERDF A way of making Europe and by the
European Union NextGenerationEU/PRTR, by the Gobierno de Aragón under Project
DGA T45-23R, by the Universidad de Zaragoza and Banco Santander, by the
Consejo Nacional de Ciencia y Tecnología (CONACYT-Mexico) with grant number
739841.
R. Aldana-López, R. Aragüés and C. Sagüés are with Departamento de Informática
e Ingeniería de Sistemas (DIIS) and Instituto de Investigación en Ingeniería
de Aragón (I3A), Universidad de Zaragoza, Zaragoza 50018, Spain. (e-mail:
<EMAIL_ADDRESS>{raragues<EMAIL_ADDRESS>
A. Macchelli and G. Notarstefano are with the Department of Electrical,
Electronic and Information Engineering (DEI), University of Bologna, 40136
Bologna, Italy (e-mail: {alessandro.macchelli,
giuseppe.notarstefano}@unibo.it).

###### Abstract

This paper introduces a novel distributed optimization technique for networked
systems, which removes the dependency on specific parameter choices, notably
the learning rate. Traditional parameter selection strategies in distributed
optimization often lead to conservative performance, characterized by slow
convergence or even divergence if parameters are not properly chosen. In this
work, we propose a systems theory tool based on the port-Hamiltonian formalism
to design algorithms for consensus optimization programs. Moreover, we propose
the Mixed Implicit Discretization (MID), which transforms the continuous-time
port-Hamiltonian system into a discrete time one, maintaining the same
convergence properties regardless of the step size parameter. The consensus
optimization algorithm enhances the convergence speed without worrying about
the relationship between parameters and stability. Numerical experiments
demonstrate the method’s superior performance in convergence speed,
outperforming other methods, especially in scenarios where conventional
methods fail due to step size parameter limitations.

##  1 Introduction

Current technological trends lean towards the extended use of machine learning
algorithms with increasing complexity and computing power requirements.
Fortunately, computation and communication capabilities are readily available
in modern devices. For this reason, distributed implementations have emerged
as scalable and robust alternatives to deploy these technologies. Since many
tasks of importance, e.g., the training of neural networks can be posed as
optimization programs, distributed optimization has become a persistently
active area of research in the last years (see [1, 2]). The problem of
consensus optimization is particularly important, where the unconstrained cost
is comprised of the sum of local costs that can be accessed at each device in
the network separately [1].

In the following, we review the literature on algorithms developed for
distributed optimization problems. First, we examine algorithms formulated in
continuous time, leveraging tools from systems theory to guarantee stability.
Subsequently, we discuss methodologies that operate in discrete time,
reflecting more practical communication scenarios. Finally, we mention works
using the Port-Hamiltonian System (PHS) formalism in multi-agent contexts and
parameter-free discretization.

Many algorithms for addressing distributed optimization problems through
continuous-time systems exist in the literature. For instance, a method for
consensus optimization over weighted digraphs was introduced in [3]. A
different approach, employing a proportional-integral protocol to tackle
constrained distributed optimization, was proposed in [4]. The strategy in [5]
aims to resolve constrained distributed convex optimization problems via non-
smooth analysis. Meanwhile, [6] presented a protocol with finite-time
convergence in solving constrained distributed optimization issues. The use of
adaptive control techniques for unconstrained distributed optimization was
explored in [7]. A passivity-based methodology for designing and analyzing
distributed protocol algorithms in the unconstrained scenario was used in [8],
with a similar approach adopted for balanced digraphs in [9]. Furthermore, a
distributed protocol leveraging Newton-Raphson iterations was detailed in
[10]. These contributions underscore the benefits of modeling distributed
optimization algorithms in continuous time, allowing for the application of
systems theory tools, such as Lyapunov theory or passivity. However, given
that communication among agents in digital platforms occurs in discrete time,
research has pivoted towards considering discrete time interactions. A leading
methodology for consensus optimization in discrete time, known as the gradient
tracking method, was introduced in [11], incorporating a dynamic average
consensus component to track the gradient of the objective function. This
method has spawned various adaptations, including those tailored for nonconvex
programs [12] and protocols integrating Newton-Raphson iterations [13].
Extensions to accommodate time-varying and stochastic networks were proposed
in [11, 14], respectively. Efforts to enhance convergence speed are used in
[15, 16]. Techniques such as push-sum and push-pull protocols for message
passing, aimed at addressing digraphs, were applied in [17, 18]. An adaptive
momentum term to boost convergence rates was introduced in [19]. Using a
systems structure perspective, most of these variations can be understood
using a general structural framework as in [20]. Despite their generally
effective performance, these algorithms constrain specific parameters, like
the learning rate. Tuning such parameters is critical, especially in multi-
agent settings where the constraints depend on global information such as the
network topology. For instance, a large learning rate may result in
instability, and a small learning rate may lead to slow convergence.

An alternative line of research for incorporating discrete time communication
involves employing continuous-time systems with event-triggered communication.
One of the first works in this context is [21]. Similar settings, including
quantized information, were explored in [22]. Extensions for quadratic
optimization programs with constraints were considered in [23]. Event-
triggered mechanisms for second-order dynamics agents were investigated in
[24]. The works [25, 26] explored edge-based approaches for the design of the
triggering conditions. A noteworthy contribution in this domain is [27], which
extends the continuous-time version of the gradient tracking algorithm to
include event-triggered communication. While adopting many advantages from
continuous-time methodologies, these approaches also need an upper limit on
the interval between communications, as indicated in [21]. Additionally, they
require the implementation of a discrete time variant for the continuous-time
algorithm, exemplified by [27], where the event-triggering condition is
evaluated at each timestep during experiments. Hence, parameter constraints
cannot be avoided either in this case.

Henceforth, this motivates the need to develop distributed optimization
methods that are not limited by design parameter constraints. An important
breakthrough for parameter-free optimization was presented in [28] where the
Port-Hamiltonian System (PHS) formalism was used to design a continuous-time
gradient flow for unconstrained optimization in a single-agent setting. This
formalism allowed the authors to propose a parameter-free symplectic
discretization. Generally, a PHS is a class of nonlinear systems describing a
wide range of systems of interest, particularly suitable for power-conserving
interconnections between subsystems [29]. This formalism has been proved
useful to describe and analyze multi-agent algorithms for consensus [30] and
distributed control in [31]. Moreover, geometric integration using the
discrete gradient has been shown to preserve the PHS structure [32, 33], with
a precise parallelism between continuous and discrete time. In addition, under
suitable conditions, this type of discretization can ensure convergence
regardless of the step size parameter. For instance, the discrete gradient
method was used in [34] to develop a parameter-free optimization algorithm.
However, despite the advantages of PHS and geometric integration methods,
these tools have yet to be used to solve distributed optimization problems.

The contributions of this work are twofold. First, we propose a PHS framework
for the first time to desig and analyze distributed optimization algorithms.
In particular, we apply it to solve an unconstrained consensus optimization
program, where the convergence proofs are based on PHS tools. Second, our
proposal allows the design in continuous time with a novel parameter-free
discretization, referred to here as Mixed Implicit Discretization (MID), that
enables the preservation of the PHS structure. In this sense, all agents in
the network choose a step size parameter without worrying about the
instability of the overall system. In addition, we discuss how having the step
size as a free design parameter can be further used to accelerate convergence.

This manuscript is organized as follows. First, in Section 2, we provide a
problem statement and some preliminaries on PHS. Section 3 presents our PHS-
based framework for consensus optimization in continuous time, including some
concrete examples. Section 4 discusses parameter-free discretization. In
Section 5, we provide several numerical examples, comparing the proposal’s
performance against other methods in the literature. Conclusions are drawn in
Section 6, and auxiliary results and proofs are placed in the appendices.

Notation: For a set $S$, let $|S|$ represent its cardinality. Given vectors
$\\{\mathbf{x}_{i}\\}_{i\in\mathcal{I}}$ with index set
$\mathcal{I}=\\{1,\dots,{\textsf{N}}\\}$, we denote with
$[\mathbf{x}_{i}]_{i\in\mathcal{I}}=[\mathbf{x}_{1}^{\top},\dots,\mathbf{x}_{\textsf{N}}^{\top}]^{\top}$.
Denote with $\nabla$ the gradient operator. Let
$\mathds{1}_{n},\mathbf{0}_{n\times m},\mathbf{I}_{n}$ denote the column
vector of $n$ ones, the zero $n\times m$ matrix and the $n\times n$ identity
matrix respectively. We write $\mathds{1},\mathbf{0},\mathbf{I}$ when there is
no ambiguity on the dimensions, according to the context. The symbol $\otimes$
represents the Kronecker product. For discrete time sequences
$\\{\mathbf{x}[k]\\}_{k\in\mathbb{N}}$, we often write $\mathbf{x}$ represents
the value of the sequence at an arbitrary time step, $\mathbf{x}^{+}$
represents its value at the next time step. Let $\|\mathbf{x}\|$ represent the
Euclidean norm of a vector $\mathbf{x}$ and $\|\mathbf{X}\|$ represent the
induced Euclidean norm of a matrix $\mathbf{X}$.

##  2 Poblem statement and preliminaries

###  2.1 Problem statement

Consider a group of ${\textsf{N}}\in\mathbb{N}$ agents labeled using an index
set $\mathcal{I}=\\{1,\dots,{\textsf{N}}\\}$, whose goal is to solve a
consensus optimization program of the form:

$\bm{\theta}^{*}=\arg\\!\min_{\bm{\theta}\in\mathbb{R}^{m}}\sum_{i=1}^{\textsf{N}}f_{i}(\bm{\theta})$
(1)

where each $f_{i}:\mathbb{R}^{m}\to\mathbb{R}$ is a function which can only be
accessed at agent $i\in\mathcal{I}$.

######  Assumption 1.

For all $i\in\mathcal{I}$, $f_{i}:\mathbb{R}^{m}\to\mathbb{R}$ is strongly
convex with constant $\mu>0$ (See Appendix A), and has Lipschitz continuous
gradient with constant $L>0$.

Note that Assumption 1 ensures that the program in (1) has a unique optimal
solution $\bm{\theta}^{*}\in\mathbb{R}^{m}$. Moreover, to obtain such optimal
solution, agents cooperate by exchanging information with their neighbors. The
communication network is modeled by a fixed undirected connected graph
$\mathcal{G}=(\mathcal{I},\mathcal{E})$. Here, $\mathcal{I}$ is the vertex
set, where each $i\in\mathcal{E}$ represents a concrete agent, and
$\mathcal{E}\subseteq\mathcal{I}\times\mathcal{I}$ is the edge set where
$(i,j)\in\mathcal{E}$ when agents $i,j\in\mathcal{I}$ can communicate.
Moreover, we denote with $\mathcal{N}_{i}\subseteq\mathcal{I}$ the neighbor
set for agent $i\in\mathcal{I}$.

Often, an algorithm to solve the program in (1) in a distributed fashion often
involves a set of design parameters. The algorithm is parameter-free if the
constraints for the parameters do not require prior information of the network
or the cost function constants $(\mu,L)$ to ensure convergence towards the
optimal solution. This work proposes a parameter-free solution to solve the
program (1) based on the PHS structure.

###  2.2 Port-Hamiltonian systems

A PHS is a particular class of dynamical system that can be described through
the interaction of ports with power-conserving features. In particular, it has
the form:

$\displaystyle\dot{\mathbf{x}}$ $\displaystyle=\mathbf{F}(\mathbf{x})\nabla
H(\mathbf{x})+\bm{B}_{I}(\mathbf{x})\mathbf{u}_{I}+\bm{B}_{C}(\mathbf{x})\mathbf{u}_{C}$
(2)

where $\mathbf{x}(t)\in\mathbb{R}^{n}$ is the state variable. The system in
(2) has four classes of ports. First, a storage port
$(-\dot{\mathbf{x}},\nabla H(\mathbf{x}))$ with scalar storage function
$H:\mathbb{R}^{n}\to\mathbb{R}$ also called the Hamiltonian energy. Second, a
resistive structure associated with the symmetric part of $\mathbf{F}$, which
satisfies
$\mathbf{F}(\mathbf{x})+\mathbf{F}(\mathbf{x})^{\top}\preceq\mathbf{0}$.
Third, an interaction port with $(\mathbf{u}_{I},\mathbf{y}_{I})$. This port
usually represents the interaction of (2) with other systems. Fourth, a
control port $(\mathbf{u}_{C},\mathbf{y}_{C})$. This port is used for feedback
control.

It can be verified that (2) is passive in the sense that:

$\dot{H}(\mathbf{x})\leq\mathbf{y}_{I}^{\top}\mathbf{u}_{I}+\mathbf{y}_{C}^{\top}\mathbf{u}_{C}$

with outputs

$\displaystyle\mathbf{y}_{I}$
$\displaystyle=\mathbf{B}_{I}(\mathbf{x})^{\top}\nabla H(\mathbf{x})$
$\displaystyle\mathbf{y}_{C}$
$\displaystyle=\mathbf{B}_{C}(\mathbf{x})^{\top}\nabla H(\mathbf{x}).$

As will become evident in subsequent sections, the passivity property helps
design control feedback laws and obtain Lyapunov functions based on $H$. We
refer the reader to [29] for more details on PHS.

##  3 PHS-based optimization in continuous time

###  3.1 System proposal

To go towards a parameter-free discrete time algorithm to solve (1), we start
by considering a continuous-time formulation based on a PHS. Assume that each
agent $i\in\mathcal{I}$ stores a state variable
$\mathbf{x}_{i}(t)\in\mathbb{R}^{n}$. The proposed dynamics for
$\mathbf{x}_{i}(t)$ are given as:

$\displaystyle\dot{\mathbf{x}}_{i}$ $\displaystyle=\mathbf{F}_{i}\nabla
H_{i}(\mathbf{x}_{i})+\bm{\phi}_{i}(\nabla
H_{i}(\mathbf{x}_{i}))+\sum_{j\in\mathcal{N}_{i}}\mathbf{M}\left(\nabla
H_{i}(\mathbf{x}_{i})-\nabla H_{j}(\mathbf{x}_{j})\right)$ (3)

Where $\mathbf{F}_{i},\mathbf{M}$ are constant matrices satisfiying
$\mathbf{F}_{i}+\mathbf{F}_{i}^{\top}\preceq\mathbf{0}$ and
$\mathbf{M}+\mathbf{M}^{\top}\preceq\mathbf{0}$. Moreover,
$H_{i}:\mathbb{R}^{n}\to\mathbb{R}$ is the storage function at agent $i$ and
$\bm{\phi}_{i}:\mathbb{R}^{n}\to\mathbb{R}^{n}$ is a nonlinear vector field.
Time dependence for the state was omitted in (3) for readability. For the same
reason, we will also omit state dependence on $H_{i}$ when there is no
ambiguity.

Note that the system in (3) is a PHS since the symmetric part of
$\mathbf{F}_{i}+|\mathcal{N}_{i}|\mathbf{M}$ is negative semi-definite.
Additionally, it includes $|\mathcal{N}_{i}|$ interaction ports, each
connected to one of the neighbors of the agent with input $\nabla
H_{j}(\mathbf{x}_{j})$ and input matrix $-\mathbf{M}$. Moreover, it features a
local feedback control port with input $\bm{\phi}_{i}(\nabla
H_{i}(\mathbf{x}_{i}))$ and input matrix $\mathbf{I}_{n}$. Note that the terms
of the form $\nabla H_{i}-\nabla H_{j}$ in (3) define an agreement error
between neighboring agents and promote consensus.

The remainder of this section aims to show that the PHS structure in (3)
naturally ensures convergence towards the equilibrium. Hence, given that
stability is provided, the remaining degrees of freedom, namely
$\mathbf{F}_{i},\mathbf{M},\bm{\phi}_{i},H_{i}$ can be used to design
equilibrium to coincide with the solution of an optimization program of
interest. To formalize these ideas, we make the following assumption.

######  Assumption 2 (Well posedness).

Let $H_{i},\mathbf{F}_{i},\mathbf{M},{\bm{\phi}}_{i}$ be such that:

  1. 1.

The PHS in (3) has a unique solution for all $t\geq
0,\mathbf{x}_{i}(0)\in\mathbb{R}^{n},i\in\mathcal{I}$.

  2. 2.

There exist unique constant vectors
$\\{\mathbf{x}_{i}^{*}\\}_{i\in\mathcal{I}}$ such that

$\displaystyle\mathbf{0}$ $\displaystyle=\mathbf{F}_{i}\nabla
H_{i}(\mathbf{x}_{i}^{*})+\bm{\phi}(\nabla
H_{i}(\mathbf{x}_{i}^{*}))+\sum_{j\in\mathcal{N}_{i}}\mathbf{M}\left(\nabla
H_{i}(\mathbf{x}_{i}^{*})-\nabla H_{j}(\mathbf{x}_{j}^{*})\right)$ (4)

To show convergence towards the unique equilibrium point satisfying (4), we
write the joint dynamics of all agents incorporating (3) for all
$i\in\mathcal{I}$. Let the joint storage function

$H(\mathbf{x})=\sum_{i\in\mathcal{I}}H_{i}(\mathbf{x}_{i}),$

the joint output $\mathbf{y}=[\nabla H_{i}(\mathbf{x}_{i})]_{i\in\mathcal{I}}$
and the joint input $\bm{\phi}(\mathbf{y})=[\bm{\phi}_{i}(\nabla
H_{i}(\mathbf{x}_{i}))]_{i\in\mathcal{I}}$ with joint state
$\mathbf{x}=[\mathbf{x}_{i}]_{i\in\mathcal{I}}$. Moreover, let
$\mathbf{F}=\text{blockdiag}(\mathbf{F}_{1},\dots,\mathbf{F}_{\textsf{N}})$.
Hence, (3) can be written in compact form as follows.

$\dot{\mathbf{x}}=(\mathbf{F}+\mathbf{L}\otimes\mathbf{M})\nabla
H(\mathbf{x})+\bm{\phi}(\nabla H(\mathbf{x}))$ (5)

where $\mathbf{L}=\text{diag}(\mathbf{A}\mathds{1})-\mathbf{A}$ is the
standard Laplacian matrix of the graph $\mathcal{G}$ associated with the $0-1$
adjacency matrix $\mathbf{A}$.

One of the most important features of (5) inherited from the PHS structure is
its passivity property concerning the input-output pair $(\bm{\phi}(\nabla
H(\mathbf{x})),\mathbf{y})$.

######  Proposition 1.

The trajectories of (5) satisfy

$\dot{H}(\mathbf{x})\leq\mathbf{y}^{\top}\bm{\phi}(\nabla H(\mathbf{x}))$

###### Proof.

Take arbitrary $\mathbf{v}\in\mathbb{R}^{{\textsf{N}}m}$ with
$\mathbf{v}\neq\mathbf{0}$. Recall
$\mathbf{M}+\mathbf{M}^{\top}\preceq\mathbf{0}$ and
$\mathbf{F}+\mathbf{F}^{\top}\preceq\mathbf{0}$ by assumption. Define
$\mathbf{w}=(\mathbf{E}^{\top}\otimes\mathbf{I}_{m})\mathbf{v}$ where
$\mathbf{E}$ is an incidence matrix of $\mathcal{G}$ satisfying
$\mathbf{L}=\mathbf{E}\mathbf{E}^{\top}$. Hence,

$\displaystyle\mathbf{v}^{\top}(\mathbf{F}+\mathbf{L}\otimes\mathbf{M})\mathbf{v}$
$\displaystyle=\mathbf{v}^{\top}\mathbf{F}\mathbf{v}+\mathbf{v}^{\top}(\mathbf{E}\mathbf{E}^{\top}\otimes\mathbf{M})\mathbf{v}$
$\displaystyle\leq\mathbf{v}^{\top}(\mathbf{E}\otimes\mathbf{I}_{m})(\mathbf{I}\otimes\mathbf{M})(\mathbf{E}^{\top}\otimes\mathbf{I}_{m})\mathbf{v}$
$\displaystyle\leq\mathbf{w}^{\top}(\mathbf{I}\otimes\mathbf{M})\mathbf{w}\leq
0$

Therefore, $\mathbf{F}+\mathbf{L}\otimes\mathbf{M}$ has symmetric part
negative semi-definite. Hence, the proposition follows by taking
$\dot{H}=\nabla H^{\top}\dot{\mathbf{x}}$ along (5). ∎

As discussed in Section 2.2, passivity can be used to construct a Lyapunov
function, ensuring convergence to the equilibrium. This idea is formalized in
the following.

######  Corollary 1.

Let Assumption 2 and consider (5). Moreover, assume that $H(\mathbf{x})$ is
strongly convex. Finally, assume that $\bm{\phi}_{i}$ satisfy the incremental
passivity property:

$(\mathbf{u}-\mathbf{v})^{\top}(\bm{\phi}_{i}(\mathbf{u})-\bm{\phi}_{i}(\mathbf{v}))\leq
0$ (6)

for arbitrary $\mathbf{u},\mathbf{v}\in\mathbb{R}^{n}$. with equality only
when $\mathbf{u}=\mathbf{v}$. Then, the unique equilibrium of (5) is globally
asymptotically stable.

###### Proof.

Let $\mathbf{x}^{*}$ denote the equilibrium of (5). To simplify the notation,
let $\nabla H^{*}=\nabla H(\mathbf{x}^{*})$. Set

$V_{\mathbf{x}^{*}}(\mathbf{x})=H(\mathbf{x})-H(\mathbf{x}^{*})-\nabla
H(\mathbf{x}^{*})^{\top}(\mathbf{x}-\mathbf{x}^{*})$ (7)

as a Lyapunov function candidate. Note that $V_{\mathbf{x}^{*}}(\mathbf{x})$
satisfy the required conditions for this purpose. First,
$V_{\mathbf{x}^{*}}(\mathbf{x}^{*})=0$. Second, since $H$ is strongly convex,
then Lemma 1 in Appendix A is used to conclude that
$V_{\mathbf{x}^{*}}(\mathbf{x})$ is always positive and radially unbounded.
Denote with $\mathbf{S}=(\mathbf{F}+\mathbf{L}\otimes\mathbf{M})$. Then,

$\displaystyle\dot{V}_{\mathbf{x}^{*}}=\left(\nabla H-\nabla
H^{*}\right)^{\top}\dot{\mathbf{x}}$ $\displaystyle=\left(\nabla H-\nabla
H^{*}\right)^{\top}(\dot{\mathbf{x}}-\mathbf{S}\nabla H^{*}-\bm{\phi}(\nabla
H^{*}))$ $\displaystyle=\left(\nabla H-\nabla
H^{*}\right)^{\top}(\mathbf{S}\nabla H+\bm{\phi}(\nabla H)-\mathbf{S}\nabla
H^{*}-\bm{\phi}(\nabla H^{*}))$ $\displaystyle=\left(\nabla H-\nabla
H^{*}\right)^{\top}\mathbf{S}\left(\nabla H-\nabla H^{*}\right)$
$\displaystyle+\left(\nabla H-\nabla H^{*}\right)^{\top}(\bm{\phi}(\nabla
H)-\bm{\phi}(\nabla H^{*}))\leq 0$

where we used the fact that $\mathbf{S}+\mathbf{S}^{\top}\preceq\mathbf{0}$
from the proof of Proposition 1. Moreover, we used the incremental passivity
property from (6) in the last step. Furtheremore, note that
$\dot{V}_{\mathbf{x}^{*}}=0$ only when $\nabla H(\mathbf{x})=\nabla
H(\mathbf{x}^{*})$, which occurs only at the equilibrium point. Henceforth,
$\dot{V}_{\mathbf{x}^{*}}<0$ for all $\mathbf{x}\neq\mathbf{x}^{*}$.
Therefore, $\mathbf{x}^{*}$ is a globally asymptotically stable equilibrium. ∎

###  3.2 Matching the equilibrium with the solution of an optimization program

Given that convergence of the PHS system in (3) is ensured by Proposition 1,
it only remains to match the equilibrium equations (4) with the solution of
the consensus optimization program in (1). This can be performed by an
appropriate design of the degrees of freedom we have in
$\mathbf{F}_{i},\mathbf{M},\bm{\phi}_{i},H_{i}$. To provide intuition, we
start with the centralized setting.

######  Example 1 (Centralized optimization).

Consider (1) with ${\textsf{N}}=1$ and $f:=f_{1}$,
$\mathbf{x}:=\mathbf{x}_{1}$. For an equilibrium $\mathbf{x}^{*}$ of (5) to
coincide with $\bm{\theta^{*}}$, the first order optimality constraints must
be met, namely, $\mathbf{0}=\nabla f(\mathbf{x}^{*})$. Comparing with the
equilibrium conditions (4), this condition can be achieved with either one of
the following choices:

  1. a)

$H_{1}(\mathbf{x})=\mathbf{x}^{\top}\mathbf{x}/2$,
$\mathbf{F}_{1}=\mathbf{M}=\mathbf{0}$, $\mathbf{\phi}_{1}(\mathbf{x})=-\nabla
f(\mathbf{x})$.

  2. b)

$H_{1}(\mathbf{x})=f(\mathbf{x})$, $\mathbf{F}_{1}=-\mathbf{I}$,
$\mathbf{M}=\mathbf{0}$, $\mathbf{\phi}_{1}(\mathbf{x})=\mathbf{0}$.

In both cases, the same system results from (3), which coincides with the
standard gradient flow,

$\dot{\mathbf{x}}=-\nabla f(\mathbf{x}).$ (8)

Note that strong convexity of $f$ imply that the incremental passivity
property (6) is satisfied for $\mathbf{\phi}_{1}(\mathbf{x})=-\nabla
f(\mathbf{x})$ due to Proposition 2 in Appendix A. Thus, convergence of
trajectories of the standard gradient flow (8) towards $\mathbf{x}^{*}$
follows from Corollary 1.

Now, we provide an example aiming to solve (1).

######  Example 2 (Unconstrained distributed optimization).

Consider (1) with arbitrary ${\textsf{N}}\in\mathbb{N}$ and connected
$\mathcal{G}$. The first order optimimality conditions for the global optimum
of the program (1) are $\mathbf{0}=\sum_{i\in\mathcal{I}}\nabla
f_{i}(\bm{\theta}^{*})$. However, this condition can be expressed in a more
appropriate distributed fashion as:

$\mathbf{0}=\sum_{i\in\mathcal{I}}\nabla f_{i}(\mathbf{q}_{i}^{*}),\ \
\mathbf{q}_{1}^{*}=\dots=\mathbf{q}_{\textsf{N}}^{*}$ (9)

where we introduce auxiliary variables
$\mathbf{q}_{i}^{*}=\bm{\theta}^{*}\in\mathbb{R}^{m}$. These include a
condition for the gradients and a condition for consensus. Setting
$\mathbf{q}^{*}=[\mathbf{q}_{i}^{*}]_{i\in\mathcal{I}}$, the conditions in (9)
at the equilibrium can be written as:

$\mathbf{0}=(\mathds{1}^{\top}_{\textsf{N}}\otimes\mathbf{I}_{m})\nabla
f(\mathbf{q}^{*}),\ \
\mathbf{0}=(\mathbf{L}\otimes\mathbf{I}_{m})\mathbf{q}^{*}$ (10)

which implies $\mathbf{q}^{*}=(\mathds{1}_{\textsf{N}}\otimes\bm{\theta}^{*})$
and with $f(\mathbf{q}^{*})=\sum_{i\in\mathcal{I}}f_{i}(\mathbf{q}_{i}^{*})$.
Having two sets of constraints motivates the addition of a second set of
auxiliary variables $\mathbf{p}_{i}\in\mathbb{R}^{m}$ and partition
$\mathbf{x}_{i}=[\mathbf{q}_{i}^{\top},\mathbf{p}_{i}^{\top}]^{\top}\in\mathbb{R}^{n}$
with $n=2m$. The PHS equilibrium equations in (4) for each partition will be
used to match either one of the two conditions (10). This can be performed by
setting

$\displaystyle\mathbf{F}_{i}=\mathbf{0}_{n\times n},$
$\displaystyle\quad\mathbf{H}_{i}(\mathbf{x}_{i})=\frac{1}{2}\mathbf{x}_{i}^{\top}\mathbf{x}_{i}$
(11) $\displaystyle\bm{\phi}_{i}(\mathbf{x}_{i})=\begin{bmatrix}-\nabla
f_{i}(\mathbf{q}_{i})\\\ \mathbf{0}_{m}\end{bmatrix},$
$\displaystyle\quad\mathbf{M}=\begin{bmatrix}-1&-1\\\
1&0\end{bmatrix}\otimes\mathbf{I}_{m}$

Henceforth, from (4), the equilibrium for the partitioned states can be
written as

$\displaystyle\mathbf{0}$
$\displaystyle=-(\mathbf{L}\otimes\mathbf{I}_{m})\mathbf{q}^{*}-(\mathbf{L}\otimes\mathbf{I}_{m})\mathbf{p}^{*}-\nabla{f}(\mathbf{q}^{*})$
(12) $\displaystyle\mathbf{0}$
$\displaystyle=(\mathbf{L}\otimes\mathbf{I}_{m})\mathbf{q}^{*}$

from which the optimally conditions in (10) follow by multiplying (12) with
$(\mathds{1}^{\top}_{\textsf{N}}\otimes\mathbf{I}_{m})$ from the left and
using $\mathds{1}^{\top}_{\textsf{N}}\mathbf{L}=\mathbf{0}$. The resulting
system at each agent from (3) is

$\displaystyle\dot{\mathbf{q}}_{i}$
$\displaystyle=-\sum_{i\in\mathcal{I}}(\mathbf{q}_{i}-\mathbf{q}_{i})-\sum_{i\in\mathcal{I}}(\mathbf{p}_{i}-\mathbf{p}_{i})-\nabla
f_{i}(\mathbf{q}_{i})$ (13) $\displaystyle\dot{\mathbf{p}}_{i}$
$\displaystyle=\sum_{i\in\mathcal{I}}(\mathbf{q}_{i}-\mathbf{q}_{i})$

######  Remark 1.

Note that in all previous examples, global convergence is ensured by Corollary
1. The only condition to be verified is existence and uniqueness of solutions
required in Assumption 2. Fortunately, this is easily verified, since
Assumption 1 for the functions $f_{i}$ ensure the Lipschitz property of the
right-hand sides of (8) and (13). Henceforth, existence and uniqueness of
solutions follow in each case.

##  4 PHS-based optimization in discrete time

###  4.1 Discrete gradient-based PHS discretization

The continuous time system in (3) can be implemented in discrete time by
taking advantage of the PHS structure one more time. For instance, we make the
following changes:

$\displaystyle i)$
$\displaystyle\quad\dot{\mathbf{x}}_{i}\to\frac{\mathbf{x}^{+}_{i}-\mathbf{x}_{i}}{\tau}$
(14) $\displaystyle ii)$ $\displaystyle\quad\nabla
H_{i}(\mathbf{x}_{i})\to\overline{\nabla}H_{i}(\mathbf{x}_{i},\mathbf{x}_{i}^{+})$

where we denote with $\mathbf{x}_{i}=\mathbf{x}_{i}[k]$ at an arbitrary time
step $k\in\mathbb{N}$ and $\mathbf{x}_{i}^{+}=\mathbf{x}_{i}[k+1]$ the updated
state at the next time step, and the object $\overline{\nabla}H_{i}$ is
defined later. Note that the change (14)-i) corresponds to the finite
difference approximation of the derivative $\dot{\mathbf{x}}_{i}$ between two
time instants $t=k\tau$ and $t=(k+1)\tau$, where $\tau>0$ represents a step
size parameter. If this was the only change to be applied, the discretization
corresponds to the standard forward Euler method. For instance, applying this
change to (8), leads to the standard gradient descent method

$\mathbf{x}^{+}=\mathbf{x}-\tau\nabla f(\mathbf{x})$

from which $\tau>0$ can be interpreted as a learning rate parameter. In
general, upper bound constraints need to be applied to $\tau$ to ensure
sufficiently a small $\tau$ which leads to convergence.

In order to circumvent any constraint to $\tau$, we introduce (14)-ii), such
that our discretization no longer corresponds to a forward Euler step, where
the symbol $\overline{\nabla}$ denotes the discrete gradient defined as
follows:

######  Definition 1.

[34] Let $g:\mathbb{R}^{n}\to\mathbb{R}$ be a scalar differentiable function.
Hence, $\overline{\nabla}g:\mathbb{R}^{n}\times\mathbb{R}^{n}\to\mathbb{R}$ is
defined by the properties:

$\displaystyle\overline{\nabla}g(\mathbf{u},\mathbf{v})^{\top}(\mathbf{v}-\mathbf{u})=g(\mathbf{u})-g(\mathbf{v})$
(15)
$\displaystyle\lim_{\mathbf{v}\to\mathbf{u}}\overline{\nabla}g(\mathbf{u},\mathbf{v})=\nabla
g(\mathbf{u})$

While there are many possible ways to construct a discrete gradient, in this
work, we consider

$\overline{\nabla}g(\mathbf{u},\mathbf{v})=\int_{0}^{1}\nabla
g((1-s)\mathbf{u}+s\mathbf{v})\text{d}s$

which for the case of quadratic functions
$H_{i}(\mathbf{x}_{i})=\mathbf{x}_{i}^{\top}\mathbf{x}_{i}/2$, results in the
average

$\overline{\nabla}H_{i}(\mathbf{x}_{i},\mathbf{x}_{i}^{+})=\frac{\mathbf{x}_{i}+\mathbf{x}_{i}^{+}}{2}$

With these changes, the discrete time equivalent of the PHS (3) corresponds
to:

$\displaystyle\frac{\mathbf{x}^{+}_{i}-\mathbf{x}_{i}}{\tau}$
$\displaystyle=\mathbf{F}_{i}\overline{\nabla}H_{i}(\mathbf{x}_{i},\mathbf{x}_{i}^{+})+\bm{\phi}_{i}(\overline{\nabla}H_{i}(\mathbf{x}_{i},\mathbf{x}_{i}^{+}))$
(16)
$\displaystyle+\sum_{j\in\mathcal{N}_{i}}\mathbf{M}\left(\overline{\nabla}H_{i}(\mathbf{x}_{i},\mathbf{x}_{i}^{+})-\overline{\nabla}H_{j}(\mathbf{x}_{j},\mathbf{x}_{j}^{+})\right)$

Note that the discrete time system in (16) imposes an implicit equation of
$\mathbf{x}_{i}^{+}$ at each time step. Hence, we require an additional
regularity condition, analogous to existence and uniqueness of solutions in
the second part of Assumption 2 which, similar to [32, 34], can be fulfilled
for $\bm{\phi}_{i},\overline{\nabla}H_{i}$ Lipschitz:

######  Assumption 3.

Given $\\{\mathbf{x}_{i}\\}_{i\in\mathcal{I}}$, there exists unique
$\\{\mathbf{x}_{i}^{+}\\}_{i\in\mathcal{I}}$ complying (16).

As discussed before, a parameter-free algorithm refers to the property of (16)
not to compromise its stability features regardless of the step size $\tau>0$.
In this sense, the advantage of (16) is twofold. First, note that since

$\overline{\nabla}H_{i}(\mathbf{x}_{i}^{*},\mathbf{x}_{i}^{+})=\nabla
H_{i}(\mathbf{x}_{i}^{*})$

by construction, then the equilibrium of (16) is the same as the one from (3),
namely the unique point $\mathbf{x}_{i}^{*}$ complying (4). Second, the
stability proof for (16) is practically the same as (3) regardless of $\tau>0$
as is formalized in the following:

######  Theorem 1.

Let Assumption 2-2), Assumption 3 and consider (16). Moreover, assume that
$H_{i},\bm{\phi}_{i}$ satisfies the conditions described in Corollary 1. Then,
the unique equilibrium of (16) is globally asymptotically stable for any
$\tau>0$.

###### Proof.

First, write the discrete time system (16) in compact form, similarly as done
before with the continuous time PHS from (5):

$\frac{\mathbf{x}^{+}-\mathbf{x}}{\tau}=(\mathbf{F}+\mathbf{L}\otimes\mathbf{M})\overline{\nabla}H(\mathbf{x},\mathbf{x}^{+})+\bm{\phi}(\overline{\nabla}H(\mathbf{x},\mathbf{x}^{+}))$
(17)

with
$\mathbf{x}=[\mathbf{x}_{i}]_{i\in\mathcal{I}},\mathbf{x}^{+}=[\mathbf{x}_{i}^{+}]_{i\in\mathcal{I}}$.
Moreover, denote $\mathbf{S}=(\mathbf{F}+\mathbf{L}\otimes\mathbf{M})$.
Similarly, as in the proof Corollary 1, set $V_{\mathbf{x}^{*}}(\mathbf{x})$
defined in (7) as a Lyapunov function candidate. Then, due to the definition
of the discrete gradient, it follows that

$\displaystyle
V(\mathbf{x}^{+})-V(\mathbf{x})=\overline{\nabla}V(\mathbf{x},\mathbf{x}^{+})^{\top}({\mathbf{x}}^{+}-{\mathbf{x}})$
$\displaystyle=\tau\overline{\nabla}V(\mathbf{x},\mathbf{x}^{+})^{\top}\left(\mathbf{S}\overline{\nabla}H(\mathbf{x},\mathbf{x}^{+})+\bm{\phi}\left(\overline{\nabla}H(\mathbf{x},\mathbf{x}^{+})\right)\right).$

Then, use

$\overline{\nabla}V(\mathbf{x},\mathbf{x}^{+})=\overline{\nabla}H(\mathbf{x},\mathbf{x}^{+})-\overline{\nabla}H(\mathbf{x}^{*},\mathbf{x}^{*})$

since $\overline{\nabla}H(\mathbf{x}^{*},\mathbf{x}^{*})=\nabla
H(\mathbf{x}^{*})$. From this point, the proof follows the same steps as those
in the proof of Corollary 1 with the discrete gradient of $H$ in place of the
gradient to show that $V(\mathbf{x}^{+})-V(\mathbf{x})\leq 0$, with equality
only when $\mathbf{x}=\mathbf{x}^{*}$, concluding the proof. ∎

The formulation presented so far from (16) is already explicitly useful for
parameter-free optimization in a non-distributed setting, as evidenced in the
following example.

######  Example 3.

Consider Example 1 and apply the two design choices of
$\mathbf{F}_{1},\mathbf{M},H_{1},\bm{\phi}_{1}$ on (16). Hence, two possible
discrete time systems for $\mathbf{x}:=\mathbf{x}_{1}$ are obtained:

$\displaystyle\frac{\mathbf{x}^{+}-\mathbf{x}}{\tau}$ $\displaystyle=-\nabla
f\left(\frac{\mathbf{x}^{+}+\mathbf{x}}{2}\right)$ (18)
$\displaystyle\frac{\mathbf{x}^{+}-\mathbf{x}}{\tau}$
$\displaystyle=-\overline{\nabla}f\left({\mathbf{x},\mathbf{x}^{+}}\right)$
(19)

Both systems are ensured to converge to the correct equilibrium regardless of
$\tau>0$. While (19) was discussed in [34], requiring computing the discrete
gradient of the cost function, which can be expensive, (18) uses the usual
gradient, which is more standard and readily available.

To implement the discrete time system in (16), each agent $i\in\mathcal{I}$
need to solve its implicit equations for $\mathbf{x}_{i}^{+}$. This can be
challenging since the implicit equations also depend on the neighbor values
$\mathbf{x}_{j}^{+}$. While a central station may solve for
$\\{\mathbf{x}_{i}^{+}\\}_{i\in\mathcal{I}}$ at each time step, distributed
equation solvers can be used at each time step. However, such an algorithm
requires additional communication overhead. Henceforth, we propose an
additional novel modification to the discretization methodology to remove the
dependence on the neighbors $\\{\mathbf{x}_{i}^{+}\\}_{i\in\mathcal{I}}$ in
the implicit equations at agent $i\in\mathcal{I}$, in a consensus optimization
problem.

###  4.2 The MID method: fully-distributed PHS in discrete time

In the following, we show a simple modification of (16), which results in a
fully distributed algorithm that does not require additional distributed
equation solvers. Aiming to obtain a solution of the program in (1), we set
(16) with the choice
$H_{i}(\mathbf{x}_{i})=\mathbf{x}_{i}^{\top}\mathbf{x}_{i}/2,\mathbf{F}_{i}=\mathbf{0}$
which coincide with the design in (13). These choices result in:

$\displaystyle\frac{\mathbf{x}^{+}_{i}-\mathbf{x}_{i}}{\tau}$
$\displaystyle=\sum_{j\in\mathcal{N}_{i}}\mathbf{M}\left(\frac{\mathbf{x}_{i}+\mathbf{x}_{i}^{+}}{2}-\frac{\mathbf{x}_{j}+\mathbf{x}_{j}^{+}}{2}\right)+\bm{\phi}_{i}\left(\frac{\mathbf{x}_{i}+\mathbf{x}_{i}^{+}}{2}\right)$
(20)

We remove the occurrence of $\mathbf{x}_{j}^{+}$ in the previous equation by
making the following modification, which we refer to as the Mixed Implicit
Discretization (MID):

$\displaystyle\frac{\mathbf{x}^{+}_{i}-\mathbf{x}_{i}}{\tau}$
$\displaystyle=\sum_{j\in\mathcal{N}_{i}}\mathbf{M}\left(\mathbf{x}_{i}^{+}-\mathbf{x}_{j}\right)+\bm{\phi}_{i}\left(\frac{\mathbf{x}_{i}+\mathbf{x}_{i}^{+}}{2}\right),$
(21)

where the disagreement difference is computed with the future state for agent
$i$ and the present state for agent $j$. This modification allows each agent
to locally solve (21) without requiring further communication between its
neighbors besides sharing $\mathbf{x}_{i}$ once at each time step. It is
important to note that this modification has the same equilibrium point as
(20). Henceforth, MID applied to the continuous time PHS in (13) results in:

$\displaystyle\frac{\mathbf{q}_{i}^{+}-\mathbf{q}_{i}}{\tau}$
$\displaystyle=-\sum_{j\in\mathcal{N}_{i}}(\mathbf{q}_{i}^{+}-\mathbf{q}_{j}+\mathbf{p}_{i}^{+}-\mathbf{p}_{j})-\nabla
f_{i}\left(\frac{\mathbf{q}_{i}^{+}+\mathbf{q}_{i}}{2}\right)$ (22)
$\displaystyle\frac{\mathbf{p}_{i}^{+}-\mathbf{p}_{i}}{\tau}$
$\displaystyle=\sum_{j\in\mathcal{N}_{i}}(\mathbf{q}_{i}^{+}-\mathbf{q}_{j})$

Some important features of (22) are described in the following result.

######  Theorem 2.

Let Assumptions 1 and 2-2) hold and consider (22). Let the matrix definitions
for $\bm{\Lambda}(\tau),\mathbf{S}_{r}(\tau)$ in (28) and (34) placed in
Appendix B for readability. Then, given $\tau>0$, the unique equilibrium of
(22) is globally asymptotically stable if the following Linear Matrix
Inequality (LMI) is satisfied:

$\displaystyle\mathbf{P}(\tau):=\begin{bmatrix}\bm{\Lambda}(\tau)&\mathbf{P}_{12}\\\
\mathbf{P}_{12}^{\top}&\mathbf{P}_{22}\end{bmatrix}\succ\mathbf{0}$ (23a)
$\displaystyle\begin{bmatrix}\mathbf{U}&\mathbf{P}_{12}\\\
\mathbf{P}_{12}^{\top}&\mathbf{I}\end{bmatrix}\succeq\mathbf{0},\quad u>0$
(23b)
$\displaystyle\mathbf{P}(\tau)\mathbf{S}_{r}(\tau)+\mathbf{S}_{r}(\tau)^{\top}\mathbf{P}(\tau)+\mathbf{S}_{\gamma}(\tau,\varepsilon)\preceq-u\begin{bmatrix}\mathbf{I},\mathbf{0}\\\
\mathbf{0},\mathbf{0}\end{bmatrix}$ (23c)

for some $\varepsilon\geq 0$ and decision variables
$\mathbf{U},\mathbf{P}_{12},\mathbf{P}_{22}\in\mathbb{R}^{{\textsf{N}}m\times{\textsf{N}}m},u\in\mathbb{R}$.

The proof for Theorem 2 is given in Appendix B.

######  Remark 2.

Note that Theorem 2 does not rely on an assumption similar to Assumption 3.
Instead, we show the existence of a unique solution to the implicit equation
in (22) directly in the proofs.

The following are some consequences of Theorem 2. While LMI formulations are
usually conservative, we provide a convergence result that holds regardless of
the value of $\tau>0$.

######  Corollary 2.

Let the conditions of Theorem 2 hold. Denote with
$\mathbf{x}^{*}=[\mathbf{x}_{i}]_{i\in\mathcal{I}}$ the the unique equilibrium
of (22). Let $\mathbf{A}$ be the $0-1$ adjacency matrix of $\mathcal{G}$ and
$\mathbf{D}=\text{diag}(\mathbf{A}\mathds{1})$ denote the degree matrix. Then,
if

$\mathbf{D}^{2}-\mathbf{A}^{2}\succeq\mathbf{0}$ (24)

it follows that $\mathbf{x}^{*}$ is globally asymptotically stable for all
$\tau>0$. In particular, (24) holds if $\mathcal{G}$ is a cycle or a fully
connected graph of arbitrary size.

From this point, the proofs for all Corollaries are given in Appendix C.

In the following, we show that even if the graph does not comply with (24),
one can find a guaranteed bound for the time step $\tau>0$, which holds for
any connected undirected graphs.

######  Corollary 3.

Let the conditions of Theorem 2 hold. Denote with
$\mathbf{x}^{*}=[\mathbf{x}_{i}]_{i\in\mathcal{I}}$ the unique equilibrium of
(22). If

$\tau<\frac{\mu}{\|\mathbf{D}^{2}-\mathbf{A}^{2}\|}$

then $\mathbf{x}^{*}$ is globally asymptotically stable for all connected
undirected graphs $\mathcal{G}$.

##  5 Numerical examples

In this section, we test the MID-based algorithm in (22) in several scenarios
and compare it with other works in the literature. For the implementation, we
used the function fsolve in MATLAB so that each agent can solve its local
implicit equation at each time step.

We start by testing the proposal’s performance in a practical data analytics
problem, similar to that in [27, Section 6]. Each agent in $i\in\mathcal{I}$
is equipped with $d_{i}\in\mathbb{N}$ points
$\mathbf{m}_{i,1},\dots,\mathbf{m}_{i,d_{i}}\in\mathbb{R}^{m-1}$ with binary
labels $l_{i,1},\dots,l_{i,d_{i}}\in\\{-1,1\\}$. The cost function is then of
the form (1) with

$f_{i}(\bm{\theta})=\sum_{\mu=1}^{d_{i}}\log(1+\exp(-l_{i,\mu}(\tilde{\bm{\theta}}^{\top}\mathbf{m}_{i,\mu}+\theta)))+\frac{C\|\bm{\theta}\|^{2}}{2{\textsf{N}}}$

where $\theta\in\mathbb{R}$ and
$\bm{\theta}=[\tilde{\bm{\theta}}^{\top},\theta]^{\top}$ is the optimization
variable, and $C>0$ is a regularization parameter, making $f_{i}(\bm{\theta})$
strongly convex. For this example, we picked ${\textsf{N}}=10$ agents,
connected in an Erdős-Réyni graph $\mathcal{G}$ with parameter 0.4 [27].
Moreover, we set $m=3,d_{i}=10$ for all agents and $C=0.1$, with randomly
chosen $\\{\mathbf{m}_{i,\mu}\\}_{\mu=1}^{d_{i}}$. In addition, initial
conditions for the state at each agent were chosen randomly.

We start by comparing the convergence speed between MID and explicit forward
Euler discretization by performing several simulations with $200$ values of
$\tau$ uniformly distributed in the interval $[0,10]$. The figure of merit for
the convergence speed is computed as the number of iterations
$K_{B}\in\mathbb{N}$ required such that the condition
$\|\mathbf{q}[k]-(\mathds{1}\otimes\mathbf{I}_{m})\bm{\theta}^{*}\|\leq B$ is
complied from that point forward. Here, we choose $B=10^{-6}$ as an accuracy
bound representing when the error between the system output and true optimal
are sufficiently close for illustrative purposes. The results are shown in
Figure 1, where it is observed that both methods achieve similar convergence
speed for small $\tau$. However, near $\tau\approx 0.77$, Euler discretization
becomes unstable, while the sampling step can be increased further for MID.
This allows us to decrease $K_{B}$ for our method, reaching a minimum near
$\tau\approx 4$. From this point, $K_{B}$ increases with $\tau$. This shows
that allowing bigger values of $\tau$ can increase convergence speed compared
to classical discretization methods, for which increasing $\tau$ leads to
divergence.

Moreover, we compare with other methods in the literature for consensus
distributed optimization. In particular, for discrete time, we compare with
the Discrete Time Gradient Tracking (DTGT) [1] and the DIGing method [11]. In
both cases, a step size parameter $\tau$ is used, which needs to be
sufficiently small. Moreover, we compare with Euler discretization of
continuous time approaches such as the Continuous Time Gradient Tracking
(CTGT)[27], and the coordination algorithm (COOR) from [21]. In both cases, we
assume communication at each time step for fairness in all methods.
Additionally, the time step for these methods is assigned to $\tau$. Finally,
we also compare with Euler discretization of the PHS system in (3).

For comparison, we used quadratic cost functions, which allows us to verify a
less conservative LMI for any graph, as stated in the following result:

######  Corollary 4.

Let the conditions of Theorem 2 hold. Denote with
$\mathbf{x}^{*}=[\mathbf{x}_{i}]_{i\in\mathcal{I}}$ the unique equilibrium of
(22). Given $\tau>0,\mathcal{G}$, as well as
$f_{i}(\mathbf{x}_{i})=\mathbf{x}_{i}^{\top}\mathbf{H}_{i}\mathbf{x}_{i}/2+\mathbf{b}^{\top}_{i}\mathbf{x}_{i}$,
$\mathbf{x}^{*}$ is globally asymptotically stable if the LMI comprised by
(23a), $u>0$ and

$\mathbf{P}(\tau)\mathbf{S}_{r}(\tau)+\mathbf{S}_{r}(\tau)^{\top}\mathbf{P}(\tau)+\mathbf{S}_{h}(\tau)\preceq-u\begin{bmatrix}\mathbf{I},\mathbf{0}\\\
\mathbf{0},\mathbf{0}\end{bmatrix}$ (25)

is satisfied, with $\mathbf{S}_{h}(\tau)$ defined in (37).

Moreover, for this comparison, we looked for the best value of $\tau>0$,
leading to the fastest behavior in each method tested. These values of $\tau$
were found by trial and error by increasing its value until divergence and
keeping the best one in terms of speed. For MID, Euler discretization, DTGT,
CTGT, DIGing, and COOR, the values of $\tau$ resulted in approximately
$3.78,0.71,0.58,0.65,0.05,0.078$ respectively. Trajectories for each case are
shown in Figure 2, where all methods converge to the global optimum, with MID
having a clear advantage in speed.

In contrast, Figure 3 shows the same scenario but now with $\tau=10$, showing
that all other methods diverge while our proposal manages to converge. For
both $\tau=3.78,10$, we verified the feasibility of the LMI in (25), which
guarantees convergence.

To further stress our proposal, we performed simulations with
$\tau=1,10,100,1000$. In all cases, the LMI in (25) was feasible. The
convergence results are shown in Figure 4, where it is observed that even in
the extreme case of $\tau=1000$, our proposal converges despite a slower
convergence rate.

Figure 1: Comparison of the number of iterations $K_{B}$ required to reach
$\|\mathbf{q}[k]-(\mathds{1}\otimes\mathbf{I}_{m})\bm{\theta}^{*}\|\leq B$
with $B=10^{-6}$. The proposed MID allows the reduction of the number of
iterations by increasing $\tau$, with an optimal value around $\tau=4$. In
contrast, increasing $\tau$ for Euler discretization does not feature such
great reduction and leads to divergence at around $\tau=0.77$. Figure 2:
Comparison of MID, Euler discretization as well as DTGT [1], CTGT [27], DIGing
[11] and COOR [21] for a quadratic optimization problem. The fastest values of
$\tau$ were found to be approximately $3.78,0.71,0.58,0.65,0.05,0.078$. The
proposal MID has a clear advantage in terms of speed of convergence. Figure 3:
Comparison of MID, Euler discretization as well as DTGT [1], CTGT [27], DIGing
[11] and COOR [21] for a quadratic optimization problem with $\tau=10$. Figure
4: Comparison of MID for a quadratic optimization problem with
$\tau=1,10,100,1000$.

##  6 Conclusions

In this work, we introduced a novel distributed optimization approach that
achieves parameter-free convergence, a significant departure from traditional
methods reliant on precise parameter selection, such as the learning rate.
Using the Mixed Implicit Discretization (MID) of a continuous-time port-
Hamiltonian system, we demonstrated how this system-based strategy enhances
convergence speeds without worrying about the relationship between parameters
and stability. Our method leveraged the inherent structure of the port-
Hamiltonian system to ensure convergence by aligning the system’s equilibrium
with the optimization problem’s solution, even when transitioning from
continuous to discrete time. Through the MID scheme, we enabled local
resolution of nonlinear implicit equations by network agents, formally proving
the algorithm’s parameter-free convergence capability. Extensive numerical
analyses highlighted our approach’s superiority in scenarios where traditional
methods fail due to parameter constraints, thereby significantly accelerating
convergence speeds.

##  Appendix A Convex functions

This section provides some useful facts used in the rest of the manuscript. We
refer the reader to [35] for standard definitions of convex functions.

######  Proposition 2.

[35, Section 9.1.2] Let $g:\mathbb{R}^{n}\to\mathbb{R}$ be a strongly convex
differentiable function. Then,

$g(\mathbf{v})\geq g(\mathbf{u})+\nabla
g(\mathbf{u})^{\top}(\mathbf{v}-\mathbf{u})+\frac{\mu}{2}\|\mathbf{v}-\mathbf{u}\|^{2}$

for arbitrary $\mathbf{u},\mathbf{v}\in\mathbb{R}^{n}$ and some $\mu>0$
referred to as the strong convexity constant. As a consequence,

$(\mathbf{v}-\mathbf{u})^{\top}(\nabla g(\mathbf{v})-\nabla
g(\mathbf{u}))\geq\mu\|\mathbf{v}-\mathbf{u}\|^{2}$

######  Lemma 1.

Let $H:\mathbb{R}^{n}\to\mathbb{R}$ be a differentiable function. Moreover,
let some fixed $\mathbf{x}^{*}\in\mathbb{R}^{n}$ and
$V_{\mathbf{x}^{*}}:\mathbb{R}^{n}\to\mathbb{R}$ be defined as

$V_{\mathbf{x}^{*}}(\mathbf{x})=H(\mathbf{x})-H(\mathbf{x}^{*})-\nabla
H(\mathbf{x}^{*})^{\top}(\mathbf{x}-\mathbf{x}^{*}).$

Then, If $H$ is strongly convex with constant $\mu>0$,
$V_{\mathbf{x}^{*}}(\mathbf{x})$ satisfy the radially unboundedness condition
given by:

$V_{\mathbf{x}^{*}}(\mathbf{x})\geq\frac{\mu}{2}\|\mathbf{x}-\mathbf{x}^{*}\|^{2}$

###### Proof.

Use Proposition 2 with $g=H$, $\mathbf{v}=\mathbf{x}$ and
$\mathbf{u}=\mathbf{x}^{*}$:

$\displaystyle
V_{\mathbf{x}^{*}}(\mathbf{x})=H(\mathbf{x})-H(\mathbf{x}^{*})-\nabla
H(\mathbf{x}^{*})^{\top}(\mathbf{x}-\mathbf{x}^{*})$
$\displaystyle\geq\left(H(\mathbf{x}^{*})+\nabla
H(\mathbf{x}^{*})^{\top}(\mathbf{x}-\mathbf{x}^{*})+\frac{\mu}{2}\|\mathbf{x}-\mathbf{x}^{*}\|^{2}\right)$
$\displaystyle-H(\mathbf{x}^{*})-\nabla
H(\mathbf{x}^{*})^{\top}(\mathbf{x}-\mathbf{x}^{*})=\frac{\mu}{2}\|\mathbf{x}-\mathbf{x}^{*}\|^{2}$

∎

##  Appendix B Proof of Theorem 2

Through this section, for any variable $\mathbf{x}$, the notations
$\overline{\bullet},\tilde{\bullet}$ act as
$\overline{\mathbf{x}}=(\mathbf{x}^{+}+\mathbf{x})/2,\tilde{\mathbf{x}}=\mathbf{x}^{+}-\mathbf{x}$.
We start by showing the well-posedness of (22).

######  Lemma 2.

Consider the conditions of Theorem 2. Then, for any
$\\{\mathbf{x}_{i}\\}_{i\in\mathcal{I}}$ there always exist a unique
$\\{\mathbf{x}_{i}^{+}\\}_{i\in\mathcal{I}}$ satisfying (22).

###### Proof.

First, recall the partition
$\mathbf{x}_{i}=[\mathbf{q}_{i}^{\top},\mathbf{p}_{i}^{\top}]^{\top}$. In this
context, and with the choices of $\mathbf{M},\bm{\phi}_{i}$ one can solve for
$\mathbf{p}_{i}^{+}$ as

$\mathbf{p}_{i}^{+}=\mathbf{p}_{i}+\tau\left(|\mathcal{N}_{i}|\mathbf{q}_{i}^{+}-\sum_{j\in\mathcal{N}_{i}}\mathbf{q}_{j}\right)$
(26)

Substitute (26) into the update equation for $\mathbf{q}_{i}$ in (22).
Rearranging:

$\mathbf{G}\mathbf{q}_{i}^{+}+\nabla
f_{i}\left(\frac{\mathbf{q}_{i}^{+}+\mathbf{q}_{i}}{2}\right)+\mathbf{c}=0$
(27)

with

$\displaystyle\mathbf{G}$
$\displaystyle=(1/\tau+|\mathcal{N}_{i}|+\tau|\mathcal{N}_{i}|^{2})\mathbf{I}_{m}$
$\displaystyle\mathbf{c}$
$\displaystyle=-\sum_{j\in\mathcal{N}_{i}}\left((1+\tau)\mathbf{q}_{j}+\mathbf{p}_{j}\right)+\mathbf{q}_{i}/\tau$

Note that (27) is equivalent to $\nabla g(\mathbf{q}_{i}^{+})=\mathbf{0}$
where $\nabla$ is taken with respect to $\mathbf{q}_{i}^{+}$ and with

$g(\mathbf{q}_{i}^{+})=\frac{1}{2}(\mathbf{q}_{i}^{+})^{\top}\mathbf{G}\mathbf{q}_{i}^{+}+2f_{i}\left(\frac{\mathbf{q}_{i}^{+}+\mathbf{q}_{i}}{2}\right)+\mathbf{c}^{\top}\mathbf{q}_{i}^{+}$

Note that $\mathbf{G}\succ\mathbf{0}$ and $f_{i}$ is a convex function. Hence,
$g$ is convex as well. Henceforth, there exist a unique solution to $\nabla
g(\mathbf{q}_{i}^{+})=\mathbf{0}$, which coincides with the solution of the
convex optimization program
$\min_{\mathbf{q}_{i}^{+}\in\mathbb{R}^{m}}g(\mathbf{q}_{i}^{+})$. This
implies the existence of a unique $\mathbf{p}_{i}^{+}$ as well from (26). ∎

Now, we write (22) in a more convenient representation.

######  Lemma 3.

Consider the conditions of Theorem 2 and let
${\mathbf{q}}=[\mathbf{q}_{i}]_{i\in\mathcal{I}}$ and
${\mathbf{p}}=[\mathbf{p}_{i}]_{i\in\mathcal{I}}$. Moreover, let
$\mathbf{A},\mathbf{D}=\text{diag}(\mathbf{A}\mathds{1})$ denote the adjacency
and degree matrices of $\mathcal{G}$ and define

$\displaystyle{\mathbf{L}_{m}}=(\mathbf{D}-\mathbf{A})\otimes\mathbf{I}_{m},\quad{\mathbf{Q}}=(\mathbf{D}+\mathbf{A})\otimes\mathbf{I}_{m}/2$
(28)
$\displaystyle\bm{\Lambda}(\tau)=\frac{1}{\tau^{2}}\mathbf{I}+\frac{1}{\tau}{\mathbf{Q}}+{\mathbf{Q}}^{2},$
$\displaystyle\mathbf{S}_{p}(\tau)=\begin{bmatrix}-\bm{\Lambda}(\tau)^{-1}(\frac{1}{\tau}\mathbf{I}+{\mathbf{Q}})\mathbf{L}_{m}&-\frac{1}{\tau}\bm{\Lambda}(\tau)^{-1}\mathbf{L}_{m}\\\
\frac{1}{\tau}\bm{\Lambda}(\tau)^{-1}\mathbf{L}_{m}&-\bm{\Lambda}(\tau)^{-1}{\mathbf{Q}}\mathbf{L}_{m}\end{bmatrix}$
$\displaystyle\mathbf{S}_{r}(\tau)=\begin{bmatrix}-\bm{\Lambda}(\tau)^{-1}(\frac{1}{\tau}\mathbf{L}_{m}+{\mathbf{Q}}\mathbf{L}_{m}+\mathbf{L}_{m}{\mathbf{Q}})&-\frac{1}{\tau}\bm{\Lambda}(\tau)^{-1}\mathbf{L}_{m}\\\
\mathbf{L}_{m}&\mathbf{0}\end{bmatrix}$

Then, $\bm{\Lambda}(\tau)$ is invertible. Moreover, using
$\mathbf{r}=\mathbf{p}-{\mathbf{Q}}\mathbf{q}$, trajectories of (22) satisfy:

$\begin{bmatrix}\tilde{\mathbf{q}}\\\
\tilde{\mathbf{r}}\end{bmatrix}=\mathbf{S}_{r}(\tau)\begin{bmatrix}\overline{\mathbf{q}}\\\
\overline{\mathbf{r}}\end{bmatrix}-\begin{bmatrix}\bm{\Lambda}(\tau)^{-1}/\tau\\\
\mathbf{0}\end{bmatrix}\nabla f(\overline{\mathbf{q}})$ (29)

###### Proof.

First, write (21) as:

$\tilde{\mathbf{x}}/\tau=(\mathbf{D}\otimes\mathbf{M})\mathbf{x}^{+}-(\mathbf{A}\otimes\mathbf{M})\mathbf{x}+\bm{\phi}(\overline{\mathbf{x}})$
(30)

Now, note that the following relations hold:

$\displaystyle\mathbf{x}^{+}$
$\displaystyle=\overline{\mathbf{x}}+\tilde{\mathbf{x}}/2$
$\displaystyle\mathbf{x}$
$\displaystyle=\overline{\mathbf{x}}-\tilde{\mathbf{x}}/2$

replacing these expressions in (30) and rearranging results in:

$(\mathbf{I}/\tau-(\mathbf{A}+\mathbf{D})\otimes\mathbf{M}/2)\tilde{\mathbf{x}}=(\mathbf{L}\otimes\mathbf{M})\overline{\mathbf{x}}+\bm{\phi}(\overline{\mathbf{x}})$

In terms of the partitions of $\mathbf{x}$ and the expression for $\mathbf{M}$
as in (11):

$\displaystyle\begin{bmatrix}\mathbf{I}/\tau+{\mathbf{Q}}&{\mathbf{Q}}\\\
-{\mathbf{Q}}&\mathbf{I}/\tau\end{bmatrix}\begin{bmatrix}\tilde{\mathbf{q}}\\\
\tilde{\mathbf{p}}\end{bmatrix}=\begin{bmatrix}-\mathbf{L}_{m}&-\mathbf{L}_{m}\\\
\mathbf{L}_{m}&\mathbf{0}\end{bmatrix}\begin{bmatrix}\overline{\mathbf{q}}\\\
\overline{\mathbf{p}}\end{bmatrix}-\begin{bmatrix}\mathbf{I}\\\
\mathbf{0}\end{bmatrix}\nabla f(\overline{\mathbf{q}})$ (31)

Note that:

$\begin{bmatrix}\mathbf{I}/\tau+{\mathbf{Q}}&{\mathbf{Q}}\\\
-{\mathbf{Q}}&\mathbf{I}/\tau\end{bmatrix}^{-1}=\bm{\Lambda}(\tau)^{-1}\begin{bmatrix}\mathbf{I}/\tau&-{\mathbf{Q}}\\\
{\mathbf{Q}}&\mathbf{I}/\tau+{\mathbf{Q}}\end{bmatrix}$

Invertibility of $\bm{\Lambda}(\tau)$ comes from the fact that
${\mathbf{Q}}=\mathbf{D}+\mathbf{A}$ is diagonally dominant and thus is
positive semi-definite [36, Section 6.1]. Hence, ${\mathbf{Q}}^{2}$ is
positive semi-definite, being $\bm{\Lambda}(\tau)$ positive definite and thus
invertible. Therefore, from (31) it follows that:

$\displaystyle\begin{bmatrix}\tilde{\mathbf{q}}\\\
\tilde{\mathbf{p}}\end{bmatrix}=\mathbf{S}_{p}(\tau)\begin{bmatrix}\overline{\mathbf{q}}\\\
\overline{\mathbf{p}}\end{bmatrix}-\begin{bmatrix}\bm{\Lambda}(\tau)^{-1}/\tau\\\
\bm{\Lambda}(\tau)^{-1}{\mathbf{Q}}\end{bmatrix}\nabla
f(\overline{\mathbf{q}})$ (32)
$\displaystyle\mathbf{S}_{p}(\tau)=\begin{bmatrix}-\bm{\Lambda}(\tau)^{-1}(\frac{1}{\tau}\mathbf{I}+{\mathbf{Q}})\mathbf{L}_{m}&-\frac{1}{\tau}\bm{\Lambda}(\tau)^{-1}\mathbf{L}_{m}\\\
\frac{1}{\tau}\bm{\Lambda}(\tau)^{-1}\mathbf{L}_{m}&-\bm{\Lambda}(\tau)^{-1}{\mathbf{Q}}\mathbf{L}_{m}\end{bmatrix}$

The rest of the lemma follows from the change of variables

$\displaystyle\begin{bmatrix}{\mathbf{q}}\\\
{\mathbf{r}}\end{bmatrix}=\mathbf{T}\begin{bmatrix}{\mathbf{q}}\\\
{\mathbf{p}}\end{bmatrix}$

with

$\displaystyle\quad\mathbf{T}=\begin{bmatrix}\mathbf{I}&\mathbf{0}\\\
-{\mathbf{Q}}/\tau&\mathbf{I}\end{bmatrix},\quad\mathbf{T}^{-1}=\begin{bmatrix}\mathbf{I}&\mathbf{0}\\\
{\mathbf{Q}}/\tau&\mathbf{I}\end{bmatrix}$

Note that $\bm{\Lambda}(\tau)$ and ${\mathbf{Q}}$ commute so that
${\mathbf{Q}}$ also commutes with $\bm{\Lambda}(\tau)^{-1}$. From this, a
straightforward calculation can verify that:

$\displaystyle\mathbf{T}\begin{bmatrix}\bm{\Lambda}(\tau)^{-1}/\tau\\\
\bm{\Lambda}(\tau)^{-1}{\mathbf{Q}}\end{bmatrix}=\begin{bmatrix}\bm{\Lambda}(\tau)^{-1}/\tau\\\
\mathbf{0}\end{bmatrix}$
$\displaystyle\mathbf{S}_{r}(\tau)=\mathbf{T}\mathbf{S}_{p}(\tau)\mathbf{T}^{-1}$

from which (29) is obtained, completing the proof. ∎

Given these results, we can show Theorem 2.

###### Proof.

(Of Theorem 2) First, use Lemma 3 to write (22) as (29). Denote with
$\mathbf{y}=[\mathbf{q}^{\top},\mathbf{r}^{\top}]^{\top}$. Let
$\mathbf{y}^{*}$ the equilibrium of (29) and set
$\mathbf{e}=\mathbf{y}-\mathbf{y}^{*}$. Hence, it follows from (29) that:

$\displaystyle\tilde{\mathbf{e}}$
$\displaystyle=\mathbf{S}_{r}(\tau)\overline{\mathbf{y}}-\mathbf{B}(\tau)\nabla
f(\overline{\mathbf{q}})$
$\displaystyle=\mathbf{S}_{r}(\tau)\overline{\mathbf{e}}-\mathbf{B}(\tau){\bm{\psi}}(\overline{\mathbf{q}})$

with $\mathbf{B}(\tau)=[\Lambda(\tau)^{-1}/\tau,\mathbf{0}]^{\top}$ and
$\bm{\psi}(\overline{\mathbf{q}})=\nabla f(\overline{\mathbf{q}})-\nabla
f({\mathbf{q^{*}}})$.

Now, assume that (23) holds and set a Lyapunov function candidate

$V(\mathbf{e})=\mathbf{e}^{\top}\mathbf{P}(\tau)\mathbf{e}/2$

with $\mathbf{P}(\tau)$ in (23a). Hence,

$\displaystyle
V(\mathbf{e}^{+})-V(\mathbf{e})=\overline{\mathbf{e}}^{\top}\mathbf{P}(\tau)\tilde{\mathbf{e}}$
$\displaystyle=\overline{\mathbf{e}}^{\top}(\mathbf{P}(\tau)\mathbf{S}_{r}(\tau)+\mathbf{S}_{r}(\tau)^{\top}\mathbf{P}(\tau))\overline{\mathbf{e}}-\overline{\mathbf{e}}^{\top}\mathbf{P}(\tau)\mathbf{B}(\tau){\bm{\psi}}(\overline{\mathbf{q}})$

Now, note that
$\mathbf{P}(\tau)\mathbf{B}(\tau)=[\mathbf{I}/\tau,\mathbf{P}_{12}\bm{\Lambda}(\tau)/\tau]^{\top}$
such that

$\displaystyle-\overline{\mathbf{e}}^{\top}\mathbf{P}(\tau)\mathbf{B}(\tau){\bm{\psi}}(\overline{\mathbf{q}})=-(\overline{\mathbf{q}}-\mathbf{q}^{*})(\nabla
f(\overline{\mathbf{q}})-\nabla f({\mathbf{q^{*}}}))/\tau$ (33)
$\displaystyle-(\overline{\mathbf{r}}-\mathbf{r}^{*})^{\top}\mathbf{P}_{12}^{\top}\bm{\Lambda}(\tau)\bm{\psi}(\overline{\mathbf{q}})/\tau$
$\displaystyle\leq-(\mu{\textsf{N}}/\tau)\|\overline{\mathbf{q}}-\mathbf{q}^{*}\|^{2}+\gamma(\tau)\|\mathbf{P}_{12}(\overline{\mathbf{r}}-\mathbf{r}^{*})\|\|\overline{\mathbf{q}}-\mathbf{q}^{*}\|$

where we used Proposition 2 due to strong convexity of $f$, used

$\|(\nabla f(\overline{\mathbf{q}})-\nabla f({\mathbf{q^{*}}})\|\leq
L\|\overline{\mathbf{q}}-\mathbf{q}^{*}\|$

due to the Lipschitz property of $\nabla f$ and setting

$\gamma(\tau)=(L/\tau)\|\bm{\Lambda}(\tau)\|$

Moreover, one can use Young’s inequality to get

$\displaystyle-(\mu{\textsf{N}}/\tau)\|\overline{\mathbf{q}}-\mathbf{q}^{*}\|^{2}+\gamma(\tau)\|\mathbf{P}_{12}(\overline{\mathbf{r}}-\mathbf{r}^{*})\|\|\overline{\mathbf{q}}-\mathbf{q}^{*}\|$
$\displaystyle\leq\left(-\frac{\mu{\textsf{N}}}{\tau}+\frac{\varepsilon\gamma(\tau)}{2}\right)\|\overline{\mathbf{q}}-\mathbf{q}^{*}\|^{2}+\frac{\gamma(\tau)}{2\varepsilon}\|\mathbf{P}_{12}(\overline{\mathbf{r}}-\mathbf{r}^{*})\|^{2}$
$\displaystyle\leq\overline{\mathbf{e}}^{\top}\mathbf{S}_{\gamma}(\tau,\varepsilon)\overline{\mathbf{e}}$

with

$\mathbf{S}_{\gamma}(\tau,\varepsilon)=\begin{bmatrix}(\gamma(\tau)\varepsilon/2-\mu{\textsf{N}}/\tau)\mathbf{I}&\mathbf{0}\\\
\mathbf{0}&\gamma(\tau)\mathbf{U}/(2\varepsilon)\end{bmatrix}$ (34)

where $\mathbf{P}_{12}^{\top}\mathbf{P}_{12}\preceq\mathbf{U}$, which follows
from (23b) due to a Schur complement. Henceforth:

$\displaystyle V(\mathbf{e}^{+})-V(\mathbf{e})$
$\displaystyle\leq\overline{\mathbf{e}}^{\top}(\mathbf{P}(\tau)\mathbf{S}_{r}(\tau)+\mathbf{S}_{r}(\tau)^{\top}\mathbf{P}(\tau)+\mathbf{S}_{\gamma}(\tau))\overline{\mathbf{e}}$
(35) $\displaystyle\leq-u\|\overline{\mathbf{q}}-\mathbf{q}^{*}\|^{2}$

due to (23c). Hence, $V(\mathbf{e})$ decrease unless
$\overline{\mathbf{q}}=\mathbf{q}^{*}$. However,
$\overline{\mathbf{q}}=\mathbf{q}^{*}$ implies
$\overline{\mathbf{q}}=(\mathds{1}\otimes\mathbf{I}_{m})\bm{\theta}^{*}$ (see
Example 2) such that
$\tilde{\mathbf{r}}=\mathbf{L}_{m}\overline{\mathbf{q}}=\mathbf{0}$ from (29),
meaning that equilibrium was reached for $\tilde{\mathbf{r}}$ as well.
Henceforth, $\mathbf{e}[k]$ converges asymptotically to the origin, completing
the proof. ∎

##  Appendix C Proof of the corollaries

###  C.1 Proof of Corollary 2

Set $\mathbf{P}_{12}=\mathbf{U}=\mathbf{0}$, $\mathbf{P}_{22}=\mathbf{I}/\tau$
and $u=\mu,\varepsilon=0$. With these choices, (23a) is complied since
$\bm{\Lambda}(\tau)$ is positive definite and (23b) is complied trivially.
Moreover,
$\mathbf{W}(\tau)=\mathbf{P}(\tau)\mathbf{S}_{r}(\tau)+\mathbf{S}_{r}(\tau)^{\top}\mathbf{P}(\tau)+\mathbf{S}_{\gamma}(\tau,0)$
for these choices which result in

$\mathbf{W}(\tau)=\begin{bmatrix}-(\frac{1}{\tau}\mathbf{L}_{m}+{\mathbf{Q}}\mathbf{L}_{m}+\mathbf{L}_{m}{\mathbf{Q}})-\mu\mathbf{I}&\mathbf{0}\\\
\mathbf{0}&\mathbf{0}\end{bmatrix}$

Since $\mathbf{L}_{m}$ is positive semi-definite and

$\mathbf{Q}\mathbf{L}_{m}+\mathbf{L}_{m}\mathbf{Q}=(\mathbf{D}^{2}-\mathbf{A}^{2})\otimes\mathbf{I}_{m}\succeq\mathbf{0}$

due to (24), therefore (23c) is complied for these choices for all $\tau>0$.
Note that if $\mathcal{G}$ is a cycle, then $\mathbf{D}=2\mathbf{I}$ and
$\mathbf{D}^{2}-\mathbf{A}^{2}=4\mathbf{I}-(2\mathbf{I}+\mathbf{A}_{2})$ where
$\mathbf{A}_{2}$ is the adjacency matrix of the graph $\mathcal{G}_{2}$ with
edges connecting agents two steps away in the original graph $\mathcal{G}$.
This means that $\mathbf{A}_{2}$ is full of zeros, except for two 1’s per row.
Henceforth,

$\mathbf{D}^{2}-\mathbf{A}^{2}=2\mathbf{I}_{\textsf{N}}-\mathbf{A}_{2}$

is diagonally dominant and thus positive semi-definite [36, Section 6.1].
Similarly, if $\mathcal{G}$ is fully connected, then
$\mathbf{D}^{2}-\mathbf{A}^{2}={\textsf{N}}\mathbf{I}_{\textsf{N}}-{\textsf{A}}^{2}$
and $\mathbf{A}=\mathds{1}\mathds{1}^{\top}-\mathbf{I}_{\textsf{N}}$ so that

$\mathbf{D}^{2}-\mathbf{A}^{2}=({\textsf{N}}-1)\mathbf{I}_{\textsf{N}}-({\textsf{N}}-2)\mathds{1}\mathds{1}^{\top}.$

Similarly as before, $\mathbf{D}^{2}-\mathbf{A}^{2}$ is diagonally dominant
and thus positive semi-definite.

###  C.2 Proof of Corollary 3

Make the same choices as in the proof of Corollary 2, but let $\mathcal{G}$ be
an arbitrary connected undirected graph and arbitrary $u>0$. Henceforth,

$\displaystyle-(\frac{1}{\tau}\mathbf{L}_{m}+{\mathbf{Q}}\mathbf{L}_{m}+\mathbf{L}_{m}{\mathbf{Q}})-\frac{\mu}{\tau}\mathbf{I}$
$\displaystyle\preceq\left(\|{\mathbf{Q}}\mathbf{L}_{m}+\mathbf{L}_{m}{\mathbf{Q}}\|-\frac{\mu}{\tau}\right)\mathbf{I}\preceq-u\mathbf{I}$

by the bound on $\tau$ for some $u>0$. Hence, (23c) is complied due to
${\mathbf{Q}}\mathbf{L}_{m}+\mathbf{L}_{m}{\mathbf{Q}}=(\mathbf{D}^{2}-\mathbf{A}^{2})\otimes\mathbf{I}_{m}$.

###  C.3 Proof of Corollary 4

First note that

$\bm{\psi}(\overline{\mathbf{q}})=\nabla f(\overline{\mathbf{q}})-\nabla
f({\mathbf{q^{*}}})=\mathbf{H}(\overline{\mathbf{q}}-\mathbf{q}^{*})$

with
$\mathbf{H}=\text{blockdiag}(\mathbf{H}_{1},\dots,\mathbf{H}_{\textsf{N}})$.
Now, we follow all steps from the proof of Theorem 2 but instead of (33), do:

$\displaystyle-\overline{\mathbf{e}}^{\top}\mathbf{P}(\tau)\mathbf{B}(\tau){\bm{\psi}}(\overline{\mathbf{q}})=\overline{\mathbf{e}}^{\top}\mathbf{S}_{h}(\tau)\overline{\mathbf{e}}$
(36)

with

$\mathbf{S}_{h}(\tau)=\begin{bmatrix}\mathbf{H}/\tau&\mathbf{0}\\\
\mathbf{P}_{12}^{\top}\bm{\Lambda}(\tau)\mathbf{H}/\tau&\mathbf{0}\end{bmatrix}$
(37)

Therefore, similar to (35):

$\displaystyle V(\mathbf{e}^{+})-V(\mathbf{e})$
$\displaystyle\leq\overline{\mathbf{e}}^{\top}(\mathbf{P}(\tau)\mathbf{S}_{r}(\tau)+\mathbf{S}_{r}(\tau)^{\top}\mathbf{P}(\tau)+\mathbf{S}_{h}(\tau))\overline{\mathbf{e}}$
(38) $\displaystyle\leq-u\|\overline{\mathbf{q}}-\mathbf{q}^{*}\|^{2}$

from which convergence follows using the same arguments at the end of the
proof of Theorem 2.

## References

  * [1] G. Notarstefano, I. Notarnicola, and A. Camisa, “Distributed optimization for smart cyber-physical networks,” _Foundations and Trends in Systems and Control_ , vol. 7, no. 3, p. 253 – 283, 2019, cited by: 36; All Open Access, Green Open Access.
  * [2] T. Yang, X. Yi, J. Wu, Y. Yuan, D. Wu, Z. Meng, Y. Hong, H. Wang, Z. Lin, and K. H. Johansson, “A survey of distributed optimization,” _Annual Reviews in Control_ , vol. 47, pp. 278–305, 2019.
  * [3] B. Gharesifard and J. Cortés, “Distributed continuous-time convex optimization on weight-balanced digraphs,” _IEEE Transactions on Automatic Control_ , vol. 59, no. 3, pp. 781–786, 2013.
  * [4] S. Yang, Q. Liu, and J. Wang, “A multi-agent system with a proportional-integral protocol for distributed constrained optimization,” _IEEE Transactions on Automatic Control_ , vol. 62, no. 7, pp. 3461–3467, 2017.
  * [5] X. Zeng, P. Yi, and Y. Hong, “Distributed continuous-time algorithm for constrained convex optimizations via nonsmooth analysis approach,” _IEEE Transactions on Automatic Control_ , vol. 62, no. 10, pp. 5227–5233, 2017.
  * [6] P. Lin, W. Ren, and J. A. Farrell, “Distributed continuous-time optimization: nonuniform gradient gains, finite-time convergence, and convex constraint set,” _IEEE Transactions on Automatic Control_ , vol. 62, no. 5, pp. 2239–2253, 2017.
  * [7] Z. Li, Z. Ding, J. Sun, and Z. Li, “Distributed adaptive convex optimization on directed graphs via continuous-time algorithms,” _IEEE Transactions on Automatic Control_ , vol. 63, no. 5, pp. 1434–1441, 2018.
  * [8] T. Hatanaka, N. Chopra, T. Ishizaki, and N. Li, “Passivity-based distributed optimization with communication delays using pi consensus algorithm,” _IEEE Transactions on Automatic Control_ , vol. 63, no. 12, pp. 4421–4428, 2018.
  * [9] M. Li, G. Chesi, and Y. Hong, “Input-feedforward-passivity-based distributed optimization over jointly connected balanced digraphs,” _IEEE Transactions on Automatic Control_ , vol. 66, no. 9, pp. 4117–4131, 2021.
  * [10] H. Moradian and S. S. Kia, “A distributed continuous-time modified newton–raphson algorithm,” _Automatica_ , vol. 136, p. 109886, 2022.
  * [11] A. Nedić, A. Olshevsky, and W. Shi, “Achieving geometric convergence for distributed optimization over time-varying graphs,” _SIAM Journal on Optimization_ , vol. 27, no. 4, pp. 2597–2633, 2017.
  * [12] P. Di Lorenzo and G. Scutari, “Next: In-network nonconvex optimization,” _IEEE Transactions on Signal and Information Processing over Networks_ , vol. 2, no. 2, pp. 120–136, 2016.
  * [13] D. Varagnolo, F. Zanella, A. Cenedese, G. Pillonetto, and L. Schenato, “Newton-raphson consensus for distributed convex optimization,” _IEEE Transactions on Automatic Control_ , vol. 61, no. 4, pp. 994–1009, 2016.
  * [14] J. Xu, S. Zhu, Y. C. Soh, and L. Xie, “Convergence of asynchronous distributed gradient methods over stochastic networks,” _IEEE Transactions on Automatic Control_ , vol. 63, no. 2, pp. 434–448, 2017.
  * [15] C. Xi, R. Xin, and U. A. Khan, “Add-opt: Accelerated distributed directed optimization,” _IEEE Transactions on Automatic Control_ , vol. 63, no. 5, pp. 1329–1339, 2017.
  * [16] G. Qu and N. Li, “Harnessing smoothness to accelerate distributed optimization,” _IEEE Transactions on Control of Network Systems_ , vol. 5, no. 3, pp. 1245–1260, 2018.
  * [17] G. Scutari and Y. Sun, “Distributed nonconvex constrained optimization over time-varying digraphs,” _Mathematical Programming_ , vol. 176, no. 1-2, pp. 497–544, 2019.
  * [18] S. Pu, W. Shi, J. Xu, and A. Nedić, “Push-pull gradient methods for distributed optimization in networks,” _IEEE Transactions on Automatic Control_ , vol. 66, no. 1, pp. 1–16, 2020.
  * [19] G. Carnevale, F. Farina, I. Notarnicola, and G. Notarstefano, “Gtadam: Gradient tracking with adaptive momentum for distributed online optimization,” pp. 1436–1448, 2023.
  * [20] B. Van Scoy and L. Lessard, “A universal decomposition for distributed optimization algorithms,” _IEEE Control Systems Letters_ , vol. 6, pp. 3044–3049, 2022.
  * [21] S. S. Kia, J. Cortés, and S. Martínez, “Distributed convex optimization via continuous-time coordination algorithms with discrete-time communication,” _Automatica_ , vol. 55, pp. 254–264, 2015.
  * [22] S. Liu, L. Xie, and D. E. Quevedo, “Event-triggered quantized communication-based distributed convex optimization,” _IEEE Transactions on Control of Network Systems_ , vol. 5, no. 1, pp. 167–178, 2016.
  * [23] Z. Zhao, G. Chen, and M. Dai, “Distributed event-triggered scheme for a convex optimization problem in multi-agent systems,” _Neurocomputing_ , vol. 284, pp. 90–98, 2018.
  * [24] X. Yi, L. Yao, T. Yang, J. George, and K. H. Johansson, “Distributed optimization for second-order multi-agent systems with dynamic event-triggered communication,” in _IEEE Conference on Decision and Control (CDC)_ , 2018, pp. 3397–3402.
  * [25] Y. Kajiyama, N. Hayashi, and S. Takai, “Distributed subgradient method with edge-based event-triggered communication,” _IEEE Transactions on Automatic Control_ , vol. 63, no. 7, pp. 2248–2255, 2018.
  * [26] T. Adachi, N. Hayashi, and S. Takai, “Distributed gradient descent method with edge-based event-driven communication for non-convex optimization,” _IET Control Theory & Applications_, 2021.
  * [27] G. Carnevale, I. Notarnicola, L. Marconi, and G. Notarstefano, “Triggered gradient tracking for asynchronous distributed optimization,” _Automatica_ , vol. 147, p. 110726, 2023.
  * [28] M. Poli, S. Massaroli, A. Yamashita, H. Asama, and J. Park, “Port-hamiltonian gradient flows,” in _ICLR 2020 Workshop on Integration of Deep Neural Models and Differential Equations_ , 2019.
  * [29] A. van der Schaft, _L2-Gain and Passivity Techniques in Nonlinear Control_ , 3rd ed., ser. Communications and Control Engineering. Springer Cham, 2016, originally published as Volume 218 in the series: Lecture Notes in Control and Information Sciences.
  * [30] J. Wei, A. R. Everts, M. K. Camlibel, and A. J. van der Schaft, “Consensus dynamics with arbitrary sign-preserving nonlinearities,” _Automatica_ , vol. 83, pp. 226–233, 2017.
  * [31] A. van der Schaft and B. Maschke, “Port-hamiltonian dynamics on graphs: Consensus and coordination control algorithms,” _IFAC Proceedings Volumes_ , vol. 43, no. 19, pp. 175–178, 2010, 2nd IFAC Workshop on Distributed Estimation and Control in Networked Systems.
  * [32] A. Macchelli, “Control design for a class of discrete-time port-hamiltonian systems,” _IEEE Transactions on Automatic Control_ , pp. 1–8, 2023.
  * [33] ——, “Trajectory tracking for discrete-time port-hamiltonian systems,” _IEEE Control Systems Letters_ , vol. 6, pp. 3146–3151, 2022.
  * [34] E. S. Riis, M. J. Ehrhardt, G. R. W. Quispel, and C.-B. Schönlieb, “A geometric integration approach to nonsmooth, nonconvex optimisation,” _Foundations of Computational Mathematics_ , vol. 22, no. 5, pp. 1351–1394, October 1 2022.
  * [35] S. Boyd and L. Vandenberghe, _Convex Optimization_. Cambridge University Press, 2004.
  * [36] R. A. Horn and C. R. Johnson, _Matrix Analysis_ , 2nd ed. USA: Cambridge University Press, 2012.