# Some remarks on segregation of $k$ species
in strongly competing systems

Flavia Lanzara
Dipartimento di Matematica G. Castelnuovo,
Sapienza Università di Roma,
p.le A. Moro 5 00185 Roma
<EMAIL_ADDRESS>
Eugenio Montefusco
Dipartimento di Matematica G. Castelnuovo,
Sapienza Università di Roma,
p.le A. Moro 5 00185 Roma
<EMAIL_ADDRESS>
###### Abstract

Spatial segregation occurs in population dynamics when $k$ species interact in
a highly competitive way. As a model for the study of this phenomenon, we
consider the competition-diffusion system of $k$ differential equations

$-\Delta u_{i}(x)=-\mu u_{i}(x)\displaystyle\sum_{j\neq i}u_{j}(x)\quad
i=1,...,k$

in a domain $D$ with appropriate boundary conditions. Any $u_{i}$ represents a
population density and the parameter $\mu$ determines the interaction strength
between the populations. The purpose of this paper is to study the geometry of
the limiting configuration as $\mu\longrightarrow+\infty$ on a planar domain
for any number of species. If $k$ is even we show that some limiting
configurations are strictly connected to the solution of a Dirichlet problem
for the Laplace equation.

_2010 Mathematics Subject Classification:_ Primary 35Bxx, 35J47; Secondary
92D25.

_Keywords:_ Pattern formation, Spatial segregation, Strong competition

##  1 Introduction and setting of the problem

When two or more species live in proximity and share the same basic
requirements, they usually compete for resources, habitat or territory. Only
the strongest prevails, driving the weaker competitors to extinction. This is
the principle of competitive exclusion (also known as Gause’s law). One
species wins because its members are more efficient at finding resources,
which leads to an increase in population. This means that a population of
competitors finds less of the same resources and cannot grow at its maximal
capacity [13].

According to the competitive exclusion principle, many competing species
cannot coexist under very strong competition, but when spatial movements are
permitted more than one species can coexist thanks to the segregation of their
habitats. For a theoretical discussion and some experimental results see [14,
18].

From a mathematical viewpoint the determination of the configuration of the
habitat segregation for some populations is an interesting problem which can
be modelled by an optimal (in a suitable sense) partition of a domain; for
example in the papers [8, 9, 10, 11, 22] the problem is studied modelling the
interspecies competition with a large interaction term in an elliptic system
of partial differential equations inspired by classical models in populations
dynamics. In [3, 4, 12] the problem is modelled as a Cauchy problem for a
parabolic system of semilinear partial differential equations describing the
dynamics of the densities of different species. In the evolutive case, in
particular see [12], it is proved that some populations can vanish under the
competition of other species; moreover, in [3, 4] the authors are able to
estimate the number of the long-term surviving populations and other
interesting qualitative properties of the spatial distributions of interacting
populations. Note that also the study of the territoriality, that is how
different groups of the same species divide an area, avoiding to effectively
fight for resources, can be viewed as a habitat segregation produced by
competition (see, for example, [5, 17, 20]); moreover, this kind of
competition is a struggle between competitors having the same features, that
is between perfect competitors.

Others approaches have been considered. In [6], for example, a phase
segregation problem is studied by minimization of integral functionals. The
authors obtain results showing mixing or separation of the two phases,
depending of the strength of the interaction between the two species (that is
depending on some relations between the parameters in the model). The main
difference with our study is that, in [6], the densities have no to satisfy a
boundary datum, so that the optimal partition of the domain is related only to
the geometry of the domain and the interaction between the parameters.

As a model for the study of the segregation phenomena, we consider the
competition-diffusion system of $k$ differential equations

