# Can mobility induce orders in active XY spins on a substrate?

Astik Haldar<EMAIL_ADDRESS>Theory Division, Saha Institute of
Nuclear Physics, HBNI, 1/AF Bidhannagar, Calcutta 700064, West Bengal, India
Apurba Sarkar<EMAIL_ADDRESS>School of Mathematical &
Computational Sciences, Indian Association for the Cultivation of Science,
Kolkata-700032, West Bengal, India  Swarnajit Chatterjee
<EMAIL_ADDRESS>Center for Biophysics & Department for
Theoretical Physics, Saarland University, 66123 Saarbrücken, Germany  Abhik
Basu<EMAIL_ADDRESS><EMAIL_ADDRESS>Theory Division, Saha
Institute of Nuclear Physics, HBNI, 1/AF Bidhannagar, Calcutta 700064, West
Bengal, India

###### Abstract

We elucidate how the interplay between diffusive mobility and number
conservation can make nearly phase-ordered active XY spins on a substrate
stable. For wide-ranging model parameters, it has stable uniform phases with
orientational order logarithmically stronger or weaker than in equilibrium,
together with miniscule (i.e., hyperuniform) or giant number fluctuations,
respectively, forming a new universality class. For other parameters, it has
no stable uniformly ordered phase. Our theory sheds light on wide-ranging
systems, e.g., active superfluids on substrates, synchronization of
oscillators, active carpets of cilia and bacterial flagella, and sandblasting.

Nonequilibrium systems can be strikingly different from their equilibrium
analogs: they can display both ordered states or instabilities not found in
their equilibrium analogs. For instance, a two-dimensional (2D) collection of
self-propelled particles can be in orientationally long-range ordered states
in the presence of finite noises Marchetti  _et al._ (2013); Toner and Tu
(1995), which is prohibited in equilibrium systems with continuous symmetries
Mermin and Wagner (1966). Driven systems can show novel instabilities as well,
e.g., driven 2D isotropic Bose systems, e.g., driven exciton polariton fluids
on a substrate Carusotto and Ciuti (2013), cannot exhibit off-diagonal
algebraic correlations (i.e., 2D superfluidity) on a substrate, but show only
short-range order (SRO) for any weak noises Altman  _et al._ (2015), unlike
quasi-long-range order (QLRO) at low enough temperature ($T$) in its
equilibrium counterpart Chaikin and Lubensky (1995).

In this Letter, we study diffusively moving nearly phase-ordered active XY
spins on a 2D solid substrate, as illustrated in Fig. 1. The “activity” of the
model stems from the propensity of the spins to rotate in response to the
local spin concentration and the magnitude of the local phase difference.
Absence of self-propulsion of the spins makes this active model distinct from
the celebrated “moving XY” Toner-Tu model Toner and Tu (1995) or the
microscopic Vicsek model Vicsek  _et al._ (1995) for flocks. We develop the
hydrodynamic theory for this system, supplemented by a lattice-gas type
microscopic agent-based model constructed by us. Surprisingly, this model
encompasses both the striking features of nonequilibrium systems vis-a-vis
their equilibrium counterparts, viz., novel order and instabilities in
different parameter ranges, and belongs to a heretofore unstudied universality
class.

Figure 1:  Schematic diagram of the model. Each lattice site contains an
arbitrary number of XY spins, with the phase at a site being the average
phases of all the spins there, marked by bluish thick arrows. The phase
changes by local “active rotations” in response to the local spin
concentration and magnitude of the nearest-neighbor phase differences, and by
simple relaxation (thick yellow arrow). The density changes via currents
having diffusive and local phase difference-dependent parts, denoted by the
red broken horizontal arrow.

Interestingly, this theory also describes the 2D Kardar-Parisi-Zhang (KPZ)
equation Kardar  _et al._ (1986); Barabási  _et al._ (1995) for surface growth
or erosion (“sandblasting”) with a conserved species on it.

Our most striking result is that, in contrast to either non-number conserving
or immobile XY spins on a substrate that admit only SRO Altman  _et al._
(2015), the mobility and number conservation together in this system can lead
to stable uniformly ordered phases; see supplemental movie MOV1 sup . These
ordered phases are distinguished by varieties of stable orientational order
for a wide range of model parameters. Although like QLRO, the variance of the
orientational fluctuations grows with the system size $L$, it does so either
logarithmically slower than QLRO (i.e. stronger order than QLRO or “SQLRO”),
or faster than QLRO (i.e. weaker order than QLRO or “WQLRO”), depending upon
the model parameters. In a surprising correspondence with phase fluctuations,
the density fluctuations are either generically miniscule, or giant,
respectively. Intriguingly, fluctuations, which are underdamped in these
ordered states, are characterized by continuously varying scaling exponents.
For other choices of the parameters, there are no stable uniform phase-ordered
states. Without mobility or number conservation, our model predicts short-
range phase-ordering; see supplemental movie MOV2 sup . In the equilibrium
limit, the model shows QLRO phase order and normal number fluctuations.

Our model provides a broad perspective concerning the impact of number
fluctuations on the orientational order in wide classes of 2D driven systems
with rotational invariance at low noises that exhibit QLRO in equilibrium. For
instance, can an isotropic driven number-conserving mobile ordered 2D
condensate exist, and if so, with what type of order? Likewise, generic
conditions for synchronization in phase-locked states of mobile oscillators,
having the same $U(1)$ internal symmetry as XY spins, is a question of
paramount significance Pikovsky  _et al._ (2003); Dörfler and Bullo (2014);
Strogatz (2018); Uriu  _et al._ (2010); Uriu and Morelli (2014); Zhou  _et
al._ (2016); Levis  _et al._ (2017); Frasca  _et al._ (2008); Peruani  _et
al._ (2010); Levis  _et al._ (2017); Banerjee and Basu (2017). Other
experimentally realizable systems motivating our theory include “active
carpets” of cilia or bacterial flagella, modeled as active rotors grafted on
the carpet and can order Uchida and Golestanian (2010), and take part in
nutrient transport Mathijssen  _et al._ (2018) or mucous in respiratory tract
D. R. Brumley and Goldstein (2012), monolayers of spinning colloids Bruot and
Cicuta (2016); Zhang  _et al._ (2021), and in vitro magnetic cilia carpets et
al (2020). Lastly, possible suppression of the thresholdless instability of 2D
KPZ surfaces is an overarching theoretical issue, given the paradigmatic
status of the KPZ equation in nonequilibrium physics Barabási  _et al._
(1995).

Below we outline our hydrodynamic theory that quantitatively explains
indispensable role of the conserved density fluctuations to sustain order
here, as dramatically manifested in movies MOV1 vis-à-vis MOV2 sup . Details
can be found in the associated long paper (ALP) Haldar  _et al._ (2021).

Due to friction from the substrate, there is no momentum conservation, so the
only conserved variable on the surface is the spin number density $c({\bf
x},t)$. In addition, for nearly phase-ordered spins, the broken symmetry phase
fluctuations $\theta({\bf x},t)$ about an arbitrary reference state are slow
variables with relaxation rates diverging in the long wavelength limit, but
the amplitude fluctuations relax fast. Therefore, $\theta({\bf x},t)$, and
$c({\bf x},t)$ are the only hydrodynamic variables. For a driven system, we
must write down the equations of motion by appealing to the general symmetries
of the underlying microscopic dynamics (i.e., translation and rotation) and
the phase-ordered state (here, arbitrariness of the reference state) and
conservation laws. Additional equilibrium requirements like detailed balance
do not apply to our nonequilibrium system. Retaining up to the lowest order
symmetry-permitted nonlinear terms in $\theta$ and spatial gradients, the
dynamical equation of $\theta({\bf x},t)$ has the form

$\partial_{t}\theta=\kappa\nabla^{2}\theta+\frac{\lambda}{2}({\bm{\nabla}}\theta)^{2}+\Omega(c)+f_{\theta},$
(1)

where $\Omega(c)$ is a general function of $c$. In case of fluctuating
surfaces, $\theta({\bf x},t)$ is the local height with respect to an arbitrary
base plane Barabási  _et al._ (1995). Number density $c$ follows a
conservation law with a current ${\bf J}_{c}$ given by

${\bf
J}_{c}=-D{\bm{\nabla}}c-\lambda_{0}\tilde{\Omega}(c){\bm{\nabla}}\theta.$ (2)

Again we have truncated up to the lowest order in $\theta$ and spatial
gradients; $\tilde{\Omega}(c)$ is yet another function of $c$. While we have
constructed (1) and (2) by writing down all the rotation-invariant leading
order terms, each term actually carries a simple physical interpretation. The
linear $\kappa\nabla^{2}\theta$ term in (1), and $D{\bm{\nabla}}c$ term in (2)
are just the equilibrium spin relaxation and particle diffusion, respectively;
$\kappa>0$ and $D>0$, respectively, are the spin stiffness and diffusivity.
The $\lambda$\- and $\Omega(c)$-terms in (1) represent, respectively, active
rotations of the spins in response to local gradients of $\theta$, and local
spin concentration. The $\lambda_{0}$-term in (2) models particle currents in
response to spatially nonuniform phases. Gaussian noises $f_{\theta}$ and
$f_{c}$ are white and conserved noises, respectively, with zero mean, and
variances

$\displaystyle\langle f_{\theta}({\bf
x},t)f_{\theta}(0,0)\rangle=2D_{\theta}\delta^{2}({\bf x})\delta(t),$ (3)
$\displaystyle\langle f_{c}({\bf
x},t)f_{c}(0,0)\rangle=2D_{c}(-\nabla^{2})\delta^{2}({\bf x})\delta(t),$ (4)

consistent with $\theta$ and $c$ being, respectively, a non-conserved and a
conserved variable. Writing $c({\bf x},t)=c_{0}+\delta c({\bf
x},t),\,\langle\delta c\rangle=0$, we obtain

$\displaystyle\frac{\partial\theta}{\partial
t}=\kappa\nabla^{2}\theta+\Omega_{1}\delta
c+\frac{\lambda}{2}({\bm{\nabla}}\theta)^{2}+\Omega_{2}(\delta
c)^{2}+f_{\theta},$ (5) $\displaystyle\frac{\partial\delta c}{\partial
t}=\lambda_{0}\tilde{\Omega}_{0}\nabla^{2}\theta+D\nabla^{2}\delta
c+\lambda_{0}\tilde{\Omega}_{1}{\bm{\nabla}}\cdot(\delta
c{\bm{\nabla}}\theta)+f_{c},$ (6)

where we have retained the most relevant nonlinear terms in fields and
gradients. Parameters $\Omega_{1}\equiv\partial\Omega/\partial
c|_{c=c_{0}},\,\Omega_{2}\equiv\partial^{2}\Omega/\partial^{2}c|_{c=c_{0}},\,\tilde{\Omega}_{0}\equiv\tilde{\Omega}|_{c=c_{0}},\,\tilde{\Omega}_{1}\equiv\partial\tilde{\Omega}/\partial
c|_{c=c_{0}}$. Couplings $\lambda,\Omega_{2},\lambda_{0}\tilde{\Omega}_{1}$
have arbitrary signs. For active carpets with grafted rotors, $c$ represents
the concentration of nutrients or other chemical substance, moved around by
the beating of cilia or bacterial flagella Uchida and Golestanian (2010);
Mathijssen  _et al._ (2018). Purely for symmetry reasons, Eqs. (5) and (6)
describe the hydrodynamics of immobile active XY spins interacting with an
incompressible binary fluid in its well-mixed phase on a substrate, where $c$
now is the binary fluid order parameter, and also chiral active hexatics on a
substrate Maitra  _et al._ (2020). One-dimensional versions of (5) and (6) are
same as the hydrodynamic equations of a one-dimensional sedimenting crystal
Lahiri and Ramaswamy (1997).

It is useful to study Eqs. (5) and (6) in their linear limit first, by setting
$\lambda=\Omega_{2}=\tilde{\Omega}_{1}=0$. When
$\Omega_{1}\tilde{\Omega}_{0}\lambda_{0}>0$, uniform states are linearly
stable. Small fluctuations around these linearly stable states travel non-
dispersively, i.e., with a wavespeed independent of wavevector $k$. The
correlators in the linear theory are exactly calculated conveniently in
Fourier space. For instance, the phase and density autocorrelation functions
show QLRO and “normal number fluctuations (NNF)”, indistinguishable from the
2D equilibrium XY model and a non-critical equilibrium system with short-range
interactions. In particular, equal-time phase correlator

$C^{0}_{\theta\theta}(k)\equiv\langle|\theta({\bf
k},t)|^{2}\rangle_{0}\approx\frac{\overline{D}}{2\Gamma k^{2}},$ (7)

in the long wavevlength limit, where $\overline{D}\equiv D_{\theta}+D_{c}$ and
$\Gamma\equiv(\kappa+D)/2$; “0” refers to a linear theory result. This in turn
gives for the variance

$\Delta_{\theta}^{0}\equiv\langle\theta^{2}({\bf
x},t)\rangle_{0}=\frac{\overline{D}}{4\pi\Gamma}\ln\left(\frac{L}{a_{0}}\right)$
(8)

in 2D; $a_{0}$ is a small-scale cutoff. Equation (8) corresponds to a
logarithmically rough Edward-Wilkinson (EW) Barabási  _et al._ (1995) surface
at 2D. Furthermore

$C_{\theta\theta}^{0}(r)\equiv\langle\left[\theta({\bf x+r},t)-\theta({\bf
x},t)\right]^{2}\rangle_{0}\approx\frac{\overline{D}}{2\pi\Gamma}\ln(r/a_{0})$
(9)

for large $r$ Chaikin and Lubensky (1995). Equations (8) and (9) imply QLRO.
The equal-time density correlator

$C_{cc}^{0}(k)\equiv\langle|\delta c({\bf
k},t)|^{2}\rangle_{0}\approx\frac{\overline{D}}{2\Gamma}$ (10)

is independent of $k$ in the long wavelength limit. In real space,
$C^{0}_{cc}(r)$ vanishes for $r\gg\zeta$, a microscopic length characterizing
the short-range interactions. This further means that the standard deviation
of the number fluctuations of $N$ spins contained in a fixed open area,
$\sqrt{\langle N^{2}\rangle-\langle N\rangle^{2}}\equiv\sigma(N_{0})$ scales
with the mean $\langle N\rangle\equiv N_{0}$ as $\sqrt{N_{0}}$, giving NNF as
expected in a non-critical equilibrium system with short-range interactions.
If $\Omega_{1}\tilde{\Omega}_{0}\lambda_{0}<0$, the uniform states are
linearly unstable with a growth rate proportional to $k$.

The Fluctuation-Dissipation-Theorem (FDT) Chaikin and Lubensky (1995) is
broken in the linearized theory. This manifests in the non-vanishing cross-
correlator $C_{\times}({\bf k})\equiv\langle\theta({\bf-k},t)\delta c({\bf
k},t)\rangle$, which is a model parameter-dependent constant in the limit of
small $k$.

It now behooves us to find whether nonlinear effects are relevant (in the RG
or renormalization group sense), and the scaling properties of any ordered
states that are robust against finite noises. To study this, we perform one-
loop perturbative RG analysis on Eqs. (5) and (6) at 2D, similar to the RG
calculations on the KPZ equation Forster  _et al._ (1977); Barabási  _et al._
(1995), or the coupled Burgers-like equation for Magnetohydrodynamics Basu
_et al._ (1999); Basu and Frey (2004, 2009). It turns out, as discussed below,
that the nonlinearities either introduce logarithmic modulations to the
scaling of the linearly stable states, or to destroy those states altogether.

As usual, the RG is done by tracing over the short wavelength Fourier modes of
the fields Forster  _et al._ (1977); Barabási  _et al._ (1995); Hohenberg and
Halperin (1977), by expanding in the dimensionless coupling constant
$g\sim\lambda^{2}\overline{D}/\Gamma^{3}$; see ALP for technical details. It
predicts that not all but some of the linearly stable states are robust
against noises. This is controlled by
$\mu_{1}\equiv\Omega_{2}/\lambda,\,\mu_{2}\equiv\lambda_{0}\tilde{\Omega}_{1}/\lambda$,
that are marginal (in the RG sense) at the one-loop order. We find that for
wide ranges of $\mu_{1},\,\mu_{2}$, renormalized scale-dependent
$g(k)=1/[\Delta_{1}(\mu_{1},\mu_{2})\ln(\Lambda/k)]$ flows to zero very
slowly, where $\Delta_{1}>0$ as $k\rightarrow 0$ at 2D, under successive
applications of the RG procedure, for stable ordered phases;
$\Lambda=2\pi/a_{0}$ is an upper wavevector cutoff. This gives, as obtained
from the RG flow equations, that renormalized, scale-dependent
$\Gamma(k)=\Gamma[\ln(\Lambda/k)]^{\eta_{2}}$ and
$\overline{D}(k)=\overline{D}[\ln(\Lambda/k)]^{\eta_{1}}$, both diverging
logarithmically when $k\rightarrow 0$;
$\eta_{1}(\mu_{1},\mu_{2}),\,\eta_{2}(\mu_{1},\mu_{2})>0$ are constants
related to $\Delta_{1}(>0)$. These log-divergences due to the “slow” vanishing
of the coupling constant are reminiscent of the logarithmic anomalous
elasticity in three-dimensional equilibrium smectics Grinstein and Pelcovits
(1981, 1982). The resulting renormalized theory, owing to $g(k)\rightarrow 0$
as $k\rightarrow 0$, is effectively linear, albeit with renormalized
parameters. Straightforward calculations of the renormalized correlation
functions show that the stable ordered states have scaling properties
essentially same as in the linear theory, modulated only by logarithmic
corrections, which either strengthen or weaken the linear theory order. For
instance, the renormalized phase correlator $C_{\theta\theta}^{R}(k)$ reads

$C^{R}_{\theta\theta}(k)\approx\frac{\overline{D}}{2\Gamma
k^{2}[\ln(\Lambda/k)]^{\eta}},$ (11)

for $k\rightarrow 0$; $\eta\equiv\eta_{2}-\eta_{1}$. Here and below, $R$
refers to renormalized quantities. Exponent $\eta$ varies continuously with
$\mu_{1},\,\mu_{2}$ and can be positive or negative; detailed calculations
show $\eta<1/3$ always Haldar  _et al._ (2021). For $\eta>(<)0$,
$C^{R}_{\theta\theta}(k)\ll(\gg)C^{0}_{\theta\theta}(k)$ as $k\rightarrow 0$,
demonstrating strong suppression (enhancement) of fluctuations in the long
wavelength limit. Next, the renormalized variance
$\Delta_{\theta}^{R}\equiv\langle\theta^{2}({\bf x},t)\rangle_{R}$ now
acquires a novel $L$-dependence:

$\Delta_{\theta}^{R}\approx\frac{\overline{D}}{4\pi\Gamma}\left[\ln\left(\frac{L}{a_{0}}\right)\right]^{1-\eta},$
(12)

This means $\Delta_{\theta}^{R}$ can grow, respectively, logarithmically
slower or faster with $\ln L$ than QLRO, representing order stronger or weaker
than QLRO, named, respectively, SQLRO, or WQLRO; for $\eta=0$, QLRO is
retrieved. Equation (12) also implies a surface logarithmically smoother or
rougher than the 2D EW surface. Likewise, the renormalized correlator
$C_{\theta\theta}^{R}(r)$, related to the inverse Fourier transform of
$C_{\theta\theta}^{R}(k)$, scales as

$C_{\theta\theta}^{R}(r)\equiv\langle[\theta({\bf x+r},t)-\theta({\bf
x},t)]^{2}\rangle_{R}\approx\frac{\overline{D}}{2\pi\Gamma}[\ln(r/a_{0})]^{1-\eta},$
(13)

for large $r$. Related to $C^{R}_{\theta\theta}(r)$, the renormalized spin
correlation function $C_{ZZ}^{R}(r)$ for large $r$ is

$C^{R}_{ZZ}(r)\equiv\langle\cos[\theta({\bf x+r},t)-\theta({\bf
x},t)]\rangle_{R}\approx(r/a_{0})^{-\tilde{\gamma}(r)}$ (14)

where the $r$-dependent exponent $\tilde{\gamma}(r)$ has a complex form:

$\tilde{\gamma}(r)\equiv\frac{\overline{D}}{4\pi\Gamma}[\ln(r/a_{0})]^{-\eta}$
(15)

for large $r$. For $\eta>0(<0)$, clearly $C_{ZZ}(r)$ decays much slower
(faster) for large $r$, giving SQLRO (WQLRO).

Next in contrast to its linear theory analog, renormalized density correlator
$C_{cc}^{R}(k)$ picks up a weak $k$-dependence:

$C_{cc}^{R}(k)\approx\frac{\overline{D}}{2\Gamma}\left[\ln\left(\frac{\Lambda}{k}\right)\right]^{-\eta}$
(16)