$\begin{array}[]{\\{}{cl}.-\Delta u_{i}(x)=-\mu
u_{i}(x)\displaystyle\sum_{j\neq i}u_{j}(x)&\hbox{ in }D\,,\\\ u_{i}(x)\geq
0&\hbox{ in }D\,,\\\ u_{i}(x)=\phi_{i}(x)&\hbox{ on }\partial
D\,.\end{array}\qquad i=1,...,k\,.$ (1.1)

Here $D\subseteq\mathbf{R}^{n}$ is an open bounded, simply connected domain
with smooth boundary $\partial D$. In this paper we consider the case $n=2$.
We assume that the function $\Phi=(\phi_{1},...,\phi_{k})$ is an admissible
datum that is $\phi_{i}\in W^{1,\infty}(\partial D)$, $\phi_{i}\geq 0$,
$i=1,...,k$, $\phi_{i}\cdot\phi_{j}=0$ a.e. in $\partial D$ for $i\neq j$, the
sets $\\{\phi_{i}>0\\}$ are nonempty, open connected arcs and the function
$\sum_{i=1}^{k}\phi_{i}$ vanishes exactly in $k$ points of $\partial D$ (the
endpoints of the $\phi_{i}$’s supports).

The system (1.1) governs the steady states of $k$ competing species coexisting
in the same domain $D$. Any $u_{i}$ represents a population density and the
parameter $\mu>0$ determines the interaction strength between the populations.
In this model the competition between two species is depicted without direct
reference to the resources they share, rather, it is assumed that the presence
of each population leads to a depression of its competitor’s growth rate.

If $\Phi$ is admissible, the existence of positive solutions of (1.1) for any
positive $\mu$ is proved in [9] using Leray-Schauder degree theory. The
uniqueness is proved in [24], using the sub- and super-solution method.

Let define the class of segregated densities

$\mathcal{U}_{\Phi}=\begin{array}[]{\\{}{rl}\\}U=(u_{1},...,u_{k})\in(H^{1}(D))^{k}\hbox{
: }&u_{i}=\phi_{i}\hbox{ on }\partial D\\\ &u_{i}\geq 0\hbox{ in }D\\\
&u_{i}\cdot u_{j}=0\hbox{ for $i\neq j$ a.e. in }D\\\ \end{array}$

and the class

$\mathcal{S}_{\Phi}=\begin{array}[]{\\{}{rl}\\}U=(u_{1},...,u_{k})\in\mathcal{U}_{\Phi}\hbox{
: }&-\Delta u_{i}\leq 0\hbox{ in }D\\\ &-\Delta\big{(}u_{i}-\sum_{j\neq
i}u_{j}\big{)}\geq 0\hbox{ in }D\\\ \end{array}\,.$

Let $U^{(\mu)}=(u_{1,\mu},...,u_{k,\mu})$ be the solution of (1.1) for every
$\mu>0$. In [9] it is proved that there exists
$\overline{U}=\\{\bar{u}_{1},...,\bar{u}_{k}\\}\in\mathcal{U}$ such that, up
to subsequences, $u_{i,\mu}\longrightarrow\bar{u}_{i}$ in $H^{1}(D)$ and
$\mathcal{S}_{\Phi}$ contains all the asymptotic limits of (1.1) that is
$\overline{U}\in\mathcal{S}_{\Phi}$.

The uniqueness of the limit solution of (1.1) as $\mu\to+\infty$ was proved in
[9] in the case $k=2$ and in [11] in the case of $k=3$ and in dimension $n=2$.
Specifically, the authors prove that the class $\mathcal{S}_{\Phi}$ consists
of one element. In [24] it is proved that $\mathcal{S}_{\Phi}$ consists of one
element also in the case of arbitrary dimension and arbitrary number of
species. A different proof of uniqueness of the limit configuration, based on
the maximum principle and on the qualitative properties of the elements of
$\mathcal{S}_{\Phi}$, is given in [2].

The description of the qualitative properties of the limiting configurations
in the planar case (i.e. $n=2$) was considered in [11] for $k=3$ and in [19]
for $k=4$. The aim of this paper is to describe the geometry of the limiting
configurations in the planar case, for any number of species.

The outline of the paper is as follows. In Section 2 we recall some basic
facts and known results of the class $\mathcal{S}_{\Phi}$, which will be used
later. In Section 3 we study the geometry of the limiting configuration in any
number of species. If $k=2s$, starting from the argument used in [19] for
$k=4$ species, we prove that some limiting configurations are strictly related
to the solution of a Dirichlet problem for the Laplace equation. Our results
rely on the construction of a harmonic function which assumes the value
$\sum_{j=1}^{2s}(-1)^{j}\phi_{j}$ on $\partial D$. This function has an even
number of nodal regions compatible with an alternate sign rule. We emphasize
that this construction cannot be done in the case of odd $k$. This will be the
object of a forthcoming paper. In Section 4 we focus on the case of $k=6$
number of species.

##  2 Basic facts

In this section we recall some basic facts that will play an important role in
our study. Suppose that $D$ is a simply connected domain in $\mathbf{R}^{2}$.
Due to the conformal invariance of the problem, with no loss of generality we
can assume

$D=B(O,1)=\\{x=(x_{1},x_{2})\in\mathbf{R}^{2}:|x|<1\\}$

and consider the class

$\mathcal{S}_{\Phi}=\begin{array}[]{\\{}{rl}\\}&U=(u_{1},...,u_{k})\in(H^{1}(D))^{k}\hbox{
: }u_{i}\geq 0\hbox{ in }D,\>u_{i}=\phi_{i}\hbox{ on }\partial D\\\
&u_{i}\cdot u_{j}=0\hbox{ for }i\neq j,-\Delta u_{i}\leq
0,-\Delta\big{(}u_{i}-\sum_{j\neq i}u_{j}\big{)}\geq 0\hbox{ in
}D\end{array}.$

The study of $\mathcal{S}_{\Phi}$ provides the understanding of the segregated
states of $k$ species induced by strong competition. If
$\Phi=(\phi_{1},...,\phi_{k})$ is an admissible datum then
$\phi_{i},i=1,...,k,$ are positive in their supports, the sets
$\\{\phi_{i}>0\\}\subset\partial D$ are open connected arcs and
$\Phi=\sum_{i=1}^{k}\phi_{i}$ vanishes at exactly $k$ points of $\partial D$,
the endpoints $p_{1},....,p_{k}$ in counter clockwise order.

In the following we will denote by $U$ both the $k-$uple $(u_{1},...,u_{k})$
and the function $\sum_{i=1}^{k}u_{i}$. For any $U\in\mathcal{S}_{\Phi}$
define the nodal regions

$\omega_{i}=\\{p\in D:u_{i}(p)>0\\}\qquad i=1,...,k,$

the multiplicity of a point $p\in\overline{D}$ with respect to $U$:

$m(p)=\\#\\{i:|\omega_{i}\cap B_{r}(p)|>0\quad\forall r>0\\}$

where $B_{r}(p)=\\{q\in\mathbf{R}^{2}:|p-q|<r\\}$, and the interfaces between
two densities

$\Gamma_{ij}=\partial\omega_{i}\cap\partial\omega_{j}\cap\\{p\in D:m(p)=2\\}.$

The element $U\in\mathcal{S}_{\Phi}$ defines exactly $k$ nodal regions and

$\overline{D}=\overline{\omega}_{1}\cup...\cup\overline{\omega}_{k}.$

Let us summarize the basic properties of the elements
$U\in\mathcal{S}_{\Phi}$:

  * (s1)

each $u_{i}\in W^{1,\infty}(\overline{D})$ ([10, Theorem 8.4]). It follows
that $u_{i}\in C(\overline{D})$, $\omega_{i}$ is open and $p\in\omega_{i}$
implies $m(p)=1$. By standard regularity theory for elliptic equations we also
have that $u_{i}\in C^{\infty}(\omega_{i})$;

  * (s2)

each $\omega_{i}$ is connected and each $\Gamma_{ij}$ is either empty or a
connected arc starting from a point $p_{i}\in\partial D$ ([11, Remark 2.1]);

  * (s3)

$u_{i}$ is harmonic in $\omega_{i}$; $u_{i}-u_{j}$ is harmonic in
$D\setminus\cup_{h\neq i,j}\overline{\omega}_{h}$ [10, Proposition 6.3];

  * (s4)

if $p\in D$ satisfies $m(p)=2$, then ([10, Remark 6.4 ])

$\lim_{\omega_{i}\ni q\to p}\nabla u_{i}(q)=-\lim_{\omega_{j}\ni q\to p}\nabla
u_{j}(q);$

  * (s5)

$U\in W^{1,\infty}(\overline{D})$ and if $p\in D$ satisfies $m(p)=2$, then

$|\nabla U(p)|=\lim_{\omega_{i}\ni q\to p}|\nabla
u_{i}(q)|=\lim_{\omega_{j}\ni q\to p}|\nabla u_{j}(q)|\neq 0$

and the set $\\{q:m(q)=2\\}$ is locally a $C^{1}$-curve through $p$ ending
either at points with higher multiplicity, or at the boundary $\partial D$
[10, Lemma 9.4];

  * (s6)

if $p\in D$ satisfies $m(p)\geq 3$, then $|\nabla U(q)|\to 0$, as $q\to p$
[10, Theorem 9.3];

  * (s7)

the set $\\{p\in D:m(p)\geq 3\\}$ consists of a finite number of points [10,
Lemma 9.11];

  * (s8)

if $p\in D$ with $m(p)=h\geq 3$ then there exists $\theta_{0}\in(-\pi,\pi]$
such that

$U(r,\theta)=r^{h/2}\left|\cos\left(\frac{h}{2}(\theta+\theta_{0})\right)\right|+o(r^{h/2})$
(2.1)

as $r\to 0$, where $(r,\theta)$ is a system of polar coordinates around $p$
[10, Theorem 9.6].

######  Remark 2.1.

The asymptotic formula (2.1) describes the behavior of $U$ in a neighborhood
of a multiple point in $D$. As a consequence, at multiple point
$U\in\mathcal{S}_{\Phi}$ shares the angle in equal parts. This property does
not hold true if $U\in\mathcal{S}_{\Phi}$ has a multiple point $p$ on the
boundary $\partial D$.

##  3 Results on $k$ species

Let $U\in\mathcal{S}_{\Phi}$, we define the set of points of multiplicity
greater than or equal to $h\in\mathbf{N}$

$\mathcal{Z}_{h}(U)=\\{p\in\overline{D}:m(p)\geq h\\}.$

The set $\mathcal{Z}_{h}(U)$ consists of a finite number of isolated points
[10, Lemma 9.11 and Theorem 9.13].

######  Proposition 3.1.

Let $k\geq 3$ and $U\in\mathcal{S}_{\Phi}$, then $\mathcal{Z}_{3}(U)$ is
nonempty and does not contain points of multiplicity higher than $k$.

Proof. If $\mathcal{Z}_{3}(U)=\emptyset$ then the interfaces $\Gamma_{ij}$
between any two densities do not intersect in $\overline{D}$. The function
$\Phi=\sum_{i=1}^{k}\phi_{i}$ vanishes in exactly $k$ points and any
$\Gamma_{ij}$ links a point $p_{i}\in\partial D$ to a point $p_{j}\in\partial
D$, $i\neq j$. Therefore, if $k$ is odd there exists at least a point
$p_{\ell}\in\partial D$ which belongs to two interfaces. Then
$p_{\ell}\in\mathcal{Z}_{3}(U)$ and we get a contradiction.

If $k$ is even, since $\phi$ vanishes in exactly $k$ points, then there are
only $k/2+1$ interfaces which are nonempty, with endpoints on the boundary and
do not intersect. This contradicts the fact that $U\in\mathcal{S}_{\Phi}$
defines $k$ nodal regions. It also implies that $\mathcal{Z}_{s}(U)=\emptyset$
for $s>k$.

Let $U\in\mathcal{S}_{\Phi}$. We associate to any $p\in\mathcal{Z}_{3}(U)$ the
number

$i(p)=m(p)-2.$

######  Proposition 3.2.

Let $k\geq 2$. The following relation holds

$k-2=\sum_{p\in\mathcal{Z}_{3}(U)}i(p).$ (3.1)

Proof. If $k=2$ then $\mathcal{Z}_{3}(U)=\emptyset$ and (3.1) is trivially
satisfied. Let $k\geq 3$. First of all we want to point out that the set
$\left(\cup\overline{\Gamma}_{ij}\right)\cup\partial D$, that is the union of
the interfaces between any two species and the boundary of the disk forms a
planar, connected graph, whose vertices are the points in $\mathcal{Z}_{3}(U)$
and the zeros of the boundary datum. From classical arguments in graph theory
(essentially the Euler polyhedral formula, see for example [16, Theorem
1.5.2]) it follows that

$n-m+f=2$

where $n$ is the number of the vertices of the graph, $m$ the number of the
edges and $f$ the number of the faces. In our case we have that
$n=\sharp\\{\mathcal{Z}_{3}(U)\cap D\\}+k$, that is the multiple points in the
disk plus the $k$ zeros of the boundary datum, $m=\sharp\\{\Gamma_{ij}\\}+k$,
the number of the arcs $\Gamma_{ij}\neq\emptyset$ plus $k$ (the number of the
arcs $\\{\phi_{i}>0\\}\subseteq\partial D$), and $f$ is the number of the
nodal regions $\omega_{i}$ plus $1$ so that $f=k+1$. Then, it follows that

$m=n+f-2=2k+\sharp\\{\mathcal{Z}_{3}(U)\cap D\\}-1.$

Now we want to point out that, for any $p\in\mathcal{Z}_{3}(U)\cap D$, the
number $m(p)$ corresponds to the number of the arcs $\Gamma_{ij}$ such that
$p\in\overline{\Gamma}_{ij}$, whileas for the others vertices $p\in\partial D$
(the zeros of $\Phi$ on the boundary) it holds that $m(p)$ is the number of
the arcs $\Gamma_{ij}$ such that $p\in\overline{\Gamma}_{ij}$ augmented of
$1$. Then, recalling that for any vertex $p\in\partial
D\setminus\mathcal{Z}_{3}(U)$ we have $m(p)=2$, we can write that

$\displaystyle\sum_{p\in\mathcal{Z}_{3}(U)}i(p)$
$\displaystyle=\sum_{p\in\mathcal{Z}_{3}(U)}\left[m(p)-2\right]=\sum_{p\in\mathcal{Z}_{3}(U)\cap
D}\left[m(p)-2\right]+\sum_{p\textrm{ vertices on }\partial
D}\left[m(p)-2\right]$ $\displaystyle=\sum_{p\in\mathcal{Z}_{3}(U)\cap
D}\left[\sharp\\{\Gamma_{ij}:p\in\overline{\Gamma}_{ij}\\}-2\right]+\sum_{p\textrm{
vertices on }\partial
D}\left[\sharp\\{\Gamma_{ij}:p\in\overline{\Gamma}_{ij}\\}+1-2\right]$
$\displaystyle=\sum_{\begin{subarray}{c}p\textrm{ vertices}\\\ \textrm{of the
graph}\end{subarray}}\left[\sharp\\{\Gamma_{ij}:p\in\overline{\Gamma}_{ij}\\}-2\right]+k$

since the vertices on the boundary are exactly $k$. Then, it follows that

$\displaystyle\sum_{p\in\mathcal{Z}_{3}(U)}i(p)$
$\displaystyle=\sum_{\begin{subarray}{c}p\textrm{ vertices}\\\ \textrm{of the
graph}\end{subarray}}\left[\sharp\\{\Gamma_{ij}:p\in\overline{\Gamma}_{ij}\\}-2\right]+k=\sum_{\begin{subarray}{c}p\textrm{
vertices}\\\ \textrm{of the
graph}\end{subarray}}\left[\sharp\\{\Gamma_{ij}:p\in\overline{\Gamma}_{ij}\\}\right]-2n+k$
$\displaystyle=2(m-k)-2n+k=2(f-2)-k=k-2$

since the edges are the union of the interfaces $\Gamma_{ij}$ and of the $k$
arcs $\\{\phi_{i}>0\\}\subseteq\partial D$, and summing over the vertices any
edge is counted twice.

######  Remark 3.3.

Note that, following [16], the set $\Gamma=\cup{\overline{\Gamma}}_{ij}$, that
is the union of the interfaces between two species, the zeros of the boundary
datum $\Phi$ and the points in $\mathcal{Z}_{3}(U)$ form a tree, since it is a
planar, acyclic, connected graph. The zeros of $\Phi$ are the leaves and the
multiple points in $\mathcal{Z}_{3}(U)$ are the other vertices. From classical
arguments in graph theory it follows that the number of the arcs composing
$\Gamma$ is $(k+\\#\\{\mathcal{Z}_{3}(U)\\}-1)$.

######  Remark 3.4.

Let $k\geq 3$. We want to point out that identity (3.1) implies that

$1\leq\\#\\{\mathcal{Z}_{3}(U)\\}\leq k-2\,.$

In the following we assume that the number of species $k$ is even, that is
$k=2\,s$, $s>1$. Then we can define a harmonic function having opposite signs
on adjacent nodal regions and strictly connected to $U\in\mathcal{S}_{\Phi}$.

Consider the boundary value problem

$\begin{array}[]{\\{}{cl}.-\Delta\psi=0&\hbox{ in }D\\\ \psi=\phi&\hbox{ on
}\partial D\,.\end{array}$ (3.2)

######  Proposition 3.5.

Let $\psi_{a}$ be the solution of (3.2) with boundary datum
$\Phi^{a}=\sum_{j=1}^{2s}(-1)^{j}\phi_{j}$ and let $U\in\mathcal{S}_{\Phi}$.
1\. If $U=|\psi_{a}|$ then for any $p\in\mathcal{Z}_{3}(U)\cap D$ we have that
$m(p)$ is even.
2\. If there exists $p\in\mathcal{Z}_{3}(U)\cap D$ such that $m(p)$ is odd
then $U\neq|\psi_{a}|$.

Proof. 1\. Let $U=|\psi_{a}|$ and $p\in\mathcal{Z}_{3}(U)$. If $p\in D$ then
$U(p)=0$ and $\nabla U(p)=(0,0)$ ([10, Theorem 9.3]). It follows that $p$ is a
critical point for $\psi_{a}$ at level $0$. By standard theory of harmonic
functions the zero set of $\psi_{a}$ around a critical point at level 0 is
made by (at least) $4$ half-lines, meeting with equal angles. We infer that
locally around $p$ the function $\psi_{a}$ defines $q$ nodal components with
$q\geq 4$ and $q$ is even because $\psi_{a}$ has alternate positive or
negative sign on adjacent sets.

2\. Let $p\in\mathcal{Z}_{3}(U)$ such that $m(p)=q$ with odd $q$. If $p\in D$
and $U=|\psi_{a}|$ then $p$ is a critical point for $\psi_{a}$ at level $0$.
We infer that locally around $p$ the function $\psi_{a}$ defines $q$ nodal
region, and $q$ is odd. This is a contradiction because $\psi_{a}$ has
alternate positive or negative sign on adjacent sets.

######  Theorem 3.6.

Let $\Phi=(\phi_{1},...,\phi_{2s})$ be an admissible datum. The harmonic
function $\psi_{a}$ which solves (3.2) with boundary datum
$\Phi^{a}=\sum_{j=1}^{2s}(-1)^{j}\phi_{j}$ possesses at most $s-1$ critical
points $q$ in $D$ such that $\psi_{a}(q)=0$.

Proof. Since $\Phi$ is an admissible datum, the solution of (3.2) with
boundary datum $\Phi^{a}$ vanishes at exactly $2s$ points on $\partial D$,
each arc $\\{\phi_{j}>0\\}\subset\partial D$ is connected and $\psi_{a}$ has
different signs on adjacent arcs.

Since a harmonic function does not admit closed level lines, the set
$\Gamma=\\{x\in D:\psi_{a}(x)=0\\}$ has no closed loop. We infer that
$\psi_{a}$ has alternate positive or negative sign on $\rho$ sets, with
$s+1\leq\rho\leq 2s$: the nodal components of $\psi_{a}$ (see figure 1 in the
case $k=6$). The zero set of the harmonic function $\psi_{a}$ around a
critical point at level 0 is made by (at least) $4$ half-lines, meeting with
equal angles. We infer that locally around each critical point at level $0$
the function $\psi_{a}$ defines $\ell$ nodal components with
$4\leq\ell\leq\rho$ and $\ell$ is even because $\psi_{a}$ has alternate
positive or negative sign on adjacent sets.

Suppose that $\psi_{a}$ has $r$ critical points $q_{1},...,q_{r}$ in $D$ such
that $\psi_{a}(q_{j})=0$, $j=1,...,r$, $r\geq 1$. Then
$q_{j}\in\mathcal{Z}_{3}(|\psi_{a}|)$ with $m(q_{j})\geq 4$, $j=1,...,r$. If
$\psi_{a}$ defines $2s$ nodal regions i.e. $\rho=2s$ then the function
$U=|\psi_{a}|$ is an element of $\mathcal{S}_{\Phi}$. From (3.1),

$2s-2=\sum_{p\in\mathcal{Z}_{3}(U)}i(p)\geq\sum_{j=1}^{r}i(q_{j})=\sum_{j=1}^{r}m(q_{j})-2r\geq
4r-2r=2r$

that is $r\leq s-1$. If $\psi_{a}$ defines $\rho$ nodal regions with $\rho<2s$
then, by repeating the same argument in the proof of formula (3.1), we obtain
that

$2\rho-2-2s=\sum_{p\in\mathcal{Z}_{3}(|\psi_{a}|)}i(p).$ (3.3)

Then,

$2\rho-2-2s=\sum_{p\in\mathcal{Z}_{3}(|\psi_{a}|)}i(p)\geq\sum_{j=1}^{r}i(q_{j})=\sum_{j=1}^{r}m(q_{j})-2r\geq
4r-2r=2r.$

We infer that $r<\rho-s-1<s-1$.

$+$$-$$+$$-$$-$$+$$a)$

$+$$-$$+$$-$$-$$+$$b)$

$+$$-$$+$$-$$-$$+$$c)$

$+$$-$$+$$-$$-$$+$$d)$

Figure 1: The level set $\Gamma=\\{x\in D:\psi_{a}(x)=0\\}$ and the nodal
components of $\psi_{a}$, where the function $\psi_{a}$ solves (3.2) with
boundary datum $\Phi^{a}=\sum_{j=1}^{6}(-1)^{j}\phi_{j}$. Three different
situations can occur: the function $\psi_{a}$ has $4$ nodal regions (figure
$a)$), $5$ nodal regions (figure $b)$) or $6$ nodal regions (figures $c)$ and
$d)$).

If the solution $\psi_{a}$ to (3.2) with alternate sign on adjacent arcs has
exactly $s-1$ critical points at level zero we can say much more.

######  Proposition 3.7.

Let $\Phi=(\phi_{1},...,\phi_{2s})$ be an admissible boundary datum and
suppose that the harmonic function $\psi_{a}$, solution to (3.2) with boundary
datum $\Phi^{a}=\sum_{j=1}^{2s}(-1)^{j}\phi_{j}$, has $q_{1},...,q_{s-1}$
critical points in $D$ such that $\psi_{a}(q_{i})=0$, $i=1,...,s-1$. Then
$q_{1},...,q_{s-1}$ are $4$-points for the function
$U=|\psi_{a}|\in\mathcal{S}_{\Phi}$.

Proof. The zero set of the harmonic function $\psi_{a}$ around a critical
point at level 0 is made by (at least) $4$ half-lines, meeting with equal
angles. We infer that locally around each $q_{i}$, $i=1,...,s-1$, the function
$\psi_{a}$ defines $\ell$ nodal components with $\ell\geq 4$. The function
$\psi_{a}$ defines $\rho\leq 2s$ nodal regions. From formula (3.3),

$2\rho-2-2s=\sum_{p\in\mathcal{Z}_{3}(|\psi_{a}|)}(m(p)-2)\geq\sum_{j=1}^{s-1}(m(q_{j})-2)\geq
2(s-1).$

We infer that $\rho\geq 2s$, hence $\rho=2s$.

The function $U=|\psi_{a}|$ is nonnegative, satisfies the boundary datum and
has exactly $2s$ nodal regions. This function generates an element of
$\mathcal{S}_{\Phi}$, with boundary datum $\Phi$ and
$\\{q_{1},...,q_{s-1}\\}\in\mathcal{Z}_{3}(U)$. Hence, we have

$\\{q_{1},...,q_{s-1}\\}\subseteq\mathcal{Z}_{3}(U),\quad m(q_{j})=m_{j}$

with $m_{j}\geq 4$, $j=1,...,s-1$. From (3.1) we get

$\sum_{j=1}^{s-1}i(q_{j})\leq 2s-2$

or equivalently

$\sum_{j=1}^{s-1}m_{j}-2(s-1)\leq 2s-2\iff\sum_{j=1}^{s-1}m_{j}\leq 4(s-1)\,.$

Since $m_{j}\geq 4$ we deduce that $m_{j}=4$, $j=1,...,s-1$.

######  Proposition 3.8.

Let $\Phi=(\phi_{1},...,\phi_{2s})$ be an admissible datum. If
$U\in\mathcal{S}_{\Phi}$ possesses a $2s$-point $a_{U}$ in $\overline{D}$ then
$U=|\psi_{a}|$, where $\psi_{a}$ is the harmonic function such that
$\psi_{a}=\Phi^{a}=\sum_{j=1}^{2s}(-1)^{j}\phi_{j}$ on $\partial D$. If
$a_{U}\in D$ then $a_{U}$ is a critical point for $\psi_{a}$ at zero level.

Proof. Let $U=(u_{1},...,u_{2s})$ be an element of $\mathcal{S}_{\Phi}$ with a
$2s$-point $a_{U}$, then there exist $2s-1$ or $2s$ arcs connecting $a_{U}$ to
any of the isolated zeros of the boundary datum (note that there are $2s-1$
arcs if and only if $a_{U}\in\partial D$). Then the function
$\psi_{a}=\sum_{j=1}^{2s}(-1)^{j}u_{j}$ is harmonic in $D\setminus\\{a_{U}\\}$
(see [10, Proposition 6.3]), moreover $\psi_{a}$ is bounded so, by Schwarz’s
removable singularity principle (see [21, Proposition 11.1]) $\psi_{a}$ is
harmonic in $D$ and, by construction, $U=|\psi_{a}|$ in $D$ and
$U=|\Phi^{a}|=\Phi$ on $\partial D$.

If $a_{U}\in\mathcal{Z}_{3}(U)\cap D$ then, for the property (s6),
$U(a_{U})=0$ and $\nabla U(a_{U})=(0,0)$. Hence $\psi_{a}(a_{U})=0$ and
$\nabla\psi_{a}(a_{U})=(0,0)$.

We can generalized Proposition 3.8 as follows.

######  Proposition 3.9.

Let $\Phi=(\phi_{1},...,\phi_{2s})$ be an admissible datum. Let
$U\in\mathcal{S}_{\Phi}$ be such that $\mathcal{Z}_{3}(U)$ contains only
points with even multiplicity in $\overline{D}$. Then $U=|\psi_{a}|$, where
$\psi_{a}$ is the solution of the Dirichlet problem (3.2) with boundary datum
$\psi_{a}=\Phi^{a}=\sum_{j=1}^{2s}(-1)^{j}\phi_{j}$ on $\partial D$. Moreover
any $p\in\mathcal{Z}_{3}(U)\cap D$ is a critical point for $\psi_{a}$ at zero
level.

Proof. Let $U$ be an element of $\mathcal{S}_{\Phi}$ such that
$\mathcal{Z}_{3}(U)=\\{q_{1},...,q_{r}\\}$ with $m(q_{i})=2\ell_{i}\geq 4$,
$i=1,...,r$. Then, for every $i=1,...,r$, there exist
$s_{1}^{(i)},...,s_{2\ell_{i}}^{(i)}\in\\{1,...,2s\\}$ such that

i)
$\overline{\omega}_{s^{(i)}_{1}}\cap...\cap\overline{\omega}_{s^{(i)}_{2\ell_{i}}}=\\{q_{i}\\}$;

ii) $\sum_{j=1}^{2\ell_{i}}(-1)^{j}u_{s_{j}}$ is harmonic in
$\dot{\Omega}_{i}\setminus\\{q_{i}\\}$ where
$\Omega_{i}=\overline{\omega}_{s^{(i)}_{1}}\cup...\cup\overline{\omega}_{s^{(i)}_{2\ell_{i}}}$
and $\dot{\Omega}_{i}$ denotes the interior of $\Omega_{i}$.

Keeping in mind that $u_{i}=0$ in $\omega_{j}$, $i\neq j$, then
$\psi_{a}=\sum_{j=1}^{2s}(-1)^{j}u_{j}$ is harmonic in
$\dot{\Omega}_{i}\setminus\\{q_{i}\\}$, $i=1,...,r$. It is easy to see that
$\dot{\Omega}_{1}\cup...\cup\dot{\Omega}_{r}=D\setminus\\{q_{1},...,q_{r}\\}$.
Then the function $\psi_{a}=\sum_{j=1}^{2s}(-1)^{j}u_{j}$ is harmonic in
$D\setminus\\{q_{1},...,q_{r}\\}$ (see [10, Proposition 6.3]), moreover
$\psi_{a}$ is bounded so, by Schwarz’s removable singularity principle (see
[21, Proposition 11.1]) $\psi_{a}$ is harmonic in $D$ and, by construction,
$U=|\psi_{a}|$ in $D$ and $U=|\Phi^{a}|=\Phi$ on $\partial D$.

If $q\in\mathcal{Z}_{3}(U)\cap D$ then $U(q)=0$ and $\nabla U(q)=(0,0)$. Hence
$\psi_{a}(q)=0$ and $\nabla\psi_{a}(q)=(0,0)$.

The next Proposition gives conditions on the admissible datum $\Phi$ such that
the limiting configuration has a point with multiplicity $2s$ in $D$.

######  Proposition 3.10.

Let $\Phi=(\phi_{1},...,\phi_{2s})$ be an admissible boundary datum and let
$\psi_{a}$ be the solution to (3.2) with boundary datum
$\Phi^{a}=\sum_{j=1}^{2s}(-1)^{j}\phi_{j}$. Let $p\in{D}$ such that $\Phi^{a}$
satisfies the conditions

$\displaystyle\int_{\partial
D}\Phi^{a}\left(\frac{\zeta+p}{\overline{p}\zeta+1}\right)T_{j}(\zeta_{1})ds_{\zeta}=0,\quad
j=0,1,...,s-1$ (3.4) $\displaystyle\int_{\partial
D}\Phi^{a}\left(\frac{\zeta+p}{\overline{p}\zeta+1}\right)\zeta_{2}U_{j-1}(\zeta_{1})ds_{\zeta}=0,\quad
j=1,...,s-1$ (3.5)

with $\zeta=(\zeta_{1},\zeta_{2})$, where $T_{j}$ and $U_{j}$ denote the
Chebyshev polynomials of the first and second kind, respectively. Then $p$ is
a $2s$-point for the function $U=|\psi_{a}|\in\mathcal{S}_{\Phi}$.

Proof. We introduce the transformation

$x=R_{p}(\zeta)=\frac{\zeta+p}{\overline{p}\zeta+1}\,.$ (3.6)

Here we identificate the complex numbers $x=x_{1}+ix_{2}$ and
$\zeta=\zeta_{1}+i\zeta_{2}$ with the points $(x_{1},x_{2})$ and
$(\zeta_{1},\zeta_{2})\in\mathbf{R}^{2}$, respectively. $R_{p}$ is a conformal
map which maps the unit disk $\overline{D}$ into itself such that
$R_{p}(\partial D)=\partial D$ and $R_{p}(0)=p$. Set

$\widetilde{\Psi}_{a}(\zeta)=\psi_{a}(R_{p}(\zeta)),\qquad\widetilde{\Phi}^{a}(\zeta)=\Phi^{a}(R_{p}(\zeta))\,.$
(3.7)

Then $\widetilde{\Psi}_{a}$ solves the problem