in the hydrodynamic limit $k\rightarrow 0$. Evidently, for $\eta>0$, the
density fluctuations are strongly suppressed in the limit $k\rightarrow 0$,
implying miniscule number fluctuations (MNF) or hyperuniformity Torquato and
Stillinger (2003) vis-a-vis for $\eta=0$, same as the equilibrium result. In
contrast, if $\eta<0$ the density fluctuations are hugely enhanced when
$k\rightarrow 0$. This is giant number fluctuations (GNF), often encountered
in orientationally ordered active fluids Marchetti  _et al._ (2013), and also
in equilibrium superfluids Chaikin and Lubensky (1995). Lastly, the equal-time
renormalized density autocorrelator $C_{cc}^{R}(r)$, the inverse Fourier
transform of $C_{cc}^{R}(k)$, is

$C^{R}_{cc}(r)\approx\frac{\overline{D}}{4\pi\Gamma}\frac{-\eta}{r^{2}[\ln(r/a_{0})]^{(1+\eta)}},$
(17)

for large $r$, $r/a_{0}\gg 1$. Thus for $\eta>(<)0$, $C^{R}_{cc}(r)$ falls off
relatively faster (slower) for MNF (GNF). By using (17) we re-express MNF
(GNF) for $\eta>(<)0$:

$\sigma(N_{0})\propto\sqrt{N_{0}/(\ln N_{0})^{\eta}}<(>)\sqrt{N_{0}},$ (18)

in the renormalized theory. For $\eta=0$, (15) and (18) reduce to the well-
known results for the 2D equilibrium XY model and an equilibrium system with
NNF having short-range interactions away from any critical point.

We thus show SQLRO (WQLRO) phase order is necessarily accompanied by MNF
(GNF).

Like $C_{cc}^{R}(k)$, renormalized cross-correlation function
$C_{\times}^{R}(k)$ also picks up a weak $k$-dependence:
$C_{\times}^{R}(k)\propto\left[\ln(\frac{\Lambda}{k})\right]^{(1-3\eta)/2}\times{\cal
O}(1)$ for small $k$. Thus FDT remains broken in the renormalized theory.

The model shows breakdown of conventional dynamic scaling in the ordered
phases: the form of $\Gamma(k)$ in the renormalized theory shows that the
diffusive scaling of time $t$ with $r$ is modulated by logarithmic corrections
as

$t\propto r^{2}/[\ln(r/a_{0})]^{(1-\eta)/2}$ (19)

for large $r$. Thus the fluctuations relax logarithmically faster than
ordinary diffusion or in the linear theory. Also, they relax faster with WQLRO
($\eta<0$) than SQLRO ($\eta>0$).

Equations (11)-(19) collectively define the new universality class,
characterized by the SQLRO-MNF and WQLRO-GNF correspondence and continuously
varying scaling exponents.

For other choices of $\mu_{1},\,\mu_{2}$, uniform ordered states get unstable
due to nonlinear effects. In the equilibrium limit, the nonlinearities in (5)
and (6) vanish, and QLRO phase-order together with density NNF follows. A
phase diagram of the model in the $\mu_{1}-\mu_{2}$ plane over limited ranges
of $\mu_{1},\,\mu_{2}$ is shown in Fig. 2(a).

Numerical results from our agent-based model indeed confirm the existence of
stable order with SQLRO/MNF and WQLRO/GNF as generic, holding beyond any
perturbative calculations, and that topological defects in the active XY model
(not included in the hydrodynamic theory) do not proliferate and destroy order
at low-enough non-zero noises. It consists of $N_{\text{tot}}$ XY spins on a
2D square lattice of size $L\times L$, with periodic boundary conditions and a
mean spin number density $c_{0}\equiv N_{\text{tot}}/L^{2}$. The microscopic
update rules for $\theta$ depend on $c$ and the lattice analogs of the phase
difference-dependent active rotations of the spins along with equilibrium
relaxation processes; $c$ changes via hopping to the nearest neighbor sites,
provided either the density or the phase in the randomly chosen neighboring
target site is different. Stochasticity is introduced in the update of
$\theta$ by adding a (small) white noise $\xi$ to it Vicsek  _et al._ (1995),
and to $c$ through hopping of particles in randomly selected sites to randomly
chosen nearest neighbors; see ALP for details. Our numerical studies reveal
stable ordered states at finite noise $\xi$; see movie MOV1 sup . Furthermore,
in these ordered states, depending upon certain parameters $g_{1},\,g_{2}$,
which are the substitutes of $\mu_{2},\,\mu_{1}$ here (see ALP for precise
correspondence) and a few other parameters that define the relative
probabilities of the various microscopic dynamical processes Haldar  _et al._
(2021), numerical analogs of $\Delta_{\theta}^{R}(L)$ can grow slower and
faster than $\ln(L/a_{0})$, implying SQLRO and WQLRO respectively. In
addition, $\sigma(N_{0})$ as a function of $\ln\,N_{0}$ in a given open area
grows logarithmically slower and faster than $\sqrt{N_{0}}$, showing MNF and
GNF, correspondingly with SQLRO and WQLRO, respectively. Generic dependence of
$\eta$ on $g_{1},\,g_{2}$ are also found; see ALP Haldar  _et al._ (2021). Two
representative plots from the simulations of the agent-based model in Fig.
2(b) and (c) show the $\ln L$-dependence of $\langle\theta^{2}\rangle$ with
SQLRO and WQLRO, and correspondingly, $\ln\,N_{0}$-dependence of
$\sigma(N_{0})$ displaying MNF and GNF, respectively. These establish the
SQLRO (WQLRO) and MNF (GNF) correspondence in the agent-based model, in
agreement with the hydrodynamic theory 111The quantitative accuracy of these
plots is restricted due to the limited range of $L$ used in the simulations..

Figure 2: (a) Schematic phase diagram in the $\mu_{1}-\mu_{2}$ plane in the
linearly stable hydrodynamic theory. SQLRO/MNF, WQLRO/GNF and disordered
regions are marked in different colors. QLRO/NNF is the (blue) line
demarcating SQLRO/MNF and WQLRO/GNF regions. Plots from agent-based
simulations: (b) $\langle\theta^{2}\rangle/\ln L$ versus $\ln L$ showing SQLRO
and WQLRO, and (c) $\sigma(N_{0})/\sqrt{N_{0}}$ versus $\ln\,N_{0}$ showing
MNF and GNF ($L=128$) with $c_{0}=5,\,\xi=0.1$ and (SQLRO/MNF):
$g_{1}=1.0,\,g_{2}=0.02$, and (WQLRO/GNF): $g_{1}=2.0,\,g_{2}=0.03$,
respectively. The red broken horizontal lines in (b) and (c) correspond to the
equilibrium results. The SQLRO-MNF and WQLRO-GNF correspondences are clearly
established in the agent-based model (see text).

The simulations also detect two kinds of disordered states, corresponding to
the two different routes to disorder (linear instability and nonlinear
effects) in the hydrodynamic theory; see ALP Haldar  _et al._ (2021) for
details.

For immobile spins with uniform density, $\delta c=0$, or if there is no
number conservation, e.g., if the spins are created and destroyed, $\delta c$
relaxes fast, and is a non-hydrodynamic variable. Then Eq. (5) reduces to the
KPZ equation with only SRO in 2D Kardar  _et al._ (1986); Barabási  _et al._
(1995); Altman  _et al._ (2015), in agreement with our numerical studies in
its immobile limit (movie MOV2 sup ). Heuristically, fluctuations of the
conserved density $c$ create an effective long-range interactions between
phases at distant locations, which for appropriately chosen parameters can
suppress the instability of the 2D KPZ equation, ultimately sustaining order.

We have, therefore, presented the hydrodynamic theory of phase-ordering in
diffusively mobile active XY model on a substrate, supplemented by numerical
studies of an equivalent agent-based model. Existence of stable SQLRO and
WQLRO phase orders, correspondingly with MNF and GNF, respectively, uncovering
a novel correspondence, are predicted in wide ranges of the parameter space.
This is our central conclusion. Absence of ordered states elsewhere in the
parameter space is also revealed.

Acknowledgement:- AH and AB thank T. Banerjee and N. Sarkar for critical
reading and helpful suggestions. AS thanks University Grants Commission (UGC),
India for a research fellowship. SC is financially supported by the German
Research Foundation (DFG) within the Collaborative Research Center SFB 1027
and Indian Association for the Cultivation of Science, Kolkata. AS and SC also
acknowledges Prof. R. Paul (IACS, Kolkata) for some invaluable suggestions and
discussions and are also thankful to him for providing the computational
resources. AB thanks J. Toner for valuable discussions in the early stages of
this work, and comments concerning derivations of Eqs. (17) and (18), and P.
K. Mohanty, A. Maitra and D. Levis for helpful comments, and the SERB, DST
(India) for partial financial support through the MATRICS scheme [file no.:
MTR/2020/000406]. AB has designed the problem, AH has contributed to the
analytical part of the work, AS and SC have contributed equally to the
numerical part of the work. All four authors jointly wrote the manuscript.