$\begin{array}[]{\\{}{cl}.-\Delta\widetilde{\Psi}_{a}=0&\hbox{ in }D\,,\\\
\widetilde{\Psi}_{a}=\widetilde{\Phi}^{a}&\hbox{ on }\partial D\,.\\\
\end{array}$ (3.8)

Introducing a system of polar coordinats $(r,\theta)$, we can write the
Fourier expansion of $\widetilde{\Psi}_{a}$

$\widetilde{\Psi}_{a}(r,\theta)=\frac{A_{0}}{2}+\sum_{j=1}^{\infty}(A_{j}\cos(j\theta)+B_{j}\sin(j\theta))r^{j}\,,\qquad\zeta=(r,\theta)\,.$
(3.9)

Let $T_{j}$ and $U_{j}$ denote the Chebychev polynomials of the first and
second kind, respectively. Keeping in mind the representation
$T_{j}(\cos(\theta))=\cos(j\theta)$ and
$U_{j-1}(\cos(\theta))=\sin(j\theta)/\sin(\theta)$ ([1, 22.3.15-16]) from the
conditions (3.4)-(3.5) we get

$\displaystyle
A_{j}:=\frac{1}{\pi}\int_{-\pi}^{\pi}\widetilde{\Phi}^{a}(e^{i\theta})\cos(j\theta)d\theta=0,\quad
j=0,1,...,s-1$ (3.10) $\displaystyle
B_{j}:=\frac{1}{\pi}\int_{-\pi}^{\pi}\widetilde{\Phi}^{a}(e^{i\theta})\sin(j\theta)d\theta=0,\quad
j=1,...,s-1.$ (3.11)

It follows that, around the origin,

$\widetilde{\Psi}_{a}(r,\theta)=\sum_{j=s}^{\infty}(A_{j}\cos(j\theta)+B_{j}\sin(j\theta))r^{j}\,,\qquad\zeta=(r,\theta)\,.$

We have $(A_{s},B_{s})\neq(0,0)$. Indeed, if not, let
$(A_{\nu},B_{\nu})\neq(0,0)$ where $\nu>s$ is the index of the first nonzero
Fourier component. Then there would be $2\nu$ arcs starting form the origin,
on which $\widetilde{\Psi}_{a}$ vanishes. Since a harmonic function does not
admit closed level lines, this contradicts the fact that
$\widetilde{\Phi}^{a}$ has exactly $2s$ zeroes.

Therefore $U(x)=|\psi_{a}(x)|=|\widetilde{\Psi}_{a}(R^{-1}_{p}(x))|$ is
nonnegative, satisfies the boundary datum $\Phi$ and has exactly $2s$ nodal
regions. This function generates an element of $\mathcal{S}_{\Phi}$, with
datum $\Phi$ and the $2s$-point $p$ (see also [11, Lemma 3.2]).

Conversely,

######  Proposition 3.11.

Suppose that $\Phi=(\phi_{1},...,\phi_{2s})$ is an admissible datum and the
function $U\in\mathcal{S}_{\Phi}$ generates a configuration with a $2s$-point
in $p\in{D}$. Then $\Phi^{a}=\sum_{j=1}^{2s}(-1)^{j}\phi_{j}$ satisfies the
conditions (3.4)-(3.5).

Proof. For Proposition 3.8, we have $U=|\psi_{a}|$ with $\psi_{a}$ solution of
(3.2) with boundary datum $\Phi^{a}=\sum_{j=1}^{2s}(-1)^{j}\phi_{j}$. Then the
function $\widetilde{\Psi}_{a}(\zeta)=\psi_{a}(R_{p}(\zeta))$, with $R_{p}$ in
(3.6), solves (3.8), it is given by

$\widetilde{\Psi}_{a}(\zeta)=\frac{1-|\zeta|^{2}}{2\pi}\int_{\partial
D}\frac{\widetilde{\Phi}^{a}(\eta)}{|\zeta-\eta|^{2}}ds_{\eta}$ (3.12)

and it belongs to $C^{2}(D)\cap C^{0}(\overline{D})$ (cf. [15, (2.27)]).

On the other hand in the Fourier expansion (3.9) of $\widetilde{\Psi}_{a}$ we
have

$A_{0}=A_{j}=B_{j}=0,\qquad j=1,...,s-1$

and $(A_{s},B_{s})\neq(0,0)$. The result follows from (3.10)-(3.11).

######  Proposition 3.12.

Let $\Phi=(\phi_{1},...,\phi_{2s})$ be an admissible boundary datum,
$\Phi^{a}=\sum_{j=1}^{2s}(-1)^{j}\phi_{j}$ and $p\in{D}$. Conditions
(3.4)-(3.5) are equivalent to

$\displaystyle\int_{\partial
D}\Phi^{a}\left(\frac{\zeta+p}{\overline{p}\zeta+1}\right)\zeta_{1}^{j-h}\zeta_{2}^{h}ds_{\zeta}=0,\qquad
h=0,...,j;j=0,...,s-1$ (3.13)

with $\zeta=(\zeta_{1},\zeta_{2})$.

Proof. From the representations

$\displaystyle\sin(j\theta)=\sum_{\begin{subarray}{c}{h=1}\\\
{h\>\>odd}\end{subarray}}^{j}(-1)^{\frac{h-1}{2}}\begin{pmatrix}j\\\
h\end{pmatrix}\cos^{j-h}\theta\sin^{h}\theta,\qquad\cos(j\theta)=\sum_{{\begin{subarray}{c}{h=0}\\\
{h\>\>even}\end{subarray}}}^{j}(-1)^{\frac{h}{2}}\begin{pmatrix}j\\\
h\end{pmatrix}\cos^{j-h}\theta\sin^{h}\theta$

we infer that, setting $\zeta_{1}=\cos(\theta),\zeta_{2}=\sin(\theta)$,

$\displaystyle\zeta_{2}U_{j-1}(\zeta_{1})=\sum_{\begin{subarray}{c}{h=1}\\\
{h\>\>odd}\end{subarray}}^{j}(-1)^{\frac{h-1}{2}}\begin{pmatrix}j\\\
h\end{pmatrix}\zeta_{1}^{j-h}\zeta_{2}^{h},\quad\zeta_{1}^{2}+\zeta_{2}^{2}=1,\quad
j\geq 1$ (3.14) $\displaystyle
T_{j}(\zeta_{1})=\sum_{{\begin{subarray}{c}{h=0}\\\
{h\>\>even}\end{subarray}}}^{j}(-1)^{\frac{h}{2}}\begin{pmatrix}j\\\
h\end{pmatrix}\zeta_{1}^{j-h}\zeta_{2}^{h},\quad\zeta_{1}^{2}+\zeta_{2}^{2}=1,j\geq
0.$ (3.15)

If conditions (3.13) are satisfied, from the representations (3.14)-(3.15), we
deduce that even conditions (3.4)-(3.5) are.

Conversely, let conditions (3.4)-(3.5) hold. For $j=0$, since $T_{0}\equiv 1$,
we get

$\int_{\partial
D}\Phi^{a}\left(\frac{\zeta+p}{\overline{p}\zeta+1}\right)ds_{\zeta}=0.$
(3.16)

For $j=1$, since $T_{1}(\zeta_{1})=\zeta_{1}$ and $U_{0}\equiv 1$, we get

$\int_{\partial
D}\Phi^{a}\left(\frac{\zeta+p}{\overline{p}\zeta+1}\right)\zeta_{1}ds_{\zeta}=\int_{\partial
D}\Phi^{a}\left(\frac{\zeta+p}{\overline{p}\zeta+1}\right)\zeta_{2}ds_{\zeta}=0.$

For $j=2$, from the relations $T_{2}(\zeta_{1})=2\zeta_{1}^{2}-1$ and
$U_{1}(\zeta_{1})=2\zeta_{1}$, and (3.16) we get

$0=\int_{\partial
D}\Phi^{a}\left(\frac{\zeta+p}{\overline{p}\zeta+1}\right)(2\zeta_{1}^{2}-1)ds_{\zeta}=2\int_{\partial
D}\Phi^{a}\left(\frac{\zeta+p}{\overline{p}\zeta+1}\right)\zeta_{1}^{2}ds_{\zeta}$
$0=\int_{\partial
D}\Phi^{a}\left(\frac{\zeta+p}{\overline{p}\zeta+1}\right)\zeta_{2}(2\zeta_{1})ds_{\zeta}=2\int_{\partial
D}\Phi^{a}\left(\frac{\zeta+p}{\overline{p}\zeta+1}\right)\zeta_{1}\zeta_{2}ds_{\zeta}$

Since $\zeta_{1}^{2}+\zeta_{2}^{2}=1$, we also get

$0=\int_{\partial
D}\Phi^{a}\left(\frac{\zeta+p}{\overline{p}\zeta+1}\right)\zeta_{2}^{2}ds_{\zeta}.$

We proceed by induction on $j$. Suppose that conditions (3.13) are valid for
$h=0,...,j$ with $j<s-1$. We prove that

$\displaystyle\int_{\partial
D}\Phi^{a}\left(\frac{\zeta+p}{\overline{p}\zeta+1}\right)\zeta_{1}^{j+1-h}\zeta_{2}^{h}ds_{\zeta}=0,\qquad
h=0,...,j+1.$ (3.17)

Using the inverse formula ([7, p.412])

$\zeta_{1}^{j+1}=2^{-j}{\sum^{j+1}_{\begin{subarray}{c}{i=0}\\\
{j+1-i\>\>even}\end{subarray}}}\\!\\!\\!\\!\\!\\!^{\prime}\begin{pmatrix}j+1\\\
i\end{pmatrix}T_{i}(\zeta_{1}),\quad j\geq 0$ (3.18)

where the prime at the sum symbol means that the first term (at $i=0$) is to
be halved unless it is skipped, and conditions (3.4) we get

$\int_{\partial
D}\Phi^{a}\left(\frac{\zeta+p}{\overline{p}\zeta+1}\right)\zeta_{1}^{j+1}ds_{\zeta}=0,$

that is (3.17) for $h=0$. If $h=2\ell$ with $0<\ell\leq(j+1)/2$, we have,
keeping in mind (3.18) and (3.4),

$\displaystyle\int_{\partial
D}\Phi^{a}\left(\frac{\zeta+p}{\overline{p}\zeta+1}\right)\zeta_{1}^{j+1-2\ell}\zeta_{2}^{2\ell}ds_{\zeta}=\int_{\partial
D}\Phi^{a}\left(\frac{\zeta+p}{\overline{p}\zeta+1}\right)\zeta_{1}^{j+1-2\ell}(1-\zeta_{1}^{2})^{\ell}ds_{\zeta}$
$\displaystyle=\sum_{i=0}^{\ell}\begin{pmatrix}\ell\\\
i\end{pmatrix}(-1)^{\ell-i}\int_{\partial
D}\Phi^{a}\left(\frac{\zeta+p}{\overline{p}\zeta+1}\right)\zeta_{1}^{j+1-2i}ds_{\zeta}=0.$

If $h=2\ell+1$, with $0\leq\ell\leq j/2$, we can write (3.17) as

$\int_{\partial
D}\Phi^{a}\left(\frac{\zeta+p}{\overline{p}\zeta+1}\right)\zeta_{1}^{j-2\ell}\zeta_{2}^{2\ell+1}ds_{\zeta}=\sum_{i=0}^{\ell}\begin{pmatrix}\ell\\\
i\end{pmatrix}(-1)^{\ell-i}\int_{\partial
D}\Phi^{a}\left(\frac{\zeta+p}{\overline{p}\zeta+1}\right)\zeta_{2}\zeta_{1}^{j-2i}ds_{\zeta}.$
(3.19)

From the relations $T_{0}(\zeta_{1})=U_{0}(\zeta_{1})$,
$T_{1}(\zeta_{1})=2^{-1}U_{1}(\zeta_{1})$,
$T_{\kappa}(\zeta_{1})=2^{-1}(U_{\kappa}(\zeta_{1})-U_{\kappa-2}(\zeta_{1}))$,
$\kappa\geq 2$, and the inverse formula (3.18), we infer that
$\zeta_{1}^{j-2i}$ can be expressed as a linear combination of the polynomials
$U_{\kappa}(\zeta_{1})$, $0\leq\kappa\leq j-2i$. Hence, conditions (3.5)
implies that the integrals in the right-hand side of (3.19) vanish.

For example, if $k=4$ the conditions (3.13) reduce to

$\int_{\partial
D}\Phi^{a}\left(\frac{\zeta+p}{\overline{p}\zeta+1}\right)ds_{\zeta}=0,\qquad\int_{\partial
D}\Phi^{a}\left(\frac{\zeta+p}{\overline{p}\zeta+1}\right)\zeta_{j}ds_{\zeta}=0,\>\>j=1,2,$

which were obtained in [19, Proposition 3.12]).

As a consequence of Propositions 3.10, 3.11 and 3.12 we deduce the following
necessary and sufficient conditions such that $p\in D$ is a point with
multipicity $2s$.

######  Theorem 3.13.

Let $\Phi=(\phi_{1},...,\phi_{2s})$ be an admissible boundary datum and let
$\psi_{a}$ be the solution to (3.2) with boundary datum
$\Phi^{a}=\sum_{j=1}^{2s}(-1)^{j}\phi_{j}$. The point $p\in{D}$ is a
$2s-$point for the function $U=|\psi_{a}|$ if and only if $\Phi^{a}$ satisfies
conditions (3.13).

##  4 Results on $6$ species

In this section we consider the case of $6$ competing species. As a
consequence of Proposition 3.2 we show that $5$ limiting configurations are
possible.

######  Proposition 4.1.

Let $U\in\mathcal{S}_{\Phi}$, then only one of the following statement is
satisfied
i. $\mathcal{Z}_{3}(U)$ consists of one point $a_{U}\in\overline{D}$ such that
$m(a_{U})=6$,
ii. $\mathcal{Z}_{3}(U)$ consists of two points $a_{U},b_{U}\in\overline{D}$,
$a_{U}\neq b_{U}$, with $m(a_{U})=m(b_{U})=4$,
iii. $\mathcal{Z}_{3}(U)$ consists of two points $a_{U},b_{U}\in\overline{D}$
with $m(a_{U})=3$ and $m(b_{U})=5$,
iv. $\mathcal{Z}_{3}(U)$ consists of three different points
$a_{U},b_{U},c_{U}\in\overline{D}$ with $m(a_{U})=4$, $m(b_{U})=m(c_{U})=3$,
v. $\mathcal{Z}_{3}(U)$ consists of four different points
$a_{U},b_{U},c_{U},d_{U}\in\overline{D}$ with
$m(a_{U})=m(b_{U})=m(c_{U})=m(d_{U})=3$.

Proof. The set $\mathcal{Z}_{3}(U)$ is nonempty and contains at most $4$
points. If $\mathcal{Z}_{3}(U)=\\{a_{U}\\}$ then, from (3.1), we deduce that
$m(a_{U})=6$ (cf. figure 2). If $\mathcal{Z}_{3}(U)=\\{a_{U},b_{U}\\}$ then,
from (3.1), we deduce that $m(a_{U})+m(b_{U})=8$, that is ii. or iii. (cf.
figure 3). If $\mathcal{Z}_{3}(U)=\\{a_{U},b_{U},c_{U}\\}$ then, from (3.1),
we deduce that $m(a_{U})+m(b_{U})+m(c_{U})=10$, that is iv. (cf. figure 4 on
the left). If $\mathcal{Z}_{3}(U)=\\{a_{U},b_{U},c_{U},d_{U}\\}$ then, from
(3.1), we deduce that $m(a_{U})+m(b_{U})+m(c_{U})+m(d_{U})=12$, that is v.
(cf. figure 4 on the right).

Figure 2: Configurations with one 6-point inside $D$ (on the left) and on
$\partial D$ (on the right).

Figure 3: Configurations when $\mathcal{Z}_{3}(U)$ consists of two points
$a_{U}\neq b_{U}$, with $m(a_{U})=m(b_{U})=4$ (on the left) and $m(a_{U})=3$
and $m(b_{U})=5$ (on the right).

Figure 4: Configurations when $\mathcal{Z}_{3}(U)$ consists of three points
(on the left) and four points (on the right).

In the case of $6$ species, Theorem 3.13 can be formulated as follows.

######  Proposition 4.2.

Let $\Phi=(\phi_{1},...,\phi_{6})$ be an admissible datum, let
$U\in\mathcal{S}_{\Phi}$ and let $\psi_{a}$ be the solution to (3.2) with
boundary datum $\Phi^{a}=\sum_{j=1}^{6}(-1)^{j}\phi_{j}$. The function $U$ has
a 6-point $p\in{D}$ if and only if the boundary datum
$\Phi=(\phi_{1},...,\phi_{6})$ satisfies the conditions

$\displaystyle\displaystyle\sum_{j=1}^{6}(-1)^{j}\int_{\partial
D}\phi_{j}\left(\frac{\zeta+p}{\overline{p}\,\zeta+1}\right)ds_{\zeta}$
$\displaystyle=0\,,$ (4.1)
$\displaystyle\displaystyle\sum_{j=1}^{6}(-1)^{j}\int_{\partial
D}\phi_{j}\left(\frac{\zeta+p}{\overline{p}\,\zeta+1}\right)\zeta_{r}ds_{\zeta}$
$\displaystyle=0\,,\qquad r=1,2$ (4.2)
$\displaystyle\displaystyle\sum_{j=1}^{6}(-1)^{j}\int_{\partial
D}\phi_{j}\left(\frac{\zeta+p}{\overline{p}\,\zeta+1}\right)\zeta_{1}^{2}ds_{\zeta}$
$\displaystyle=0\,,$ (4.3)
$\displaystyle\displaystyle\sum_{j=1}^{6}(-1)^{j}\int_{\partial
D}\phi_{j}\left(\frac{\zeta+p}{\overline{p}\,\zeta+1}\right)\zeta_{1}\zeta_{2}ds_{\zeta}$
$\displaystyle=0.$ (4.4)

Moreover the function $U=|\psi_{a}|$.

Proof. Conditions (4.1)-(4.4) are easily obtained by assuming $s=3$ in (3.13).

######  Proposition 4.3.

Let $\Phi$ be an admissible datum and let $\psi_{a}$ be the solution to (3.2)
with boundary datum $\Phi^{a}=\sum_{j=1}^{6}(-1)^{j}\phi_{j}$. Conditions
(4.1)-(4.4) are equivalent to

$\psi_{a}(p)=0,\quad\nabla\psi_{a}(p)=(0,0),\quad H\psi_{a}(p)=\bf 0$ (4.5)

where $H\psi_{a}=\\{\partial^{2}_{x_{i}x_{j}}\psi_{a}\\}_{i,j=1,2}$ denotes
the hessian matrix of the function $\psi_{a}$.

Proof. We introduce the transformation (3.6) and define $\widetilde{\Psi}_{a}$
and $\widetilde{\Phi}^{a}$ according to (3.7). We have

$\psi_{a}(p)=\psi_{a}(R_{p}(0))=\widetilde{\Psi}_{a}(0);$
$\nabla_{x}\psi_{a}(p)=\begin{pmatrix}(1-|p|^{2})^{-1}&0\\\
0&(1-|p|^{2})^{-1}\end{pmatrix}\nabla_{\zeta}\widetilde{\Psi}_{a}(0)$

and, assuming $p=(p_{1},p_{2})$,

$H_{x}\psi_{a}(p)=\frac{1}{(1-|p|^{2})^{2}}H_{\zeta}\widetilde{\Psi}_{a}(0)+\frac{2}{1-|p|^{2}}\left(\partial_{\zeta_{1}}\widetilde{\Psi}_{a}(0)\begin{pmatrix}p_{1}&p_{2}\\\
p_{2}&-p_{1}\end{pmatrix}+\partial_{\zeta_{2}}\widetilde{\Psi}_{a}(0)\begin{pmatrix}-p_{2}&p_{1}\\\
p_{1}&p_{2}\end{pmatrix}\right)\,.$

By construction, $\widetilde{\Psi}_{a}$ solves (3.8). By the Poisson integral
formula (3.12) we deduce that

$\widetilde{\Psi}_{a}(0)=\frac{1}{2\pi}\int_{\partial
D}\frac{\widetilde{\Phi}^{a}(\eta)}{|\eta|^{2}}ds_{\eta}=\frac{1}{2\pi}\int_{\partial
D}\widetilde{\Phi}^{a}(\eta)ds_{\eta}=\frac{1}{2\pi}\int_{\partial
D}\Phi^{a}(R_{p}(\eta))ds_{\eta}.$

On the other hand, by direct differentation of the Poisson integral, we get

$\frac{\partial}{\partial\zeta_{j}}\widetilde{\Psi}_{a}(0)=\frac{1}{\pi}\int_{\partial
D}\widetilde{\Phi}^{a}(\eta)\eta_{j}ds_{\eta}=\frac{1}{\pi}\int_{\partial
D}\Phi^{a}(R_{p}(\eta))\eta_{j}ds_{\eta},\quad j=1,2$
$\frac{\partial^{2}}{\partial^{2}\zeta_{j}}\widetilde{\Psi}_{a}(0)=-\frac{2}{\pi}\int_{\partial
D}\widetilde{\Phi}^{a}(\eta)ds_{\eta}+\frac{4}{\pi}\int_{\partial
D}\widetilde{\Phi}^{a}(\eta)\eta_{j}^{2}ds_{\eta},\qquad j=1,2$
$\frac{\partial^{2}}{\partial\zeta_{1}\partial\zeta_{2}}\widetilde{\Psi}_{a}(0)=\frac{4}{\pi}\int_{\partial
D}\widetilde{\Phi}^{a}(\eta)\eta_{1}\eta_{2}ds_{\eta}.$

The equivalence between (4.1)-(4.4) and (4.5) easily follows.

######  Proposition 4.4.

Let $\Phi=(\phi_{1},...,\phi_{6})$ be an admissible datum and suppose that the
related harmonic function $\psi_{a}$, solution to (3.2) with boundary datum
$\phi^{a}$, has two critical points $p,q\in D$ such that
$\psi_{a}(p)=\psi_{a}(q)=0$, then $p,q$ are 4-points for the function
$U=|\psi_{a}|\in\mathcal{S}_{\Phi}$.

Proof. By standard theory of harmonic functions the zero set of $\psi_{a}$
around a critical point at level 0 is made by (at least) $4$ half-lines,
meeting with equal angles. We infer that locally around $p$ the function
$\psi_{a}$ defines $k_{p}$ nodal components with $k_{p}\geq 4$ and $k_{p}$ is
even because $\psi_{a}$ has alternate positive or negative sign on adjacent
sets. Analogously around $q$ the function $\psi_{a}$ defines $k_{q}$ nodal
components with $k_{q}\geq 4$ and $k_{q}$ is even. By Proposition 4.1 we infer
that $k_{p}=k_{q}=4$. Hence the function $U=|\psi_{a}|$ is nonnegative,
satisfies the boundary datum $\Phi$, has exactly $6$ nodal regions and
generates an element of $\mathcal{S}_{\Phi}$, with $p,q\in\mathcal{Z}_{3}(U)$.

Conversely, as a direct application of Proposition 3.9 we have the following
result.

######  Proposition 4.5.

Let $\Phi=(\phi_{1},...,\phi_{6})$ be an admissible datum. Let
$U\in\mathcal{S}_{\Phi}$ such that $\mathcal{Z}_{3}(U)$ contains two points
$p,q$ in ${D}$ with multiplicity $4$. Then $U=|\psi_{a}|$, where $\psi_{a}$ is
the harmonic function such that
$\psi_{a}=\Phi^{a}=\sum_{j=1}^{6}(-1)^{j}\phi_{j}$ on $\partial D$. Moreover
$p,q$ are critical points for $\psi_{a}$ at zero level.

##  5 Conclusion and open problems

Propositions 4.4 and 4.5 show that the limiting configurations, which elements
of $\mathcal{Z}_{3}(U)$ have only even multiplicity, are closely connected to
harmonic solutions of (3.2) with alternate boundary datum. Since such
solutions have to satisfy some integral conditions (see (3.13)), it follows
that the most probable segregation configurations have only points in
$\mathcal{Z}_{3}(U)$ with odd multiplicity (see also Remark 3.13 in [19]).

In fact we think that the most probable configurations in nature are those
with points with multiplicity three. But this means that odd multiple points
(with multiplicity greater than 3) have to verify some conditions in order to
belong to subsets of admissible data with nontrivial codimension. Now this is
an open problem, since our techniques are strongly related to the properties
of the harmonic functions, and these functions are related only to critical
points with even multiplicity.

An open problem we wish to handle is to compare results by reaction-diffusion
systems in populations competition, as in [3, 4, 9, 11, 12, 19, 23], and
results obtained by differential games theory contained in [20, 17] (see also
the references therein), in order to reach a better understanding of the
territoriality of the competing species or groups (see also [5]).

### Acknowledgment

The authors are very grateful to the anonymous referees for their insightful
comments and helpful suggestions.

## References

  * [1] Abramowitz, M. & Stegun, I., _Handbook of Mathematical functions_ , Dover Publ. 1970.
  * [2] Arakelyan, V. & Bozorgnia, F., Uniqueness of Limiting Solution to a Strongly Competing System, _Elect. J. Diff. Eqs._ 2017 (2017).
  * [3] Berestycki, H. & Zilio, A., Predators-Prey models with Competition. Part I: Existence, Bifurcation and Qualitative Properites, _Comm. Contemp. Math._ 20 (2018), 1850010.
  * [4] Berestycki, H. & Zilio, A., Predators-prey models with competition: the emergence of packs and territoriality, _The American Naturalist_ 193 (2019), 436-446.
  * [5] Börger, L., Dalziel, B.D. & Fryxell, J.M., Are there general mechanisms of animal home range behaviour? A review and prospects for future research, _Ecol. Lett._ 11 (2008), 637–650.
  * [6] Cicalese M., De Luca L., Novaga M., & Ponsiglione M., Ground States of a Two Phase Model with Cross and Self Attractive Interactions, _SIAM J. Math. Anal._ 48 (2016), 3412-3443.
  * [7] Cody W.J., A survey of practical rational and polynomial approximation of functions, _SIAM Rev._ 12 (1970), 400-423.
  * [8] Conti, M., Terracini, S. & Verzini, G., Nehari’s problem and competing species systems, _Ann. I.H. Poincaré A.N._ 19 (2002), 871-888.
  * [9] Conti, M., Terracini, S. & Verzini, G., Asymptotic estimates for the spatial segregation of competitive systems, _Adv. Math._ 195 (2005), 524-560.
  * [10] Conti, M., Terracini, S. & Verzini, G., A Variational Problem for the Spatial Segregation of Reaction-diffusion Systems, _Indiana Univ. Math. J._ 54 (2005), 779-815.
  * [11] Conti, M., Terracini, S. & Verzini, G., Uniqueness and least energy property for solutions to strongly competing systems, _Interf. Free Bound._ 8 (2006), 437-446.
  * [12] Dancer, E.N., Wang, K. & Zhang, Z., Dynamics of Strongly Competing Systems with many Species, _Trans. A.M.S._ 364 (2012), 961-1005.
  * [13] Edelstein-Keshet, L., _Mathematical Models in Biology_ , SIAM 2005.
  * [14] Gause, G.F., Experimental studies on the struggle for existence: 1. Mixed population of two species of yeast, _J. Exper. Biol._ 9 (1932), 389-402.
  * [15] Gilbarg, D. & Trudinger, N.S., _Elliptic Partial Differential Equations of Second Order_ , Springer-Verlag 2001.
  * [16] Jungnickel, D., _Graphs, Networks and Algorithms_ , Springer-Verlag 2008.
  * [17] Hamelin, F.M. & Lewis, M.A., _A differential game theoretical analysis of mechanistic models for territoriality_ , J. Math. Biol. 61 (2010), 665-694.
  * [18] Hardin, G., The Competive Exclusion Principle, _Science, New Series_ 131 (1960), 1292-1297.
  * [19] Lanzara, F. & Montefusco, E., On the limit configuration of four species strongly competing systems, _Nonlinear Differ. Equ. Appl._ 26 (2019), 19.
  * [20] Morrell, L.J. & Kokko, H., Bridging the gap between mechanistic and adaptive explanations of territory formation, _Behav. Ecol. Sociobiol._ 57 (2005), 381-390.
  * [21] Ponce, A.C., _Elliptic PDEs, Measures and Capacities_ , EMS Tracts in Mathematics 23, EMS 2016.
  * [22] Tavares, H. & Terracini, S., Regularity of the nodal set of segregated critical configurations under weak reflection law, _Calc. Var._ 45 (2012), 273-317.
  * [23] Terracini, S., Verzini, G. & Zilio, A., Spiraling asymptotic profiles of competition-diffusion systems, _Comm. Pure Appl. Math._ 72 (2019), 2578-2620.
  * [24] Wang, K, & Zhang, Z., Some new results in competing systems with many species, _Ann. I.H. Poincaré A.N._ 27 (2010), 739-761.


11institutetext:  Moscow Institute of Physics and Technology, Moscow, Russia
22institutetext: Saint Petersburg State University, Saint Petersburg, Russia
33institutetext: Caucasus Mathematical Center, Adyghe State University,
Maikop, Russia 44institutetext: IITP RAS, Moscow, Russia

# Decentralized Strongly-Convex Optimization with Affine Constraints: Primal
and Dual Approaches††thanks: The work of D. Yarmoshik in Sections 1, 6, 7 was
supported by the program “Leading Scientific Schools” (grant no.
NSh-775.2022.1.1). The work of A. Rogozin and A. Gasnikov in Sections 2–5 was
supported by Russian Science Foundation (project No. 21-71- 30005).

Alexander Rogozin 11 Demyan Yarmoshik 11 Ksenia Kopylova 22 Alexander Gasnikov
11 3 3 4 4

###### Abstract

Decentralized optimization is a common paradigm used in distributed signal
processing and sensing as well as privacy-preserving and large-scale machine
learning. It is assumed that several computational entities locally hold
objective functions and are connected by a network. The agents aim to commonly
minimize the sum of the local objectives subject by making gradient updates
and exchanging information with their immediate neighbors. Theory of
decentralized optimization is pretty well-developed in the literature. In
particular, it includes lower bounds and optimal algorithms. In this paper, we
assume that along with an objective, each node also holds affine constraints.
We discuss several primal and dual approaches to decentralized optimization
problem with affine constraints.

###### Keywords:

distributed optimization, convex optimization, constrained optimization

##  1 Introduction

Many distributed systems such as distributed sensor networks, systems for
power flow control and large-scale architectures for machine learning use
decentralized optimization as a basic mathematical tool. Several applications
such as power systems control [11, 17] lead to problems where the agents
locally hold optimization objectives and aim to cooperatively minimize the sum
of the objectives. Moreover, every node locally holds affine constraints for
its decision variable.