## References

  * Marchetti  _et al._ (2013) M. C. Marchetti, J. F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao, and R. A. Simha, “Hydrodynamics of soft active matter,” Rev. Mod. Phys. 85, 1143–1189 (2013).
  * Toner and Tu (1995) J. Toner and Y. Tu, “Long-Range Order in a Two-Dimensional Dynamical $\mathrm{XY}$ Model: How Birds Fly Together,” Phys. Rev. Lett. 75, 4326–4329 (1995).
  * Mermin and Wagner (1966) N. D. Mermin and H. Wagner, “Absence of ferromagnetism or antiferromagnetism in one-or two-dimensional isotropic Heisenberg models,” Phys. Rev. Lett. 17, 1133 (1966).
  * Carusotto and Ciuti (2013) I. Carusotto and C. Ciuti, “Quantum fluids of light,” Rev. Mod. Phys. 85, 299–366 (2013).
  * Altman  _et al._ (2015) E. Altman, L. M. Sieberer, L. Chen, S. Diehl, and J. Toner, “Two-dimensional superfluidity of exciton polaritons requires strong anisotropy,” Phys. Rev. X 5, 011017 (2015).
  * Chaikin and Lubensky (1995) P. M. Chaikin and T. C. Lubensky,  _Principles of condensed matter physics_ (Cambridge university press Cambridge, 1995).
  * Vicsek  _et al._ (1995) T. Vicsek, A. Czirók, E. Ben-Jacob, I. Cohen, and O. Shochet, “Novel type of phase transition in a system of self-driven particles,” Phys. Rev. Lett. 75, 1226 (1995).
  * Kardar  _et al._ (1986) M. Kardar, G. Parisi, and Y.-C. Zhang, “Dynamic Scaling of Growing Interfaces,” Phys. Rev. Lett. 56, 889–892 (1986).
  * Barabási  _et al._ (1995) A-L Barabási, H. E Stanley, _et al._ ,  _Fractal concepts in surface growth_ (Cambridge university press, 1995).
  * (10) “See Supplemental Material for movies showing the time evolution of the spin configurations in a typical ordered phase (MOV1) and disordered phase without mobility (MOV2), respectively, of the model,” .
  * Pikovsky  _et al._ (2003) A. Pikovsky, J. Kurths, M. Rosenblum, and J. Kurths,  _Synchronization: a universal concept in nonlinear sciences_ , 12 (Cambridge university press, 2003).
  * Dörfler and Bullo (2014) F. Dörfler and F. Bullo, “Synchronization in complex networks of phase oscillators: A survey,” Automatica 50, 1539–1564 (2014).
  * Strogatz (2018) S. H Strogatz,  _Nonlinear dynamics and chaos with student solutions manual: With applications to physics, biology, chemistry, and engineering_ (CRC press, 2018).
  * Uriu  _et al._ (2010) K. Uriu, Y. Morishita, and Y. Iwasa, “Random cell movement promotes synchronization of the segmentation clock,” Proc. Natl. Acad. Sci. U.S.A. 107, 4979–4984 (2010).
  * Uriu and Morelli (2014) K. Uriu and L. G. Morelli, “Collective cell movement promotes synchronization of coupled genetic oscillators,” Biophys. J. 107, 514–526 (2014).
  * Zhou  _et al._ (2016) J. Zhou, Y. Zou, S. Guan, Z. Liu, and S. Boccaletti, “Synchronization in slowly switching networks of coupled oscillators,” Sc. Rep. 6, 35979 (2016).
  * Levis  _et al._ (2017) Demian Levis, Ignacio Pagonabarraga, and Albert Díaz-Guilera, “Synchronization in Dynamical Networks of Locally Coupled Self-Propelled Oscillators,” Phys. Rev. X 7, 011028 (2017).
  * Frasca  _et al._ (2008) M. Frasca, A. Buscarino, A. Rizzo, L. Fortuna, and S. Boccaletti, “Synchronization of moving chaotic agents,” Phys. Rev. Lett. 100, 044102 (2008).
  * Peruani  _et al._ (2010) F. Peruani, E. M Nicola, and L. G. Morelli, “Mobility induces global synchronization of oscillators in periodic extended systems,” New J. Phys. 12, 093029 (2010).
  * Banerjee and Basu (2017) T. Banerjee and A. Basu, “Active hydrodynamics of synchronization and ordering in moving oscillators,” Phys. Rev. E 96, 022201 (2017).
  * Uchida and Golestanian (2010) N. Uchida and R. Golestanian, “Synchronization and Collective Dynamics in a Carpet of Microfluidic Rotors,” Phys. Rev. Lett. 104, 178103 (2010).
  * Mathijssen  _et al._ (2018) A. J. T. M. Mathijssen, F. Guzmán-Lastra, A. Kaiser, and H. Löwen, “Nutrient Transport Driven by Microbial Active Carpets,” Phys. Rev. Lett. 121, 248101 (2018).
  * D. R. Brumley and Goldstein (2012) T. J. Pedley D. R. Brumley, M. Polin and R. E. Goldstein, “Hydrodynamic Synchronization and Metachronal Waves on the Surface of the Colonial Alga Volvox carteri,” Phys. Rev. Lett. 109, 268102 (2012).
  * Bruot and Cicuta (2016) N. Bruot and P. Cicuta, “Realizing the Physics of Motile Cilia Synchronization with Driven Colloids,” Annu. Rev. Condens. Matter Phys 7, 323–348 (2016).
  * Zhang  _et al._ (2021) Z. Zhang, H. Yuan, Y. Dou, M. O. de la Cruz, and K. J. M. Bishop, “Quincke Oscillations of Colloids at Planar Electrodes,” Phys. Rev. Lett. 126, 258001 (2021).
  * et al (2020) H. Gu et al, “Magnetic cilia carpets with programmable metachronal waves,” Nat. Comm. 11, 2637 (2020).
  * Haldar  _et al._ (2021) A. Haldar, A. Sarkar, S. Chatterjee, and A. Basu, “Active XY model on a substrate: Density fluctuations and phase ordering,” companion long paper; ArXiv: 2105.03919 (2021).
  * Maitra  _et al._ (2020) A. Maitra, M. Lenz, and R. Voituriez, “Chiral Active Hexatics: Giant Number Fluctuations, Waves, and Destruction of Order,” Phys. Rev. Lett. 125, 238005 (2020).
  * Lahiri and Ramaswamy (1997) R. Lahiri and S. Ramaswamy, “Are Steadily Moving Crystals Unstable?” Phys. Rev. Lett. 79, 1150–1153 (1997).
  * Forster  _et al._ (1977) D. Forster, D. R Nelson, and M. J Stephen, “Large-distance and long-time properties of a randomly stirred fluid,” Phys. Rev. A 16, 732 (1977).
  * Basu  _et al._ (1999) A. Basu, J. K Bhattacharjee, and S. Ramaswamy, “Mean magnetic field and noise cross-correlation in magnetohydrodynamic turbulence: results from a one-dimensional model,” Eur. Phys. J. B 9, 725–730 (1999).
  * Basu and Frey (2004) A. Basu and E. Frey, “Novel universality classes of coupled driven diffusive systems,” Phys. Rev. E 69, 015101 (2004).
  * Basu and Frey (2009) A. Basu and E. Frey, “Scaling and universality in coupled driven diffusive models,” J. Stat. Mech.: Theory Exp. 2009, P08013 (2009).
  * Hohenberg and Halperin (1977) P. C. Hohenberg and B. I. Halperin, “Theory of dynamic critical phenomena,” Rev. Mod. Phys. 49, 435 (1977).
  * Grinstein and Pelcovits (1981) G. Grinstein and R. A. Pelcovits, “Anharmonic Effects in Bulk Smectic Liquid Crystals and Other ”One-Dimensional Solids”,” Phys. Rev. Lett. 47, 856–859 (1981).
  * Grinstein and Pelcovits (1982) G. Grinstein and R. A. Pelcovits, “Smectic-$A-C$ transition in three dimensions,” Phys. Rev. A 26, 2196–2217 (1982).
  * Torquato and Stillinger (2003) S. Torquato and F. H. Stillinger, “Local density fluctuations, hyperuniformity, and order metrics,” Phys. Rev. E 68, 041113 (2003).
  * Note (1) The quantitative accuracy of these plots is restricted due to the limited range of $L$ used in the simulations.


# A Zeroth-Order Momentum Method for Risk-Averse
Online Convex Games