Decentralized optimization without affine constraints can be called a well-
examined area of research. It is known that the performance of optimization
algorithms executed over strongly-convex smooth objectives is lower bounded by
a multiple of the graph condition number and objective condition number (up to
a logarithmic factor) [19]. Both primal [8] and dual [19] algorithms that
reach the lower bounds have been proposed. The algorithms are based on
reformulating network communication constraints as affine constraints via a
communication matrix associated with the network (i.e. Laplacian matrix).
Introduction of affine constraints at the nodes leads to new classes of
algorithms that can be divided into two main types. The first type are
consensus-based methods that can be either primal or dual [2, 10, 23, 9, 12,
13, 14]. The second type are ADMM-based methods [1, 18, 3, 21]. Let us briefly
review some of the closely related papers.

The paper [12] is dedicated to constrained distributed optimization and
consider only separable objective functions (each agent has its own
independent variable). Moreover, affine constraints are supposed to be
network-compatiable (constraint matrix can have a non-zero element on position
$(i,j)$ only if there is an edge in communication graph between agents $i$ and
$j$). We do not impose such limitations: in our case each term in the
objective functions depends on the same shared variable (formulation in [12]
is obviously a special case of this) and matrix of constraints can have
arbitrary structure.

In [14] the authors present various formulations of distributed optimization
problems with different types of interconnections between constraints and
objectives, including the case, when the objective (cost) cannot be
represented as sum of cost functions of each agent. However, their algorithms
for problems with coupled affine constraints require to solve a “master
problem” on central node at each iteration and thus are not decentralized.

The authors of [20] consider multi-cluster distributed problem formulation
which is a generalization of multi-agent approach. In multi-cluster case
agents within one cluster have the same decision variable while different
clusters corresponds to different decision variables. All variables are
subject to a coupled affine constraint. By incorporating consensus constraints
into dual problem with Lagrangian multipliers the author comes to solving a
saddle point problem and prove asymptotic $O(1/N)$ ergodic convergence rate
for their method. Dependency of convergence rate on problem parameters in
saddle point approach was studied in [22].

Our paper studies the application of different techniques to decentralized
problems with affine constraints. We obtain linear convergence rates with
(explicitly specified) accelerated dependencies on function properties,
constraint matrix spectrum and communication graph properties.

The paper outline is as follows. In Section 4 we discuss a primal approach,
that is based on reformulation the initial distributed problem as a saddle-
point problem and applying algorithm of paper [7] afterwards. In Section 5, we
describe a method that allows to incorporate both affine and communication
constraints to the dual function. We refer the approach in Section 5 as a
globally dual approach. Finally, in Section 6 we describe a slightly different
dual approach that firstly takes dual functions locally at the nodes and
incorporates consensus constraints afterwards. We refer to the latter method
as a locally dual approach.

##  2 Preliminaries

Let $\operatorname{col}(x_{1},\ldots,x_{m})$ define a column vector of
$x_{1},\ldots,x_{m}\in\mathbb{R}^{d}$, i.e.
$\operatorname{col}(x_{1},\ldots,x_{m})=[x_{1}^{\top}\ldots
x_{m}^{\top}]^{\top}$. For matrices $P$ and $Q$, their Kronecker product is
defined as $P\otimes Q$. Identity matrix of size $p\times p$ is denoted ${\bf
I}_{p}$. Moreover, given a symmetric positive semi-definite matrix, we denote
$\lambda_{\max}(\cdot),\leavevmode\nobreak\
\lambda_{\min}(\cdot),\leavevmode\nobreak\ \lambda_{\min}^{+}(\cdot)$ its
maximal, minimal and minimal nonzero eigenvalues, respectively. We also let
$\sigma_{\max}(\cdot),\leavevmode\nobreak\ \sigma_{\min}(\cdot)$ and
$\sigma_{\min}^{+}(\cdot)$ be the maximal, minimal and minimal nonzero
singular values of a matrix, respectively.

In the forthcoming analysis, we will need the following basic lemma concerning
Kronecker product properties.

###### Lemma 1

Given two matrices $P$ and $Q$ such that
$\sigma_{\min}(P)=\sigma_{\min}(Q)=0$, we have

$\displaystyle\sigma_{\max}(P\otimes{\bf I}+{\bf I}\otimes Q)$
$\displaystyle=\sigma_{\max}(P)+\sigma_{\max}(Q),$
$\displaystyle\sigma_{\min}^{+}(P\otimes{\bf I}+{\bf I}\otimes Q)$
$\displaystyle=\min\left\\{\sigma_{\min}^{+}(P),\sigma_{\min}^{+}(Q)\right\\}$

###### Proof

Consider decompositions $P=U_{P}\Sigma_{P}V_{P}^{\top}$ and
$Q=U_{Q}\Sigma_{Q}V_{Q}^{\top}$, where $U_{P},V_{P},U_{Q},V_{Q}$ are
orthogonal matrices and $\Sigma_{P}$ and $\Sigma_{Q}$ are diagonal matrices
with corresponding eigenvalues at the diagonal. We have

$\displaystyle(U_{P}^{\top}\otimes U_{Q}^{\top})(P\otimes{\bf I}+{\bf
I}\otimes Q)(V_{P}\otimes V_{Q})=\Sigma_{P}\otimes{\bf I}+{\bf
I}\otimes\Sigma_{Q}.$

Denote singular values of $P$ as $\alpha_{1},\ldots,\alpha_{n}$ and the
singular values of $Q$ as $\beta_{1},\ldots,\beta_{m}$. Singular values of
$P\otimes{\bf I}+{\bf I}\otimes Q$ have form

$\displaystyle\lambda(\alpha_{i},\beta_{j})=\alpha_{i}+\beta_{j},\leavevmode\nobreak\
i=1,\ldots,n,\leavevmode\nobreak\ j=1,\ldots,m.$

Therefore, $\sigma_{\max}(P\otimes{\bf I}+{\bf I}\otimes
Q)=\sigma_{\max}(P)+\sigma_{\max}(Q)$. For the minimal nonzero singular values
we obtain

$\displaystyle\sigma_{\min}^{+}(P\otimes{\bf I}+{\bf I}\otimes
Q)=\min\left\\{\sigma_{\min}^{+}(P),\sigma_{\min}^{+}(Q)\right\\}.$

##  3 Problem Statement

Consider minimization problem with affine constraints.

$\displaystyle\min_{x\in\mathbb{R}^{d}}\leavevmode\nobreak\
\sum_{i=1}^{m}f_{i}(x)\leavevmode\nobreak\ \leavevmode\nobreak\ \text{s.t.
}Bx=0.$ (1)

We assume that each $f_{i}$ is held by a separate agent, and the agents can
exchange information through some communication network. Each agent also
locally holds affine optimization constraints $Bx=0$, where
$B\in\mathbb{R}^{p\times d}$. Further we assume that
$\operatorname{Ker}B\neq\\{0\\}$, because otherwise the constraints $Bx=0$
define a set consisting of only $\\{0\\}$, which is not an interesting case.

We make assumptions on the optimization objectives that are standard for
optimization literature [16].

###### Assumption 3.1

Each $f_{i}$ ($i=1,\ldots,m$) is differentiable, $\mu$-strongly convex and
$L$-smooth, i.e.

$\displaystyle f(y)$ $\displaystyle\geq f(x)+\left\langle\nabla
f(x),y-x\right\rangle+\frac{\mu}{2}\left\|y-x\right\|_{2}^{2},$ $\displaystyle
f(y)$ $\displaystyle\leq f(x)+\left\langle\nabla
f(x),y-x\right\rangle+\frac{L}{2}\left\|y-x\right\|_{2}^{2}.$

The communication network is represented by an undirected connected graph
${\mathcal{G}}=({\mathcal{V}},{\mathcal{E}})$. The communication constraints
are represented by a specific matrix $W$ associated with the graph
${\mathcal{G}}$.

###### Assumption 3.2

  1. 1.

${W}$ is a symmetric positive semi-definite matrix.

  2. 2.

(Network compatibility) For all $i,j=1,\dots,m$ it holds $[{W}]_{ij}=0$ if
$(i,j)\notin{\mathcal{E}}$ and $i\neq j$.

  3. 3.

(Kernel property) For any $v=[v_{1},\ldots,v_{m}]^{\top}\in\mathbb{R}^{m}$,
${W}v=0$ if and only if $v_{1}=\ldots=v_{m}$, i.e.
$\operatorname{Ker}{W}=\operatorname{span}\left\\{{\mathbf{1}}\right\\}$.

An explicit example of a matrix that satisfies Assumption 3.2 is the Graph
Laplacian $W\in\mathbb{R}^{m\times m}$:

$\displaystyle[W]_{ij}\triangleq\begin{cases}-1,&\text{if }(i,j)\in E,\\\
\text{deg}(i),&\text{if }i=j,\\\ 0,&\text{otherwise.}\end{cases}$ (2)

Let us introduce ${\bf x}=\operatorname{col}\left(x_{1}\ldots x_{m}\right)$
and ${\bf W}=W\otimes{\bf I}$. According to Assumption 3.2, communication
constraints $x_{1}=\ldots=x_{m}$ can be equivalently rewritten as ${\bf W}{\bf
x}=0$. Also introduce ${\bf B}={\bf I}\otimes B$ and $F({\bf
x})=\sum_{i=1}^{m}f_{i}(x_{i})$. That allows to rewrite problem (1) in the
following way.

$\displaystyle\min_{{\bf x}\in\mathbb{R}^{md}}\leavevmode\nobreak\ $
$\displaystyle F({\bf x})$ (3) s.t. $\displaystyle{\bf W}{\bf
x}=0,\leavevmode\nobreak\ {\bf B}{\bf x}=0.$

Reformulation 3 admits implementation of optimization methods for affinely
constrained minimization. The iterations of such methods become automatically
decentralized in the following sense. Let the optimization algorithm use
primal or dual oracle calls of the objective function and use multiplications
by the matrices representing affine constraints. In the case of problem (3)
the gradient $\nabla F({\bf x})=\operatorname{col}\left[\nabla
f_{1}(x_{1})\ldots\nabla f_{m}(x_{m})\right]$ is computed locally on the nodes
and stored in a distributed manner across the network. Multiplication by ${\bf
B}$ is also performed locally due to its definition (i.e. the $i$-th node
computes $Bx_{i}$), and the multiplication by ${\bf W}$ is performed in a
decentralized manner due to the network compatibility property of $W$ (see
Assumption 3.2).

##  4 Primal Approach

In this section, we discuss the solution of problem (3) by an algorithm APDG
[7] that only uses primal oracle calls. The algorithm is designed for saddle-
point problems, so we reformulate (3) as a saddle-point problem.

We add dual multipliers for the constraints and get a saddle-point problem

$\displaystyle\min_{{\bf x}\in\mathbb{R}^{md}}\max_{{\bf
u}\in\mathbb{R}^{mp},{\bf v}\in\mathbb{R}^{md}}\leavevmode\nobreak\ F({\bf
x})+\left\langle{\bf u},{\bf B}{\bf x}\right\rangle+\gamma\left\langle{\bf
v},{\bf W}{\bf x}\right\rangle=F({\bf x})+\left\langle\begin{pmatrix}{\bf
u}\\\ {\bf v}\end{pmatrix},\begin{pmatrix}{\bf B}\\\ \gamma{\bf
W}\end{pmatrix}{\bf x}\right\rangle.$ (4)

Algorithm 1 APDG: Accelerated Primal-Dual Gradient Method

1: Input: ${\bf x}^{0}\in\operatorname{Range}{\bf A}^{\top},{\bf
y}^{0}\in\operatorname{Range}{\bf A},$
$\eta_{x},\eta_{y},\alpha_{x},\beta_{x},\beta_{y}>0$,
$\tau_{x},\tau_{y},\sigma_{x},\sigma_{y}\in(0,1]$, $\theta\in(0,1)$

2: ${\bf x}_{f}^{0}={\bf x}^{0}$

3: ${\bf y}_{f}^{0}={\bf y}^{-1}={\bf y}^{0}$

4: for $k=0,1,2,\ldots$ do

5: ${\bf y}_{m}^{k}={\bf y}^{k}+\theta({\bf y}^{k}-{\bf y}^{k-1})$

6: ${\bf x}_{g}^{k}=\tau_{x}{\bf x}^{k}+(1-\tau_{x}){\bf x}_{f}^{k}$

7: ${\bf y}_{g}^{k}=\tau_{y}{\bf y}^{k}+(1-\tau_{y}){\bf y}_{f}^{k}$

8: ${\bf x}^{k+1}={\bf x}^{k}+\eta_{x}\alpha_{x}({\bf x}_{g}^{k}-{\bf
x}^{k})-\eta_{x}\beta_{x}{\bf A}^{\top}{\bf A}{\bf x}^{k}-\eta_{x}\left(\nabla
F({\bf x}_{g}^{k})+{\bf A}^{\top}{\bf y}_{m}^{k}\right)$

9: ${\bf y}^{k+1}={\bf y}^{k}-\eta_{y}\beta_{y}{\bf A}({\bf A}^{\top}{\bf
y}^{k}+\nabla F({\bf x}_{g}^{k}))+\eta_{y}{\bf A}{\bf x}^{k+1}$

10: ${\bf x}_{f}^{k+1}={\bf x}_{g}^{k}+\sigma_{x}({\bf x}^{k+1}-{\bf x}^{k})$

11: ${\bf y}_{f}^{k+1}={\bf y}_{g}^{k}+\sigma_{y}({\bf y}^{k+1}-{\bf y}^{k})$

12: end for

Denote ${\bf A}=\begin{pmatrix}{\bf B}\\\ \gamma{\bf W}\end{pmatrix}$. In
order to get complexity bounds for APDG applied to problem (4), we need to
bound the spectrum of ${\bf A}$. Note that ${\bf A}^{\top}{\bf A}={\bf
B}^{\top}{\bf B}+\gamma^{2}{\bf W}^{2}={\bf
I}_{m}\otimes(B^{\top}B)+\gamma^{2}W^{2}\otimes{\bf I}_{d}$. By Lemma 1 we
have

$\displaystyle\lambda_{\max}({\bf A}^{\top}{\bf A})$
$\displaystyle=\lambda_{\max}(B^{\top}B)+\gamma^{2}\lambda_{\max}^{2}(W),$
$\displaystyle\lambda_{\min}^{+}({\bf A}^{\top}{\bf A})$
$\displaystyle=\min\left\\{\lambda_{\min}^{+}(B^{\top}B),\gamma^{2}(\lambda_{\min}^{+}(W))^{2}\right\\}.$

We can also compute the condition number of ${\bf A}^{\top}{\bf A}$:

$\displaystyle\chi({\bf A}^{\top}{\bf A})=\frac{\lambda_{\max}({\bf
A}^{\top}{\bf A})}{\lambda_{\min}^{+}({\bf A}^{\top}{\bf
A})}=\frac{\lambda_{\max}(B^{\top}B)+\gamma^{2}\lambda_{\max}^{2}(W)}{\min\left\\{\lambda_{\min}^{+}(B^{\top}B),\gamma^{2}(\lambda_{\min}^{+}(W))^{2}\right\\}}.$

By accurately choosing factor $\gamma$, we can control the condition number
$\chi({\bf A}^{\top}{\bf A})$. The minimal value of $\chi({\bf A}^{\top}{\bf
A})$ is attained at
$\gamma^{2}=\frac{\lambda_{\min}^{+}(B^{\top}B)}{(\lambda_{\min}^{+}(W))^{2}}$
and equals $\chi({\bf A}^{\top}{\bf A})=\chi(B^{\top}B)+\chi^{2}(W)$.
Therefore, if we apply APDG directly to problem (3), the complexity would be

$\displaystyle
O\left(\max\left(\sqrt{\chi^{2}(W)+\chi(B^{\top}B)}\sqrt{\frac{L}{\mu}},\chi^{2}(W)+\chi(B^{\top}B)\right)\log\frac{1}{\varepsilon_{\bf
x}}\right)$

calls of $\nabla f_{i}(\cdot)$ at each node and communication rounds, with
$\varepsilon_{\bf x}$ being the desired distance to the solution: $\|{\bf
x}^{N}-{\bf x}^{*}\|\leq\varepsilon_{\bf x}$. In the smooth, strongly convex
case it is also the complexity for satisfying $F({\bf x}^{N})-F({\bf
x})\leq\varepsilon_{F}$ or $\|{\bf A}{\bf x}^{N}\|\leq\varepsilon_{\bf A}$ (up
to logarithmic dependencies on the problem parameters). Indeed, from Lipschitz
smoothness we have $F({\bf x}^{N})-F({\bf x})\leq L\varepsilon^{2}_{\bf x}/2$
and $\|{\bf A}{\bf x}^{N}\|=\|{\bf A}{\bf x}^{N}-{\bf A}{\bf
x}^{*}\|\leq\sigma_{max}({\bf A})\varepsilon_{\bf x}$. By that means, in the
following inequalities $\varepsilon$ can be replaced by any of
$\varepsilon_{\bf x},\leavevmode\nobreak\ \varepsilon_{f},\leavevmode\nobreak\
\varepsilon_{\bf A}$.

The dependence on network parameters $W$ and affine constraints parameters $B$
can be enhanced by using Chebyshev acceleration [19]. Let us replace $W$ by a
Chebyshev polynomial $P_{K}(W)$ such that it has degree
$K=O\left(\sqrt{\chi(W)}\right)$ and condition number $\chi(P_{K}(W))=O(1)$.
Multiplication by $P_{K}(W)$ is equivalent to making $K$ communication rounds.
Analogically, let us replace $B^{\top}B$ by a Chebyshev polynomial
$P_{M}(B^{\top}B)$ with degree $M=O\left(\sqrt{\chi(B^{\top}B)}\right)$ and
condition number $\chi\left(P_{M}(B^{\top}B)\right)=O(1)$. As a result, we
obtain

$\displaystyle
N=O\left(\sqrt{\frac{L}{\mu}}\log\frac{1}{\varepsilon}\right)\text{ oracle
calls at each node},$ $\displaystyle O\left(N\sqrt{\chi(W)}\right)\text{
communications},$ $\displaystyle O\left(N\sqrt{\chi(B^{\top}B)}\right)\text{
multiplications by $B$, $B^{\top}$ at each node}.$

##  5 Globally Dual Approach

In this section, we describe an approach to solving (3) that is based on
passing to the dual problem. We call this approach “global” since both
constraints, that is, affine constraints ${\bf B}{\bf x}=0$ and communication
constraints ${\bf W}{\bf x}=0$ are used in the dual reformulation.

Let $\gamma$ be a positive scalar and ${\bf A}^{\top}=[{\bf
B}^{\top}\leavevmode\nobreak\ \gamma{\bf W}]$ and introduce dual function

$\displaystyle\Phi({\bf y})=\max_{{\bf x}\in\mathbb{R}^{md}}\left[-F({\bf
x})+\left\langle{\bf y},{\bf A}{\bf x}\right\rangle\right]=F^{*}({\bf
A}^{\top}{\bf y}).$

We have $\displaystyle\nabla\Phi({\bf y})={\bf A}\nabla F^{*}({\bf
A}^{\top}{\bf y})={\bf A}\cdot\operatorname*{arg\,min}_{{\bf
x}\in\mathbb{R}^{md}}\left[-F({\bf x})+\left\langle{\bf y},{\bf A}{\bf
x}\right\rangle\right]$. Note that multiplication by ${\bf A}$ is performed in
a distributed manner: indeed, it includes local multiplications by $B$ and a
consensus round, which is a multiplication by ${\bf W}$. Moreover, the
$\operatorname*{arg\,min}$ operation is computed locally, which is standard
for decentralized optimization [19]. Finally, dual function $\Phi$ is
$\frac{\lambda_{\max}({\bf A}^{\top}{\bf A})}{\mu}$-smooth on
$\mathbb{R}^{m(p+d)}$ and $L_{\Phi}=\frac{\lambda_{\min}^{+}({\bf
A}^{\top}{\bf A})}{L}$-strongly convex on $(\operatorname{Ker}{\bf
A}^{\top})^{\bot}$. Solving dual problem

$\displaystyle\min_{{\bf y}\in\mathbb{R}^{m(p+d)}}\leavevmode\nobreak\
\Phi({\bf y})$

by a fast gradient method (see i.e. accelerated Nesterov method in Section 2.2
of [16]) until accuracy $\Phi({\bf y}^{N})-\Phi({\bf
y})\leq\varepsilon_{\Phi}$ requires
$N=O\left(\sqrt{\frac{L}{\mu}}{\sqrt{\chi({\bf A}^{\top}{\bf
A})}}\log\frac{1}{\varepsilon_{\Phi}}\right)$ iterations.

Following the same arguments as in Section 4, we compute the condition number
$\chi({\bf A}^{\top}{\bf A})$:

$\displaystyle\chi({\bf A}^{\top}{\bf
A})=\frac{\lambda_{\max}(B^{\top}B)+\gamma^{2}\lambda_{\max}^{2}(W)}{\min\left\\{\lambda_{\min}^{+}(B^{\top}B),\gamma^{2}(\lambda_{\min}^{+}(W)^{2})\right\\}}.$

The minimal value of $\chi({\bf A}^{\top}{\bf A})$ is attained at
$\gamma^{2}=\frac{\lambda_{\min}^{+}(B^{\top}B)}{(\lambda_{\min}^{+}(W))^{2}}$
and equals $\chi({\bf A}^{\top}{\bf A})=\chi(B^{\top}B)+\chi^{2}(W)$.
Communication and computation complexities of fast dual method equal

$\displaystyle
O\left(\sqrt{\frac{L}{\mu}}\left(\chi(B^{\top}B)+\chi^{2}(W)\right)^{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\frac{1}{2}}\log\frac{1}{\varepsilon_{\Phi}}\right).$

To obtain desired complexity estimates for the algorithm to find the
approximate solution ${\bf x}^{N}$ satisfying $F({\bf x}^{N})-F({\bf
x})\leq\varepsilon$ and $\|{\bf A}{\bf x}^{N}\|\leq\varepsilon$, we refer to
the following properties of dual function (see, e.g. Theorem 5.2 from [4]):

$\displaystyle\|\nabla\Phi({\bf y})\|\leq\epsilon/R_{\bf y}$
$\displaystyle\Rightarrow F({\bf x}({\bf y}))-F({\bf x}^{*})\leq\epsilon,$
$\displaystyle\|\nabla\Phi({\bf y})\|\leq\epsilon$
$\displaystyle\Rightarrow\|{\bf A}{\bf x}({\bf y})\|\leq\epsilon,$

where $\|{\bf y}\|\leq 2R_{\bf y}$, and ${\bf x}({\bf
y})=\operatorname*{arg\,min}_{{\bf x}\in\mathbb{R}^{md}}\left[-F({\bf
x})+\left\langle{\bf y},{\bf A}{\bf x}\right\rangle\right]$. Combining it with
$\Phi({\bf y}^{N})-\Phi({\bf y})\geq\|\Phi({\bf y}^{N})\|^{2}/2L_{\Phi}$,
which is true for a smooth convex function, we justify substitution of
$\varepsilon_{\Phi}$ by $\varepsilon$ in the complexity estimate. This
transition will only change the constant hidden by big-O notation (by the
factor of two), and affect omitted logarithmic dependencies on the problem
parameters.

To employ Chebyshev acceleration in this case we do substitution ${\bf
A}^{\top}{\bf y}\to{\bf p}$. In this variables accelerated Nesterov method
turns into Algorithm 2, where ${\bf x}({\bf q})=\nabla F^{*}({\bf
q})=\operatorname*{arg\,min}\left[-F({\bf x})+\left\langle{\bf
q},x\right\rangle\right]$:

Algorithm 2 Globally Dual Method

1: Input: ${\bf p}^{0}\in\operatorname{Range}{\bf A}^{\top}$, $\eta>0$,
$\beta\in(0,1)$

2: ${\bf p}^{-1}={\bf p}^{0}$

3: for $k=0,1,2,\ldots$ do

4: ${\bf q}={\bf p}^{k}+\beta\left({\bf p}^{k}-{\bf p}^{k-1}\right)$

5: ${\bf p}^{k+1}={\bf q}-\eta{\bf A}^{\top}{\bf A}{\bf x}({\bf q})$

6: end for

For the algorithm in this form we can replace ${\bf A}^{\top}{\bf A}$ with
Chebyshev polynomial of it, as we did in Section 4, and obtain the same
complexity estimates as for APDG:

$\displaystyle
N=O\left(\sqrt{\frac{L}{\mu}}\log\frac{1}{\varepsilon}\right)\text{ oracle
calls at each node},$ $\displaystyle O\left(N\sqrt{\chi(W)}\right)\text{
communications},$ $\displaystyle O\left(N\sqrt{\chi(B^{\top}B)}\right)\text{
multiplications by $B$, $B^{\top}$ at each node}.$

##  6 Locally Dual approach

In Section 5 we discussed a dual reformulation of (3) where both constraints
${\bf B}{\bf x}=0$ and ${\bf W}{\bf x}=0$ are used simultaneously. This
section describes a dual approach, as well, but the difference is that we
firstly pass to dual functions locally at the nodes and impose the
communication constraints only afterwards.

###  6.1 Utilizing locality on ${\bf u}$

One can note that in the above approaches optimization over ${\bf u}$ could be
done locally at each node. This is equivalent to including affine constraints
into the objective (as an indicator function) instead of handling them with
Lagrangian multipliers. In settings there the “cost” of communication is
limiting or comparable to that of local computations, we can find the solution
faster by going this way. It may be the case when $x$ has a small dimension
and decentralization is desirable due to privacy constraints.

Dual problem in this approach will be

$\max_{\bf v}\min_{{\bf B}{\bf x}=0}\\{F({\bf x})+\langle{\bf v},{\bf W}{\bf
x}\rangle\\}=-\min_{\bf v}F^{*}_{[{\bf B}{\bf x}=0]}({\bf W}^{\top}{\bf v}),$

where $\displaystyle F^{*}_{[{\bf B}{\bf x}=0]}({\bf v})=\max_{{\bf B}{\bf
x}=0}\\{\langle{\bf v},{\bf x}\rangle-F({\bf x})\\}$ denotes a convex
conjugate under affine constraints.