Zifan Wang, Yi Shen, Zachary I. Bell, Scott Nivison, Michael M. Zavlanos, and
Karl H. Johansson *This work is supported in part by AFOSR under award
#FA9550-19-1-0169 and by NSF under award CNS-1932011.Zifan Wang and Karl H.
Johansson are with Division of Decision and Control Systems, School of
Electrical Enginnering and Computer Science, KTH Royal Institute of
Technology, and also with Digital Futures, SE-10044 Stockholm, Sweden. Email:
<EMAIL_ADDRESS>Shen and Michael M. Zavlanos are with the Department
of Mechanical Engineering and Materials Science, Duke University, Durham, NC,
USA. Email: {yi.shen478<EMAIL_ADDRESS>I. Bell and Scott
Nivison are with the Air Force Research Laboratory, Eglin AFB, FL, USA. Email:
{scott.nivison<EMAIL_ADDRESS>
###### Abstract

We consider risk-averse learning in repeated unknown games where the goal of
the agents is to minimize their individual risk of incurring significantly
high cost. Specifically, the agents use the conditional value at risk (CVaR)
as a risk measure and rely on bandit feedback in the form of the cost values
of the selected actions at every episode to estimate their CVaR values and
update their actions. A major challenge in using bandit feedback to estimate
CVaR is that the agents can only access their own cost values, which, however,
depend on the actions of all agents. To address this challenge, we propose a
new risk-averse learning algorithm with momentum that utilizes the full
historical information on the cost values. We show that this algorithm
achieves sub-linear regret and matches the best known algorithms in the
literature. We provide numerical experiments for a Cournot game that show that
our method outperforms existing methods.

##  I Introduction

In online convex games [1, 2], agents interact with each other in the same
environment in order to minimize their individual cost functions. Many
applications, including traffic routing [3] and online marketing [4], among
many others [5], can be modeled as online convex games. The individual cost
functions depend on the joint decisions of all agents and are typically
unknown but can be accessed via querying. Using this information, the agents
can sequentially select the best actions that minimize their expected
accumulated costs. This can be done using online learning algorithms whose
performance is typically measured using the notions of regret [5] that
quantify the gap between the agents’ online decisions and the best decisions
in hindsight. Many such online learning algorithms have been recently
proposed, and have been shown to achieve sub-linear regret, which indicates
that optimal decisions are eventually selected; see e.g., [1, 2, 6]. In this
paper, here we focus on online convex games for high-stakes applications,
where decisions that minimize the expected cost functions are not necessarily
desirable. For example, in portfolio management, a selection of assets with
the highest expected return is not necessarily desirable since it may have
high volatility that can lead to catastrophic losses. To capture unexpected
and catastrophic outcomes, risk-averse criteria have been widely proposed in
place of the expectation, such as the Sharpe Ratio [7] and Conditional Value
at Risk (CVaR) [8].

In this paper, we consider online convex games with risk-averse agents that
employ CVaR as risk measure. We assume that only cost function evaluations at
selected decisions are accessible. Despite the practical utility of CVaR as a
measure of risk in high-stakes applications, the theoretical analysis of risk-
averse learning methods that employ the CVaR is limited, mainly since the CVaR
values must be estimated from the distributions of the unknown cost functions.
To avoid this problem, [9] proposes a variational definition of CVaR that
allows to reformulate the computation of CVaR values into an optimization
problem by introducing additional variables. Building on this reformulation,
many risk-averse learning methods have been proposed, including [10, 11, 12],
and their convergence rates have been analysed. Specifically, [10] provides
performance guarantees for stochastic gradient descent-type algorithms in
risk-averse statistical learning problems; [11] addresses risk-averse online
convex optimization problems with bandit feedback; and [12] proposes an
adaptive sampling strategy for CVaR learning, and transforms the supervised
learning problem into a zero-sum game between the sampler and the learner.
However, all the above works focus on single-agent learning problems.

To the best of our knowledge, multi-agent risk-averse learning problems have
not been analyzed in the literature, with only a few exceptions [13], [14].
The work in [13] solves multi-armed bandit problems with finite actions, which
are different from the convex games with continuous actions considered here.
Most closely related to this paper is our recent work in [14] proposing a new
risk-averse learning algorithm for online convex games, where the agents
utilize the cost evaluations at each time step to update their actions and
achieve no-regret learning with high probability. To improve the learning rate
of this algorithm, [14] also proposes two variance reduction methods on the
CVaR estimates and gradient estimates, respectively. In this paper, we build
on this work to significantly improve the performance of the algorithm using a
notion of momentum that combines these variance reduction techniques.
Specifically, the proposed new algorithm reuses samples of the cost functions
from the full history of data while appropriately decreasing the weights of
out-of-date samples. The idea is that, using more samples, the unknown
distribution of the cost functions needed to estimate the CVaR values can be
better estimated. To do so, at each time step, the proposed algorithm
estimates the distribution of the cost functions using both the immediate
empirical distribution estimate and the distribution estimate from the
previous time step, which summarizes all past cost values. Then, to construct
the gradient estimates, more weight is allocated to more recent information,
similar to momentum gradient descent methods [15]. To reduce the variance of
the gradient estimate, we employ residual feedback [16], which uses previous
gradient estimates to update the current one.

Although [14] showed that reusing samples from the last iteration or residual
feedback can individually improve the performance of the algorithm, it also
identified that the combination of sample reuse and residual feedback is an
unexplored and theoretically nontrivial problem. In this work, we combine
sample reuse and residual feedback by proposing a novel zeroth-order momentum
method that reuses all historical samples, compared to samples from only the
last iteration. An additional important contribution of this work is that the
proposed zeroth-order momentum method subsumes Algorithm 3 in [14] as a
special case by appropriately choosing the momentum parameter. We show that
our proposed method theoretically matches the best known result in [14] but
empirically outperforms other existing methods.

The rest of the paper is organized as follows. In Section II, we formally
define the problem and provide some preliminary results. In Section III, we
summarize results in [14] that are needed to develop our proposed method. In
Section IV, we present our proposed method and analyze its regret. In Section
V, we numerically verify our method using a Cournot game example. Finally, in
Section VI, we conclude the paper.

##  II Problem definition

We consider a repeated game $\mathcal{G}$ with $N$ risk-averse agents. Each
agent selects an action $x_{i}$ from the convex set
$\mathcal{X}_{i}\subseteq\mathbb{R}^{d_{i}}$, and then receives a stochastic
cost
$J_{i}(x_{i},x_{-i},\xi_{i}):\mathcal{X}\times\Xi_{i}\rightarrow\mathbb{R}$,
where $x_{-i}$ represents all agents’ actions except for agent $i$, and
$\mathcal{X}=\Pi_{i=1}^{N}\mathcal{X}_{i}$ is the joint action space. We
sometimes instead write $J_{i}(x,\xi_{i})$ for ease of notation, where
$x=(x_{i},x_{-i})$ is the concatenated vector of all agents’ actions. The cost
$J_{i}(x_{i},x_{-i},\xi_{i})$ is stochastic, as $\xi_{i}\in\Xi_{i}$ is a
stochastic variable. We assume that the diameter of the convex set
$\mathcal{X}_{i}$ is bounded by $D_{x}$ for all $i=1,\ldots,N$. In what
follows, we make the following assumptions on the class of games we consider
in this paper.

######  Assumption 1.

The function $J_{i}(x_{i},x_{-i},\xi_{i})$ is convex in $x_{i}$ for every
$\xi_{i}\in\Xi_{i}$ and bounded by $U$, i.e., $|J_{i}(x,\xi_{i})|\leq U$, for
all $i=1,\ldots,N$.

######  Assumption 2.

$J_{i}(x,\xi_{i})$ is $L_{0}$-Lipschitz continuous in $x$ for every
$\xi_{i}\in\Xi_{i}$, for all $i=1,\ldots,N$.

Assumptions 1 and 2 are common assumptions in the analysis of online learning
in games [17].

The goal of the risk-averse agents is to minimize the risk of incurring
significantly high costs, for possibly different risk levels. In this work, we
utilize CVaR as the risk measure. For a given risk level $\alpha_{i}\in[0,1]$,
CVaR is defined as the average of the $\alpha_{i}$ percent worst-case cost.
Specifically, we denote by $F_{i}^{x}(y)=\mathbb{P}\\{J_{i}(x,\xi_{i})\leq
y\\}$ as the cumulative distribution function (CDF) of the random cost
$J_{i}(x,\xi_{i})$ for agent $i$; and by $J^{\alpha_{i}}$ the cost value at
the $1-\alpha_{i}$ quantile of the distribution, also called the Value at Risk
(VaR). Then, for a given risk level $\alpha_{i}\in[0,1]$, CVaR of the cost
function $J_{i}(x,\xi_{i})$ of agent $i$ is defined as

$\displaystyle C_{i}(x):$
$\displaystyle={\rm{CVaR}}_{\alpha_{i}}[J_{i}(x,\xi_{i})]$ $\displaystyle:$
$\displaystyle=\mathbb{E}_{F}[J_{i}(x,\xi_{i})|J_{i}(x,\xi_{i})\geq
J^{\alpha_{i}}].$

Notice that the CVaR value is determined by the distribution function
$F_{i}^{x}(y)$ for a given $\alpha_{i}$, so we sometimes write CVaR as a
function of the distribution function, i.e.,
${\rm{CVaR}}_{\alpha_{i}}[F_{i}^{x}]:={\rm{CVaR}}_{\alpha_{i}}[J_{i}(x,\xi_{i})]$.
In addition, we assume that the agents cannot observe other agents’ actions,
but the only information they can observe is the cost evaluations of the
jointly selected actions at each time step.

Given Assumptions 1 and 2, the following lemma lists properties of CVaR, that
will be important in the analysis that follows. The proof can be found in
[11].

######  Lemma 1.

Given Assumptions 1 and 2, we have that $C_{i}(x_{i},x_{-i})$ is convex in
$x_{i}$ and $L_{0}$-Lipschitz continuous in $x$, for all $i=1,\ldots,N$.

The following additional assumption on the variation of the CDFs is needed for
the analysis that follows.

######  Assumption 3.

Let $F^{w}_{i}(y)=\mathbb{P}\\{J_{i}(w,\xi_{i})\leq y\\}$ and
$F^{v}_{i}(y)=\mathbb{P}\\{J_{i}(v,\xi_{i})\leq y\\}$. There exist a constant
$C_{0}>0$ such that

$\displaystyle\mathop{\rm{sup}}_{y}|F^{w}_{i}(y)-F^{v}_{i}(y)|\leq
C_{0}\left\|w-v\right\|.$

Assumption 3 states that the difference between two CDFs can be bounded by the
difference between the corresponding two action profiles. It is less
restrictive than Assumption 3 in [14] since it does not depend on the
algorithm but only on the cost functions.

A common measure of the ability of the agents to learn online optimal
decisions that minimize their individual risk-averse objective functions
$C_{i}(x)$ is the algorithm regret, which is defined as the difference between
the actual rewards and best rewards that the agent could have achieved by
playing the single best action in hindsight. Suppose the action sequences of
agent $i$ and the other agents in the team are $\\{\hat{x}_{i,t}\\}_{t=1}^{T}$
and $\\{\hat{x}_{-i,t}\\}_{t=1}^{T}$, respectively. Then, we define the regret
of agent $i$ as

$\displaystyle{\rm{R}}_{C_{i}}(T)=\sum_{t=1}^{T}C_{i}(\hat{x}_{i,t},\hat{x}_{-i,t})-\mathop{\rm{min}}_{\tilde{x}_{i}\in\mathcal{X}_{i}}\sum_{t=1}^{T}C_{i}(\tilde{x}_{i},\hat{x}_{-i,t}).$

An algorithm is said to be no-regret if its regret grows sub-linearly with the
number of episodes $T$. In this paper, we propose a no-regret learning
algorithm, so that $\lim_{T\rightarrow\infty}\frac{{\rm{R}}_{C_{i}}(T)}{T}=0$,
$i=1,\ldots,N$.

##  III Preliminary Results

In this section, we summarize some results from [14] that lay the foundation
for the subsequent analysis. Note that the CVaR values depend on the
distributions of the cost functions, which are generally unknown. To estimate
the distribution of the cost functions, and subsequently the CVaR values with
only a few samples, a sampling strategy is proposed in [14] that uses a
decreasing number of samples with the number of iterations. Specifically, at
episode $t$, the sampling strategy used by the agents is designed as follows

$\displaystyle n_{t}=\lceil bU^{2}(T-t+1)^{a}\rceil,$ (1)

where $\lceil\cdot\rceil$ is the regularized ceiling function, $T$ is the
number of episodes, $U$ is the bound of $J_{i}$ as in Assumption 1, and
$a,b\in(0,1)$ are parameters to be selected later.

Using the cost evaluations at each episode, zeroth-order methods can be
employed to estimate the CVaR gradients and then update the agents’ actions.
Specifically, at episode $t$, the agents calculate the number of samples
$n_{t}$ according to (1). Then, each agent perturbs the current action
$x_{i,t}$ by an amount $\delta u_{i,t}$, where $u_{i,t}\in\mathbb{S}^{d_{i}}$
is a random perturbation direction sampled from the unit sphere
$\mathbb{S}^{d_{i}}\subset\mathbb{R}^{d_{i}}$ and $\delta$ is the size of this
perturbation. Next the agents play their perturbed actions
$\hat{x}_{i,t}=x_{i,t}+\delta u_{i,t}$ for $n_{t}$ times, and obtain $n_{t}$
cost evaluations which are utilized to update their actions. To facilitate the
theoretical analysis, we define the $\delta$-smoothed function
$C_{i}^{\delta}(x)=\mathbb{E}_{w_{i}\sim\mathbb{B}_{i},u_{-i}\sim\mathbb{S}_{-i}}[C_{i}(x_{i}+\delta
w_{i},x_{-i}+\delta u_{-i})]$, where $\mathbb{S}_{-i}=\Pi_{j\neq
i}\mathbb{S}_{j}$, and $\mathbb{B}_{i}$, $\mathbb{S}_{i}$ denote the unit ball
and unit sphere in $\mathbb{R}^{d_{i}}$, respectively. The size of the
perturbation $\delta$ here is related to a smoothing parameter that controls
how well $C_{i}^{\delta}(x)$ approximates $C_{i}(x)$. As shown in [14], the
function $C_{i}^{\delta}(x)$ satisfies the following properties.

######  Lemma 2.

Let Assumptions 1 and 2 hold. Then we have that

  1. 1.

$C_{i}^{\delta}(x_{i},x_{-i})$ is convex in $x_{i}$,

  2. 2.

$C_{i}^{\delta}(x)$ is $L_{0}$-Lipschitz continuous in $x$,

  3. 3.

$|C_{i}^{\delta}(x)$-$C_{i}(x)|\leq\delta L_{0}\sqrt{N}$,

  4. 4.

$\mathbb{E}[\frac{d_{i}}{\delta}C_{i}(\hat{x}_{t})u_{i,t}]=\nabla_{i}C_{i}^{\delta}(x_{t})$.

The last property in Lemma 2 shows that the term
$\frac{d_{i}}{\delta}C_{i}(\hat{x}_{t})u_{i,t}$ is an unbiased estimate of the
gradient of the smoothed function $C_{i}^{\delta}(x)$. Next we present a lemma
that helps bound the distance between two CVaR values by the distance between
the two corresponding CDFs. The proof can be found in [14].

######  Lemma 3.

Let $F$ and $G$ be two CDFs of two random variables and suppose the random
variables are bounded by $U$. Then we have that

$\displaystyle|{\rm{CVaR}}_{\alpha}[F]-{\rm{CVaR}}_{\alpha}[G]|\leq\frac{U}{\alpha}\mathop{\rm{sup}}_{y}|F(y)-G(y)|.$

##  IV A zeroth-order momentum method

In this section, we present a new one-point zeroth-order momentum method that
combines sample reuse as in Algorithm 2 in [14] and residual feedback as in
Algorithm 3 in [14], but unlike Algorithm 2 in [14] that reuses samples only
from the last iteration, it uses all past samples to update the agents’
actions. The new algorithm is given as Algorithm 1.

Algorithm 1 Risk-averse learning with momentum

0: Initial value $x_{0}$, step size $\eta$, parameters $a$, $b$, $\delta$,
$T$, risk level $\alpha_{i}$, $i=1,\cdots,N$.

1: for ${\rm{episode}}\;t=1,\ldots,T$ do

2: Select $n_{t}=\lceil bU^{2}(T-t+1)^{a}\rceil$

3: Each agent samples $u_{i,t}\in\mathbb{S}^{d_{i}}$, $i=1,\ldots,N$

4: Each agent play $\hat{x}_{i,t}=x_{i,t}+\delta u_{i,t}$, $i=1,\ldots,N$

5: for $j=1,\ldots,n_{t}$ do

6: Let all agents play $\hat{x}_{i,t}$

7: Obtain $J_{i}(\hat{x}_{i,t},\hat{x}_{-i,t},\xi_{i}^{j})$

8: end for

9: for agent $i=1,\ldots,N$ do

10: Build EDF $\bar{F}_{i,t}(y)$

11: Estimate CVaR: ${\rm{CVaR}}_{\alpha_{i}}[\bar{F}_{i,t}]$

12: Construct gradient
estimate$\bar{g}_{i,t}=\frac{d_{i}}{\delta}\left({\rm{CVaR}}_{\alpha_{i}}[\bar{F}_{i,t}]-{\rm{CVaR}}_{\alpha_{i}}[\bar{F}_{i,t-1}]\right)u_{i,t}$

13: Update $x$:
$x_{i,t+1}\leftarrow\mathcal{P}_{\mathcal{X}_{i}^{\delta}}(x_{i,t}-\eta\bar{g}_{i,t})$

14: end for

15: end for

Using the sampling strategy in (1), each agent plays the perturbed action
$\hat{x}_{i,t}$ for $n_{t}$ times and obtains $n_{t}$ samples at episode $t$.
For agent $i$, we denote the CDF of the random cost
$J_{i}(\hat{x}_{t},\xi_{i})$ that is returned by the perturbed action
$\hat{x}_{t}$ as
$F_{i,t}^{\hat{x}_{t}}(y)=\mathbb{P}\\{J_{i}(\hat{x}_{t},\xi_{i})\leq y\\}$.
Since $\hat{x}_{t}$ depends on $t$, we write $F_{i,t}$ for
$F_{i,t}^{\hat{x}_{t}}$ for ease of notation. Using bandit feedback in the
form of finitely many cost evaluations, the agents cannot obtain the accurate
CDF $F_{i,t}$. Instead, they construct an empirical distribution function
(EDF) $\hat{F}_{i,t}$ of the cost $J_{i}(\hat{x}_{t},\xi_{i})$ using $n_{t}$
cost evaluations by

$\displaystyle\hat{F}_{i,t}(y)=\frac{1}{n_{t}}\sum_{j=1}^{n_{t}}\mathbf{1}\\{J_{i}(\hat{x}_{t},\xi_{i}^{j})\leq
y\\}.$ (2)

To improve estimate of the CDF, we utilize past samples for all episodes
$t\geq 2$, and construct a modified distribution estimate $\bar{F}_{i,t}$ by
adding a momentum term:

$\displaystyle\bar{F}_{i,t}(y)=\beta\bar{F}_{i,t-1}(y)+(1-\beta)\hat{F}_{i,t}(y),$
(3)

where $\beta\in[0,1)$ is the momentum parameter. For $t=1$, we set
$\bar{F}_{i,t}(y)=\hat{F}_{i,t}(y)$. The agents utilize the distribution
estimate $\bar{F}_{i,t}$ to calculate the CVaR estimates and further construct
the gradient estimate using residual feedback as

$\displaystyle\bar{g}_{i,t}=\frac{d_{i}}{\delta}\left({\rm{CVaR}}_{\alpha_{i}}[\bar{F}_{i,t}]-{\rm{CVaR}}_{\alpha_{i}}[\bar{F}_{i,t-1}]\right)u_{i,t},$
(4)

where $\delta$ is the size of the perturbation on the action $x_{i,t}$ defined
above. Depending on the size of the perturbation, the played action
$\hat{x}_{i,t}$ may be infeasible, i.e.,
$\hat{x}_{i,t}\not\in\mathcal{X}_{i}$. To handle this issue, we define the
projection set
$\mathcal{X}_{i}^{\delta}=\\{x_{i}\in\mathcal{X}_{i}|{\rm{dist}}(x_{i},\partial\mathcal{X}_{i})\geq\delta\\}$
that we use in the projected gradient-descent update

$\displaystyle
x_{i,t+1}=\mathcal{P}_{\mathcal{X}_{i}^{\delta}}(x_{i,t}-\eta\bar{g}_{i,t}).$
(5)

Note that the agents are not able to obtain accurate values of
$C_{i}(\hat{x}_{t})$ using finite samples. In fact, if we use the distribution
estimate $\bar{F}_{i,t}$ in (3) to calculate the CVaR values, there will be a
CVaR estimation error, which is defined as

$\displaystyle\bar{\varepsilon}_{i,t}:={\rm{CVaR}}_{\alpha_{i}}[\bar{F}_{i,t}]-{\rm{CVaR}}_{\alpha_{i}}[F_{i,t}].$

As a result, the gradient estimate $\bar{g}_{i,t}$ in (4) is biased since
$\mathbb{E}[\bar{g}_{i,t}]=\mathbb{E}\left[\frac{d_{i}}{\delta}(C_{i}(\hat{x}_{t})+\bar{\varepsilon}_{i,t})u_{i,t}\right]=\nabla_{i}C_{i}^{\delta}(x_{t})+\mathbb{E}\left[\frac{d_{i}}{\delta}\bar{\varepsilon}_{i,t}u_{i,t}\right]$,
with the bias captured by the last term. In what follows, we present a lemma
that bounds the sum of the CVaR estimation errors.

######  Lemma 4.

Suppose that Assumption 3 holds. Then, the sum of the CVaR estimation errors
satisfies

$\displaystyle\sum_{t=1}^{T}\left\|\bar{\varepsilon}_{i,t}\right\|\leq\frac{U}{\alpha_{i}}\Big{(}\frac{\beta}{1-\beta}e_{i,1}+\frac{\beta\Sigma_{1}}{1-\beta}T+\sum_{t=1}^{T}r_{t}\Big{)}$
(6)

with probability at least $1-\gamma$, where
$e_{i,t}:=\mathop{\rm{sup}}_{y}|\bar{F}_{i,t}(y)-F_{i,t}(y)|$,
$r_{t}:=\sqrt{\frac{{\rm{ln}}(2T/\gamma)}{2bU^{2}(T-t+1)^{a}}}$, and
$\Sigma_{1}:=C_{0}\eta\frac{d_{i}}{\delta}U\sqrt{N}+2C_{0}\delta$.

###### Proof.

To bound the CVaR estimation error, it suffices to bound the CDF differences
using Lemma 3. Recall the definition of $\bar{F}_{i,t}$ in (3). Then, we have
that

$\displaystyle\mathop{\rm{sup}}_{y}|\bar{F}_{i,t}(y)-F_{i,t}(y)|$
$\displaystyle=$
$\displaystyle\mathop{\rm{sup}}_{y}|\beta\bar{F}_{i,t-1}(y)+(1-\beta)\hat{F}_{i,t}(y)-F_{i,t}(y)|$
$\displaystyle=$
$\displaystyle\mathop{\rm{sup}}_{y}|\beta(\bar{F}_{i,t-1}(y)-{F}_{i,t-1}(y))+\beta({F}_{i,t-1}(y)-{F}_{i,t}(y))$
$\displaystyle+(1-\beta)(\hat{F}_{i,t}(y)-{F}_{i,t}(y))|$ $\displaystyle\leq$
$\displaystyle\beta\mathop{\rm{sup}}_{y}|\bar{F}_{i,t-1}(y)-{F}_{i,t-1}(y)|+\beta\mathop{\rm{sup}}_{y}|{F}_{i,t-1}(y)-{F}_{i,t}(y)|$
$\displaystyle+(1-\beta)\mathop{\rm{sup}}_{y}|\hat{F}_{i,t}(y)-{F}_{i,t}(y)|.$
(7)

We now bound the latter two terms in the right-hand-side of (IV) separately.
Using Assumption 3 and Lemma 2, we have that

$\displaystyle\mathop{\rm{sup}}_{y}|F_{i,t}(y)-F_{i,t-1}(y)|$
$\displaystyle\leq$ $\displaystyle
C_{0}\left\|\hat{x}_{t}-\hat{x}_{t-1}\right\|\leq
C_{0}\left\|x_{t}-x_{t-1}\right\|+2C_{0}\delta$ $\displaystyle\leq$
$\displaystyle C_{0}\eta\left\|\bar{g}_{t}\right\|+2C_{0}\delta\leq
C_{0}\eta\frac{d_{i}}{\delta}U\sqrt{N}+2C_{0}\delta:=\Sigma_{1},$ (8)

where the last inequality holds since
$\left\|\bar{g}_{t}\right\|=\sqrt{\sum_{i=1}^{N}\left\|\bar{g}_{i,t}\right\|^{2}}\leq\frac{\sqrt{N}d_{i}U}{\delta}$.
Then, applying the Dvoretzky–Kiefer–Wolfowitz (DKW) inequality, we obtain that

$\displaystyle\mathbb{P}\left\\{\mathop{\rm{sup}}_{y}|F_{i,t}(y)-\hat{F}_{i,t}(y)|\geq\sqrt{\frac{{\rm{ln}}(2/\bar{\gamma})}{2n_{t}}}\right\\}\leq\bar{\gamma}.$
(9)

Define the events in (9) as $A_{t}$ and denote $\gamma=\bar{\gamma}T$. Then,
the following inequality holds

$\displaystyle\mathop{\rm{sup}}_{y}|F_{i,t}(y)-\hat{F}_{i,t}(y)|\leq\sqrt{\frac{{\rm{ln}}(2T/\gamma)}{2n_{t}}},\forall
t=1,\ldots,T,$ (10)

with probability at least $1-\gamma$ since
$1-\mathbb{P}\\{\bigcup_{t=1}^{T}A_{t}\\}\geq
1-\sum_{t=1}^{T}\mathbb{P}\\{A_{t}\\}\geq 1-T\frac{\gamma}{T}\geq 1-\gamma$.
Substituting the bounds in (IV) and (10) into (IV), and using the definition
of $e_{i,t}$, we get that

$\displaystyle e_{i,t}\leq\beta e_{i,t-1}+\beta\Sigma_{1}+(1-\beta)r_{t}.$
(11)

Note that (11) holds for all $t\in\\{2,\cdots,T\\}$. By iteratively using this
inequality, we have that

$\displaystyle
e_{i,t}\leq\frac{\beta}{1-\beta}\Sigma_{1}+(1-\beta)\sum_{k=0}^{k-2}\beta^{k}r_{t-k}+\beta^{t-1}e_{1,i}.$
(12)

Taking the sum over $t=1,\ldots,T$ of both sides of (12), we have that

$\displaystyle\sum_{t=1}^{T}e_{i,t}\leq\frac{\beta}{1-\beta}e_{i,1}+\frac{\beta\Sigma_{1}}{1-\beta}T+\sum_{t=1}^{T}r_{t},$
(13)

with probability at least $1-\gamma$. Finally, using Lemma 3, we can obtain
the desired result. ∎

Lemma 4 bounds the CVaR estimation error caused by using past samples. The
next lemma quantifies the error due to residual feedback in the gradient
estimate (4). Specifically, it bounds the second moment of the gradient
estimate.

######  Lemma 5.

Let Assumptions 1, 2 and 3 hold. Suppose that the action is updated as in (5).
Then, the second moment of the gradient estimate satisfies

$\displaystyle\sum_{t=1}^{T}\left\|\bar{g}_{t}\right\|^{2}\leq\frac{1}{1-\sigma}\left\|\bar{g}_{1}\right\|^{2}+\frac{1}{1-\sigma}\Sigma_{2}T,$
(14)

where $\bar{g}_{t}:=(\bar{g}_{i,t},\bar{g}_{-i,t})$ is the concatenated vector
of all agents’ gradient estimates,
$\sigma=\frac{4d_{i}^{2}L_{0}^{2}N\eta^{2}}{\delta^{2}}$,
$\Sigma_{2}=16d_{i}^{2}N^{2}L_{0}^{2}+\frac{48C_{0}^{2}d_{i}^{4}U^{4}N\eta^{2}\beta^{2}}{\delta^{4}(1-\beta)^{2}}(\sum_{i}\frac{1}{\alpha_{i}^{2}})+\frac{192C_{0}^{2}d_{i}^{2}U^{2}\beta^{2}}{(1-\beta)^{2}}(\sum_{i}\frac{1}{\alpha_{i}^{2}})+\frac{12\ln(2T/\gamma)d_{i}^{2}}{b\delta^{2}}(\sum_{i}\frac{1}{\alpha_{i}^{2}})+\frac{24d_{i}^{2}U^{2}\beta^{2}e_{m}^{2}}{\delta^{2}}(\sum_{i}\frac{1}{\alpha_{i}^{2}})$,
and $e_{m}:=\mathop{\max}_{i}\\{e_{i,1}\\}$.

###### Proof.

Using the fact that $r_{t}\leq\sqrt{\frac{{\rm{ln}}(2T/\gamma)}{2bU^{2}}}:=r$
and (12), we have that

$\displaystyle\left\|e_{i,t}\right\|^{2}$
$\displaystyle\leq\left\|\frac{\beta}{1-\beta}\Sigma_{1}+r+\beta
e_{m}\right\|^{2}$ $\displaystyle\leq
3\left\|\frac{\beta}{1-\beta}\Sigma_{1}\right\|^{2}+3\left\|r\right\|^{2}+3\left\|\beta
e_{m}\right\|^{2},$ (15)

where the last inequality holds since $\left\|a+b+c\right\|^{2}\leq
3\left\|a\right\|^{2}+3\left\|b\right\|^{2}+3\left\|c\right\|^{2}$. Then,
applying Lemma 3 to (IV), we obtain that

$\displaystyle\left\|\bar{\varepsilon}_{i,t}\right\|^{2}\leq\frac{U^{2}}{\alpha_{i}^{2}}\left(\frac{3\beta^{2}\Sigma_{1}^{2}}{(1-\beta)^{2}}+3\left\|r\right\|^{2}+3\beta^{2}\left\|e_{m}\right\|^{2}\right).$
(16)

By the definition of $\bar{g}_{i,t}$ in (4), we have that

$\displaystyle\left\|\bar{g}_{i,t}\right\|^{2}=\frac{d_{i}^{2}}{\delta^{2}}\left(({\rm{CVaR}}_{\alpha_{i}}[\bar{F}_{i,t}]-{\rm{CVaR}}_{\alpha_{i}}[\bar{F}_{i,t-1}])u_{i,t}\right)^{2}$
$\displaystyle\leq\frac{d_{i}^{2}}{\delta^{2}}\left({\rm{CVaR}}_{\alpha_{i}}[\bar{F}_{i,t}]-{\rm{CVaR}}_{\alpha_{i}}[\bar{F}_{i,t-1}]\right)^{2}\left\|u_{i,t}\right\|^{2}$
$\displaystyle\leq\frac{d_{i}^{2}}{\delta^{2}}\Big{(}2({\rm{CVaR}}_{\alpha_{i}}[F_{i,t}]-{\rm{CVaR}}_{\alpha_{i}}[F_{i,t-1}])^{2}$
$\displaystyle\quad+2(\bar{\varepsilon}_{i,t}-\bar{\varepsilon}_{i,t-1})^{2}\Big{)}\left\|u_{i,t}\right\|^{2}$
$\displaystyle\leq\frac{d_{i}^{2}}{\delta^{2}}\Big{(}2L_{0}^{2}\left\|\hat{x}_{t}-\hat{x}_{t-1}\right\|^{2}+4|\bar{\varepsilon}_{i,t}|^{2}+4|\bar{\varepsilon}_{i,t-1}|^{2}\Big{)},$
(17)

where the last inequality is due to the fact that the CVaR function is
Lipschitz continuous as shown in Lemma 1. We now bound the term
$\left\|\hat{x}_{t}-\hat{x}_{t-1}\right\|^{2}$ in (IV). Recalling the update
rule
$x_{i,t}=\mathcal{P}_{\mathcal{X}_{i}^{\delta}}(x_{i,t-1}-\eta_{i}\bar{g}_{i,t-1})$,
we have that

$\displaystyle\left\|\hat{x}_{t}-\hat{x}_{t-1}\right\|^{2}=\left\|x_{t}-x_{t-1}+\delta
u_{t}-\delta u_{t-1}\right\|^{2}$ $\displaystyle\leq$ $\displaystyle
2\left\|x_{t}-x_{t-1}\right\|^{2}+2\delta^{2}(2\left\|u_{t}\right\|^{2}+2\left\|u_{t-1}\right\|^{2})$
$\displaystyle\leq$ $\displaystyle
2\left\|x_{t}-x_{t-1}\right\|^{2}+8N\delta^{2}\leq
2\eta^{2}\left\|\bar{g}_{t-1}\right\|^{2}+8N\delta^{2}.$

Substituting the above inequality and the inequality in (16) into the right
hand side of (IV) and summing both sides over all agents $i=1,\ldots,N$, we
have that

$\displaystyle\left\|\bar{g}_{t}\right\|^{2}=\sum_{i=1}^{N}\left\|\bar{g}_{i,t}\right\|^{2}$
$\displaystyle\leq\frac{d_{i}^{2}N}{\delta^{2}}\Big{(}4L_{0}^{2}\eta^{2}\left\|\bar{g}_{t-1}\right\|^{2}+16L_{0}^{2}N\delta^{2}\Big{)}$
$\displaystyle\quad+\frac{8d_{i}^{2}}{\delta^{2}}\sum_{i=1}^{N}\frac{U^{2}}{\alpha_{i}^{2}}\left(\frac{3\beta^{2}\Sigma_{1}^{2}}{(1-\beta)^{2}}+3\left\|r\right\|^{2}+3\beta^{2}\left\|e_{m}\right\|^{2}\right)$
$\displaystyle\leq\sigma\left\|\bar{g}_{t-1}\right\|^{2}+16d_{i}^{2}N^{2}L_{0}^{2}+\frac{12\ln(2T/\gamma)d_{i}^{2}}{b\delta^{2}}(\sum_{i}\frac{1}{\alpha_{i}^{2}})$
$\displaystyle\quad+\frac{24d_{i}^{2}U^{2}\beta^{2}}{\delta^{2}(1-\beta)^{2}}(C_{0}\eta\frac{d_{i}}{\delta}U\sqrt{N}+2C_{0}\delta)^{2}(\sum_{i}\frac{1}{\alpha_{i}^{2}})$
$\displaystyle\quad+\frac{24d_{i}^{2}U^{2}\beta^{2}e_{m}^{2}}{\delta^{2}}(\sum_{i}\frac{1}{\alpha_{i}^{2}})$
$\displaystyle\leq\sigma\left\|\bar{g}_{t-1}\right\|^{2}+\Sigma_{2},$ (18)

where the last inequality is due to the fact that
$(C_{0}\eta\frac{d_{i}}{\delta}U\sqrt{N}+2C_{0}\delta)^{2}\leq
2C_{0}^{2}\eta^{2}\frac{d_{i}^{2}}{\delta^{2}}U^{2}N+8C_{0}^{2}\delta^{2}$.
Summing (IV) over all episode $t=1,\ldots,T$ completes the proof. ∎

Before proving the main theorem, we present a regret decomposition lemma that
links the regret to the errors bounds on the CVaR estimates and gradient
estimates as provided in Lemmas 4 and 5.

######  Lemma 6.

Let Assumptions 1, 2 and 3 hold. Then, the regret of Algorithm 1 satisfies

$\displaystyle{\rm{R}}_{C_{i}}(T)\leq$
$\displaystyle\frac{D_{x}^{2}}{2\eta}+\frac{\eta}{2}\mathbb{E}\big{[}\sum_{t=1}^{T}\left\|\bar{g}_{i,t}\right\|^{2}\big{]}+(4\sqrt{N}+\Omega)L_{0}\delta
T$
$\displaystyle+\frac{d_{i}D_{x}}{\delta}\sum_{t=1}^{T}|\bar{\varepsilon}_{i,t}|,$
(19)

where $\Omega>0$ is a constant that represents the error from projection
$\mathcal{P}$.

###### Proof.

The proof can be adapted from Lemma 5 in [14] and is omitted due to the space
limitations. ∎

######  Theorem 1.

Let Assumptions 1, 2 and 3 hold and select
$\eta=\frac{D_{x}}{d_{i}L_{0}N}T^{-\frac{3a}{4}}$,
$\delta=\frac{D_{x}}{N^{\frac{1}{6}}}T^{-\frac{a}{4}}$,
$\beta=\frac{1}{U^{2}T^{\frac{a}{4}}}$. Suppose that $n_{t}$ is chosen as in
(1) with $a\in(0,1)$, and the EDF and the gradient estimate are defined as in
(2) and (4), respectively. Then, when $T\geq(8N^{\frac{2}{3}})^{\frac{1}{a}}$,
Algorithm 1 achieves regret

${\rm{R}}_{C_{i}}(T)=\mathcal{O}(D_{x}d_{i}L_{0}NS(\alpha)\ln(T/\gamma)T^{1-\frac{a}{4}}),$

with probability at least $1-\gamma$ and where
$S(\alpha):=\sum_{i=1}^{N}\frac{1}{\alpha_{i}^{2}}$.

###### Proof.

Substituting the bounds in Lemmas 4 and 5 into Lemma 6, we have that

$\displaystyle{\rm{R}}_{C_{i}}(T)\leq$
$\displaystyle\frac{D_{x}^{2}}{2\eta}+\frac{\eta}{2}\left(\frac{1}{1-\sigma}\left\|\bar{g}_{1}\right\|^{2}+\frac{1}{1-\sigma}\Sigma_{2}T\right)$
$\displaystyle+\frac{d_{i}D_{x}U}{\delta\alpha_{i}}\Big{(}\frac{\beta}{1-\beta}e_{i,1}+\frac{\beta\Sigma_{1}}{1-\beta}T+\sum_{t=1}^{T}r_{t}\Big{)}$
$\displaystyle+(4\sqrt{N}+\Omega)L_{0}\delta T.$ (20)

Substituting $\delta$, $\eta$ and $\beta$ into the above inequality, we can
obtain the desired result. The detailed proof is straightforward and is
omitted due to the space limitations. ∎

Note that Algorithm 1 in fact achieves the regret
${\rm{R}}_{C_{i}}(T)=\tilde{\mathcal{O}}(T^{1-\frac{a}{4}})$, in which we use
the notation $\tilde{\mathcal{O}}$ to hide constant factors and poly-
logarithmic factors of $T$. This result matches the best result in [14].
Compared to [14], the analysis of Algorithm 1 that combines historical
information and residual feedback is nontrivial and requires an appropriate
selection of the momentum parameter to show convergence. Note that Algorithm 3
in [14] is a special case of Algorithm 1 for the momentum parameter $\beta=0$.

##  V Numerical Experiments

In this section, we illustrate the proposed algorithm on a Cournot game.
Specifically, we consider two risk-averse agents $i=1,2$ that compete with
each other. Each agent determines the production level $x_{i}$ and receives an
individual cost feedback which is given by
$J_{i}=1-(2-\sum_{j}x_{j})x_{i}+0.2x_{i}+\xi_{i}x_{i}$, where $\xi_{i}\sim
U(0,1)$ is a uniform random variable. Here we utilize the term $\xi_{i}x_{i}$
to represent the uncertainty occurred in the market, which is proportional to
the production level $x_{i}$. The agents have their own risk levels
$\alpha_{i}$ and they aim to minimize the risk of incurring high costs, i.e.,
the CVaR of their cost functions.

Note that the regret for agent $i$ depends on the sequence of the other
agent’s actions $\\{x_{-i,t}\\}_{t=1}^{T}$, and this sequence depends on the
algorithm. Therefore, it is not appropriate to compare different algorithms in
terms of the regret. Instead, we use the empirical performance to evaluate the
algorithms. Specifically, we compute the CVaR values at each episode, and
observe how fast the agents can minimize the CVaR values during the learning
process.

To compute the distribution function $\bar{F}_{i,t}$, we divide the interval
$[0,U]$ into $n$ bins of equal width, and we approximate the expectation by
the sum of finite terms. We compare our zeroth-order momentum method with the
two algorithms in [14], which we term as the algorithm with sample reuse and
the algorithm with residual feedback, respectively. The number of samples
$n_{t}$ at each episode for Algorithm 1 is shown in Figure 1. We set the
momentum parameter $\beta=0.5$. Each algorithm is run for 20 trials and the
parameters of these algorithms are separately optimally tuned. The CVaR values
of these algorithms are presented in Figure 2. We observe that all the
algorithms finally converge to the same CVaR values, but the zeroth-order
momentum method proposed here converges faster. This is because our method
uses all past samples to get a better estimate of the CVaR values and thus
allows for a larger step size. Moreover, we also observe that the standard
deviation of both our proposed algorithm and the algorithm with residual
feedback is sufficiently small, which verifies the variance reduction effect
of using residual feedback.

Figure 1: The number of samples of Algorithm 1.

Figure 2: Comparative results among our proposed zeroth-order momentum method,
and the three algorithms in [14]. Shaded areas represent $\pm$ one standard
deviation over 20 runs.

Figure 3: CVaR values achieved by Algorithm 1 with different $\beta$ values.
Shaded areas represent $\pm$ one standard deviation over 20 runs.

Moreover, we explore the effect of different momentum parameters $\beta$. In
Figure 3, we present the achieved CVaR values of our method for various
$\beta$ values and all other parameters unchanged. Recalling the definition of
the distribution estimate in (3), large values of $\beta$ place more weight on
outdated cost values and thus reduce the convergence speed.

##  VI Conclusion

In this work, we proposed a zeroth-order momentum method for online convex
games with risk-averse agents. The use of momentum that employs past samples
allowed to improve the ability of the algorithm to estimate the CVaR values.
Zeroth-order estimation of the CVaR gradients using residual feedback allowed
us to reduce the variance of the CVaR gradient estimates. We showed that the
proposed algorithm theoretically achieves no-regret learning with high
probability, matching the best known result in literature. Moreover, we
provided numerical simulations that demonstrated the superior performance of
our method in practice. While sample reuse and residual feedback had been both
separately shown to improve the performance of online learning, here we showed
that their combination yields yet better performance.

## References

  * [1] S. Shalev-Shwartz and Y. Singer, “Convex repeated games and fenchel duality,” in _NIPS_ , vol. 6. Citeseer, 2006, pp. 1265–1272.
  * [2] G. J. Gordon, A. Greenwald, and C. Marks, “No-regret learning in convex games,” in _Proceedings of the 25th international conference on Machine learning_ , 2008, pp. 360–367.
  * [3] P. G. Sessa, I. Bogunovic, M. Kamgarpour, and A. Krause, “No-regret learning in unknown games with correlated payoffs,” _Advances in Neural Information Processing Systems_ , vol. 32, pp. 13 624–13 633, 2019.
  * [4] Y. Shi and B. Zhang, “No-regret learning in cournot games,” _arXiv preprint arXiv:1906.06612_ , 2019.
  * [5] E. Hazan, “Introduction to online convex optimization,” _arXiv preprint arXiv:1909.05207_ , 2019.
  * [6] S. Shalev-Shwartz _et al._ , “Online learning and online convex optimization,” _Foundations and trends in Machine Learning_ , vol. 4, no. 2, pp. 107–194, 2011.
  * [7] W. F. Sharpe, “The sharpe ratio,” _Journal of portfolio management_ , vol. 21, no. 1, pp. 49–58, 1994.
  * [8] P. Artzner, F. Delbaen, J.-M. Eber, and D. Heath, “Coherent measures of risk,” _Mathematical finance_ , vol. 9, no. 3, pp. 203–228, 1999.
  * [9] R. T. Rockafellar, S. Uryasev _et al._ , “Optimization of conditional value-at-risk,” _Journal of risk_ , vol. 2, pp. 21–42, 2000.
  * [10] T. Soma and Y. Yoshida, “Statistical learning with conditional value at risk,” _arXiv preprint arXiv:2002.05826_ , 2020.
  * [11] A. R. Cardoso and H. Xu, “Risk-averse stochastic convex bandit,” in _The 22nd International Conference on Artificial Intelligence and Statistics_. PMLR, 2019, pp. 39–47.
  * [12] S. Curi, K. Levy, S. Jegelka, A. Krause _et al._ , “Adaptive sampling for stochastic risk-averse learning,” _arXiv preprint arXiv:1910.12511_ , 2019\.
  * [13] A. Tamkin, R. Keramati, C. Dann, and E. Brunskill, “Distributionally-aware exploration for cvar bandits,” in _NeurIPS 2019 Workshop on Safety and Robustness on Decision Making_ , 2019.
  * [14] Z. Wang, Y. Shen, and M. Zavlanos, “Risk-averse no-regret learning in online convex games,” in _International Conference on Machine Learning_. PMLR, 2022, pp. 22 999–23 017.
  * [15] N. Qian, “On the momentum term in gradient descent learning algorithms,” _Neural networks_ , vol. 12, no. 1, pp. 145–151, 1999.
  * [16] Y. Zhang, Y. Zhou, K. Ji, and M. M. Zavlanos, “Boosting one-point derivative-free online optimization via residual feedback,” _arXiv preprint arXiv:2010.07378_ , 2020.
  * [17] M. Bravo, D. Leslie, and P. Mertikopoulos, “Bandit learning in concave n-person games,” _Advances in Neural Information Processing Systems_ , vol. 31, 2018.