We can reduce the problem of computing the gradient of such a modified
conjugate function to calling conventional dual oracle. Let $E$ be a matrix,
the rows of which constitute an orthogonal basis in the null space of $B$
(matrix $E$ can be computed at the preprocessing stage of an optimization
algorithm). Then instead of working with functions $f_{i}(x)$ we can optimize
the sum of functions $h_{i}(t)=f_{i}(Et)$.

Denote ${\bf t}=\operatorname{col}(t_{1},\ldots,t_{m})$, $H({\bf
t})=\sum_{i=1}^{m}h_{i}(t_{i})$. Then problem (1) could be written in
decentralized way as follows

$\displaystyle\min_{\bf t}$ $\displaystyle\sum_{i=1}^{m}h_{i}(t_{i})$ s.t.
$\displaystyle{\bf W}_{\bf t}{\bf t}=0.$

Its dual form is

$\displaystyle\max_{{\bf t}}\\{\langle{\bf z},{\bf W}_{\bf t}{\bf
t}\rangle-H({\bf t})\\}=-\min_{{\bf z}}H^{*}({\bf W}_{\bf t}^{\top}{\bf z}),$

and the gradient of the objective can be computed using Demyanov–Danskin’s
theorem:

$\displaystyle\nabla H^{*}({\bf z})=\arg\max_{\bf t}\\{\langle{\bf z},{\bf
t}\rangle-F({\bf E}{\bf t})\\}.$

From smaller dimension of $t$ comparing to $x$ we can expect that computation
of $\nabla H^{*}({\bf z})$ is easier than calling conventional first-order
dual oracle, the only drawback is the necessity of storing matrix $E$ and
performing multiplications by $E$.

Let $\mu_{t}$ and $L_{t}$ be the constants of strong convexity and Lipschitz
smoothness of $h_{i}$ respectively for all $i=1,\ldots,m$. Then, obviously,
$\mu_{t}\geq\mu$ and $L_{t}\leq L$. For example, if $f_{i}(x)$ is twice
continuously differentiable, then its smoothness constant can be computed as
$\displaystyle
L_{x,i}=\sup_{x\in\mathbb{R}^{n}}\lambda_{\max}(\nabla^{2}f_{i}(x))$. The
smoothness constant of $h_{i}(t)$ is given by $\displaystyle
L_{t,i}=\sup_{t\in\mathbb{R}^{d_{t}}}\lambda_{\max}(E^{\top}\nabla^{2}f_{i}(Et)E)$.
Note that the dimension of $t$ can be computed as $d_{t}=d-\text{rank}(B)$. In
the latter variant the maximum is taken over a smaller set of points, and
multiplication by $E$ is likely to further reduce the smoothness constant (and
increase strong convexity constant).

Since $H({\bf t})$ is $L_{t}$-smooth and $\mu_{t}$-strongly convex, we have
that $F^{*}_{[{\bf B}{\bf x}=0]}({\bf z})=H^{*}({\bf z})$ is
$\frac{1}{\mu_{t}}$-smooth and $\frac{1}{L_{t}}$-strongly convex [6].

Thus, the fast gradient method [15] applied to the dual problem requires

$\displaystyle
O\left(\sqrt{\frac{L_{t}}{\mu_{t}}}\chi(W)\log\frac{1}{\varepsilon}\right),$

dual-oracle calls and communication rounds to ensure $F({\bf x}^{N})-F({\bf
x})\leq\varepsilon$ and $\|{\bf A}{\bf x}^{N}\|\leq\varepsilon$ (see Section 5
for details). And using Chebyshev acceleration as described in Section 4 we
can reduce the complexities to

$\displaystyle
N=O\left(\sqrt{\frac{L_{t}}{\mu_{t}}}\log\frac{1}{\varepsilon}\right)\text{
oracle calls at each node},$ $\displaystyle
O\left(N\sqrt{\chi(W)}\right)\text{ communications}.$

##  7 Numerical Experiments

In the simulation we consider the following smooth, strongly convex objective
function:

$\displaystyle f_{i}(x)$
$\displaystyle=\frac{1}{2}\|C_{i}x-d_{i}\|_{2}^{2}+\frac{\theta}{2}\|x\|_{2}^{2},$
$\displaystyle F({\bf x})$ $\displaystyle=\frac{1}{2}\|{\bf C}{\bf x}-{\bf
d}\|_{2}^{2}+\frac{\theta}{2}\|{\bf x}\|_{2}^{2},$ $\displaystyle{\bf C}$
$\displaystyle=\text{diag}(C_{1},\ldots,C_{m}),\leavevmode\nobreak\ {\bf
d}=\operatorname{col}{(d_{1},\ldots,d_{m})}.$

We consider different parameters of the problem such as the dimension of $x$,
the rank of $B\in\mathbb{R}^{\dim(x)\times\dim(x)}$ and the number of nodes.
For each case we plot function error and constraints violation norm at each
iteration for all our algorithms: APDG, Locally and Globally Dual approaches.
The Chebyshev acceleration is not applied in the experiments, so each
iteration corresponds to one gradient computation (gradient of primal function
in case of APDG, and gradient of dual function is case of dual approaches). We
also provide tables with comparison of time and number of iterations required
to achieve given accuracy. Time is measured with our Python/NumPy [5]
implementation of the algorithms, which is available on GitHub111Source code:
https://github.com/niquepolice/decentr_constr_dual.

  1. 1.

For the first case we consider the ring network with $m=5$ nodes,
$x\in\mathbb{R}^{40}$ and $\operatorname{rank}B=1$. Typical convergence plot
is shown on Fig. 1. One can see that all algorithms converge linearly, with
the fastest one in terms of iterations number being Locally Dual, and the
slowest one being APDG. However, computing the gradient of a dual function
might be an arithmetically more expensive operation than computing primal
gradient in the black-box scenario. In our implementation we compute the
gradient of dual function by numerically solving the system of linear
equations with its right-hand part being changed between iterations. It means
that one iteration of the Dual methods is more time-consuming than one
iteration of APDG. In the Table 1, we compare computational time and number of
iterations required to achieve given accuracy. The results are averaged for
100 randomly generated problems.

Figure 1: 5 nodes, $\dim(x)=40$, $\operatorname{rank}B=1$. Table 1: Time and iterations for achieving $\|{\bf A}{\bf x}^{k}\|<10^{-2}$. Averaged over 100 experiments. Problem parameters: 5 nodes, $\dim(x)=40$, $\operatorname{rank}B=1$. | APDG | Globally Dual | Locally Dual
---|---|---|---
Iterations | 875.3 | 502.7 | 276.7
Time (s) | 0.193 | 0.510 | 0.233
  2. 2.

Next we use the same number of nodes and the dimension of $x$, but increase
the rank of $B$. Even for $\operatorname{rank}B=3$ the condition number of the
locally dual problem usually is about two orders of magnitude smaller than the
condition number of the globally dual problem, therefore the globally dual
approach has a significant advantage in that case. Typical convergence plots
are shown in Figure 2, averaged iteration and time complexities for satisfying
stopping criteria are shown in Table 2.

Figure 2: 5 nodes, $\dim(x)=40$, $\operatorname{rank}B=3$. Table 2: Time and iterations for achieving $\|{\bf A}{\bf x}^{k}\|<10^{-1}$. Averaged over 100 experiments. Problem parameters: 5 nodes, $\dim(x)=40$, $\operatorname{rank}B=3$. | APDG | Globally Dual | Locally Dual
---|---|---|---
Iterations | 1555.5 | 1551.7 | 123.1
Time (s) | 0.337 | 1.577 | 0.127
  3. 3.

In the case of higher dimension ($10$ nodes, $\dim(x)=100$,
$\operatorname{rank}B=1$) we used Erdős-Rényi random communication graphs with
edge probability = $0.3$. APDG seems to converge much faster by constraints
violation norm at first iterations then other methods (Fig. 3), and its
convergence rate is close to other methods. See also Table 3 for averaged
results of multiple experiments.

Figure 3: Erdős-Rényi graph on 10 nodes, average degree = 3.6. $\dim(x)=100$, $\operatorname{rank}B=1$. Table 3: Time and iterations for achieving accuracy $\|{\bf A}{\bf x}^{k}\|<10^{1}$. Averaged over 10 experiments. Problem parameters: 10 nodes, edge probability = $0.3$, $\dim(x)=100$, $\operatorname{rank}B=1$. | APDG | Globally Dual | Locally Dual
---|---|---|---
Iterations | 404.3 | 2227.9 | 1425.5
Time (s) | 2.561 | 54.024 | 16.544

## References

  * [1] T. Erseghe.  Distributed optimal power flow using admm.  IEEE Transactions on Power Systems, 29(5):2370–2380, Sep. 2014\.
  * [2] A. Falsone, K. Margellos, S. Garatti, and M. Prandini.  Dual decomposition for multi-agent distributed optimization with coupling constraints.  Automatica, 84:149–158, 2017.
  * [3] A. Falsone, I. Notarnicola, G. Notarstefano, and M. Prandini.  Tracking-admm for distributed constraint-coupled optimization.  Automatica, 117:1–13, 202.
  * [4] E. Gorbunov, D. Dvinskikh, and A. Gasnikov.  Optimal decentralized distributed algorithms for stochastic convex optimization.  arXiv:1911.07363, 2019.
  * [5] C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane, J. F. del Río, M. Wiebe, P. Peterson, P. Gérard-Marchant, K. Sheppard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke, and T. E. Oliphant.  Array programming with NumPy.  Nature, 585(7825):357–362, Sept. 2020.
  * [6] S. M. Kakade, S. Shalev-Shwartz, and A. Tewari.  On the duality of strong convexity and strong smoothness: learning applications and matrix regularization.  Technical report, Toyota Technological Institute, 2009.
  * [7] D. Kovalev, A. Gasnikov, and P. Richtárik.  Accelerated primal-dual gradient method for smooth and convex-concave saddle-point problems with bilinear coupling.  arXiv preprint arXiv:2112.15199, 2021.
  * [8] D. Kovalev, A. Salim, and P. Richtárik.  Optimal and practical algorithms for smooth and strongly convex decentralized optimization.  Advances in Neural Information Processing Systems, 33:18342–18352, 2020.
  * [9] S. Liang, L. Y. Wang, and G. Yin.  Distributed smooth convex optimization with coupled constraints.  IEEE Transactions on Automatic Control, 65:347–353, Jan 2020.
  * [10] S. Liang, X. Zheng, and Y.Hong.  Distributed ninsmooth optimization with coupled inequality constraints via modified lagrangian function.  IEEE Transaction on Automatic Control, 63:1753–1759, 2018.
  * [11] D. K. Molzahn, F. Dörfler, H. Sandberg, S. H. Low, S. Chakrabarti, R. Baldick, and J. Lavaei.  A survey of distributed optimization and control algorithms for electric power systems.  IEEE Transactions on Smart Grid, 8(6):2941–2962, Nov 2017.
  * [12] I. Necoara and V. Nedelcu.  Distributed dual gradient methods and error bound conditions.  arXiv:1401.4398, 2014.
  * [13] I. Necoara and V. Nedelcu.  On linear convergence of a distributed dual gradient algorithm for linearly constrained separable convex problems.  Automatica, 55:209–216, 2015.
  * [14] I. Necoara, V. Nedelcu, and I. Dumitrache.  Parallel and distributed optimization methods for estimation and control in networks.  Journal of Process Control, 21(5):756–766, 2011.  Special Issue on Hierarchical and Distributed Model Predictive Control.
  * [15] Y. Nesterov.  A method of solving a convex programming problem with convergence rate $o(1/k^{2})$.  Soviet Mathematics Doklady, 27(2):372–376, 1983.
  * [16] Y. Nesterov.  Introductory Lectures on Convex Optimization: a basic course.  Kluwer Academic Publishers, Massachusetts, 2004.
  * [17] N. Patari, V. Venkataramanan, A. Srivastava, D. K. Molzahn, N. Li, and A. Annaswamy.  Distributed optimization in distribution systems: Use cases, limitations, and research needs.  IEEE Transactions on Power Systems, pages 1–1, 2021.
  * [18] V. Rostampour, O. t. Haar, and T. Keviczky.  Distributed stochastic reserve scheduling in ac power systems with uncertain generation.  IEEE Transactions on Power Systems, 34(2):1005–1020, 2019.
  * [19] K. Scaman, F. Bach, S. Bubeck, Y. T. Lee, and L. Massoulié.  Optimal algorithms for smooth and strongly convex distributed optimization in networks.  In D. Precup and Y. W. Teh, editors, Proceedings of the 34th International Conference on Machine Learning, volume 70 of Proceedings of Machine Learning Research, pages 3027–3036, International Convention Centre, Sydney, Australia, 06–11 Aug 2017. PMLR.
  * [20] J. Wang and G. Hu.  Distributed optimization with coupling constraints in multi-cluster networks based on dual proximal gradient method.  arXiv preprint arXiv:2203.00956, 2022.
  * [21] Z. Wnag and C. J. Ong.  Distributed model predictive control of linear descrete-times systems with local and global cosntraints.  Automatica, 81:184–195, 2017.
  * [22] D. Yarmoshik, A. Rogozin, O. Khamisov, P. Dvurechensky, A. Gasnikov, et al.  Decentralized convex optimization under affine constraints for power systems control.  arXiv preprint arXiv:2203.16686, 2022.
  * [23] D. Yuan, D. W. C. Ho, and G. P. Jiang.  An adaptive primal-dual subgradient algorithm for online distributed constrained optimization.  IEEE Transactions on Cybernetics, 48:3045–3055, 2018.