tcb@breakable

# Second-Order Mirror Descent: Convergence in Games Beyond Averaging and
Discounting

Bolin Gao and Lacra Pavel This work was supported by a grant from NSERC and
Huawei Technologies Canada.B. Gao and L. Pavel are with the Department of
Electrical and Computer Engineering, University of Toronto, Toronto, ON, M5S
3G4, Canada. Emails<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

In this paper, we propose a second-order extension of the continuous-time
game-theoretic mirror descent (MD) dynamics, referred to as MD2, which
converges to mere (but not necessarily strict) variationally stable states
(VSS) without using common auxiliary techniques such as averaging or
discounting. We show that MD2 enjoys no-regret as well as an exponential rate
of convergence towards a strong VSS upon a slight modification. Furthermore,
MD2 can be used to derive many novel primal-space dynamics. Lastly, we use
stochastic approximation techniques to provide a convergence guarantee of
discrete-time MD2 with noisy observations towards interior mere VSS. Selected
simulations are provided to illustrate our results.

##  1 Introduction

A central problem of Nash equilibrium (NE) seeking in games is the design of
dynamics (in continuous-time) or algorithms (in discrete-time), which a set of
subscribing players can utilize to arrive at a desired collective outcome.
These methods often involve an information processing step, whereby game-
relevant information is made available to each player and then collected for
individual processing, followed by a decision step, whereby the gathered
information is converted back into a strategy for each player. A unified
mathematical codification of this two-step process that emerged in recent
decades is the mirror descent algorithm, first studied in the context of
optimization by [4]. Here, the information aggregation typically involves a
linear combination of all the past gradients, followed by a Bregman-type
projection that sends this combination into the next iterate while obeying the
constraints of the problem. Mirror descent has since been investigated under a
variety of perspectives in optimization, see [7, 8, 6, 10, 9, 11, 5, 12],
among others.

In games, the most popular manifestation of the mirror descent algorithm is
the continuous-time mirror descent dynamics (MD), which models after [5] in a
single-player/optimization setup. MD is variously referred to as dual
averaging (DA) in continuous games [2, 20] and Follow-the-Regularized-Leader
(FTRL) or exponential learning in finite games [21, 19, 22]. The core appeal
of MD lies in its flexibility in adapting to a wide variety of setups as well
as its rich theoretical properties. In particular, MD is known to converge in
terms of its generated strategies under a trio of strict assumptions:

  * (i)

the game is strictly monotone, i.e., the pseudo-gradient of the game is a
strictly monotone operator[51]. The class of strictly monotone games captures
games that admit strictly concave potential functions, as well as saddle-point
problems with strictly convex, strictly concave saddle functions, a special
case of diagonally strict concavity of [52].

  * (ii)

the NE is strictly variationally stable (in the sense of [2]). The class of
games with a strictly stable NE (or strictly variationally stable games [18])
contains the class of strictly monotone games, as well as strictly coherent
saddle point problems [47].

  * (iii)

the NE is strict [15], which coincides with a locally strictly stable NE in
finite games [2].

Despite these theoretical guarantees, the applicability of MD and its variants
remains limited as many games do not conveniently exhibit these strict
properties. For example, MD does not converge in many games, most notably in
zero-sum (ZS) finite games that feature a unique interior NE [22]. [13] showed
that no two-player ZS finite game is strictly monotone. Furthermore,
strictness is also intimately related to the uniqueness of the solution. In
practice, many games exhibit a convex set of solutions as opposed to a unique
one. This means strictness imposes serious constraints on the problem
parameters or problem types. Moreover, many algorithms that are directly
derived from MD (through discretization or otherwise) share similar
limitations in that some type of strictness needs to be assumed to ensure
convergence, e.g., [2, 19, 3, 17, 18, 20, 21, 23].

To achieve convergence beyond strictness in a continuous-time setup, the
existing literature typically relies on two primary approaches: averaging and
discounting. The first method utilizes the time-averaged (ergodic) strategies
generated by MD as opposed to the actual strategies [19]. However, this
approach has two drawbacks. Firstly, the players need to combine their time-
averaged strategies in order to recover the NE, which is unrealistic in non-
cooperative settings. Secondly, time-averaging could fail outside of ZS games
[58]. Another method is via a discounting procedure [28, 29, 16], which is
conceptually related to the weight decay method of [42]. Discounting can be
seen as a regularization technique that makes use of a strictly convex
function to offset poor game properties, such as the lack of strict
monotonicity. However, in general, discounting cannot yield exact convergence
[16, 28, 29]. Since averaging and discounting represent two of the most widely
studied method for improving the convergence properties of continuous-time MD,
yet at the same time both have been met with challenges, therefore it is
natural for us to set our sights on alternative approaches.

Contributions: In this work, in a game-theoretic setup, we propose a second-
order variant of the continuous-time MD (the first-order of which was studied
in [5, 6, 20, 19, 22, 21]), which we refer to as MD2. Second-order means that
the set of ODEs is second-order in time, thus MD2 can be seen as a dual-space
formulation of the heavy-ball method [33] and is closely related to the class
of $n$th order discounted game dynamics of [28], adapted here in a more
general game setting. Like MD, we show that MD2 satisfies the basic assumption
of no-regret. However, unlike MD, MD2 converges without using averaging even
when the NE is not strictly variationally stable. Moreover, unlike the
discounted MD [29], which partially overcomes the convergence issue of MD at
the cost of inexact convergence, MD2 converges exactly. Exact convergence via
a primal-space, second-order pseudo-gradient-type dynamics has been achieved
in (merely) monotone games [27]. We provide a dual-space generalization of
this result, beyond the monotone game setting, into the case where the
equilibria are merely variationally stable. Furthermore, MD2 can exactly
recover the unconstrained, primal-space dynamics of [27] in the full-
information case. Finally, we use our continuous-time convergence results to
guide the proof of convergence of a discrete-time MD2. In gist, we show that
higher-order dynamics can converge towards Nash solutions with more general
stability properties.

Related Works: This paper incorporates several ideas from the variational
stability, higher-order dynamics, and stochastic (semi-bandit) Nash-seeking
literature. The concept of a variationally stable state (VSS) has roots in
evolutionary game theory, and it is motivated by and analogous to that of an
evolutionarily stable state (ESS) [38]. Namely, an ESS is a locally strict VSS
for games with one or more populations of agents, also known as population
games (similarly, a global ESS or GESS is a globally strict VSS) [13]. Many
evolutionary game dynamics such as the replicator and projection dynamics have
been shown to converge towards NEs of strictly stable games, which are GESS
[13, 3]. In contrast to ESS, VSS is defined for the more general class of
continuous games [2]. The literature on strict VSS seeking without structural
property such as strict monotonicity is relatively recent [2, 20, 17, 18, 21].
[20] showed that DA converges towards a globally strict VSS. [18, 17] studied
online gradient descent (discrete MD with projection map) where the feedback
is received asynchronously between agents.

In contrast to a strict VSS, the class of mere VSS is motivated by the more
general notion of a neurally stable state (NSS) and global NSS (GNSS) [39].
Comparatively speaking, the set of literature that deals with convergence
towards a mere VSS is fewer, especially in the absence of structural
assumptions such as the game being merely monotone (also known as stable game)
[13]. In a population game setup, [13] showed that the best response,
integrable excess payoff and impartial pairwise comparison dynamics converge
globally to the set of NEs in (merely) stable games, which are GNSS. With a
few exceptions such as [47], the current set of literature on the (non-
ergodic) convergence towards a mere VSS outside of a finite/population game
setup typically also assumes monotonicity of the game [25, 59, 27]. In
contrast to the above references, in this work we generalize the convergence
behavior of MD to non-strict NE, i.e,. a mere VSS without any structural
assumption, through the use of a higher-order augmentation.

Higher-order dynamics was pioneered by Polyak in [33], whose algorithm later
became widely known as the heavy-ball method (HB) [33]. HB was recognized as
an analogue of a second-order ODE upon its inception and has since been
heavily investigated in the optimization literature [40] and arguably forms
the backbone of deep learning, where non-convexity and local minima dominate
[41]. Inspired by [33], the authors of [35] studied a second-order dynamics
for continuous concave games, which can be seen as a second-order extension of
the projection dynamics of [37] and showed that their dynamics converge
whenever there exists a Lipschitz and bounded potential function. In a finite
game setup, [34] studied $n$th-order variant of exponential learning [19],
whereby the payoff is processed via successive integration, and showed that
these higher-order learning schemes can achieve faster rate of elimination of
iteratively dominated solutions and convergence towards strict NE. An
important connection between higher-order dynamics and the stability of a NE
was established by [36], who showed that the addition of an anticipatory
process to first-order gradient-play and replicator dynamics can result in
local convergence towards an interior (non-strict) NE, despite the dynamics
themselves being uncoupled. For continuous-time equilibrium seeking via a
higher-order augmentation in the dual-space, the closest work related to ours
is the second-order discounted exponential learning scheme in [28], which
overcomes the problem with non-convergence towards interior NE of [34] at the
cost of inexact convergence. It was found empirically in [28] that such
higher-order augmentation improves the convergence property of its first-order
counterpart and converges in non-strictly monotone games where the first-order
fails. We improve [28] by generalizing these results from finite games into a
continuous-game setup and go beyond a monotone game setup while simultaneously
achieving exact convergence. Exact convergence (also known as last iterate
convergence) via a primal-space, second-order pseudo-gradient dynamics was
achieved in merely monotone games [27], which we also generalize, by stripping
away monotonicity assumptions and showing that MD2 encapsulates the dynamics
of [27] as a special case.

In the discrete-time, semi-bandit setup, whereby each player receives a
(possibly) noise-corrupted version of the pseudo-gradient, the closest work
related to ours is [2]. [2] studied the first-order MD, known as DA therein,
and showed that under imperfect gradient setup, DA converges to a globally
variationally stable NE under diminishing step-sizes and its ergodic average
converges to the set of NE in 2-player ZS games. Several variations to DA have
been studied in the semi-bandit literature, whereby their convergence is
towards either strict (strong) VSS or strict NE. For instance, in finite
games, [24] studied a Hedge-variant of exponential weight algorithm (HEDGE)
and provided convergence when the NE is strongly VS with respect to the
$L^{1}$ norm. [16] studied a similar setup for a discounted variant of HEDGE
for potential games. A similar analysis was also performed by [23] for FTRL in
finite games. [26] provided primal-space algorithms for merely monotone games.
[47] studied the optimistic MD algorithm, which converges in coherent games,
i.e., two-player games that possess mere VSS. Compared to optimistic MD,
discrete-time MD2 only requires one mirror projection instead of two prox
projections. Furthermore, optimistic MD cannot converge in the presence of
noise [47] unless the game is strictly coherent (a stronger assumption),
whereas the discrete-time MD2 can converge under noisy conditions.

Organization: We provide the background materials in Section II. Section III
reviews several basic properties of the first-order MD. We discuss the
convergence properties of MD2 in Section IV and provide the rate of
convergence as well as regret minimization guarantee in Section V. We derive
associated primal-space dynamics in Section VI. Section VII discusses MD2 in
discrete-time with noisy observations. Simulations are presented in Section
VIII. Section IX presents our conclusions. For readability, all proofs are
relegated to Section X. A table of main notations used in this paper (Table 1)
is provided in the Appendix.

##  2 Background

The following material for convexity and duality are drawn from [49, 63].
Those for game theory are drawn from [50, 51, 60].

###  2.1 Sets and Vectors

Given a convex set $\mathcal{C}\subseteq\mathbb{R}^{n}$, the (relative)
interior of the set is denoted as ($\operatorname{rint}(\mathcal{C})$)
$\operatorname{int}(\mathcal{C})$. $\operatorname{rint}(\mathcal{C})$
coincides with $\operatorname{int}(\mathcal{C})$ whenever
$\operatorname{int}(\mathcal{C})$ is non-empty. The closure of a set
$\mathcal{C}$ is denoted as $\operatorname{cl}(\mathcal{C})$. We denote the
non-negative orthant of $\mathbb{R}^{n}$ as $\mathbb{R}^{n}_{\geq 0}$, and the
positive orthant as $\mathbb{R}^{n}_{>0}$. A column vector in $\mathbb{R}^{n}$
is denoted as
$x=(x_{1},\ldots,x_{n})=\begin{bmatrix}x_{1},\ldots,x_{n}\end{bmatrix}^{\top}$.
$\mathbf{1}$ and $\mathbf{0}$ denote the column vector all ones and all zeros.
$\mathpzc{I}_{n\times n}$ and $\mathpzc{O}_{n\times n}$ denote the $n\times n$
identity and zero matrices. Subscript is omitted when the dimensionality is
unambiguous. Suppose $n$ is a natural number, then
$[n]\coloneqq\\{1,\ldots,n\\}$.

###  2.2 Convex Functions and Duality

Let $\mathcal{M}=\mathbb{R}^{n}$ be endowed with norm $\|\cdot\|$ and dot
product $\langle\cdot,\cdot\rangle$. An extended real-valued function is a
function $f$ that maps from $\mathcal{M}$ to $[-\infty,\infty]$. The
(effective) domain of $f$ is
$\operatorname{dom}(f)=\\{x\in\mathcal{M}:f(x)<\infty\\}$. A function
$f:\mathcal{M}\to[-\infty,\infty]$ is proper if $f(x)\neq-\infty,\forall x$
and there exists at least one $x\in\mathcal{M}$ such that $f(x)<\infty$; $f$
is closed if its epigraph is closed. A function
$f:\mathcal{M}\to[-\infty,\infty]$ is supercoercive if
$\lim_{\|x\|\to\infty}f(x)/\|x\|\to\infty$. The indicator function over
$\mathcal{C}$ is denoted by $\delta_{\mathcal{C}}$. The normal cone of
$\mathcal{C}$ is defined as,
$N_{\mathcal{C}}(x)=\\{v\in\mathbb{R}^{n}|v^{\top}(y-x)\leq 0,\forall
y\in\mathcal{C}\\}$. The tangent cone of a non-empty convex set $\mathcal{C}$
at $x\in\mathcal{C}$ is
$T_{\mathcal{C}}(x)=\overline{\bigcup_{\lambda\in\mathbb{R}_{>0}}\lambda(\mathcal{C}-x)}$
and equals $\emptyset$ for all $x\notin\mathcal{C}$. $Id$ denotes the identity
function.

Recall that
$\pi_{\mathcal{C}}(x)={\text{argmin}}_{y\in\mathcal{C}}\|y-x\|_{2}^{2}$ is the
Euclidean projection of $x$ onto $\mathcal{C}$, where we refer to
$\pi_{\mathcal{C}}$ as the Euclidean projector. Let $\partial f(x)$ denote a
the set of subgradients of $f$ at $x$ and $\nabla f(x)$ the gradient of $f$ at
$x$. For differentiation on the boundary of a closed set $\mathcal{C}$, in
lieu of the subgradient, we can also assume $f$ is defined and differentiable
on an open set containing $\mathcal{C}$. Given $f$, the function
$f^{\star}\\!:\\!\mathcal{M}^{\star}\\!\to\\![-\infty,\infty]$ defined by
$f^{\star}(z)\\!=\\!\sup_{x\in\mathcal{M}}\\!\big{[}x^{\top}z\\!-\\!f(x)\big{]}$,
is called the conjugate function of $f$, where $\mathcal{M}^{\star}\\!$ is the
dual space of $\mathcal{M}$, endowed with the dual norm $\|\cdot\|_{\star}$.
$f^{\star}$ is closed and convex if $f$ is proper. By the conjugate
subgradient theorem [49], suppose $f:\mathbb{R}^{n}\to(-\infty,+\infty]$ is
proper, closed and convex and $f^{\star}$ is its Fenchel conjugate, then for
any $x,z\in\mathbb{R}^{n}$, $x^{\top}z=f(x)+f^{\star}(z)\Longleftrightarrow
z\in\partial f(x)\Longleftrightarrow x\in\partial f^{\star}(z).$ The Bregman
divergence of a proper, closed, convex function $f$ is
$D_{f}:\operatorname{dom}(f)\times\operatorname{dom}(\partial
f)\to\mathbb{R},D_{f}(x,y)=f(x)-f(y)-g^{\top}(x-y),g\in\partial f(y)$. Given a
vector-valued function $F$, the Jacobian of $F$ is denoted as
$\mathbf{J}_{F}$.

###  2.3 N-Player Concave Games

Let
$\mathcal{G}=(\mathcal{N},\\{\Omega^{p}\\}_{p\in\mathcal{N}},\\{\mathcal{U}^{p}\\}_{p\in\mathcal{N}})$
be a game, where $\mathcal{N}=\\{1,\ldots,N\\}$ is the set of players,
$\Omega^{p}\subseteq\mathbb{R}^{n_{p}}$ is the set of player $p$’s strategies
(actions). We denote the strategy (action) set of player $p$’s opponents as
$\Omega^{-p}\subseteq\prod_{q\in\mathcal{N},q\neq p}\mathbb{R}^{n_{q}}$. We
denote the set of all the players strategies as
$\Omega=\prod_{p\in\mathcal{N}}\Omega^{p}\subseteq\prod_{p\in\mathcal{N}}\mathbb{R}^{n_{p}}=\mathbb{R}^{n},n=\sum_{p\in\mathcal{N}}n_{p}$.
We refer to $\mathcal{U}^{p}:\Omega\to\mathbb{R},x\mapsto\mathcal{U}^{p}(x)$
as player $p$’s real-valued payoff function, where
$x=(x^{p})_{p\in\mathcal{N}}\in\Omega$ is the action profile of all players,
and $x^{p}\in\Omega^{p}$ is the action of player $p$. We also denote $x$ as
$x=(x^{p};x^{-p})$ where $x^{-p}\in\Omega^{-p}$ is the action profile of all
players except $p$. For differentiability purposes, we make the implicit
assumption that there exists some open set, on which $\mathcal{U}^{p}$ is
defined and continuously differentiable, such that it contains $\Omega$.

######  Assumption 1.

For all $p\in\mathcal{N}$, $\Omega^{p}$ is a non-empty, closed, convex, subset
of $\mathbb{R}^{n_{p}}$. $\mathcal{U}^{p}(x^{p};x^{-p})$ is (jointly)
continuous in $x=(x^{p};x^{-p})$. $\mathcal{U}^{p}(x^{p};x^{-p})$ is concave
and continuously differentiable ($\mathcal{C}^{1}$) in $x^{p}$ for all
$x^{-p}\in\Omega^{-p}$.

Under 1, we refer to $\mathcal{G}$ as a (continuous) concave game. Given
$x^{-p}\in\Omega^{p}$, each agent $p\in\mathcal{N}$ aims to find the solution
of the following optimization problem,

$\displaystyle\underset{x^{p}}{\text{maximize}}$
$\displaystyle\mathcal{U}^{p}(x^{p};x^{-p})$ subject to $\displaystyle
x^{p}\in\Omega^{p}.$ (1)

A profile ${{x}}^{\star}=({{x}^{p}}^{\star})_{p\in\mathcal{N}}$ is a Nash
equilibrium (NE) if,

$\mathcal{U}^{p}({x^{p}}^{\star};{x^{-p}}^{\star})\geq\mathcal{U}^{p}(x^{p};{x^{-p}}^{\star}),\forall
x^{p}\in\Omega^{p},\forall p\in\mathcal{N}.$ (2)

At a NE, no player can increase its payoff by unilateral deviation.

######  Remark 1.

Under Assumption 1, existence of a NE is guaranteed for bounded $\Omega^{p}$
[50, Theorem 4.4]. When $\Omega^{p}$ is closed but not bounded, the existence
of a NE is guaranteed under the additional assumption that $-\mathcal{U}^{p}$
is coercive in $x^{p}$, that is,
$\textstyle\lim_{\|x^{p}\|\to\infty}-\mathcal{U}^{p}(x^{p};x^{-p})=+\infty,$
for all $x^{-p}\in\Omega^{-p},p\in\mathcal{N}$ [50, Corollary 4.2]. A NE is
said to be strict when the inequality in (2) is strict [15].

A useful characterization of a NE of a concave game $\mathcal{G}$ is given in
terms of the pseudo-gradient [52], which is defined as,

$U:\Omega\to\mathbb{R}^{n},U(x)=(U^{p}(x))_{p\in\mathcal{N}},$ (3)

where $U^{p}(x)=\nabla_{x^{p}}\mathcal{U}^{p}(x^{p};x^{-p})$ is the partial-
gradient. By [51, Proposition 1.4.2], $x^{\star}\in\Omega$ is a NE iff,

$(x-x^{\star})^{\top}U(x^{\star})\leq 0,\forall x\in\Omega.$ (4)

Equivalently, $x^{\star}$ is a solution of the Stampacchia variational
inequality $\text{VI}(\Omega,-U)$, [51]. Following [3], we say that a NE is
globally strict if the inequality of (4) is held strictly for all
$x^{\star}\neq x$.

###  2.4 Monotonicity

A general class of games in which many dynamics are guaranteed to converge is
the class of monotone games, also known as stable games [13, 14] or
dissipative games [45]. We contrast some known definitions associated with
monotone games in the literature.

######  Definition 1.

The game $\mathcal{G}$ is,

  * (i)

$\eta$-strongly monotone if
$\\!(U(x)\\!-\\!U(x^{\prime}))\\!^{\top}\\!(x\\!-\\!x^{\prime})\\!\leq\\!-\eta\|x\\!-\\!x^{\prime}\|^{2}_{2}$,
$\\!\forall x,x^{\prime}\\!\in\\!\Omega$, for some $\eta\\!>\\!0$.

  * (ii)

strictly monotone if
$(U(x)\\!-\\!U(x^{\prime}))^{\\!\top}\\!(x\\!-\\!x^{\prime})\\!<\\!0,\forall
x\in\Omega\backslash\\{x^{\prime}\\}$, with equality if and only if
$x=x^{\prime}$.

  * (iii)

(merely) monotone if $(U(x)-U(x^{\prime}))^{\top}(x-x^{\prime})\leq 0,\forall
x,x^{\prime}\in\Omega.$

  * (iv)

$\mu$-weakly monotone if
$\\!(U(x)\\!-\\!U(x^{\prime}))\\!^{\top}\\!(x\\!-\\!x^{\prime})\\!\leq\\!\mu\|x\\!-\\!x^{\prime}\|^{2}_{2}$,
$\\!\forall x,x^{\prime}\\!\in\\!\Omega$, for some $\mu\\!>\\!0$.

Strongly and strictly monotone games have been extensively investigated in the
literature, at least dating from [52]. Merely monotone games have been studied
in [25, 26, 59, 29, 28, 27, 13]. Weakly monotone games were considered in [28,
44, 31]. When $U$ is $\mathcal{C}^{1}$, there exists a natural
characterization in terms of definiteness of its Jacobian [51, p. 155, Prop.
2.3.2]. Note that strictly monotone games can have at most one NE, whereas NEs
can form a non-singleton convex set in merely monotone games [51].

###  2.5 Variational Stablility

Although many classical examples of games satisfy monotonicity properties
[13], these conditions may not hold in more complex scenarios. A recent line
of research has started to relax monotonicity notions from a game to that of
an equilibrium via the notion of a variationally stable (VS) equilibrium or
state [2].

######  Definition 2.

A Nash equilibrium $x^{\star}\in\Omega$ is,

  * (i)

$\eta$-strongly VS if
$U(x)^{\top}(x-x^{\star})\leq-\eta\|x-x^{\star}\|_{2}^{2},\forall x\in\Omega,$
for some $\eta\\!>\\!0$.

  * (ii)

strictly VS if $U(x)^{\top}(x-x^{\star})\leq 0,\forall x\in\Omega$, with
equality if and only if $x=x^{\star}$.

  * (iii)

(merely) VS if $U(x)^{\top}(x-x^{\star})\leq 0,\forall x\in\Omega.$

  * (iv)

$\mu$-weakly VS if
$U(x)^{\top}(x-x^{\star})\leq\mu\|x-x^{\star}\|_{2}^{2},\forall x\in\Omega,$
for some $\mu\\!>\\!0$.

If a condition (i) - (iv) hold on $\mathcal{D}\subset\Omega$, then the
associated definition is said to hold locally. Otherwise the definition is
said to hold globally.

######  Remark 2.

We refer to $x^{\star}$ as a global (strong/strict/mere/weak) variationally
stable state (VSS) if it satisfies one of the corresponding VS notions in 2 on
all of $\Omega$; $x^{\star}$ is a local (strong/strict/mere/weak) VSS
otherwise.

Globally mere VSS are the solutions to Minty variational inequality [51]. In a
population game context, the set of globally mere VSS is called GNSS and (the
unique) globally strict VSS is called GESS [3, 13, 14]. Strict VSS was
extended to a local/global set-wise definition in [2]. [17] introduced a
slight variation of strict VSS called $\lambda$-VS. [48] studied a local
version of strict VSS under the name locally asymptotically stable
differential NE, which requires twice continuous-differentiability of the
payoff functions. Strong VSS (2(i)) was also studied in [2]. [24] studied a
variant of the strong VSS with $L^{1}$ norm.

It can be seen that any NE of the strongly/strictly/merely/weak monotone game
is a strong/strict/mere/weak VSS. In particular, the non-empty set of NE
coincides with the set of mere VSS for merely monotone concave games [3,
Theorem 2]. When $\operatorname{int}(\Omega)\neq\emptyset$ and $x^{\star}$ is
an interior NE, i.e., $U(x^{\star})=\mathbf{0}$, the condition for VS recovers
the condition associated with monotone at $x^{\prime}=x^{\star}$. Hence VS can
be seen as a type of point-wise monotonicity and it is known to have a second-
order characterization similar to that for monotone games, but specifically at
the NE, see [2]. Throughout this work, we make the following assumption.

######  Assumption 2.

$\mathcal{G}$ admits an interior mere VSS.

###  2.6 Second-Order Characterization of VSS

In practice, outside of monotone games, it is often difficult to verify that a
NE is a particular type of VSS directly through 2. A more common approach to
characterize the VSS conditions of a solution is by looking at the symmetric
game Jacobian whenever the pseudo-gradient $U$ is $\mathcal{C}^{1}$, which is
given by,

$\operatorname{\boldsymbol{\mathfrak{J}}}_{U}(x)\coloneqq\frac{1}{2}(\operatorname{\mathbf{J}}_{U}(x)+\operatorname{\mathbf{J}}^{\top}_{U}(x)),x\in\Omega.$
(5)

The matrix
$\operatorname{\boldsymbol{\mathfrak{J}}}_{U}(x)\in\mathbb{R}^{n\times n}$ is
symmetric and thus guaranteed to have real eigenvalues, which makes it
amenable to analysis.

We now provide several sufficient conditions for verifying whether a VSS
$x^{\star}$ is strict, mere or weak. This is performed by checking the
definiteness of $\operatorname{\boldsymbol{\mathfrak{J}}}_{U}(x)$ or
$\operatorname{\boldsymbol{\mathfrak{J}}}_{U}(x^{\star})$. Our condition for
strict VSS is the same as those in [2]. Before proceeding, we say that the
(symmetric) game Jacobian
$\operatorname{\boldsymbol{\mathfrak{J}}}_{U}(x)\in\mathbb{R}^{n\times n}$ (5)
is negative definite on $T_{\Omega}(x)$ if,

$y^{\top}\operatorname{\boldsymbol{\mathfrak{J}}}_{U}(x)y<0,\forall
x\in\Omega,\forall y\in T_{\Omega}(x),y\neq\mathbf{0},$ (6)

and negative semi-definite on $T_{\Omega}(x)$ if the preceding inequality is
non-strict. We use the shorthand notations
$\operatorname{\boldsymbol{\mathfrak{J}}}_{U}(x)\prec\mathpzc{O}$ and
$\operatorname{\boldsymbol{\mathfrak{J}}}_{U}(x)\preceq\mathpzc{O}$
respectively. For $\mu>0$,
$\operatorname{\boldsymbol{\mathfrak{J}}}_{U}(x)-\mu\mathpzc{I}\prec\mathpzc{O}$
is written as
$\operatorname{\boldsymbol{\mathfrak{J}}}_{U}(x)\preceq\mu\mathpzc{I}$.
Similar conventions for when $x=x^{\star}$.

######  Proposition 1.

Let $x^{\star}\in\Omega$ be a NE of $\mathcal{G}$. Suppose $U$ is continuously
differentiable, and,

  1. (i)

$\operatorname{\boldsymbol{\mathfrak{J}}}_{U}(x)\prec\mathpzc{O}$ on
$T_{\Omega}(x),\forall x\in\Omega$, then $x^{\star}$ is globally strictly VS
and isolated.

  2. (ii)

$\operatorname{\boldsymbol{\mathfrak{J}}}_{U}(x)\preceq\mathpzc{O}$ on
$T_{\Omega}(x),\forall x\in\Omega$, then $x^{\star}$ is globally merely VS.

  3. (iii)

$\operatorname{\boldsymbol{\mathfrak{J}}}_{U}(x)\preceq\mu\mathpzc{I}$ on
$T_{\Omega}(x),\forall x\in\Omega$, then $x^{\star}$ is globally $\mu$-weakly
VS.

Suppose instead,

  1. (i’)

$\operatorname{\boldsymbol{\mathfrak{J}}}_{U}(x^{\star})\prec\mathpzc{O}$ on
$T_{\Omega}(x^{\star})$, then $x^{\star}$ is locally strictly VS and isolated.

  2. (ii’)

$\operatorname{\boldsymbol{\mathfrak{J}}}_{U}(x^{\star})\preceq\mathpzc{O}$ on
$T_{\Omega}(x^{\star})$, then $x^{\star}$ is locally merely VS.

  3. (iii’)

$\operatorname{\boldsymbol{\mathfrak{J}}}_{U}(x^{\star})\preceq\mu\mathpzc{I}$
on $T_{\Omega}(x^{\star})$, then $x^{\star}$ is locally $\mu$-weakly VS.

######  Remark 3.

These conditions can be verified by calculating the maximum eigenvalue of
$\operatorname{\boldsymbol{\mathfrak{J}}}_{U}(x)$ (resp.
$\operatorname{\boldsymbol{\mathfrak{J}}}_{U}(x^{\star})$), where
$\lambda_{\text{max}}(\operatorname{\boldsymbol{\mathfrak{J}}}_{U}(x))$ (resp.
$\lambda_{\text{max}}(\operatorname{\boldsymbol{\mathfrak{J}}}_{U}(x^{\star}))$)
denotes the (real) eigenvalue with the largest magnitude associated with an
eigenvector in $T_{\Omega}(x)$ (resp. $T_{\Omega}(x^{\star})$). When
$\mathbf{J}_{U}(x)$ (resp. $\mathbf{J}_{U}(x^{\star})$) is symmetric, then
analysis can be performed directly without resorting to calculating
$\operatorname{\boldsymbol{\mathfrak{J}}}_{U}$.

Note that, for instance, even if $x^{\star}$ is shown to be $\mu$-weak through
1(iii) or (iii’), it does not preclude the possibility that $x^{\star}$ could
in fact be strict or mere due to the looseness of the bound in the proof of 1.
Furthermore, unlike strict VSS, mere VSS need not be isolated. The next
examples illustrate various notions of variational stability.

######  Example 1.

(Monotone game with a mere VSS) Every NE of a merely monotone game is globally
mere VSS. Perhaps the simplest example of a merely monotone game is the so-
called bilinear saddle point problem, whereby, suppose we have a saddle
function, $f(x^{1},x^{2})=x^{1}x^{2}$, $x^{p}\in\mathbb{R}$ for which we want
to minimize in $x^{1}$ and maximize in $x^{2}$. We can cast this saddle point
problem as a game by utilizing two payoff functions
$\mathcal{U}^{1}(x^{1};x^{2})=-x^{1}x^{2}$ and
$\mathcal{U}^{2}(x^{2};x^{1})=x^{1}x^{2}$. This game has a pseudo-gradient of
$U(x)=(-x^{2},x^{1})$ and a game Jacobian
$\operatorname{\boldsymbol{\mathfrak{J}}}_{U}(x)=\mathpzc{O}$. By 1(ii), the
unique interior NE $x^{\star}=(0,0)$ is the globally mere VSS.

######  Example 2.

(Non-monotone potential game with a mere VSS) This example shows that the mere
monotonicity enjoyed by 1 can be destroyed by the addition of another player,
even when all the player’s payoff functions remain linear in its own argument.
Considered a three-player game, whereby,

$\textstyle\mathcal{U}^{p}(x^{p};x^{-p})=-x^{p}\prod_{q\neq p}x^{q},\quad
p\in\\{1,2,3\\},$ (7)

and each $x^{p}\in[-1,1]$, i.e., $\Omega=[-1,1]^{3}$. For this game, the
pseudo-gradient is $U(x)=-(x^{2}x^{3},x^{1}x^{3},x^{1}x^{2})$ which shows that
there exists an interior NE at $x^{\star}=(0,0,0)$. The Jacobian, which is
symmetric, is,

$\operatorname{\mathbf{J}}_{U}(x)=-\begin{bmatrix}0&x^{3}&x^{2}\\\
x^{3}&0&x^{1}\\\ x^{2}&x^{1}&0\end{bmatrix},$

hence this game is non-monotone in general and in fact possesses a non-convex
and non-concave potential function $P(x)=-x^{1}x^{2}x^{3}$, i.e., $\nabla
P=U$. Since $\operatorname{\mathbf{J}}_{U}(x^{\star})=\mathpzc{O}$ is negative
semi-definite, hence by 1(ii’), $x^{\star}$ is locally merely VS. Note that if
we restrict the strategy set to $\Omega=[0,1]^{3}$ instead, every NE
$x^{\star}\in\\{x\in\Omega|x^{p}=x^{q}=0,p,q\in\mathcal{N},p\neq q\\}$ is a
globally mere VSS, as $U(x)^{\top}(x-x^{\star})=-3x^{1}x^{2}x^{3}\leq
0,\forall x\in[0,1]^{3}$.

######  Example 3.

(RPS with non-negative payoff for ties) Consider a two-player Rock-Paper-
Scissors (RPS) game with $\mathpzc{A}$ and $\mathpzc{B}$ being the payoff
matrices for player $1$ and $2$ respectively,

$\mathpzc{A}=\begin{bmatrix}\varsigma&-\mathpzc{l}&\mathpzc{w}\\\
\mathpzc{w}&\varsigma&-\mathpzc{l}\\\
-\mathpzc{l}&\mathpzc{w}&\varsigma\end{bmatrix},\mathpzc{B}=\mathpzc{A}^{\top},$
(8)

where $\mathpzc{l},\mathpzc{w}\geq 0$ are the values associated with a loss or
a win and $\varsigma\in\mathbb{R}$ is the payoff of a tie. The strategy set
associated with this game is the simplex
$\Omega^{p}=\\{x^{p}\in\mathbb{R}^{n_{p}}|\sum_{i=1}^{n_{p}}x^{p}_{i}=1,x^{p}_{i}\geq
0\\}$. The pseudo-gradient (as known as payoff vector [13, 14]) and the game
Jacobian are as follows,

$U(x)=\begin{bmatrix}\mathpzc{O}&\mathpzc{A}\\\
\mathpzc{A}&\mathpzc{O}\end{bmatrix}\begin{bmatrix}x^{1}\\\
x^{2}\end{bmatrix}\quad\operatorname{\boldsymbol{\mathfrak{J}}}_{U}(x)=\dfrac{1}{2}\begin{bmatrix}\mathpzc{O}&\mathpzc{A}+\mathpzc{A}^{\top}\\\
\mathpzc{A}+\mathpzc{A}^{\top}&\mathpzc{O}\end{bmatrix}\\!.$ (9)

To simplify our analysis, consider an example where $l=0$. In this game,
$x^{\star}=({x^{p}}^{\star})_{p\in\mathcal{N}},{x^{p}}^{\star}=(1/3,1/3,1/3)$
is the unique interior NE for any $\varsigma\in[0,w)$, $\mathpzc{w}>0$. The
eigenvalues of $\operatorname{\boldsymbol{\mathfrak{J}}}_{U}(x)$ on
$T_{\Omega}(x)$ are
$\\{-(\varsigma+\mathpzc{w}),\pm(\varsigma-\mathpzc{w}/2),\pm(\varsigma-\mathpzc{w}/2)\\}$,
hence this game is $\mu$-weakly merely monotone with
$\mu=\varsigma-\mathpzc{w}/2$ whenever
$\varsigma\in(\mathpzc{w}/2,\mathpzc{w})$. When $\varsigma=\mathpzc{w}/2$, the
game is merely monotone (1). Taken together, this means for
$\varsigma\in(\mathpzc{w}/2,\mathpzc{w})$, $x^{\star}$ is a unique, locally
$\mu$-weak VSS where $\mu=\varsigma-\mathpzc{w}/2$, whereas for
$\varsigma=\mathpzc{w}/2$, $x^{\star}$ is a unique, globally mere VSS.
Finally, note that when $\mathpzc{l}=\mathpzc{w}=\varsigma=0$ (and hence both
$\mathpzc{A},\mathpzc{B}$ are zero matrices), every strategy is globally
merely VS and forms a convex set (the strategy set itself).

##  3 Review of First-Order MD Dynamics

We now describe the dynamic process that a group of agents utilizes in order
to arrive at one of the equilibrium notions discussed in previous sections.
One such general model is the mirror descent (MD) dynamics. Intuitively, the
family of MD dynamics ascribes an abstract behavior model to each player,
which states that, upon receiving the partial-gradient of its payoff function,
each player processes the partial-gradient information (typically via an
aggregation), then converts the processed information into the next strategy.
This process can be described by the following set of ODEs [20, 19, 22, 21],

$\vspace{-0.1cm}\dot{z}^{p}=\gamma\nabla_{x^{p}}\mathcal{U}^{p}(x^{p};x^{-p})=\gamma
U^{p}(x),\,\,\,x^{p}=C^{p}_{\epsilon}(z^{p}),$ (MD)

where $z^{p}$ is referred to as score or dual variable, $\gamma>0$ is a rate
parameter, and $C^{p}_{\epsilon}:\mathbb{R}^{n_{p}}\to\Omega^{p}$ is referred
to as the mirror map,

$C^{p}_{\epsilon}(z^{p})=\text{argmax}_{y^{p}\in\Omega^{p}}\left[{y^{p}}^{\top}z^{p}-\epsilon\vartheta^{p}(y^{p})\right],\epsilon>0.$
(10)

Here, $\vartheta^{p}:\mathbb{R}^{n_{p}}\to\mathbb{R}\cup\\{\infty\\}$ is
assumed to be a closed, proper and (at least) strictly convex function,
referred to as a regularizer, where $\operatorname{dom}(\vartheta^{p})$ is
assumed be a non-empty, closed and convex set, which agrees with the strategy
set $\Omega^{p}$, $\epsilon>0$ is referred to as the regularization constant.
The regularizer often satisfies one of the following distinct assumptions that
ensure the existence of a unique solution of MD and (10):

######  Assumption 3.

$\vartheta^{p}\\!:\\!\mathbb{R}^{n_{p}}\\!\\!\to\\!\mathbb{R}\\!\cup\\!\\{\\!\infty\\!\\}$
is closed, proper, convex, with $\operatorname{dom}(\vartheta^{p})=\Omega^{p}$
non-empty, closed and convex. In addition, $\vartheta^{p}$ is, (i) Legendre
(i.e., strictly convex, steep,
$\operatorname{int}(\operatorname{dom}(\vartheta^{p}))\\!\neq\\!\emptyset$)
and supercoercive, or (ii) $\rho$-strongly convex, $\rho>0$.

We note that for $\vartheta^{p}$ to be steep, it means that
$\|\nabla\vartheta^{p}(x_{k}^{p})\|_{\star}\to+\infty$, $\|\cdot\|_{\star}$ is
the dual norm, whenever $\\{x_{k}^{p}\\}_{k=1}^{\infty}$ is a sequence in
$\operatorname{rint}(\operatorname{dom}(\vartheta^{p}))=\operatorname{rint}(\Omega^{p})$
converging to a point in the (relative) boundary. Let
${\psi^{p}_{\epsilon}}^{\star}$ be the convex conjugate of
${\psi^{p}_{\epsilon}}$. Then by 3 and 2 (in the Appendix), under 3,
$C^{p}_{\epsilon}=\nabla{\psi^{p}_{\epsilon}}^{\star}$. When $\vartheta^{p}$
is steep, $C^{p}_{\epsilon}$ maps from $\mathbb{R}^{n_{p}}$ to all values in
$\Omega^{p}$ except those at the boundary; $C^{p}_{\epsilon}$ could map to
boundary points otherwise. We refer to $C^{p}_{\epsilon}$ as the mirror map
induced by $\psi^{p}_{\epsilon}\coloneqq\epsilon\vartheta^{p}\\!$ and Legendre
or strongly-induced (by $\psi^{p}_{\epsilon}$) when $\vartheta^{p}$ satisfies
Assumption 2(i) or (ii). Note that the notions of Legendre and strongly convex
need not be dichotomous. It is possible for a $C^{p}_{\epsilon}$ to be both
Legendre and strongly-induced. The following examples capture the three types
of mirror maps that are extensively used in the literature [5, 6, 16, 20, 19,
21, 22, 28, 29, 18].

######  Example 4.

Let $\Omega^{p}=\mathbb{R}^{n_{p}}$,
$\vartheta^{p}(x^{p})=\frac{1}{2\epsilon}\|x^{p}\|_{2}^{2}$, hence it
satisfies 3(i), (ii). In this case, $C^{p}_{\epsilon}=Id$ and it is both
Legendre and strongly-induced.

######  Example 5.

Let
$\Omega^{p}=\\{x^{p}\in\mathbb{R}^{n_{p}}|\sum_{i=1}^{n_{p}}x_{i}^{p}=1,x_{i}^{p}\geq
0\\}$ (the unit simplex),
$\vartheta^{p}(x^{p})=\sum_{i=1}^{n_{p}}{x^{p}_{i}}^{\top}\log(x_{i}^{p})$,
where we assume $0\log(0)=0$. $\vartheta^{p}$ can be shown to be $1$-strongly
convex in $\|\cdot\|_{2}$ (or $\|\cdot\|_{1}$), steep, hence it satisfies
3(ii). $\textstyle
C^{p}_{\epsilon}(z^{p})=(\exp(\epsilon^{-1}z^{p}_{i})(\sum_{j=1}^{n_{p}}\exp(\epsilon^{-1}z^{p}_{j}))^{-1})_{i\in[n_{p}]}$
is the softmax function.

######  Example 6.

Let $\Omega^{p}\subset\mathbb{R}^{n_{p}}$ be a non-empty, convex, compact set,
assume $\vartheta^{p}(x^{p})=\frac{1}{2}\|x^{p}\|_{2}^{2}$, hence
$\vartheta^{p}$ satisfies 3(ii) and is non-steep. Here,
$C_{\epsilon}^{p}=\pi_{\Omega^{p}}$ is the Euclidean projection onto
$\Omega^{p}$.

We can represent MD in a more compact stacked notation,

$\dot{z}=\gamma U(x),\quad x=C_{\epsilon}(z),$ (11)

where
$x=(x^{p})_{p\in\mathcal{N}},z=(z^{p})_{p\in\mathcal{N}},C_{\epsilon}=(C^{p}_{\epsilon})_{p\in\mathcal{N}},U=(U^{p})_{p\in\mathcal{N}}$.
Observe that at rest, MD (11) satisfies,

$\mathbf{0}=U(\overline{x}),\quad\overline{x}=C_{\epsilon}(\overline{z}),\quad\overline{z}\in
C_{\epsilon}^{-1}(\overline{x}),$ (12)

where
$\overline{x}=({\overline{x}^{p}})_{p\in\mathcal{N}},\overline{z}=({\overline{z}^{p}})_{p\in\mathcal{N}}$
are the rest points of MD and
$C_{\epsilon}^{-1}=((C^{p}_{\epsilon})^{-1})_{p\in\mathcal{N}}=(\partial\psi^{p}_{\epsilon})_{p\in\mathcal{N}}$
is the inverse (or the pre-image) of $C_{\epsilon}$. The rest condition (12)
implies that $\overline{x}$ is an interior NE. Hence if a trajectory $z(t)$ of
MD comes to a rest, $x(t)=C_{\epsilon}(z(t))$ converges to an interior NE. We
note that the uniqueness of $\overline{x}$ does not imply the uniqueness of
$\overline{z}$ unless $C_{\epsilon}$ is Legendre-induced (for which
$C_{\epsilon}$ is one-to-one; this follows from Legendre theorem [53]). Hence,
in general, the convergence of MD is to a set and not to an equilibrium and
$z(t)$ may continue to evolve even after $x(t)$ has reached an equilibrium.
Note that rest points are not the only game relevant solutions for which MD
may converge to. As pointed out in [19], there are non-rest points that occur
on the boundary of $\Omega$ that are asymptotically stable under MD. The
following lemma broadly summarizes the key convergence properties of MD (refer
to [20, 19] for proofs).

######  Lemma 1.

Let $\mathcal{G}$ be a concave game with a globally strict VSS $x^{\star}$.
Suppose that all players choose strategies according to MD (11). Let
$x(t)=(x^{p}(t))_{p\in\mathcal{N}}=C_{\epsilon}(z(t))$ be generated by (11)
and $C_{\epsilon}=(C_{\epsilon}^{p})_{p\in\mathcal{N}}$ be induced by
$\psi^{p}_{\epsilon}\coloneqq\epsilon\vartheta^{p}$. For any
$\gamma,\epsilon>0$ and any
$x(0)=C_{\epsilon}(z(0))\in\Omega,z(0)\in\mathbb{R}^{n}$,

  * (i)

suppose $\vartheta^{p}$ satisfies 3(i), $\forall p$, then
$x(t)=C_{\epsilon}(z(t))$ converges to $x^{\star}$, whenever $x^{\star}$ is
interior.

  * (ii)

suppose $\Omega^{p}$ is compact and $\vartheta^{p}$ satisfies 3(ii), $\forall
p$, then $x(t)=C_{\epsilon}(z(t))$ converges to $x^{\star}$.

In general, MD does not converge beyond strict VSS, i.e., a mere VSS. This
poses considerable limitations in practice, for instance, mere VSS are
commonly found in ZS games. [13, 14] showed that every ZS finite game is
merely (but not strictly) monotone, hence all of its NEs are mere (non-strict)
VSS. A standard method to overcome non-convergence to the NE in ZS games is to
calculate a time-averaged (ergodic) trajectory
$x_{\text{avg}}(t)=t^{-1}\int\nolimits_{0}^{t}x(\tau)\mathrm{d}\tau$ in tandem
with MD, that is,

$\dot{z}=\gamma U(x),\quad x=C_{\epsilon}(z),\quad
x_{\text{avg}}(t)=t^{-1}\int\nolimits_{0}^{t}x(\tau)\mathrm{d}\tau,$ (MDA)

for which the time-averaged strategy $x_{\text{avg}}$ has been shown in many
contexts to converge, e.g., [19]. The main critique of using MDA is that the
actual strategies do not arrive at the NE in the long-run, thereby making it
unsuitable for on-line equilibrium seeking in the absence of a central planner
or coordination between players. Furthermore, averaging may fail to converge
outside of ZS games [58], which makes this approach vulnerable to parameter
perturbation.

Another method for overcoming non-convergence is through discounting, which
was studied in [29],

$\dot{z}=\gamma(-z+U(x)),\quad x=C_{\epsilon}(z),$ (DMD)

where DMD stands for discounted mirror descent. Compared to MD, an extra $-z$
term is inserted in the $\dot{z}$ system, which translates into an exponential
weighted decay (or discounting) term in the closed-form solution $z(t)$. DMD
has a connection with the so-called weight decay method in the machine
learning literature [42], as $z\in C_{\epsilon}^{-1}(x)$ can be shown to be
equivalent to a (usually non-Euclidean) regularization term, which directly
interacts with the monotonicity property of $U$. While it is known that DMD
could converge exactly in subclasses of finite games that exhibit symmetric
interior NEs, whereby $C_{\epsilon}$ is also chosen to enforce symmetry (see
[28]), in general, it cannot converge exactly to a NE, which also means that
it cannot converge exactly to a VSS. This directs our attention to
alternatives methods, such as higher-order augmentation of game dynamics [28,
34, 36, 35].

##  4 Second-Order Mirror Descent Dynamics

We now propose the second-order mirror descent, which in terms of each player
$p$ appears as,

$\begin{cases}\dot{z}^{p}&\\!=\gamma(U^{p}(x)\\!-\\!\alpha(x^{p}-\xi^{p})),\,\,\quad\dot{\xi}^{p}\\!=\beta(x^{p}\\!-\\!\xi^{p}),\\\
x^{p}&\\!=C^{p}_{\epsilon}(z^{p}),\end{cases}$ (MD2)

and in stacked notation,

$\dot{z}=\gamma(U(x)\\!-\\!\alpha(x\\!-\\!\xi)),\,\,\quad\dot{\xi}=\beta(x-\xi),\quad
x=C_{\epsilon}(z),$ (13)

where $x=(x^{p})_{p\in\mathcal{N}}$, $z=(z^{p})_{p\in\mathcal{N}}$,
$\xi=(\xi^{p})_{p\in\mathcal{N}}$,
$C_{\epsilon}=(C^{p}_{\epsilon})_{p\in\mathcal{N}}$,
$U=(U^{p})_{p\in\mathcal{N}}$. The rest point condition for MD2 is,

$\mathbf{0}=U(\overline{x}),\quad\mathbf{0}=\overline{x}-\overline{\xi},\quad\overline{x}=C_{\epsilon}(\overline{z}),\quad\overline{z}\in
C^{-1}_{\epsilon}(\overline{x}),$ (14)

which coincides with that of MD, i.e., $\overline{x}=x^{\star}$ is an interior
NE.

######  Remark 4.

(Intuitive learning interpretation of MD2) We can explicitly write
$\xi^{p}(t)=(\xi^{p}_{i}(t))_{i\in[n_{p}]}$ as,

$\xi^{p}_{i}(t)=e^{{-\beta}t}\xi_{i}^{p}(0)+\beta\smallint\nolimits_{0}^{t}e^{{-\beta}(t-\tau)}x^{p}_{i}(\tau)\mathrm{d}\tau,$
(15)

which represents an “exponential weighting” of the strategies. Hence, we refer
to $\xi$ as the primal aggregate, and $z$ can be referred as the dual
aggregate. We note that $\xi^{p}$ resides in the unconstrained (ambient) space
$\mathbb{R}^{n_{p}}$ containing the action set. Note that, whenever
$\xi^{p}(0)$ is initialized in $\Omega^{p}$, then $\xi^{p}(t)\in\Omega^{p}$
for all times. Furthermore, noting that we can re-write $\dot{z}^{p}$
subsystem in MD2 as
$\dot{z}^{p}\\!=\gamma(U^{p}(x)\\!-\\!\alpha\beta^{-1}\dot{\xi}^{p})$. By
assuming $z^{p}_{i}(0)=0,\forall i,p$,

$\displaystyle
z^{p}_{i}(t)=\gamma(\smallint\nolimits_{0}^{t}U^{p}_{i}(x(\tau))\mathrm{d}\tau-\alpha\beta^{-1}\smallint\nolimits_{0}^{t}\dot{\xi}^{p}_{i}(\tau)\mathrm{d}\tau)$
(16)
$\displaystyle=\gamma(\smallint\nolimits_{0}^{t}U^{p}_{i}(x)\mathrm{d}\tau\\!-\\!\alpha
e^{{-\beta}t}(\beta^{-1}(1\\!-\\!e^{{\beta}t})\xi_{i}^{p}(0)\\!+\\!\smallint\nolimits_{0}^{t}e^{\beta\tau}x^{p}_{i}\mathrm{d}\tau)).$

It can be seen, when $t$ is small, the effects of $\xi^{p}_{i}(0)$ and
$x^{p}_{i}$ are at their largest, so the presence of $\dot{\xi}$ term takes
into consideration the uncertainty during the initial periods of play. Hence,
the incorporation of $\dot{\xi}$ has an exploratory effect on the play, which
can result in strategies being played more conservatively, preventing the
players from immediately reaching a deadlock.

######  Remark 5.

(Comparison with Existing Dynamics) The closest dynamics related to MD2 is the
higher-order exponentially discounted dynamics (H-EXPD-RL) for finite games
[28]. Specifically, when the higher-order augmentation is taken to be a high-
pass filter, we obtain,

$\begin{cases}\dot{z}&=\gamma(U(x)-z-\alpha(x+\xi)),\,\,\quad\dot{\xi}=\beta(-x-\xi),\\\
x&=C_{\epsilon}(z).\end{cases}$ (17)

(17) can be seen as the second-order extension of DMD. In [28] it was observed
that (17) is robust to a greater degree of parameter perturbation in monotone
games as compared to DMD.

Let’s now consider MD2 as a purely second-order ODE in the dual-space. To
simplify our calculations, we assume $x(t)=C_{\epsilon}(z(t))$ is
differentiable. By using MD2 and noting that
$\dot{z}\\!=\gamma(U(x)\\!-\\!\alpha\beta^{-1}\dot{\xi}))$, taking a time
derivative of $\dot{z}$ and making the assumption that $U$ and $C_{\epsilon}$
are $\mathcal{C}^{1}$, we have,
$\ddot{z}=\gamma(\operatorname{\mathbf{J}}_{U}(x)\dot{x}-\alpha\beta^{-1}\ddot{\xi})$,
which, along with $\ddot{\xi}=\beta(\dot{x}-\dot{\xi})$, $x=C_{\epsilon}(z)$,
$\dot{x}=\operatorname{\mathbf{J}}_{C_{\epsilon}(z)}\dot{z}$,
$\dot{\xi}=\beta\alpha^{-1}(U(x)-\gamma^{-1}\dot{z})$, we obtain,

$\ddot{z}=\gamma\left[\operatorname{\mathbf{J}}_{U}(C_{\epsilon}(z))\operatorname{\mathbf{J}}_{C_{\epsilon}}(z)-\alpha\operatorname{\mathbf{J}}_{C_{\epsilon}}(z)-\beta\gamma^{-1}\mathpzc{I}\right]\dot{z}+\gamma\beta
U(x),$ (18)

or by using the chain-rule for Jacobian,

$\displaystyle\begin{cases}\ddot{z}&=\gamma\left[\operatorname{\mathbf{J}}_{U\circ
C_{\epsilon}}(z)-\alpha\operatorname{\mathbf{J}}_{C_{\epsilon}}(z)-\beta\gamma^{-1}\mathpzc{I}\right]\dot{z}+\gamma\beta
U(x),\\\ x&=C_{\epsilon}(z),\end{cases}$ (19)

where $U\circ C_{\epsilon}\coloneqq U(C_{\epsilon})$. To the best of our
knowledge, most of the existing second-order continuous gradient-type dynamics
(see [40] and references therein) cannot recover MD2 due to the presence of
the Jacobian of $C_{\epsilon}$. This is because in the unconstrained, primal
setting (as studied in [40, 27]), $C_{\epsilon}$ is the identity and hence,
$\operatorname{\mathbf{J}}_{C_{\epsilon}}(z)=\operatorname{\mathpzc{I}},\forall
z$ and its contribution is lumped together with
$\beta{\gamma}^{-1}\operatorname{\mathpzc{I}}$.

We proceed to demonstrate that the primal aggregate $\xi$ contributes to the
convergence of MD2 beyond strict VSS and such convergence depends on both the
property of the regularizer as well as the topological properties of the
underlying strategy set.

######  Theorem 1.

Let $\mathcal{G}$ be a concave game and assume that every interior NE is
globally merely VS. Suppose that all players choose strategies according to
MD2 (13). Let $x=(x^{p}(t))_{p\in\mathcal{N}}=C_{\epsilon}(z(t))$ be generated
by (13) and $C_{\epsilon}=(C_{\epsilon}^{p})_{p\in\mathcal{N}}$ be induced by
$\psi^{p}_{\epsilon}\coloneqq\epsilon\vartheta^{p}$. Let $x^{\star}$ denote an
interior mere VSS (possibly from a set of them). For any
$\alpha,\beta,\gamma,\epsilon>0$,
$x(0)=C_{\epsilon}(z(0))\in\Omega,z(0),\xi(0)\in\mathbb{R}^{n}$,

  * (i)

suppose $\vartheta^{p}$ is twice-continuously differentiable
($\mathcal{C}^{2}$) and satisfies 3(i), $\forall p$, then
$x(t)=C_{\epsilon}(z(t))$ converges to $x^{\star}$.

  * (ii)

assume $\Omega^{p}$ is compact and $\vartheta^{p}$ satisfies 3(ii), $\forall
p$, then $x(t)=C_{\epsilon}(z(t))$ converges to $x^{\star}$.

######  Remark 6.

In the appendix, we show that the convergence guarantee of 1 can in fact go
beyond a mere VSS to that of a nearly mere VSS, i.e., a weak VSS that is
approximately a mere VSS. This shows that MD2 is robust around the mere VSS.

We now sharpen 1 by considering several other classes of games in which MD2 is
guaranteed to converge. From [3, 51],

######  Definition 3.

The game $\mathcal{G}$ is,

  1. (i)

pseudo-monotone (resp. strictly pseudo-monotone), if
$U(x^{\prime})^{\top}(x-x^{\prime})\leq 0\implies
U(x)^{\top}(x-x^{\prime})\leq 0,\forall x,x^{\prime}\in\Omega$ (resp.
$U(x)^{\top}(x-x^{\prime})<0$, with equality iff $x=x^{\prime}$.)

  2. (ii)

quasi-monotone (resp. strictly quasi-monotone),
$U(x^{\prime})^{\top}(x-x^{\prime})<0\implies U(x)^{\top}(x-x^{\prime})\leq
0,\forall x,x^{\prime}\in\Omega$ (resp. $U(x)^{\top}(x-x^{\prime})<0$, with
equality iff $x=x^{\prime}$).

Observe that when $x^{\prime}$ in the above definitions is NE $x^{\star}$ of a
pseudo-monotone game (or a globally strict NE of a quasi-monotone game), it
automatically implies the associated NE is globally mere VS. Moreover,
$x^{\star}$ is globally strictly VS when the above games are strict. Since MD
converges to globally strict VSS [2], therefore it converges to NE in strictly
pseudo-monotone games and strict NE in strictly quasi-monotone games. Our next
corollary, which immediately follows from 1, partially generalizes these
results to the non-strict setting under MD2.

######  Corollary 1.

Let $\mathcal{G}$ be a concave game and assume that every NE is interior.
Suppose that all players choose strategies according to MD2 (13). Let
$x=(x^{p}(t))_{p\in\mathcal{N}}=C_{\epsilon}(z(t))$ be generated by (13) and
$C_{\epsilon}=(C_{\epsilon}^{p})_{p\in\mathcal{N}}$ be induced by
$\psi^{p}_{\epsilon}\coloneqq\epsilon\vartheta^{p}$, $\vartheta^{p}$ is
$\mathcal{C}^{2}$ and satisfies 3(i), $\forall p$, then for any
$\alpha,\beta,\gamma,\epsilon>0$ and any $x(0)=C_{\epsilon}(z(0))$,

  * (i)

$x(t)$ converges to an interior NE $x^{\star}$ whenever $\mathcal{G}$ is
merely monotone or pseudo-monotone, and,

  * (ii)

$x(t)$ converges to an interior strict NE $x^{\star}$ if $\mathcal{G}$ is
quasi-monotone, whenever $x^{\star}$ exists.

The same conclusions hold whenever $\vartheta^{p}$ satisfies 3(ii) and
$\Omega^{p}$ is compact $\forall p$.

######  Remark 7.

We note that the existence of a NE in pseudo-monotone games is guaranteed
under our continuous game assumption [51, Theorem 2.3.5] (possibly requiring
$-\mathcal{U}^{p}$ to be coercive, see 1), and the non-empty set of NEs
coincides with the set of globally mere VSS [3, Theorem 2]. However, [3]
pointed out that the existence of a globally merely VS for quasi-monotone
games is not guaranteed in general. Instead, a globally merely VS exists under
the stronger property of properly quasi-monotone [3, Theorem 3].

##  5 Additional Convergence Properties of MD2

In this section, we investigate two additional properties of MD2, namely, that
of rate of convergence and regret minimization. To account for the geometry of
the problem, we provide a relative extension to the strong VSS.

######  Definition 4.

Let $h:\mathbb{R}^{n}\to\mathbb{R}\cup\\{\infty\\}$ be any differentiable,
convex function with domain $\operatorname{dom}(h)=\Omega$. Then
$x^{\star}\in\Omega$ is $\eta$-relatively strongly VS (with respect to $h$) if
for all $x\in\Omega$,
$U(x)^{\top}(x-x^{\star})\leq-\eta(D_{h}(x,x^{\star})+D_{h}(x^{\star},x)),$
for some $\eta\\!>\\!0$, where $D_{h}$ is the Bregman divergence of $h$.

We note that 4 is analogous to that of a relatively strongly monotone game,
which was previously introduced in [30]. Following [30, Theorem 4.4], it can
be shown that MD converges to a relatively strongly VSS in
$\mathcal{O}(e^{-\gamma\eta\epsilon^{-1}t})$ given a strongly-induced
$C_{\epsilon}$ that is adapted to the geometry of this VSS. We now wish to
provide a similar result for the convergence of MD2 towards a relatively
strong VSS. However, exponential convergence does not follow from 1. Instead,
we propose an augmented version of MD2 that exhibits exponential convergence,

$\begin{cases}\dot{z}&=(U(x)-\gamma\alpha(x-\xi)),\qquad\dot{\xi}=\beta(x-\xi),\\\
x&=C_{\epsilon}(z),\qquad\qquad\qquad\qquad\dot{\gamma}=-\eta\epsilon^{-1}\gamma,\end{cases}$
(MD2$\gamma$)

where $\gamma(0)>0$. Observe that MD2γ is equivalent to the non-autonomous
system,

$\begin{cases}\dot{z}&=(U(x)-e^{-\eta\epsilon^{-1}t}\gamma(0)\alpha(x-\xi)),\quad\dot{\xi}=\beta(x-\xi),\\\
x&=C_{\epsilon}(z),\end{cases}$ (20)

whose rest point condition coincides with that of MD2.

######  Remark 8.

From (20), MD2$\gamma$ can be seen as MD with a vanishing perturbation
$g(t,z,\xi)=e^{-\eta\epsilon^{-1}t}\alpha(0)(x-\xi)$, which allows for the
following simple learning interpretation: when the players are aware that the
game being played has a $\eta$-strong VSS, they no longer bother with
exploring the strategy space during the initial stages and instead discard the
extra information represented by $x-\xi$ exponentially fast.

######  Theorem 2.

Let $\mathcal{G}$ be a concave game with a unique interior $\eta$-strongly VSS
relative to $h(x)=\sum_{p\in\mathcal{N}}\vartheta^{p}(x^{p})$, denoted by
$x^{\star}$, where $x=(x^{p}(t))_{p\in\mathcal{N}}\\!=\\!C_{\epsilon}(z(t))$
is generated by MD2$\gamma$ and
$C_{\epsilon}\\!=\\!(C_{\epsilon}^{p})_{p\in\mathcal{N}}$ is induced by
$\psi^{p}_{\epsilon}\\!\coloneqq\\!\epsilon\vartheta^{p}$, $\vartheta^{p}$
satisfies 3(ii). Then for any
$\alpha,\beta,\epsilon,\eta,\gamma_{0}\\!=\\!\gamma(0)\\!>\\!0$ and any
$x_{0}\\!=\\!(x^{p}(0))_{p\in\mathcal{N}}\\!=\\!C_{\epsilon}(z(0)),z(0),\xi_{0}\\!=\\!\xi(0)\in\mathbb{R}^{n}$,
$x$ converges to $x^{\star}$ with the rate,

$D_{h}(x^{\star},x)\leq
e^{-\eta{\epsilon}^{-1}t}(\dfrac{\alpha\gamma_{0}}{2\beta\epsilon}\|\xi_{0}-x^{\star}\|_{2}^{2}+D_{h}(x^{\star},x_{0})).$
(21)

Furthermore, since $\vartheta^{p}$ satisfies 3(ii), i.e., $\vartheta^{p}$ is
$\rho$-strongly convex, therefore,

$\|x^{\star}-x\|_{2}^{2}\leq
2\rho^{-1}e^{-\eta{\epsilon}^{-1}t}(\dfrac{\alpha\gamma_{0}}{2\beta\epsilon}\|\xi_{0}-x^{\star}\|_{2}^{2}+D_{h}(x^{\star},x_{0})).$
(22)

######  Remark 9.

The above rate results confirm our intuition: the upper-bound rate of
convergence diminishes faster when an equilibrium is more strongly VSS
$(\eta\uparrow)$, when $C_{\epsilon}$ is induced by an even more strongly
convex regularizer ($\rho\uparrow$), when the regularization parameter goes
down ($\epsilon\downarrow$), so that $C_{\epsilon}$ approximates a best-
response map, or when the extra information $x-\xi$ diminishes more rapidly
($\beta\uparrow$).

Next, we turn our attention towards the question of regret minimization. Let
us define the time-averaged (external) regret function of player
$p\in\mathcal{N}$ as, $\mathpzc{R}^{p}:[0,\infty)\to\mathbb{R}$,
$\mathpzc{R}^{p}(0)=0$,

$\mathpzc{R}^{p}(t)=\max\limits_{y^{p}\in\Omega^{p}}\dfrac{1}{t}\smallint\nolimits_{0}^{t}\mathcal{U}^{p}(y^{p};x^{-p}(\tau))-\mathcal{U}^{p}(x(\tau))\mathrm{d}\tau.$
(23)

Intuitively, the value of the regret at time $t$ represents the time-averaged
sum total of the payoff difference between the actual strategy
$x(t)=(x^{p}(t);x^{-p}(t))$ versus the best strategy $(y^{p};x^{-p}(t))$ that
the player $p$ could have played in hindsight, given any opponents’ strategy
$x^{-p}(t)\in\Omega^{-p}$. Player $p$’s dynamics that generate $x^{p}(t)$ is
said to achieve no-regret with respect to (23) if
$\textstyle\limsup_{t\to\infty}\mathpzc{R}^{p}(t)\leq 0$.

While MD was shown to achieve no-regret in finite games, it cannot converge in
ZS finite games with an interior mixed equilibrium [22]. In contrast to [22,
21], our next result, coupled with 1, shows that MD2 achieves both no-regret
as well as exact convergence whenever the (interior) equilibrium is a mere
VSS.

######  Theorem 3.

Let $\mathcal{G}$ be a concave game, where $\Omega^{p}$ is compact for all
$p\in\mathcal{N}$. Suppose that each player $p\in\mathcal{N}$ is using MD2
(13), and $C_{\epsilon}\\!=\\!(C_{\epsilon}^{p})_{p\in\mathcal{N}}$ is induced
by $\psi^{p}_{\epsilon}\\!\coloneqq\\!\epsilon\vartheta^{p}$, $\vartheta^{p}$
satisfies 3(i) or 3(ii). Then for every continuous trajectory of $x^{-p}(t)$
of opponents of player $p$, each player achieves no-regret independently from
the rest of the players.

##  6 Construction of Second-Order Primal-Space Dynamics

In this section, we show that MD2 can either be used to create new primal-
space dynamics (that is, dynamics that solely evolve on $\Omega$) or recover
existing ones. To simplify our presentation and to avoid the technicality of
non-differentiable mirror maps $C_{\epsilon}$ (for which examples of induced
primal-space dynamics can still be generated through a technical treatment,
see the so-called projection dynamics discussed in [19]), we impose the
following general assumption,

######  Assumption 4.

$\vartheta^{p}\\!:\\!\mathbb{R}^{n_{p}}\\!\\!\to\\!\mathbb{R}\\!\cup\\!\\{\\!\infty\\!\\}$
is steep and induces a $\mathcal{C}^{1}$ mirror map
$C_{\epsilon}=(C^{p}_{\epsilon})_{p\in\mathcal{N}}=(\nabla\psi^{p\star}_{\epsilon})_{p\in\mathcal{N}}$.

Under this assumption, taking the time-derivative of $x=C_{\epsilon}(z)$ shows
that MD2 can be written as a pair of differential inclusions,

$\dot{\xi}\in\beta(x-\xi),\quad\dot{x}\in\gamma\operatorname{\mathbf{J}}_{C_{\epsilon}}(z)(U(x)-\alpha(x-\xi)),$
(24)

where $\gamma>0,\alpha,\beta\geq 0$ and the Jacobian of $C_{\epsilon}$ is,

$\operatorname{\mathbf{J}}_{C_{\epsilon}}(z)\coloneqq\text{blkdiag}(\operatorname{\mathbf{J}}_{C^{p}_{\epsilon}}(z^{p}))=\begin{bmatrix}\operatorname{\mathbf{J}}_{C^{1}_{\epsilon}}(z^{1})&&\\\
&\ddots&\\\
&&\operatorname{\mathbf{J}}_{C^{N}_{\epsilon}}(z^{N})\end{bmatrix}\in\mathbb{R}^{n\times
n},$

and each
$\operatorname{\mathbf{J}}_{C^{p}_{\epsilon}}(z^{p})=\operatorname{\mathbf{J}}_{{\nabla\psi^{p}_{\epsilon}}^{\star}}(z^{p})=\nabla^{2}{\psi_{\epsilon}^{p}}^{\star}(z^{p})\in\mathbb{R}^{n_{p}\times
n_{p}}$. We note that the inclusions in (24) come from $z\in
C^{-1}_{\epsilon}(x)$. We can further express (24) in each $p$ as,

$\dot{\xi}^{p}\in\beta(x^{p}-\xi^{p}),\,\,\dot{x}^{p}\in\gamma\nabla^{2}{\psi_{\epsilon}^{p}}^{\star}(z^{p})(U^{p}(x)-\alpha(x^{p}-\xi^{p})).$
(25)

In the following, we offer two general ways of re-writing (25) as a set of
ODEs in the primal-space.

General Case: Our next result shows that, in the general case, finding the
primal-space dynamics associated with MD2 generally amounts to evaluating the
operator
$\nabla^{2}{\psi^{p}_{\epsilon}}^{\star}\circ\nabla\psi^{p}_{\epsilon}\coloneqq\nabla^{2}{\psi^{p}_{\epsilon}}^{\star}(\nabla\psi^{p}_{\epsilon})$.

######  Proposition 2.

Suppose 4 holds, then (25) reduces to,

$\begin{cases}\dot{\xi}^{p}&=\beta(x^{p}-\xi^{p})\\\
\dot{x}^{p}&=\gamma\nabla^{2}{\psi^{p}_{\epsilon}}^{\star}(\nabla\psi^{p}_{\epsilon}(x^{p}))(U^{p}(x)-\alpha(x^{p}-\xi^{p})).\end{cases}$
(26)

######  Example 7.

(Simplex) Let $\Omega^{p}=\\{x^{p}\in\mathbb{R}^{n_{p}}_{\geq
0}|\|x^{p}\|_{1}=1\\}$ and consider
$\vartheta^{p}(x^{p})=\sum_{i=1}^{n_{p}}x_{i}^{p}\log(x_{i}^{p})$,
$\psi_{\epsilon}^{p}(x^{p})=\epsilon\vartheta^{p}(x^{p})$ and
${\psi^{p}_{\epsilon}}^{\star}(z^{p})=\epsilon\log(\sum_{i=1}^{n_{p}}\exp(\epsilon^{-1}z^{p}_{i}))$.
It can be shown that,
$\nabla^{2}{\psi_{\epsilon}^{p}}^{\star}(z^{p})=\epsilon^{-1}(\operatorname{diag}(C^{p}_{\epsilon}(z^{p}))-C^{p}_{\epsilon}(z^{p})C^{p}_{\epsilon}(z^{p})^{\top})$
where $\textstyle
C^{p}_{\epsilon}(z^{p})=(\exp(\epsilon^{-1}z^{p}_{i})(\sum_{j=1}^{n_{p}}\exp(\epsilon^{-1}z^{p}_{j}))^{-1})_{i\in[n_{p}]}$
is the softmax function (5). Let $x^{p}\in\operatorname{rint}(\Omega^{p})$. We
can show $\nabla\psi_{\epsilon}^{p}(x^{p})=\epsilon(\log(x^{p})+\mathbf{1})$
($\log$ is applied component-wise) and
$C^{p}_{\epsilon}(\nabla\psi^{p}_{\epsilon}(x^{p}))=x^{p}$. Taken together, we
have,
$\nabla^{2}{\psi^{p}_{\epsilon}}^{\star}\circ\nabla\psi^{p}_{\epsilon}=\nabla^{2}{\psi^{p}_{\epsilon}}^{\star}(\nabla\psi^{p}_{\epsilon}(x^{p}))=\epsilon^{-1}(\operatorname{diag}(x^{p})-x^{p}{x^{p}}^{\top})$.
By applying (26) and defining
$\mathpzc{H}(x^{p})\coloneqq\operatorname{diag}(x^{p})\\!-\\!x^{p}{x^{p}}^{\top}$,
we obtain (27) as the induced primal-space dynamics of MD2 on the simplex,

$\dot{\xi}^{p}\\!=\\!\beta(x^{p}\\!-\\!\xi^{p}),\quad\dot{x}^{p}\\!=\\!\gamma\epsilon^{-1}\mathpzc{H}(x^{p})(U^{p}(x)\\!-\\!\alpha(x^{p}\\!-\\!\xi^{p})),$
(27)

We note that (27) represents a version of second-order replicator dynamics
(RD), which is different from the ones introduced in [34, 36]. In particular,
when $\alpha=\beta=0$, we recover the (multi-population) replicator dynamics
[3, 13, 14, 19, 22, 45].

We can also reduce (27) to a single population model ($p=1$), which takes the
following form:

$\textstyle\begin{cases}\dot{\xi}_{i}&=\beta(x_{i}-\xi_{i}),\\\
\dot{x}_{i}&=x_{i}(U_{i}(x)-\dfrac{\alpha}{\beta}\dot{\xi}_{i}-\sum_{j\in\mathcal{A}}x_{j}(U_{j}(x)-\dfrac{\alpha}{\beta}\dot{\xi}_{j})),\end{cases}$
(28)

where $\mathcal{A}=\\{1,\ldots,n\\}$ represents $n$ subpopulations of non-
atomic players, and $U_{i}$ represents the fitness of subpopulation $i$. We
note that structurally (28) shares similarity with the anticipatory replicator
dynamics studied in [36], which can be written as,

$\textstyle\begin{cases}\dot{\xi}_{i}&=\beta(x_{i}-\xi_{i}),\\\
\dot{x}_{i}&=x_{i}(U_{i}(x+\alpha\dot{\xi})-\sum_{j\in\mathcal{A}}x_{j}(U_{j}(x+\alpha\dot{\xi}))).\end{cases}$
(29)

Note that the additional variable $\xi_{j}$ in (28) do not directly interacts
with the pseudo-gradient/payoff vector $U(x)$, as is the case for the
anticipatory RD.

Legendre and Supercoercive: A more structured scenario is when
$\psi^{p}_{\epsilon}\coloneqq\epsilon\vartheta^{p}$ is Legendre and
supercoercive for all $p$ (3(i)). Since $\psi^{p}_{\epsilon}$ is finite on
$\operatorname{int}(\Omega^{p})$, coercive, strictly convex and
$\mathcal{C}^{2}$, then applying [63, Example 11.9, p. 480],
${\psi^{p}_{\epsilon}}^{\star}$ shares the same properties and we have that
$[\nabla^{2}\psi^{p}_{\epsilon}(x^{p})]^{-1}=\nabla^{2}{\psi^{p}_{\epsilon}}^{\star}(z^{p})=({\nabla^{2}\psi_{\epsilon}^{p}}^{\star}\circ\nabla\psi_{\epsilon}^{p})(x^{p}),\forall
x^{p}=\nabla{\psi^{p}_{\epsilon}}^{\star}(z^{p})\in\operatorname{int}(\Omega^{p})$,
which allows us to write (25) or (26) as,

$\begin{cases}\dot{\xi}^{p}&=\beta(x^{p}-\xi^{p}),\\\
\dot{x}^{p}&=\gamma[\nabla^{2}\psi^{p}_{\epsilon}(x^{p})]^{-1}(U^{p}(x)-\alpha(x^{p}-\xi^{p})),\end{cases}$
(NG2)

which we refer to as the second-order natural gradient descent (NG2). NG2 is
so named because in the optimization setup ($p=1$), for $U=-\nabla f$, where
$f:\mathbb{R}^{n}\to\mathbb{R}\cup\\{\infty\\}$ is some objective function and
$\alpha,\beta=0$, is analogous to the so-called natural gradient descent [32],
$\dot{x}=-[\nabla^{2}\psi_{\epsilon}(x)]^{-1}\nabla f(x).$

######  Example 8.

(Unconstrained) Let $\Omega^{p}=\mathbb{R}^{n_{p}}$,
$\vartheta^{p}(x^{p})=\frac{1}{2}\|x^{p}\|_{2}^{2}$, then
$\nabla^{2}{\psi^{p}_{\epsilon}}(x^{p})=\epsilon\mathpzc{I}$ where
${\psi^{p}_{\epsilon}}=\epsilon\vartheta^{p}$. By comparing with NG2, we
obtain,

$\dot{\xi}^{p}=\beta(x^{p}-\xi^{p}),\quad\dot{x}^{p}=\gamma\epsilon^{-1}(U^{p}(x)-\alpha(x^{p}-\xi^{p})),$
(30)

which recovers the one recently introduced in [27]. (30) can be further shown
via discretization to recover algorithms such as optimistic gradient-
descent/ascent, Polyak’s heavy-ball method, among others (see [27]). While
(30) is by far the most standard choice for unconstrained action sets, more
than one type of dynamics can reside on the same action sets. To witness, let
$\Omega^{p}=\mathbb{R}^{n_{p}}$,
$\vartheta^{p}(x^{p})=\sum_{i=1}^{n_{p}}x_{i}^{p}\operatorname{arcsinh}(x^{p}_{i}/\wp)-\mathchoice{{\hbox{$\displaystyle\sqrt[\
]{x_{i}^{p2}+\wp^{2}\,}$}\lower 0.4pt\hbox{\vrule
height=6.10999pt,depth=-4.88802pt}}}{{\hbox{$\textstyle\sqrt[\
]{x_{i}^{p2}+\wp^{2}\,}$}\lower 0.4pt\hbox{\vrule
height=6.10999pt,depth=-4.88802pt}}}{{\hbox{$\scriptstyle\sqrt[\
]{x_{i}^{p2}+\wp^{2}\,}$}\lower 0.4pt\hbox{\vrule
height=4.30276pt,depth=-3.44223pt}}}{{\hbox{$\scriptscriptstyle\sqrt[\
]{x_{i}^{p2}+\wp^{2}\,}$}\lower 0.4pt\hbox{\vrule
height=3.44165pt,depth=-2.75334pt}}}\kern 1.00006pt,\wp>0$, which represents
an interpolation between entropic and Euclidean terms [62]. $\vartheta^{p}$ is
Legendre and supercoercive (its convex conjugate is
$\sum_{i=1}^{n_{p}}\wp\cosh(x^{p}_{i})$, which is also Legendre). Let
$\psi^{p}_{\epsilon}=\epsilon\vartheta^{p}$, then
$\nabla\psi^{p}_{\epsilon}=\epsilon\sum_{i=1}^{n_{p}}\operatorname{arcsinh}(x_{i}^{p}/\wp)$
and
$\nabla^{2}\psi^{p}_{\epsilon}=\epsilon\left[\operatorname{diag}(\begin{bmatrix}\mathchoice{{\hbox{$\displaystyle\sqrt[\
]{{x^{p}_{1}}^{2}+\wp^{2}\,}$}\lower 0.4pt\hbox{\vrule
height=7.31554pt,depth=-5.85246pt}}}{{\hbox{$\textstyle\sqrt[\
]{{x^{p}_{1}}^{2}+\wp^{2}\,}$}\lower 0.4pt\hbox{\vrule
height=7.31554pt,depth=-5.85246pt}}}{{\hbox{$\scriptstyle\sqrt[\
]{{x^{p}_{1}}^{2}+\wp^{2}\,}$}\lower 0.4pt\hbox{\vrule
height=5.16388pt,depth=-4.13113pt}}}{{\hbox{$\scriptscriptstyle\sqrt[\
]{{x^{p}_{1}}^{2}+\wp^{2}\,}$}\lower 0.4pt\hbox{\vrule
height=4.30276pt,depth=-3.44223pt}}}\kern
1.00006pt,\ldots,\mathchoice{{\hbox{$\displaystyle\sqrt[\
]{x^{p2}_{n^{p}}+\wp^{2}\,}$}\lower 0.4pt\hbox{\vrule
height=6.10999pt,depth=-4.88802pt}}}{{\hbox{$\textstyle\sqrt[\
]{x^{p2}_{n^{p}}+\wp^{2}\,}$}\lower 0.4pt\hbox{\vrule
height=6.10999pt,depth=-4.88802pt}}}{{\hbox{$\scriptstyle\sqrt[\
]{x^{p2}_{n^{p}}+\wp^{2}\,}$}\lower 0.4pt\hbox{\vrule
height=4.30276pt,depth=-3.44223pt}}}{{\hbox{$\scriptscriptstyle\sqrt[\
]{x^{p2}_{n^{p}}+\wp^{2}\,}$}\lower 0.4pt\hbox{\vrule
height=3.44165pt,depth=-2.75334pt}}}\kern
1.00006pt)\end{bmatrix}^{\top})\right]^{-1}.$ By comparing with NG2, we
obtain,

$\dot{\xi}^{p}=\beta(x^{p}-\xi^{p}),\quad\dot{x}^{p}=\gamma\epsilon^{-1}\mathpzc{S}(x^{p})(U^{p}(x)-\alpha(x^{p}-\xi^{p})).$
(31)

where
$\mathpzc{S}(x^{p})\coloneqq\operatorname{diag}(\begin{bmatrix}\mathchoice{{\hbox{$\displaystyle\sqrt[\
]{{x^{p}_{1}}^{2}+\wp^{2}\,}$}\lower 0.4pt\hbox{\vrule
height=7.31554pt,depth=-5.85246pt}}}{{\hbox{$\textstyle\sqrt[\
]{{x^{p}_{1}}^{2}+\wp^{2}\,}$}\lower 0.4pt\hbox{\vrule
height=7.31554pt,depth=-5.85246pt}}}{{\hbox{$\scriptstyle\sqrt[\
]{{x^{p}_{1}}^{2}+\wp^{2}\,}$}\lower 0.4pt\hbox{\vrule
height=5.16388pt,depth=-4.13113pt}}}{{\hbox{$\scriptscriptstyle\sqrt[\
]{{x^{p}_{1}}^{2}+\wp^{2}\,}$}\lower 0.4pt\hbox{\vrule
height=4.30276pt,depth=-3.44223pt}}}\kern
1.00006pt,\ldots,\mathchoice{{\hbox{$\displaystyle\sqrt[\
]{x^{p2}_{n^{p}}+\wp^{2}\,}$}\lower 0.4pt\hbox{\vrule
height=6.10999pt,depth=-4.88802pt}}}{{\hbox{$\textstyle\sqrt[\
]{x^{p2}_{n^{p}}+\wp^{2}\,}$}\lower 0.4pt\hbox{\vrule
height=6.10999pt,depth=-4.88802pt}}}{{\hbox{$\scriptstyle\sqrt[\
]{x^{p2}_{n^{p}}+\wp^{2}\,}$}\lower 0.4pt\hbox{\vrule
height=4.30276pt,depth=-3.44223pt}}}{{\hbox{$\scriptscriptstyle\sqrt[\
]{x^{p2}_{n^{p}}+\wp^{2}\,}$}\lower 0.4pt\hbox{\vrule
height=3.44165pt,depth=-2.75334pt}}}\kern 1.00006pt\end{bmatrix}^{\top})$,
which translates into a version of (30) with a strategy-dependent coefficient.

######  Example 9.

(Orthant) Let $\Omega^{p}=\mathbb{R}^{n_{p}}_{\geq 0}$ and consider
$\vartheta^{p}(x^{p})=\sum_{i=1}^{n_{p}}x_{i}^{p}\log(x_{i}^{p})-x_{i}^{p}$
with $0\log(0)=0$. We can show
$\nabla^{2}\psi_{\epsilon}^{p}(x^{p})=\epsilon[\operatorname{diag}({x^{p}})]^{-1}$
where $\psi_{\epsilon}^{p}=\epsilon\vartheta^{p}$. Using NG2, we immediately
arrive at,

$\dot{\xi}^{p}=\beta(x^{p}-\xi^{p})\quad\dot{x}^{p}=\gamma\epsilon^{-1}\operatorname{diag}(x^{p})(U^{p}(x)-\alpha(x^{p}-\xi^{p})).$
(32)

We note that the first-order variant $(\alpha=\beta=0)$ of (32) was recently
studied by [61] as the mean-dynamics to a payoff-based learning dynamics in
continuous games with non-negative orthant action sets. Suppose instead
$\vartheta^{p}(x^{p})=-\sum_{i=1}^{n}\ln(x^{p}_{i})$ (the log-barrier
function), which leads to
$\nabla^{2}\psi_{\epsilon}^{p}(x^{p})=\epsilon[\operatorname{diag}(x^{p2})]^{-1}$
and,

$\\!\dot{\xi}^{p}=\beta(x^{p}\\!-\\!\xi^{p})\quad\\!\dot{x}^{p}=\gamma\epsilon^{-1}\operatorname{diag}(x^{p2})(U^{p}(x)\\!-\\!\alpha(x^{p}\\!-\\!\xi^{p})).$
(33)

We stress however that the log-barrier does not fall under our definition of a
regularizer, as $\operatorname{dom}(\vartheta^{p})$ is a proper subset of
$\Omega^{p}$. (33) with $\alpha=\beta=0$ has been previously studied in
optimization, see [64].

##  7 Discrete-Time Second-Order Dual Averaging with Noisy Observations

So far we have considered a continuous-time setup whereby each player is able
to acquire a partial-gradient $U^{p}$ at each time instance. In practical
scenarios where the games are played in discrete-time, the acquired pseudo-
gradient information could be corrupted due to a multitude of reasons, such as
a noisy communication channel. This leads us to consider the so-called “noise-
corrupted” pseudo-gradient scenario (also known as “semi-bandit” learning). In
this setting, each player receives a realization of a so-called noise-
corrupted version of the true pseudo-gradient,

$\widehat{U}^{p}_{k+1}\coloneqq U^{p}(x_{k})+\zeta^{p}_{k+1},$

where $\zeta^{p}_{k}$ is some noise process. A special case of the semi-bandit
scenario is when the payoff is obtained as the expectation of the true payoff,
i.e.,
$\mathcal{U}^{p}(x_{k}^{p};x^{-p}_{k})=\mathbb{E}[\mathcal{U}^{p}(x^{p}_{k};x^{-p}_{k},v^{p}_{k})]$,
$v^{p}_{k}$ some random vector, where $\mathbb{E}$ denotes the expectation
operator. Here $\widehat{U}^{p}_{k}$ is an estimate of the expected partial-
gradient
$\nabla_{x^{p}}\mathbb{E}[\mathcal{U}^{p}(x^{p}_{k};x^{-p}_{k},v^{p}_{k})]$.

Let’s consider the convergence of the discrete-time MD2 with noisy
observations which we refer to as second-order dual averaging or DA2,

$\begin{cases}X^{p}_{k}&=C^{p}_{\epsilon}(Z^{p}_{k})\\\
Z^{p}_{k+1}&=Z^{p}_{k}+\gamma
t_{k+1}(\widehat{U}^{p}_{k+1}-\alpha(X^{p}_{k}-\Xi^{p}_{k}))\\\
\Xi^{p}_{k+1}&=\Xi^{p}_{k}+\tau_{k+1}\beta(X^{p}_{k}-\Xi^{p}_{k}),\end{cases}$
(DA2)

where $X^{p}_{k},Z^{p}_{k},\Xi^{p}_{k}$ are the stochastic counter-part of
(resp.) $x^{p},z^{p},\xi^{p}$ at time $k$. $\\{t_{k}\\}_{k\in\mathbb{N}}$,
$\\{\tau_{k}\\}_{k\in\mathbb{N}}$ denote deterministic sequences of non-
increasing step-sizes, assumed to be common for all players.
$\alpha,\beta,\gamma>0$ are the auxiliary parameters as from before. When
$\alpha=\beta=0$, we refer to the resulting expression as the dual averaging
with noisy observations,

$\begin{cases}X^{p}_{k}&=C^{p}_{\epsilon}(Z^{p}_{k})\\\
Z^{p}_{k+1}&=Z^{p}_{k}+\gamma t_{k+1}\widehat{U}^{p}_{k+1},\\\ \end{cases}$
(DA)

which coincides with the dual averaging scheme studied in [2] for $\gamma=1$.

To analyze the convergence behavior of DA2, we impose the following set of
regularity assumptions [55, 56].

######  Assumption 5.

($\ell^{2}\backslash\ell^{1}$-Summability and Diminishing Step-sizes)

$\begin{split}&\textstyle\sum_{k\in\mathbb{N}}t_{k}=\infty\quad\textstyle\sum_{k\in\mathbb{N}}t^{2}_{k}<\infty\quad\textstyle\lim_{k\to\infty}t_{k}=0,\\\
&\textstyle\sum_{k\in\mathbb{N}}\tau_{k}=\infty\quad\sum_{k\in\mathbb{N}}\tau^{2}_{k}<\infty\quad\lim_{k\to\infty}\tau_{k}=0.\end{split}$
(34)

######  Assumption 6.

($L^{2}$-Bounded Martingale Difference Noise) Assume
$\\{\zeta_{k}\\}_{k\in\mathbb{N}},\zeta_{k}=(\zeta^{p}_{k})_{p\in\mathcal{N}}$
is a $L^{2}$-bounded marginal difference process adapted to the filtration
$\\{\mathcal{F}_{k}\\}_{k\in\mathbb{N}}$: each $\zeta_{k}$ is a random vector
that is measurable with respect to $\mathcal{F}_{k}$ for each $k$, where each
$\mathcal{F}_{k}$ is the $\sigma$-field, i.e.,
$\mathcal{F}_{k}=\sigma(\Xi_{0},Z_{0},\zeta_{0},\ldots,\zeta_{k})$ and
$\mathcal{F}_{k}\subseteq\mathcal{F}_{k+1}$. In particular,
$\\{\zeta_{k}\\}_{k\in\mathbb{N}}$ satisfies,

$\operatorname{\mathbb{E}}[\zeta_{k+1}|\mathcal{F}_{k}]=0,\forall
k\in\mathbb{N}\quad a.s.,$ (35)

and for some $\sigma\geq 0$,

$\operatorname{\mathbb{E}}[\|\zeta_{k+1}\|_{\star}^{2}|\mathcal{F}_{k}]\leq\sigma^{2},\forall
k\in\mathbb{N}\quad a.s.$ (36)

######  Assumption 7.

(Bounded Iterates)

$\textstyle\sup_{k}\|Z^{p}_{k}\|<\infty\quad\sup_{k}\|\Xi^{p}_{k}\|<\infty,\forall
p.$ (37)

Finally, we impose “global integrability” on MD2, i.e., MD2 has a complete
vector field. To do so, we need to convert MD2 into a first-order system by
defining: $\omega^{p}\coloneqq(\xi^{p},z^{p})$, which generates the following
stacked system on $\mathbb{R}^{2n}\cong\mathbb{R}^{n}\times\mathbb{R}^{n}$,

$\begin{split}\dot{\omega}^{p}&=\begin{bmatrix}\dot{\xi}^{p}\\\
\dot{z}^{p}\end{bmatrix}=\begin{bmatrix}\beta(C_{\epsilon}^{p}(z^{p})-\xi^{p})\\\
\gamma(U^{p}(C_{\epsilon}(z))-\alpha(C_{\epsilon}^{p}(z^{p})-\xi^{p}))\end{bmatrix},\\\
&=\begin{bmatrix}\beta C_{\epsilon}^{p}(z^{p})\\\
\gamma(U^{p}(C_{\epsilon}(z))-\alpha
C_{\epsilon}^{p}(z^{p}))\end{bmatrix}+\begin{bmatrix}-\beta\xi^{p}\\\
\gamma\alpha\xi^{p}\end{bmatrix},\end{split}$ (38)

or equivalently, $\dot{\omega}^{p}=F^{p}(\omega)=g^{p}(z)+A^{p}(\xi^{p}),$
where, $A^{p}(\xi^{p})=\begin{bmatrix}-\beta\xi^{p}\\\
\gamma\alpha\xi^{p}\end{bmatrix}$ and $g^{p}(z)=\begin{bmatrix}\beta
C_{\epsilon}^{p}(z^{p})\\\ \gamma(U^{p}(C_{\epsilon}(z))-\alpha
C_{\epsilon}^{p}(z^{p}))\end{bmatrix}.$

We then proceed to impose the following assumption on the overall system,

$\dot{\omega}=F(\omega),\omega=(\omega^{p})_{p\in\mathcal{N}},F=(F^{p})_{p\in\mathcal{N}}.$
(39)

######  Assumption 8.

(Global Integrability) The vector field $F:\mathbb{R}^{2n}\to\mathbb{R}^{2n}$
of (39) is continuous globally integrable, that is, for every initial
condition $(\xi(0),z(0))\in\mathbb{R}^{2n}$, the unique solution of (39) is
defined for all $t\in\mathbb{R}$.

######  Remark 10.

We note that global integrability of MD2 is satisfied whenever the following
holds,

  * (i)

Suppose $C^{p}_{\epsilon}$ is bounded on $\mathbb{R}^{n_{p}}$ and $U^{p}$ is
bounded continuous locally Lipschitz on $\operatorname{ran}(C^{p}_{\epsilon})$
for all $p\in\mathcal{N}$, then $F$ (39) is bounded locally Lipschitz and
hence continuous globally integrable. This follows from [55].

  * (ii)

Suppose $U^{p}$ is continuous locally Lipschitz and both $C^{p}_{\epsilon}$
and $U^{p}\circ C_{\epsilon}$ are sublinear for all $p$, that is,

$\limsup_{\|z\|_{2}\to\infty}\dfrac{\|C^{p}_{\epsilon}(z)\|_{2}}{\|z\|_{2}}<\infty,\quad\limsup_{\|z\|_{2}\to\infty}\dfrac{\|U^{p}\circ
C_{\epsilon}(z)\|_{2}}{\|z\|_{2}}<\infty,$ (40)

then $F$ is sublinear and hence continuous and globally integrable. This
follows from [57].

Under 5 – 8, DA2 can be shown to track the continuous trajectories generated
by MD2 via stochastic approximation arguments [55, 56].

######  Theorem 4.

Let $\mathcal{G}$ be a concave game and assume that every interior NE is
globally merely VS. Suppose that all players choose strategies according to
DA2 and $C_{\epsilon}=(C_{\epsilon}^{p})_{p\in\mathcal{N}}$ is induced by
$\psi^{p}_{\epsilon}\coloneqq\epsilon\vartheta^{p}$, where $\vartheta^{p}$ is
$\mathcal{C}^{2}$ and satisfies 3(i). Suppose that the following assumptions
hold,

  * •

5 ($\ell^{2}\backslash\ell^{1}$ Summability and diminishing step-sizes)

  * •

6 ($L^{2}$-bounded Martingale difference noise),

  * •

7 (Bounded iterates)

  * •

8 (Global integrability)

then $X_{k}$ converges to an interior mere VSS of $\mathcal{G}$ almost surely.
The same conclusions hold whenever $\vartheta^{p}$ satisfies 3(ii) and
$\Omega^{p}$ is compact $\forall p$.

######  Remark 11.

The closest result to ours is [2, Theorem 4.7], but for DA (i.e.,
$\alpha,\beta=0$). 5, 6 are identical to theirs and 8 encapsulates the
Lipschitz continuous requirement for $U$. The major departure is 7, which does
not hold under certain circumstances, such as convergence towards (boundary)
strict VSS in finite games. Hence our result only deals with interior mere VSS
in general. To account for these boundary VSS possibly involves an extension
of the inductive shadowing argument used by [2] which we leave for future
work.

######  Remark 12.

In addition to DA [2], there are several other algorithms that converge in
similar settings. The most notable example is the mirror-prox/optimistic
mirror descent algorithm, which converges to mere VSS in perfect-gradient
feedback setting [47]. However, the authors of [47] noted that convergence in
null-coherent saddle point problems (a game with a mere VSS) fails in the
presence of noise. This is a key advantage of DA2 over mirror-prox. Another
algorithm is the stochastic iterative Tikhonov method of [26], which converges
in the semi-bandit setting. However, [26] requires the game to be strictly
monotone, whereas DA2 does not require monotonicity.

##  8 Simulations

In this section, we consider three illustrative examples. The first example is
the RPS game with a non-negative payoff for ties as in 3. We provide
convergence behavior of MD and MD2 as well as MDA towards mere and weak VSS.
We then provide the convergence behavior of DA2. The second example concerns a
wireless power control game previously studied in [59]. We show that DA2
converges in this game under different noise assumptions and study the effect
of the number of players. Finally, we provide an example involving a
generative adversarial network (GAN), which admits a locally mere VSS. We show
that DA2 also converges in this game under different noise assumptions. In
lieu of exact estimation of the basin of attraction for locally mere VSS,
which is difficult, we deal with all such cases through appropriate
initialization.

######  Example 10.

(RPS with non-negative payoff for ties) Consider the RPS game as discussed in
3. We set the game parameters to be $\mathpzc{l}=0,\mathpzc{w}=2$ and vary the
tie payoff parameter $\varsigma\in[1,\mathpzc{w})=[1,2)$ to induce different
VS properties on the NE
$x^{\star}=(x^{p})_{p\in\mathcal{N}},x^{p}=(1/3,1/3,1/3)$. We simulate both MD
and MD2 at initial conditions
$z(0)=\begin{bmatrix}3&2&1\end{bmatrix}^{\top},\xi(0)=\mathbf{0}$.
$\gamma,\beta,\alpha$ are kept as $1$.

We begin by contrasting the continuous-time MD and MD2. For $\varsigma=1$,
this game is merely monotone, $x^{\star}$ is globally merely VS, MD2 converges
by 1 while MD diverges (Figure 1). For $\varsigma=1.1$, the game is
$0.1$-weakly monotone, but since $x^{\star}$ is only $0.1$-weak VSS therefore
it can be considered a “nearly” mere VSS. By 6 and 13, MD2 still converges
while MD approaches a heteroclinic orbit (Figure 2).

Figure 1: $\varsigma=1$, $\mathcal{G}$ merely monotone, $x^{\star}$ is merely
VS. MD2 converges (Embedded: MD, cycling). Figure 2: $\varsigma=1.1$,
$\mathcal{G}$ $0.1$-weakly monotone, $x^{\star}$ is $0.1$-weakly VS. MD2
converges (Embedded: MD, cycling).

It is useful to also examine the advantage of MD2 over time-averaged MD (MDA).
While MDA does converge for $\varsigma=1$ towards the mere VSS, when the
equilibrium becomes just slightly weak ($\varsigma=1.1$), it no longer
converges and instead approaches a Shapley triangle [58] (Figure 3). This
shows time-averaging is in general not a panacea to non-convergence.

Figure 3: $\varsigma=1.1$, $\mathcal{G}$ $0.1$-weakly monotone, $x^{\star}$ is
$0.1$-weakly VS. MD2 converges (Embedded: MDA, cycling).

Next, we perform a set of experiments for the globally merely VS case
($\varsigma=1$) using DA2, where we assume the same initial condition as
before. For each of the following simulations, we separately perturb the
pseudo-gradient $U^{p}$ (9) with zero-mean Gaussian noise with the same
variance $\sigma^{2}_{\zeta}>0$ across all players. Such a perturbed gradient
setting could describe a game involving multiple users interacting over a
fully connected network with noisy communication channels. Figure 4 shows that
DA2 easily converge in the low variance regimes $\sigma^{2}_{\zeta}=1$. As we
employ larger variance, e.g., $\sigma^{2}_{\zeta}=10$, the standard step-sizes
no longer leads to convergence: a larger $t_{k}$ will amplify the additive
noise. Instead, we utilize step-size sequences of the form
$\dfrac{c}{k^{\iota}}$ and search over both $\iota\in(0,1),c>0$ for optimal
sets of parameters. The simulation with tuned step-size is shown in Figure 5.
We see that despite the larger noise injected into $U^{p}$, the strategies
$X=(X^{p})_{p\in\mathcal{N}}$ still converge toward the vicinity of the mere
VSS.

Figure 4: $\varsigma=1,\sigma_{\zeta}^{2}=1$, $t_{k}=4/k,\tau_{k}=1/k$. Figure
5: $\varsigma=1,\sigma^{2}_{\zeta}=10$,
$t_{k}=0.23/k^{0.48},\tau_{k}=0.34/k^{0.88}$.

######  Example 11.

(Wireless power control game with a non-concave potential) Consider the
wireless power control game in [59], where $\mathcal{N}=\\{1,\ldots,N\\}$
network users decide on intensities $x^{p}\in\mathbb{R}$ of power flow to send
over a wireless network. These intensities are converted to the transmitted
power through an exponential function, i.e., $\exp(x^{p})\in\mathbb{R}_{>0}$.
The payoff function for each user $p\in\mathcal{N}$ is modeled as,

$\textstyle\mathcal{U}^{p}(x^{p};x^{-p})=\log\left(1+\dfrac{a^{p}\exp(x^{p})}{1+\sum_{p\neq
r}a^{r}\exp(x^{r})}\right)-\mathpzc{K}^{p}(x^{p})$ (41)

where $a^{p}\in(0,1],\forall p$ and
$\mathpzc{K}^{p}(x^{p})=b^{p}\log(1+\exp(x^{p}))-c^{p}x^{p}$ is the cost of
user $p$ for transmission, $b^{p}>0,c^{p}\geq 0$. The partial-gradient $U^{p}$
is calculated to be,

$\textstyle
U^{p}(x)=\dfrac{a^{p}\exp(x^{p})}{1+\sum_{r\in\mathcal{N}}a^{r}\exp(x^{r})}-\dfrac{b^{p}\exp(x^{p})}{1+\exp(x^{p})}+c^{p}.$
(42)

The second-order partial derivative of $\mathcal{U}^{p}$ can be calculated to
be,

$\dfrac{\partial^{2}\mathcal{U}^{p}(x^{p};x^{-p})}{\partial x^{q}\partial
x^{p}}=\begin{cases}\dfrac{-a^{p}a^{q}\exp(x^{p}+x^{q})}{(1+\sum_{r\in\mathcal{N}}a^{r}\exp(x^{r}))^{2}}&p\neq
q\\\\[12.91663pt] \begin{aligned} &\dfrac{a^{p}\exp(x^{p})(1+\sum_{r\neq
p}a^{r}\exp(x^{r}))}{(1+\sum_{r\in\mathcal{N}}a^{r}\exp(x^{r}))^{2}}\\\
&-\dfrac{b^{p}\exp(x^{p})}{(1+\exp(x^{p}))^{2}}\end{aligned}&p=q\end{cases}$
(43)

which generates a symmetric Jacobian $\operatorname{\mathbf{J}}_{U}(x),\forall
x$,

${\mathbf{J}_{U}(x)}_{p,p}=\dfrac{a^{p}\exp(x^{p})(1+\sum_{r\neq
p}a^{r}\exp(x^{r}))}{(1+\sum_{r\in\mathcal{N}}a^{r}\exp(x^{r}))^{2}}-\dfrac{b^{p}\exp(x^{p})}{(1+\exp(x^{p}))^{2}}$
(44)

and

${\mathbf{J}_{U}(x)}_{p,q}=\dfrac{-a^{p}a^{q}\exp(x^{p}+x^{q})}{(1+\sum_{r\in\mathcal{N}}a^{r}\exp(x^{r}))^{2}},p\neq
q,$ (45)

hence the game is a potential game with a non-concave potential function $P$
such that $\nabla P=U$.

Consider an example with $N=2$ and problem parameters,
$a=(a^{1},a^{2})=(1,1),b=(b^{1},b^{2})=(4,4),c=(c^{1},c^{2})=(3,3)$. For all
players, $Z^{p}(0)$ is sampled uniformly from $[0,10]$, $\Xi^{p}(0)=0$. The NE
of this game can be found at $x^{\star}=(1.8663,1.8663)$. The Jacobian at the
NE is
$\operatorname{\mathbf{J}}(x^{\star})=\operatorname{diag}(\begin{bmatrix}-4.308&0\end{bmatrix}^{\top})$,
hence by 1, the NE is locally mere VSS. By 4, DA2 converges towards this NE.
For each of the following experiments, we simulate for $T=10^{4}$ steps and
report the strategies at termination. We restrict $\Omega^{p}$ to be a large
but compact set $[-1000,1000]^{2}$ and set $C^{p}_{\epsilon}$ to be the
projection operator,
$C^{p}_{\epsilon}(Z^{p})=\text{argmin}_{y^{p}\in\Omega^{p}}\|\epsilon^{-1}Z^{p}-y^{p}\|_{2}^{2}$
with $\epsilon=1$. Unless specified otherwise, all other parameters
$\gamma,\alpha,\beta$ are kept as $1$. All additive noise are zero-mean
Gaussian with the same variance $\sigma^{2}_{\zeta}$.

We start our experiment with a small variance $\sigma^{2}_{\zeta}=0.1$ and set
the step-sizes to be $t_{k}=0.39/k^{0.26},\tau_{k}=0.12/k^{0.64}$. Figure 6
shows that DA2 converges to the NE. We increase the variance to
$\sigma^{2}_{\zeta}=10$ and Figure 7 shows similar observation.

Figure 6: $\sigma^{2}_{\zeta}=0.1$,
$t_{k}=0.39/k^{0.26},\tau_{k}=0.12/k^{0.64}$ Figure 7:
$\sigma^{2}_{\zeta}=10$, $t_{k}=0.38/k^{0.47},\tau_{k}=0.13/k^{0.71}$

Next, we investigate some possible computational issues can arise when we
increase the number of players in this game. Consider an example with $N=10$
players. The game parameters are:

$\displaystyle a\\!=(12,\\!4,\\!2,\\!5,\\!20,\\!20,\\!12,\\!3,\\!1,\\!2),$
$\displaystyle
b\\!=(15,\\!20,\\!20,\\!14,\\!21,\\!20,\\!16,\\!19,\\!17,\\!17),$
$\displaystyle c\\!=(13,\\!12,\\!14,\\!8,\\!10,\\!19,\\!10,\\!12,\\!15,\\!1).$

A variance of $\sigma^{2}_{\zeta}=1$ is applied for the following experiments.
We employ step-size sequences
$t_{k}=\nicefrac{{0.4}}{{k^{0.033}}},\tau_{k}=\nicefrac{{0.78}}{{k^{0.001}}}$
across all players. The NE of this game is located at
$x^{\star}=(1.87,0.41,0.85,0.28,-0.10,16,0.51,0.54,2.01,-2.77)$ (rounded to
two decimal places) and was confirmed to be a nearly mere VSS (13). We note
that an exact mere VSS is difficult to produce for this game when there are a
large number of players.

Figure 8: $\sigma^{2}_{\zeta}=1$,
$t_{k}=0.40/k^{0.033},\tau_{k}=0.78/k^{0.001}$ ($X^{2}$ omitted).

Figure 8 shows that the strategy $X^{1}$ converges to
$(1.87,0.41,0.85,0.28,-0.10,11.29,0.51,0.54,2.01,-2.77)$ which is exactly the
same as $x^{\star}$ in all entries except for the sixth entry ($16.51$ vs
$11.29$), which is likely due to the small rate parameters employed in the
step-sizes. This highlights the difficulty of employing the same step-size
sequences across all players. We adjust all the parameters $\beta,\alpha$ and
step-sizes $ck^{-\iota},\iota\in(0,1),c>0$ on a per-player basis, which
resulted in a much closer convergence to the NE as opposed to applying them
across all players uniformly (Figure 9).

Figure 9: $\sigma^{2}_{\zeta}=1$, per-player tuned parameters ($X^{2}$
omitted).

######  Example 12.

(Learning to Generate a Gaussian)

Let $\mathpzc{Z}\sim\mathpzc{P}(\mathpzc{z})$ and
$\mathpzc{X}\sim\mathpzc{Q}(\mathpzc{x})$ be two random variables. We wish to
construct a model $\mathpzc{G}_{\theta}:\mathbb{R}{r}\to\mathbb{R}{m}$,
$\mathpzc{r},\mathpzc{m}\geq 1$, with an unknown, continuous parameter
$\theta$ such that $\mathpzc{G}_{\theta}(\mathpzc{Z})$ recovers the statistics
(mean, variance, etc.) of $\mathpzc{X}$. Following [54],
$\mathpzc{G}_{\theta}$ can be constructed through solving,

${\text{min}_{\theta}}\thinspace{\text{max}_{w}}\thinspace\mathbb{E}_{\mathpzc{X}}(\mathpzc{D}_{w}(\mathpzc{X}))-\mathbb{E}_{\mathpzc{Z}}(\mathpzc{D}_{w}(\mathpzc{G}_{\theta}(\mathpzc{Z}))),$
($\star$)

where the model $\mathpzc{D}_{w}:\mathbb{R}{m}\to\mathbb{R}$ is parametrized
by an unknown continuous parameter $w$. This problem can be thought of as a
game between the owners of the models $\mathpzc{G}_{\theta}$ and
$\mathpzc{D}_{w}$ whereby $\theta,w$ are their respective strategies. The
owner of $\mathpzc{G}_{\theta}$ (or the designer) varies $\theta$ and uses
$G_{\theta}(\mathpzc{Z})$ to estimate the statistics of $\mathpzc{X}$, whereas
the owner of $\mathpzc{D}_{w}$ (or the tester) varies $w$ in an attempt to
incur the largest penalty possible for the discrepancy between
$G_{\theta}(\mathpzc{Z})$ and $\mathpzc{X}$.

Figure 10: Diagram representation of $\mathpzc{G}_{\theta}$ and
$\mathpzc{D}_{w}$ models

We construct a more complicated variant of the examples in [54, 29], whereby
we assume that the $\mathpzc{G}_{\theta}$ owner wishes to learn the mean and
variance of a one-dimensional Gaussian at the same time: let
$\mathpzc{Z}\sim\mathpzc{N}(0,1)$, and
$\mathpzc{X}\sim\mathpzc{N}(v,\sigma^{2})$, where
$v\in\mathbb{R},\sigma\in\mathbb{R}_{>0}$ are the mean and standard deviation
of a Gaussian distribution. Let
$\mathpzc{D}_{w}(\mathpzc{X})=w_{1}\mathpzc{X}^{2}+w_{2}\mathpzc{X},w=(w_{1},w_{2})\in\mathbb{R}\times\mathbb{R}$
and
$\mathpzc{G}_{\theta}(\mathpzc{Z})=\theta_{1}\mathpzc{Z}+\theta_{2},\theta=(\theta_{1},\theta_{2})\in\mathbb{R}\times\mathbb{R}$,
then a brief calculation shows that

$(\star)=w_{1}((\sigma^{2}+v^{2})-(\theta^{2}_{1}+\theta_{2}^{2}))+w_{2}(v-\theta_{2}).$

Let $x^{1}=\theta,x^{2}=w$, we obtain a two-player ZS game,

$\textstyle\mathcal{U}^{1}(x^{1};x^{2})=x^{2}_{1}(\sigma^{2}+v^{2}-\sum_{i=1}^{2}(x_{i}^{1})^{2})+x^{2}_{2}(v-x_{2}^{1})$
(46)

and $\mathcal{U}^{2}(x^{2};x^{1})=-\mathcal{U}^{1}(x^{1};x^{2})$ where the
action sets are, $\Omega^{1}=\mathbb{R}\times\mathbb{R}=\Omega^{2}$. The
pseudo-gradient is,
$U(x)=(-2{x^{2}_{1}}x^{1}_{1},-2{x^{2}_{1}}x^{1}_{2}-{x^{2}_{2}},\sum_{i=1}^{2}(x_{i}^{1})^{2}-(\sigma^{2}+v^{2}),x^{1}_{2}-v),$
which implies an interior NE at
$x^{\star}=({x^{1}}^{\star},{x^{2}}^{\star})=((\sigma,v),(0,0))$. The Jacobian
of $U$ is shown to be,

$\mathbf{J}_{U}(x)=\begin{bmatrix}-2x_{1}^{2}&0&-2x^{1}_{1}&0\\\
0&-2x^{2}_{1}&-2x^{1}_{2}&-1\\\ 2x^{1}_{1}&2x^{1}_{2}&0&0\\\
0&1&0&0\end{bmatrix}.$ (47)

Since
$y^{\top}\operatorname{\boldsymbol{\mathfrak{J}}}_{U}(x)y=\frac{1}{2}y^{\top}(\mathbf{J}_{U}(x)+\mathbf{J}^{\top}_{U}(x))y=-2((y_{1})^{2}+(y_{2})^{2})x^{2}_{1}$
and $x^{2}_{1}$ is not restricted to be positive, therefore, this game is not
monotone [51, Prop. 2.3.2]. However, at the interior NE,
${x^{2}}^{\star}=({x_{1}^{2}}^{\star},{x_{2}^{2}}^{\star})=(0,0)\implies{x_{1}^{2}}^{\star}=0$,
and hence
$y^{\top}\operatorname{\boldsymbol{\mathfrak{J}}}_{U}(x^{\star})y=\frac{1}{2}y^{\top}(\mathbf{J}_{U}(x^{\star})+\mathbf{J}^{\top}_{U}(x^{\star}))y=-2((y_{1})^{2}+(y_{2})^{2})(0)=0$.
By 1(ii’), this implies that $x^{\star}$ is locally merely VS.

Consider an example with $v=10$, $\sigma=5$. We perform two simulations of the
discrete-time MD2 with or without perturbation on the pseudo-gradient. We
assume $Z(0)$ is uniformly sampled as follows:
$Z^{1}(0)\sim\text{unif}([10,50]\times[20,40])$,
$Z^{2}(0)\sim\text{unif}([0,5]\times[0,5])$, and $\Xi(0)=\mathbf{0}$. We
restrict $\Omega^{p}$ to be a compact set $[-1000,1000]^{2}$ and set
$C^{p}_{\epsilon}$ to be the projection operator as in 11. We set
$\lambda=\beta=\alpha=1$. Each of the following simulations is ran for
$T=4\times 10^{4}$ steps, with initial conditions randomly sampled as stated
above. Note that our simulations has experimentally shown to accommodate
larger sets of initial conditions than the ones provided above.

First, we consider the noiseless case, i.e., $\sigma^{2}_{\zeta}=0$. From
Figure 11, we see that $X^{1}$ converges to $X^{1\star}$ (highlighted using a
green star).

Figure 11: $\sigma^{2}_{\zeta}=0$, $t_{k}=\dfrac{0.0108}{k^{0.3971}}$,
$\tau_{k}=\dfrac{0.0485}{k^{0.7877}}$ ($x^{2}$ omitted).

Next, we consider the harder case where $U^{p}$ is subjected to zero-mean
Gaussian noise with a very large variance $\sigma^{2}_{\zeta}=30^{2}$ for each
$p$. The trajectories of $X^{1}$ are shown in Figure 12. Our previous
observation continues to hold.

Figure 12: $\sigma^{2}_{\zeta}=30^{2}$, $t_{k}=\dfrac{0.0108}{k^{0.3971}}$,
$\tau_{k}=\dfrac{0.0485}{k^{0.7877}}$ ($X^{2}$ omitted).

##  9 Conclusions

In this paper, we showed that MD2 overcomes many shortcomings of the first-
order MD as well as other members of this class of dynamics, such as MD with
time-averaging or discounting. To summarize our findings: first, we showed
that the trajectories of MD2 converge exactly to an interior mere (and not
necessarily strictly) VSS (1). As such, it converges to interior mere NEs of
monotone and pseudo-monotone games, and interior strict NE of quasi-monotone
games (1). We also found that MD2 is robust around interior mere VSS (6, 13,
10). Next, we provided a modification to MD2 based on vanishing perturbation
which results in an exponential rate of convergence towards a strongly stable
VSS (2). We then showed MD2 can achieve no-regret while converging beyond a
strict VSS, such as to an interior/fully-mixed NEs of zero-sum finite games
(3). This result improves our current understandings between no-regret
learning and convergence [22, 21]. We then provided several ways of deriving
primal-space dynamics from MD2. Finally, we moved to the scenario where the
game evolves in discrete-time and the pseudo-gradient is subjected to noise
perturbation. Using stochastic approximation techniques, we were able to
relate the limiting behavior of MD2 with that of its discrete-time, noisy
counterpart, and showed that the discrete MD2 with noisy observations or
second-order dual averaging (DA2) can converge to the interior mere VSS under
standard regularity assumptions (4).

There are several outstanding issues that could be better understood. First,
our work did not provide a thorough convergence proof for local VSS, which are
dealt with appropriate initialization. Moreover, we did not address
convergence towards boundary point for MD2, as our proofs relied on an
interiority condition associated with interior NE (14). One possible approach
for dealing with these boundary solutions is to adopt the approach of [19], by
restricting ourselves to the context of finite games, where these boundary
solutions are quite relevant. Another basic issue is that it is still not
entirely clear to us at this time how MD2 relates to other existing
continuous-time dynamics for optimization or games. Aside from the cases that
we know of, e.g., [27], it is possible that there are other existing
algorithms that appear as specific instances of MD2.

Our work opens up an extensive array of directions on the interface between
dynamical systems and games. As a starting point, in the continuous-time, an
interesting direction would be to explain the difference between MD and MD2 in
zero-sum games through a volume compressibility perspective as in [23, 16]. We
can also try to craft third or even higher-order versions of MD based on the
technique in [28]. In the discrete-time case, we can further reduce the
information, to that of the zero-order information case (or full-bandit
feedback). Moreover, we can tackle the case whereby the game itself has
parameters that are time-varying. Finally, it is also worthwhile to examine
MD2 in a discrete-time setting whereby the mirror map is replaced with a
proximal operator.

##  10 Appendix

Symbol |  Definition
---|---
$x/x(t)/X_{k}$ |  Strategy/at time $t$/at time $k$
$z/z(t)/Z_{k}$ |  Dual aggregate/at time $t$/at time $k$
$\xi/\xi(t)/\Xi_{k}$ |  Primal aggregate/at time $t$/at time $k$
$\mathcal{U}^{p}/U^{p}/U$ |  Payoff/partial-gradient/pseudo-gradient
$\mathbf{J}_{U}/\operatorname{\boldsymbol{\mathfrak{J}}}_{U}$ |  Jacobian of $U$/Symmetric game Jacobian of $U$
$\vartheta^{p}/\psi^{p}_{\epsilon}/\psi^{p\star}_{\epsilon}$ |  Regularizer/$\psi^{p}_{\epsilon}=\epsilon\vartheta^{p}$/convex conjugate of $\psi^{p}_{\epsilon}$
$C^{p}_{\epsilon}/C_{\epsilon}/\mathbf{J}_{C_{\epsilon}}$ |  Player $p$’s/overall mirror map/Jacobian of $C_{\epsilon}$
$\gamma,\epsilon,\alpha,\beta$ |  Parameters associated with MD/MD2/DA2
$\mathpzc{I}/\mathpzc{O}$ |  Identity matrix/zero matrix
Table 1: A list of main notations used in this paper.

The following two lemmas are standard in this area of literature, see [19, 23,
29]. In particular, 2 follows directly from the Legendre theorem [53].

######  Lemma 2.

Let $\psi^{p}_{\epsilon}\\!\\!\coloneqq\\!\epsilon\vartheta^{p}\\!$,
$\epsilon\\!>\\!\\!0$, where $\vartheta^{p}$ satisfies 3(ii), and let
${\psi^{p}_{\epsilon}}^{\star}$ be the convex conjugate of
$\psi^{p}_{\epsilon}$. Then,

  * (i)

${\psi^{p}_{\epsilon}}^{\star}\\!:\\!\mathbb{R}^{n_{p}}\\!\to\\!\mathbb{R}\\!\cup\\!\\{\infty\\}$
is closed, proper, strictly convex, steep, and finite-valued over
$\mathbb{R}^{n_{p}}$.

  * (ii)

$\nabla\psi^{p}_{\epsilon}\\!:\\!\operatorname{int}(\operatorname{dom}(\psi^{p}_{\epsilon}))\\!\to\\!\operatorname{int}(\operatorname{dom}({\psi^{p}_{\epsilon}}^{\star}))$
is a homeomorphism with inverse mapping
$(\nabla\psi^{p}_{\epsilon})^{-1}\\!=\\!\nabla{\psi^{p}_{\epsilon}}^{\star}\\!=\\!C^{p}_{\epsilon}.$

  * (iii)

$C^{p}_{\epsilon}$ is strictly monotone on
$\operatorname{int}(\operatorname{dom}({\psi^{p}_{\epsilon}}^{\star}))$.

######  Lemma 3.

Let $\psi^{p}_{\epsilon}\\!\\!\coloneqq\\!\epsilon\vartheta^{p}\\!$,
$\epsilon\\!>\\!\\!0$, where $\vartheta^{p}$ satisfies 3(i), and let
${\psi^{p}_{\epsilon}}^{\star}$ be the convex conjugate of
$\psi^{p}_{\epsilon}$. Then,

  1. (i)

${\psi^{p}_{\epsilon}}^{\star}\\!:\\!\mathbb{R}^{n_{p}}\\!\to\\!\mathbb{R}\cup\\{\infty\\}$
is closed, proper, convex and finite-valued over $\mathbb{R}^{n_{p}}$, i.e.,
$\operatorname{dom}({\psi^{p}_{\epsilon}}^{\star})=\mathbb{R}^{n_{p}}$.

  2. (ii)

${\psi^{p}_{\epsilon}}^{\star}$ is continuously differentiable and
$\nabla{\psi^{p}_{\epsilon}}^{\star}\\!=\\!C^{p}_{\epsilon}$.

  3. (iii)

$C^{p}_{\epsilon}\\!$ is surjective from $\mathbb{R}^{n_{p}}\\!$ onto
$\operatorname{rint}(\Omega^{p})$ whenever $\psi^{p}_{\epsilon}\\!$ is steep,
and onto $\Omega^{p}$ whenever $\psi^{p}_{\epsilon}\\!$ is non-steep.

###### Proof.

(Proof of 1) Let $x,x^{\prime}\in\Omega$ be two arbitrarily chosen strategies.
Consider $\overline{x}=\theta x+(1-\theta)x^{\prime},\theta\in[0,1]$. Then the
proof for all these claims boil down to the equality condition found in [60,
Prop. 1.44], in which it be can shown,

$\textstyle(x-x^{\prime})^{\top}(U(x)-U(x^{\prime}))=\textstyle\int_{0}^{1}(x^{\prime}-x)^{\top}\operatorname{\boldsymbol{\mathfrak{J}}}_{U}(\overline{x})(x^{\prime}-x)\,\,\mathrm{d}\theta$
(48)

Rearranging (48), we have,

$\begin{split}\textstyle(x-x^{\prime})^{\top}U(x)=&\textstyle\int_{0}^{1}(x^{\prime}-x)^{\top}\operatorname{\boldsymbol{\mathfrak{J}}}_{U}(\overline{x})(x^{\prime}-x)\,\,\mathrm{d}\theta\\\
&+(x-x^{\prime})^{\top}U(x^{\prime}).\end{split}$ (49)

Let $x^{\prime}=x^{\star}$, then $(x-x^{\star})^{\top}U(x^{\star})\leq 0$ and
we have,

$\textstyle(x-x^{\star})^{\top}U(x)\leq\int_{0}^{1}(x^{\star}-x)^{\top}\operatorname{\boldsymbol{\mathfrak{J}}}_{U}(\overline{x})(x^{\star}-x)\,\,\mathrm{d}\theta.$
(50)

By the definiteness assumptions of
$\operatorname{\boldsymbol{\mathfrak{J}}}_{U}(x)$ on $T_{\Omega}(x)$ for all
$x\in\Omega$ and $\overline{x}=\theta x+(1-\theta)x^{\star}$ is contained in
$\Omega$ for any $\theta\in[0,1]$, (50) implies statements $(i),(ii),(iii)$.

Next, suppose the definiteness condition on
$\operatorname{\boldsymbol{\mathfrak{J}}}_{U}(x^{\star})$ holds for all $y\in
T_{\Omega}(x^{\star})$. Since $U$ is assumed to be continuously
differentiable, therefore
$y^{\top}\operatorname{\boldsymbol{\mathfrak{J}}}_{U}(x)y$ is continuous for
all $x\in\Omega$. Using the standard $\epsilon-\delta$ definition of
continuity, this means for every $\epsilon>0$ there exists a
$\delta(\epsilon)>0$ such that for all $x\in\Omega$,
$\|x-x^{\star}\|<\delta(\epsilon)\implies|y^{\top}\operatorname{\boldsymbol{\mathfrak{J}}}_{U}(x)y-y^{\top}\operatorname{\boldsymbol{\mathfrak{J}}}_{U}(x^{\star})y|<\epsilon$.
The latter statement,
$|y^{\top}\operatorname{\boldsymbol{\mathfrak{J}}}_{U}(x)y-y^{\top}\operatorname{\boldsymbol{\mathfrak{J}}}_{U}(x^{\star})y|<\epsilon$,
implies
$y^{\top}\operatorname{\boldsymbol{\mathfrak{J}}}_{U}(x)y<\epsilon+y^{\top}\operatorname{\boldsymbol{\mathfrak{J}}}_{U}(x^{\star})y$.
Since this holds for every $\epsilon>0$, therefore this is equivalent to
$y^{\top}\operatorname{\boldsymbol{\mathfrak{J}}}_{U}(x)y\leq
y^{\top}\operatorname{\boldsymbol{\mathfrak{J}}}_{U}(x^{\star})y$, and hence
$\operatorname{\boldsymbol{\mathfrak{J}}}_{U}(x)$ shares the same definiteness
property as $\operatorname{\boldsymbol{\mathfrak{J}}}_{U}(x^{\star})$ for all
$x$ in some ball
$\mathcal{D}=\\{x\in\Omega|\|x-x^{\star}\|<\delta(\epsilon)\\}$. Let
$\overline{x}=\theta x+(1-\theta)x^{\star}$, $\theta\in[0,1]$,
$x\in\mathcal{D}$. This implies (50) holds locally on the set of all
$\overline{x}$, and $(i^{\prime}),(ii^{\prime}),(iii^{\prime})$ follows. The
isolation condition of $(i),(i^{\prime})$ are proven in [2, Prop. 2.7]. ∎

###### Proof.

(Proof of 1) (i) Consider the Lyapunov function,

$\begin{split}\textstyle
V(\xi,z)=&\dfrac{\alpha}{2\beta}\|\xi(t)-x^{\star}\|_{2}^{2}\\\
&+\textstyle\gamma^{-1}\sum_{p\in\mathcal{N}}\psi^{p}_{\epsilon}({x^{p}}^{\star})-{z^{p}}^{\top}{x^{p}}^{\star}+\psi_{\epsilon}^{p\star}(z^{p}).\end{split}$
(51)

Here, $V$ is composed of a quadratic distance term, which measures the
progress of $\xi(t)$ towards $\xi^{\star}=x^{\star}$, added onto a function
which consists the collection of all the terms associated with the Fenchel-
Young inequality[49, p. 88], also known as the Fenchel coupling in [2]. Since
$\vartheta^{p}$ satisfies 3(i), $V(\xi,z)$ is continuous, positive definite
(due to the Legendre theorem), $V(\xi,z)=0$ iff $\xi=x^{\star}$ and
$z=C^{-1}_{\epsilon}(x^{\star})$. Taking the time-derivative of $V(\xi,z)$
along the solutions of MD2, using $x=C_{\epsilon}(z)$,
$\nabla{\psi^{p}_{\epsilon}}^{\star}=C^{p}_{\epsilon}$ and
$\dot{z}=\gamma(U(x)-\alpha\beta^{-1}\dot{\xi})$,

$\displaystyle\dot{V}(\xi,z)$
$\displaystyle=\\!\dfrac{\alpha}{\beta}(\xi-x^{\star})^{\top}\dot{\xi}+\gamma^{-1}(C_{\epsilon}(z)-x^{\star})^{\top}\dot{z}$
$\displaystyle=\\!\dfrac{\alpha}{\beta}(\xi-x^{\star})^{\top}\dot{\xi}+(C_{\epsilon}(z)-x^{\star})^{\top}(U(x)-\dfrac{\alpha}{\beta}\dot{\xi})$
$\displaystyle=\\!\dfrac{\alpha}{\beta}(\xi-x^{\star}-x+x^{\star})^{\top}\dot{\xi}+(C_{\epsilon}(z)-x^{\star})^{\top}U(x)$
$\displaystyle=\\!\dfrac{\alpha}{\beta}(\xi-x)^{\top}\dot{\xi}+(x-x^{\star})^{\top}U(x).$

Substituting in $\dot{\xi}=\beta(x-\xi)$, we have,

$\displaystyle\dot{V}(\xi,z)$
$\displaystyle=-\alpha\beta^{-2}\|\dot{\xi}\|^{2}_{2}+(x-x^{\star})^{\top}U(x)$
(52) $\displaystyle\leq-\alpha\beta^{-2}\|\dot{\xi}\|^{2}_{2}\leq 0,$

where we have used $x^{\star}$ is a mere VSS. Observe that $\dot{V}(\xi,z)=0$
only if $\dot{\xi}=\mathbf{0}\Longleftrightarrow x=\xi$. By the Legendre
theorem [53], ${\psi^{p}_{\epsilon}}^{\star}$ is coercive, which implies
$V(\xi,z)$ is coercive (radially unbounded) on
$\mathbb{R}^{n}\times\mathbb{R}^{n}$ and hence all of its sublevel sets are
compact. Let
$\mathcal{D}_{c}=\\{(\xi,z)\in\mathbb{R}^{n}\times\mathbb{R}^{n}|V(\xi,z)\leq
c\\}$ be a compact sublevel set for some $c>0$. Consider
$\mathcal{E}=\\{(\xi,z)\in\mathcal{D}_{c}|\dot{V}(\xi,z)=0\\}=\\{(\xi,z)\in\mathcal{D}_{c}|\xi=x=C_{\epsilon}(z)\\}$.
Let $(\xi(t),z(t))$ be some solution starting in $\mathcal{E}$, then
$\xi=x=C_{\epsilon}(z)\implies\dot{\xi}=\dot{x}=\operatorname{\mathbf{J}}_{C_{\epsilon}}(z)\dot{z}\implies\mathbf{0}=\operatorname{\mathbf{J}}_{C_{\epsilon}}(z)\dot{z}$.
By Legendre theorem and the $\mathcal{C}^{2}$ assumption,
$\operatorname{\mathbf{J}}_{C_{\epsilon}}(z)$ exists and is invertible,
therefore
$\dot{z}=\operatorname{\mathbf{J}}_{C_{\epsilon}}(z)^{-1}\mathbf{0}=\mathbf{0}$.
This means the largest invariant set contained in $\mathcal{E}$ is
$\mathcal{S}=\\{(\xi,z)\in\mathcal{E}|\dot{z}=\mathbf{0}\\}$. By LaSalle’s
invariance principle [43, Theorem 3.3], any solution starting from
$\mathcal{D}_{c}$ converges to $\mathcal{S}$ as $t\to\infty$. On
$\mathcal{S}$, $\dot{z}=\mathbf{0}\implies
U(x)=U(C_{\epsilon}(z))=\mathbf{0}$, which means $x=C_{\epsilon}(z)$ converges
to some interior NE (which we may denote as $x^{\star}$, since $x^{\star}$ is
arbitrary), which by assumption is a globally mere VSS. Global convergence
follows from coercivity of $V(\xi,z)$.

(ii) Suppose instead $\Omega^{p}$ is compact for all players $p$. Let
$x^{\omega}$ be an $\omega$-limit point of $x(t)=C_{\epsilon}(z(t))$ of MD2.
Since $x(t)$ is bounded for all $t\geq 0$, the existence of $x^{\omega}$
follows from Bolzano-Weierstrass theorem. Suppose that $x^{\omega}$ is not a
globally mere VSS $x^{\star}$. Take $\mathcal{O}_{\omega}$ be an open ball
around $x^{\omega}$. By definition, there exists a sequence
$\\{x(t_{k})\\}_{k\in\mathbb{N}}$, $x(t_{k})\in\mathcal{O}_{\omega}$,
converging towards $x^{\omega}$, where $\\{t_{k}\\}_{k\in\mathbb{N}}$ is an
increasing sequence of times. Following the technique in [20], we wish to
build an auxiliary sequence in $\mathcal{O}_{\omega}$ and show that $V$ is
unbounded below along this sequence, thereby obtaining a contradiction.

First, note that since $\vartheta^{p}$ is $\rho$-strongly convex,
$\psi^{p}_{\epsilon}\coloneqq\epsilon\vartheta^{p}$ is $\epsilon\rho$-strong
convex, and hence by the conjugate correspondence theorem [49, Theorem 5.26],
${\psi^{p}_{\epsilon}}^{\star}$ is $(\epsilon\rho)^{-1}$-smooth, i.e.,
$C^{p}_{\epsilon}=\nabla{\psi^{p}_{\epsilon}}^{\star}$ is
$(\epsilon\rho)^{-1}$-Lipschitz. Let $\tau>0$ be a time increment, then, by
$z(t_{k}+\tau)-z(t_{k})=\int_{t_{k}}^{t_{k}+\tau}\dot{z}(s)\mathrm{d}s=\gamma\int_{t_{k}}^{t_{k}+\tau}U(x(s))-\alpha(x(s)-\xi(s))\mathrm{d}s$,
we have,

$\displaystyle\|x(t_{k}+\tau)-x(t_{k})\|=\|C_{\epsilon}(z(t_{k}+\tau)-C_{\epsilon}(z(t_{k}))\|$
(53) $\displaystyle\leq(\epsilon\rho)^{-1}\|z(t_{k}+\tau)-z(t_{k})\|_{\star}$
(54)
$\displaystyle\textstyle\leq(\epsilon\rho)^{-1}\gamma\|\smallint_{t_{k}}^{t_{k}+\tau}U(x(s))-\alpha(x(s)-\xi(s))\mathrm{d}s\|_{\star}$
(55)
$\displaystyle\textstyle\leq(\epsilon\rho)^{-1}\gamma\smallint_{t_{k}}^{t_{k}+\tau}\|U(x(s)\|_{\star}+\alpha\|x(s)-\xi(s)\|_{\star}\mathrm{d}s$
(56)
$\displaystyle\leq(\epsilon\rho)^{-1}\tau\gamma\max_{x\in\Omega}(\|U(x)\|_{\star}+\alpha\|x(s)-\xi(s)\|_{\star})\leq
b(\tau),$ (57)

where $b(\tau)$ is a $\tau$-dependent upper-bound. Let
$\mathcal{B}_{\tau}(x(t_{k}))=\\{x\in\Omega|\|x-x(t_{k})\|\leq
b(\tau),x(t_{k})\in\mathcal{O}_{\omega}\\}$ be a ball around $x(t_{k})$. Since
$\mathcal{O}_{\omega}$ is open, therefore there exists some $\delta>0$, such
that $\mathcal{B}_{\tau}(x(t_{k}))\subseteq\mathcal{O}_{\omega}$ for all
$\tau\in[0,\delta]$ and hence $x(t_{k}+\tau)\in\mathcal{O}_{\omega}$.

Integrating (52), we have,

$\displaystyle V(\xi(t),z(t))$
$\displaystyle\textstyle=V_{0}-\smallint_{0}^{t}\alpha\|x-\xi\|^{2}_{2}+(x-x^{\star})^{\top}U(x)\mathrm{d}\tau,$
(58)

where $V_{0}\coloneqq V(\xi(0),z(0))$ is some constant. Since $x^{\star}$ is
assumed to be a globally mere VSS, therefore $(x-x^{\star})^{\top}U(x)\leq
0,\forall x$. Tossing out $(x-x^{\star})^{\top}U(x)\leq 0,\forall x$, (58), we
obtain,

$\displaystyle V(\xi(t),z(t))$ $\displaystyle\textstyle\leq
V_{0}-\smallint_{0}^{t}\alpha\|x(\tau)-\xi(\tau)\|^{2}_{2}\mathrm{d}\tau.$

Expressing this inequality in terms of the sequence
$\\{x(t_{k})\\}_{k\in\mathbb{N}}$ and $\\{\xi(t_{k})\\}_{k\in\mathbb{N}}$,
$-\|x(\tau)-\xi(\tau)\|_{2}\leq 0$, and using the $\delta$ which we have
previous found, we have,

$\displaystyle V(\xi(t_{k}+\delta),z(t_{k}+\delta))$
$\displaystyle\textstyle\leq
V_{0}-\sum_{\ell=1}^{k}\smallint_{t_{\ell}}^{t_{\ell}+\delta}\alpha\|x(\tau)-\xi(\tau)\|^{2}_{2}\mathrm{d}\tau$

for some $k$. Since $\Omega$ is compact, $x$ is continuous, $\xi$ is a
continuous function of $x$ (15), and $\xi$ is not eventually identically equal
to $x$ (otherwise, using the rest point condition for MD2 (14) with our
assumption that all interior NE are globally mere, we arrive at
$x^{\omega}=x^{\star}$ hence a contradiction), $\|x-\xi\|_{2}^{2}$ achieves a
lower-bound on an interval $[t_{\ell},t_{\ell}+\delta]$. Let
$c=\max\\{l>0|\|x(\tau)-\xi(\tau)\|_{2}^{2}>l,\tau\in[t_{\ell},t_{\ell}+\delta],\ell=1,\ldots,k\\}$.
This means,

$\displaystyle V(\xi(t_{k}+\delta),z(t_{k}+\delta))\textstyle\leq
V_{0}-\sum_{\ell=1}^{k}\smallint_{t_{\ell}}^{t_{\ell}+\delta}\alpha
c\mathrm{d}\tau\leq V_{0}-\alpha c\delta k,$

which shows that $V$ is unbounded below as $k\to\infty$, a contradiction. ∎

######  Remark 13.

(Exact convergence towards globally weak VSS) Consider (52) in the proof of 1,

$\displaystyle\dot{V}(\xi,z)$
$\displaystyle=-\alpha\beta^{-2}\|\dot{\xi}\|^{2}_{2}+(x-x^{\star})^{\top}U(x).$

Suppose that $x^{\star}$ is a globally weak VSS, then for fixed
$\alpha,\beta,\mu>0$,

$\displaystyle\dot{V}(\xi,z)$
$\displaystyle\leq-\alpha\beta^{-2}\|\dot{\xi}\|^{2}_{2}+\mu\|x-x^{\star}\|_{2}^{2}.$

Observe that the right-hand side of this inequality is negative as long as
$\alpha\beta^{-2}\|\dot{\xi}\|^{2}_{2}>\mu\|x-x^{\star}\|_{2}^{2}$, or
equivalently, $\alpha\|x-\xi\|^{2}_{2}>\mu\|x-x^{\star}\|_{2}^{2}$ for all $t$
prior to convergence. Suppose this preceding inequality holds, then MD2 can
converge to equilibria that are weaker than mere VSS. This is likely to happen
when one or more of the following holds: (i) $\alpha$ is large, (ii) the
difference between $x,\xi$ is kept reasonably large, which can be achieved by
decreasing $\beta$, or (iii) $\mu$ is small. We can informally refer
$\mu$-weak VSS with a small $\mu$ as a type of nearly mere VSS, i.e., a weak
VSS within a small distance away from becoming a mere VSS. The class of nearly
mere VSS is practically significant, as in realistic applications, it is often
unlikely for a NE to be an exactly mere VSS. For instance, by 1, the mere VSS
condition amounts to requiring the largest eigenvalue of the game Jacobian at
the NE, $\operatorname{\boldsymbol{\mathfrak{J}}}_{U}(x^{\star})$, to be
exactly located at $0$, which is often difficult to achieve due to numerical
precision, especially in the cases where there are many players in the game
(11). The preceding observation also shows that MD2 enjoys a certain
robustness property as compared to MD: for MD, if the problem is perturbed so
that $x^{\star}$ goes from strict to mere, then MD can cease to converge,
whereas for MD2, when $x^{\star}$ goes from mere to weak, MD2 can still
converge.

######  Lemma 4.

Suppose $\vartheta^{p}:\mathbb{R}^{n_{p}}\to\mathbb{R}\cup\\{\infty\\}$
satisfies 3(i) or (ii). Let $h=\sum_{p\in\mathcal{N}}\vartheta^{p}$,
${x^{p}}^{\star}\in\operatorname{dom}(\vartheta^{p}),x^{p}\in\operatorname{dom}(\partial\vartheta^{p}),z^{p}\in\partial\psi_{\epsilon}^{p}(x^{p}),\psi_{\epsilon}^{p}=\epsilon\vartheta^{p},\forall
p$, then,

$\textstyle
D_{h}(x^{\star},x)=\epsilon^{-1}\sum_{p\in\mathcal{N}}\psi_{\epsilon}^{p}({x^{p}}^{\star})\\!-\\!{\psi_{\epsilon}^{p}}^{\star}(z^{p})\\!-\\!{z^{p}}^{\top}{x^{p}}^{\star},$
(59)

where
$x^{\star}=({x^{p}}^{\star})_{p\in\mathcal{N}},x=({x^{p}})_{p\in\mathcal{N}}$.

###### Proof.

(Proof of 4) By definition,

$\displaystyle D_{h}(x^{\star},x)$
$\displaystyle\textstyle=\sum_{p\in\mathcal{N}}\vartheta^{p}({x^{p}}^{\star})\\!-\vartheta^{p}({x^{p}})\\!-\epsilon^{-1}{z^{p}}^{\top}({x^{p}}^{\star}\\!-x^{p})$
(60)
$\displaystyle\textstyle=\sum_{p\in\mathcal{N}}\vartheta^{p}({x^{p}}^{\star})+{\vartheta^{p}}^{\star}(\epsilon^{-1}z^{p})-\epsilon^{-1}{z^{p}}^{\top}{x^{p}}^{\star}$
(61)
$\displaystyle\textstyle=\epsilon^{-1}\sum_{p\in\mathcal{N}}\epsilon\vartheta^{p}({x^{p}}^{\star})+\epsilon{\vartheta^{p}}^{\star}(\epsilon^{-1}z^{p})-{z^{p}}^{\top}{x^{p}}^{\star}$
(62)
$\displaystyle\textstyle=\epsilon^{-1}\sum_{p\in\mathcal{N}}\psi_{\epsilon}^{p}({x^{p}}^{\star})+{\psi_{\epsilon}^{p}}^{\star}(z^{p})-{z^{p}}^{\top}{x^{p}}^{\star},$
(63)

where we used
$\psi_{\epsilon}^{p}({x^{p}}^{\star})=\epsilon\vartheta^{p}({x^{p}}^{\star}),z^{p}\in\partial\psi_{\epsilon}^{p}(x^{p})=\epsilon\partial\vartheta^{p}(x^{p}),{\psi_{\epsilon}^{p}}^{\star}(z^{p})=\epsilon{\vartheta^{p}}^{\star}(\epsilon^{-1}z^{p}),{x^{p}}\in\partial{\psi_{\epsilon}^{p}}^{\star}(z^{p})$.
We note that
${\psi_{\epsilon}^{p}}^{\star}(z^{p})=\epsilon{\vartheta^{p}}^{\star}(\epsilon^{-1}z^{p})$
follows from the conjugate property
$\epsilon\vartheta^{p}(x^{p})\Longleftrightarrow\epsilon\vartheta^{p\star}(\epsilon^{-1}z^{p})$
[49, p. 93]. ∎

###### Proof.

(Proof of 2) Consider the following Lyapunov function,

$\displaystyle\textstyle V(\xi,\\!z,\\!\gamma)$
$\displaystyle\textstyle\\!=\\!\dfrac{\alpha\gamma}{2\beta}\|\xi-x^{\star}\|_{2}^{2}\\!+\\!\sum_{p\in\mathcal{N}}\psi^{p}_{\epsilon}({x^{p}}^{\star})-{z^{p}}^{\top}{x^{p}}^{\star}\\!+\\!\psi_{\epsilon}^{p\star}(z^{p}),$

Then taking the time-derivative of $V(\xi,z,\gamma)$ along MD2$\gamma$, and
using $\nabla{\psi^{p}_{\epsilon}}^{\star}={C^{p}_{\epsilon}}$,
$x^{p}={C^{p}_{\epsilon}}(z^{p})$,
$\dot{z}=U(x)-\gamma\alpha\beta^{-1}\dot{\xi}$, we have,

$\displaystyle\dot{V}(\xi,z,\gamma)=\dfrac{\alpha\dot{\gamma}}{2\beta}\|\xi\\!-\\!x^{\star}\|_{2}^{2}\\!+\\!\dfrac{\alpha\gamma}{\beta}(\xi-x^{\star})^{\top}\dot{\xi}\\!+\\!(C_{\epsilon}(z)\\!-\\!x^{\star})^{\top}\dot{z}$
$\displaystyle=\\!\dfrac{\alpha\dot{\gamma}}{2\beta}\|\xi\\!-\\!x^{\star}\|_{2}^{2}\\!+\\!\dfrac{\alpha\gamma}{\beta}(\xi\\!-\\!x^{\star}\\!-\\!x\\!+\\!x^{\star})^{\top}\dot{\xi}\\!+\\!(x\\!-\\!x^{\star})^{\top}U(x)$
$\displaystyle=\\!\dfrac{-\alpha\eta\gamma}{2\beta\epsilon}\|\xi-x^{\star}\|_{2}^{2}-\dfrac{\alpha\gamma}{\beta^{2}}\|\dot{\xi}\|_{2}^{2}+(x-x^{\star})^{\top}U(x)$
$\displaystyle\leq\\!\dfrac{-\alpha\eta\gamma}{2\beta\epsilon}\|\xi-x^{\star}\|_{2}^{2}+(x-x^{\star})^{\top}U(x)$
$\displaystyle\textstyle\leq\\!\dfrac{-\alpha\eta\gamma}{2\beta\epsilon}\|\xi-x^{\star}\|_{2}^{2}-\eta(D_{h}(x^{\star},x)+D_{h}(x,x^{\star}))$
$\displaystyle\textstyle\leq\\!\dfrac{-\alpha\eta\gamma}{2\beta\epsilon}\|\xi-x^{\star}\|_{2}^{2}-\eta
D_{h}(x^{\star},x)$
$\displaystyle\textstyle=\\!-\eta\epsilon^{-1}(\dfrac{\alpha\gamma}{2\beta}\|\xi-x^{\star}\|_{2}^{2}+\sum_{p\in\mathcal{N}}\psi_{\epsilon}^{p}({x^{p}}^{\star})\\!-\\!{\psi_{\epsilon}^{p}}^{\star}(z^{p})\\!-\\!{z^{p}}^{\top}{x^{p}}^{\star})$
$\displaystyle=\textstyle-\eta\epsilon^{-1}V(\xi,z,\gamma),$ (64)

where we have used
$D_{h}(x^{\star},x)=\epsilon^{-1}\sum_{p\in\mathcal{N}}\psi_{\epsilon}^{p}({x^{p}}^{\star})\\!-\\!{\psi_{\epsilon}^{p}}^{\star}(z^{p})\\!-\\!{z^{p}}^{\top}{x^{p}}^{\star}),\forall
x=C_{\epsilon}(z),z\in C_{\epsilon}^{-1}(x)$ (4). (64) implies
$V(\xi,z,\gamma)\leq e^{-\eta\epsilon^{-1}t}V(\xi_{0},z_{0},\gamma_{0})$.
Substituting in the exact expression for $V(\xi,z,\gamma)$ and
$V(\xi_{0},z_{0},\gamma_{0})$, then rearranging the resulting expression, we
obtain (21). (22) follows from $D_{h}(x^{\star},x)\geq
2^{-1}\rho\|x^{\star}-x\|_{2}^{2}$ whenever $\vartheta^{p}$ is $\rho$-strongly
convex. ∎

###### Proof.

(Proof of 3) Without loss of generality, let $\gamma=1$. Let
$y^{p}\in\Omega^{p}$ be an arbitrary strategy. Since
$\mathcal{U}^{p}(y^{p};x^{-p}(t))$ is concave in $y^{p}$,
$\mathcal{U}^{p}(y^{p};x^{-p}(t))\leq\mathcal{U}^{p}(x^{p}(t);x^{-p}(t))+\nabla_{x^{p}}\mathcal{U}^{p}(x^{p}(t);x^{-p}(t))^{\top}(y^{p}-x^{p}(t))$
= $\mathcal{U}^{p}(x^{p}(t);x^{-p}(t))+U^{p}(x)^{\top}(y^{p}-x^{p}(t))$,
$\forall t\geq 0$, therefore, the integral term of $\eqref{eqn:regret}$, can
be written as (suppressing time-index for brevity),

$\displaystyle\textstyle\int_{0}^{t}\mathcal{U}^{p}(y^{p};x^{-p})-\mathcal{U}^{p}(x)\mathrm{d}\tau\leq\int_{0}^{t}({y^{p}}-{x^{p}})^{\top}U^{p}(x)\mathrm{d}\tau,$
(65)

Next, define,

$\mathpzc{Q}^{p}(t)\coloneqq{\psi^{p}_{\epsilon}}^{\star}(z^{p}(t))-{y^{p}}^{\top}z^{p}(t)+\dfrac{\alpha}{\beta}(\dfrac{1}{2}\|\xi^{p}(t)\|_{2}^{2}-{y^{p}}^{\top}{\xi^{p}(t)}).$
(66)

Taking the time-derivative of $\mathpzc{Q}^{p}$ along the solutions of MD2 and
using $\nabla{\psi^{p}_{\epsilon}}^{\star}={C^{p}_{\epsilon}}$,
$x^{p}={C^{p}_{\epsilon}}(z^{p})$,
$\dot{z}^{p}=U^{p}(x)-\beta^{-1}\alpha\dot{\xi}^{p}$, yields,

$\displaystyle\dot{\mathpzc{Q}}^{p}(t)$
$\displaystyle\\!=\\!\nabla{\psi^{p}_{\epsilon}}^{\star}(z^{p})^{\top}\dot{z}^{p}\\!-\\!{y^{p}}^{\top}{\dot{z}^{p}}+\dfrac{\alpha}{\beta}({\xi^{p}}^{\top}\dot{\xi}^{p}\\!-\\!{y^{p}}^{\top}\dot{\xi}^{p})$
(67)
$\displaystyle={C^{p}_{\epsilon}}(z^{p})^{\top}\dot{z}^{p}-{y^{p}}^{\top}{\dot{z}^{p}}+\dfrac{\alpha}{\beta}({\xi^{p}}^{\top}\dot{\xi}^{p}-{y^{p}}^{\top}\dot{\xi}^{p})$
(68)
$\displaystyle=({x^{p}}-y^{p})^{\top}\dot{z}^{p}+\dfrac{\alpha}{\beta}(\xi^{p}-y^{p})^{\top}\dot{\xi}^{p}$
(69)
$\displaystyle=({x^{p}}\\!-\\!y^{p})^{\top}(U^{p}(x)\\!-\\!\dfrac{\alpha}{\beta}\dot{\xi}^{p})+\dfrac{\alpha}{\beta}(\xi^{p}\\!-\\!y^{p})^{\top}\dot{\xi}^{p}$
(70)
$\displaystyle=({x^{p}}-y^{p})^{\top}U^{p}(x)+\dfrac{\alpha}{\beta}(-{x^{p}}+\xi^{p})^{\top}\dot{\xi}^{p}$
(71)
$\displaystyle=({x^{p}}-y^{p})^{\top}U^{p}(x)-\dfrac{\alpha}{\beta^{2}}\|\dot{\xi}^{p}\|_{2}^{2}.$
(72)

Carrying on from (65) and using (72), we have,

$\displaystyle\textstyle\int_{0}^{t}\mathcal{U}^{p}(y^{p};x^{-p})-\mathcal{U}^{p}(x)\mathrm{d}\tau$
$\displaystyle\textstyle=\int_{0}^{t}-{\dot{\mathpzc{Q}}}^{p}(\tau)-\dfrac{\alpha}{\beta^{2}}\|\dot{\xi}^{p}\|_{2}^{2}\mathrm{d}\tau$
(73) $\displaystyle\leq\mathpzc{Q}^{p}(0)-\mathpzc{Q}^{p}(t).$ (74)

Applying Fenchel’s inequality [53] to the function $\mathpzc{Q}^{p}$, we
obtain,
$\mathpzc{Q}^{p}(t)\geq-\psi^{p}_{\epsilon}(y^{p})-\dfrac{\alpha}{2\beta}\|y^{p}\|_{2}^{2}.$
Hence
$\mathpzc{Q}^{p}(0)-\mathpzc{Q}^{p}(t)\leq\mathpzc{Q}^{p}(0)+\psi^{p}_{\epsilon}(y^{p})+\dfrac{\alpha}{2\beta}\|y^{p}\|_{2}^{2}$,
where $\mathpzc{Q}^{p}(0)$ is some finite constant (due to 3(i)). Therefore we
obtain the following inequality,

$\displaystyle\textstyle\int_{0}^{t}\\!\mathcal{U}^{p}(y^{p};x^{-p})\\!-\\!\mathcal{U}^{p}(x)\mathrm{d}\tau\
$
$\displaystyle\\!\leq\\!\mathpzc{Q}^{p}(0)\\!+\\!\psi^{p}_{\epsilon}(y^{p})\\!+\\!\dfrac{\alpha}{2\beta}\|y^{p}\|_{2}^{2}.$
(75)

Plugging into (75) the definition of regret, we obtain,

$\displaystyle\mathpzc{R}^{p}(t)$
$\displaystyle\textstyle=\max\limits_{y^{p}\in\Omega^{p}}t^{-1}\int_{0}^{t}\mathcal{U}^{p}(y^{p};x^{-p}(\tau))-\mathcal{U}^{p}(x(\tau))\mathrm{d}\tau$
$\displaystyle\leq\max\limits_{y^{p}\in\Omega^{p}}t^{-1}(\mathpzc{Q}^{p}(0)+\psi^{p}_{\epsilon}(y^{p})+\dfrac{\alpha}{2\beta}\|y^{p}\|_{2}^{2}).$
(76)

Since $\Omega^{p}$ is assumed to be compact, and the numerator is continuous
on all of $\Omega^{p}$, therefore by Weierstrass theorem [49, Theorem 2.12] a
maximum is achieved. Taking the limsup yields the desired result. When
$\alpha=0$, this proof recovers no-regret bound of MD [22]. ∎

###### Proof.

(Proof of 2) By 3(ii),
$x^{p}=C_{\epsilon}^{p}(z^{p})=\nabla{\psi_{\epsilon}^{p}}^{\star}(z^{p})$,
therefore $x^{p}\in\partial{\psi^{p}_{\epsilon}}^{\star}(z^{p})$. Using the
fact that $\partial{\psi_{\epsilon}^{p}}$ is the inverse map of
$\partial{\psi^{p}_{\epsilon}}^{\star}$, for all
$x^{p}\in\operatorname{dom}(\partial\psi^{p\star}_{\epsilon})=\operatorname{rint}(\Omega^{p})$,

$z^{p}\in\partial\psi^{p}_{\epsilon}(x^{p})\\!=\\!\\{\nabla\psi^{p}_{\epsilon}(x^{p})\\}+N_{\Omega^{p}}(x^{p}),$
(77)

From (77), we can write $z^{p}=\nabla\psi^{p}_{\epsilon}(x^{p})+n^{p}(x^{p})$
where $n^{p}(x^{p})\in N_{\Omega^{p}}(x^{p})$. Plugging this expression into
(25), we have,

$\begin{cases}\dot{\xi}^{p}&\\!\in\beta(x^{p}-\xi^{p})\\\
\dot{x}^{p}&\\!\in\gamma\nabla^{2}{\psi^{p}_{\epsilon}}^{\star}(\nabla\psi^{p}_{\epsilon}(x^{p})\\!+\\!n^{p}(x^{p}))(U^{p}(x)\\!-\\!\alpha(x^{p}\\!-\\!\xi^{p})).\end{cases}$
(78)

To reduce (78) from a set of differential inclusions to a set of ODEs, we will
show, for all $x^{p}\in\operatorname{rint}(\Omega^{p})$,

$\nabla^{2}{\psi^{p}_{\epsilon}}^{\star}(\nabla\psi^{p}_{\epsilon}(x^{p})+n^{p}(x^{p}))=\nabla^{2}{\psi^{p}_{\epsilon}}^{\star}(\nabla\psi^{p}_{\epsilon}(x^{p})).$
(79)

Indeed, $x^{p}\in\partial\psi^{p\star}_{\epsilon}(z^{p})\Longleftrightarrow
z^{p}\in\partial\psi^{p}_{\epsilon}(x^{p})\overset{\eqref{eqn:sub_diff_psi_p}}{\Longleftrightarrow}\nabla\psi^{p}_{\epsilon}(x^{p})+n^{p}(x^{p})\in\partial\psi^{p}_{\epsilon}(x^{p})$.
The last equation,
$\nabla\psi^{p}_{\epsilon}(x^{p})+n^{p}(x^{p})\in\partial\psi^{p}_{\epsilon}(x^{p})$
means, $\forall y^{p}\in\Omega^{p}$
$\psi^{p}_{\epsilon}(y^{p})\geq\psi^{p}_{\epsilon}(x^{p})+(\nabla\psi^{p}_{\epsilon}(x^{p})+n^{p}(x^{p}))^{\top}(y^{p}-x^{p})$.
Since $n^{p}(x^{p})\in N_{\Omega^{p}}(x^{p})$, therefore
$n^{p}(x^{p})^{\top}(y^{p}-x^{p})\leq 0$, hence it is also true that
$\psi^{p}_{\epsilon}(y^{p})\geq\psi^{p}_{\epsilon}(x^{p})+\nabla\psi^{p}_{\epsilon}(x^{p})^{\top}(y^{p}-x^{p})$
or equivalently,
$\nabla\psi^{p}_{\epsilon}(x^{p})\in\partial\psi^{p}_{\epsilon}(x^{p})$, or
$\partial\psi^{p\star}_{\epsilon}(\nabla\psi^{p}_{\epsilon}(x^{p}))=\partial\psi^{p\star}_{\epsilon}(\partial\psi^{p}_{\epsilon}(x^{p}))$.
Hence by re-substituting in (77) and using the single-valuedness of
$\partial\psi^{p\star}_{\epsilon}$, we obtain,
$\nabla{\psi_{\epsilon}^{p}}^{\star}(\nabla{\psi^{p}_{\epsilon}}(x^{p})+n^{p}(x^{p}))=\nabla{\psi_{\epsilon}^{p}}^{\star}(\nabla{\psi^{p}_{\epsilon}}(x^{p}))$.
Taking the Jacobian of the preceding equality yields our desired result. ∎

The proof of 4 requires the following definitions from [55, 57, 56]. Recall
that a semiflow $\Phi$ on a metric space $(\mathcal{M},d)$ is a continuous map
$\Phi:\mathbb{T}\times\mathcal{M}\to\mathcal{M},(t,x)\mapsto\Phi(t,x)=\Phi_{t}(x)$,
where $\mathbb{T}=\mathbb{R}_{\geq 0}$, such that, $\Phi_{0}=Id$ and
$\Phi_{t+s}=\Phi_{t}\circ\Phi_{s}$, for all $(t,s)\in\mathbb{R}_{\geq
0}\times\mathbb{R}_{\geq 0}$. If $\mathbb{T}=\mathbb{R}$, then $\Phi$ defines
a flow. We say that a continuous function $f:\mathbb{R}_{>0}\to\mathcal{M}$ is
an asymptotic pseudotrajectory of $\Phi$ if
$\lim\limits_{t\to\infty}\sup_{0\leq h\leq T}d(f(t+h),\Phi_{h}(f(t)))=0$ for
any $T>0$. A set $\mathcal{A}\subset\mathcal{M}$ is positively invariant if
$\Phi_{t}(\mathcal{A})\subset\mathcal{A}$ for all $t\geq 0$ and invariant if
$\Phi_{t}(\mathcal{A})=\mathcal{A}$ for all $t\in\mathbb{T}$.

###### Proof.

(Proof of 4) For simplicity, we assume that all constants associated with MD2
are set to $1$. We break up our proof into the following steps:

  1. 1.

We first need to show that DA2 is the stochastic approximation of MD2.

Rearranging the $Z^{p}_{k}$ update for DA2,
$\dfrac{Z^{p}_{k+1}-Z^{p}_{k}}{t_{k+1}}=\widehat{U}^{p}_{k+1}=U^{p}(X_{k})+\zeta^{p}_{k+1}$
and taking the expectation with respect to the filtration
$\mathcal{F}_{k}=\sigma(\Xi_{0},Z_{0},\zeta_{0},\ldots,\zeta_{k})$ yields,

$\operatorname{\mathbb{E}}(\dfrac{Z^{p}_{k+1}-Z^{p}_{k}}{t_{k+1}}|\mathcal{F}_{k})=\operatorname{\mathbb{E}}(U^{p}(X_{k})+\zeta^{p}_{k+1}-(X^{p}_{k}-\Xi^{p}_{k})|\mathcal{F}_{k}),$
(80)

Since
$\operatorname{\mathbb{E}}(U^{p}(X_{k})|\mathcal{F}_{k})=U^{p}(X_{k})-(X^{p}_{k}-\Xi^{p}_{k})$
and $\operatorname{\mathbb{E}}(\zeta^{p}_{k+1}|\mathcal{F}_{k})=0$ by 6,
therefore,

$\operatorname{\mathbb{E}}(\dfrac{Z^{p}_{k+1}-Z^{p}_{k}}{t_{k+1}}|\mathcal{F}_{k})=U^{p}(X_{k})-(X^{p}_{k}-\Xi^{p}_{k}),$
(81)

And similarly for $\Xi_{k+1}^{p}$, we have,

$\operatorname{\mathbb{E}}(\dfrac{\Xi^{p}_{k+1}-\Xi^{p}_{k}}{\tau_{k+1}}|\mathcal{F}_{k})=X^{p}_{k}-\Xi^{p}_{k},$
(82)

which shows that MD2 is the mean ODE of DA2.

  2. 2.

Next, we need to build interpolated processes for $Z_{k}$ and $\Xi_{k}$ and
show that the interpolated processes converge to a rest point of MD2. The full
construction process is provided as follows: let
$\\{\ell_{k}\\}_{k\in\mathbb{N}}$ be a sequence such that,

$\textstyle\ell_{0}=0\quad\ell_{k}=\sum_{i=1}^{k}t_{i},k\geq 1,$ (83)

and define the interpolated process,

$\textstyle\underline{Z}^{p}(\ell_{k}+s)=Z^{p}_{k}+s\frac{Z^{p}_{k+1}-Z^{p}_{k}}{\ell_{k+1}-\ell_{k}},$
(84)

where $0\leq s<t_{k+1}$. Similarly, let
$\\{\ell^{\prime}_{k}\\}_{k\in\mathbb{N}}$ be a sequence such that,

$\textstyle\ell^{\prime}_{0}=0\quad\ell^{\prime}_{k}=\sum_{i=1}^{k}\tau_{i},k\geq
1,$ (85)

and define the interpolated process associated with $\Xi^{p}$ as,

$\textstyle\underline{\Xi}^{p}(\ell^{\prime}_{k}+s^{\prime})=\Xi^{p}_{k}+s^{\prime}\frac{\Xi^{p}_{k+1}-\Xi^{p}_{k}}{\ell^{\prime}_{k+1}-\ell^{\prime}_{k}},$
(86)

where $0\leq s^{\prime}<\tau_{k+1}$.

  3. 3.

Let $\underline{Z}^{p}:\mathbb{R}_{\geq 0}\to\mathbb{R}^{n_{p}}$ and
$\underline{\Xi}^{p}:\mathbb{R}_{\geq 0}\to\mathbb{R}^{n_{p}}$ denote the
continuous functions (as functions of $s,s^{\prime}$ respectively) associated
with the above processes and let
$\underline{Z}=(\underline{Z}^{p})_{p\in\mathcal{N}},\underline{\Xi}=(\underline{\Xi}^{p})_{p\in\mathcal{N}}$
be the stacked-vector of all the individual interpolated processes.

  4. 4.

Next, define the overall process,

$\underline{W}:\mathbb{R}_{\geq
0}\to\mathbb{R}^{2n}\quad\underline{W}=(\underline{\Xi},\underline{Z}),$ (87)

Define the interpolated process $\underline{W}$. By Proposition 4.1 and
Proposition 4.2 of [55], under 5 \- 8 $\underline{W}$ is an asymptotic pseudo-
trajectory of the semiflow $\Phi:\mathbb{R}_{\geq
0}\times\mathbb{R}^{2n}\to\mathbb{R}^{2n}$ induced by $\dot{\omega}=F(\omega)$
(39), that is,

$\textstyle\lim_{t\to\infty}\sup_{0\leq h\leq
T}\|\underline{W}(t+h)-\Phi(h,\underline{W}(t))\|_{2}=0,$ (88)

for any $T>0$.

  5. 5.

By 7, $\underline{W}$ has compact closure, i.e., is pre-compact. Since we have
shown that $\underline{W}$ is a pre-compact APT of $\Phi$ induced by
$\dot{\omega}$, by [55, Theorem 5.7(i)] the limit set

$\textstyle L(\underline{W})=\bigcap_{t\geq
0}\operatorname{cl}({\underline{W}([t,\infty))}),$ (89)

is internally chain transitive, which by [55, Prop 5.3], $L(\underline{W})$ is
a compact invariant set.

  6. 6.

Recall that $V:\mathbb{R}^{2n}\to\mathbb{R}$ was used as the Lyapunov function
associated with $\dot{\omega}=F(\omega)$ in 1. Consider the set of critical
points of $V$,
$\mathcal{E}=\\{(\xi,z)\in\mathbb{R}^{2n}|\dot{V}=0\\}=\\{(\xi,z)\in\mathbb{R}^{n}\times\mathbb{R}^{n}|\xi=x=C(z)\\}$.
By the arguments in the proof of 1, since
$V(\xi^{\prime},z^{\prime})<V(\xi,z)$ for all
$(\xi,z)\in\mathbb{R}^{2n}\backslash\mathcal{E}$,
$(\xi^{\prime},z^{\prime})=\Phi(t,(\xi,z))$ and
$V(\xi^{\prime},z^{\prime})\leq V(\xi,z)$ for all
$(\xi,z)\in\mathcal{E},(\xi^{\prime},z^{\prime})=\Phi(t,(\xi,z))$, therefore
$V$ is a Lyapunov function for $\mathcal{E}$ [56].

  7. 7.

Since $V(\xi,z)$ is a Lyapunov function for
$\mathcal{E}\subset\mathbb{R}^{2n}$, and $V(\mathcal{E})$ is a constant, which
implies $\operatorname{int}(V(\mathcal{E}))=\emptyset$, therefore by [56,
Prop. 3.27], $L(\underline{W})$ is contained in $\mathcal{E}$. By our
assumption that every NE is a mere VSS, hence every point in $\mathcal{E}$
corresponds to an interior mere VSS and $L(\underline{W})\subset\mathcal{E}$,
therefore $L(\underline{W})$ contains a compact subset of the rest points of
$\dot{\omega}$.

  8. 8.

By the definition of a limit set, for any $\underline{W}(0)$, the interpolated
process $\underline{W}(t)$ converges as $t\to\infty$.

  9. 9.

From our construction of the interpolated process and the diminishing step-
size assumption, i.e., $0\leq s<\tau_{k+1}$, $0\leq s^{\prime}<t_{k+1}$ and
$\lim_{k\to\infty}t_{k}=\lim_{k\to\infty}\tau_{k}=0$, the convergence of the
interpolated processes $\underline{\Xi}$ and $\underline{Z}$ implies the
convergence of $\Xi_{k}$, $Z_{k}$ respectively.

  10. 10.

By continuity of $C_{\epsilon}$, it follows that $X_{k}=C_{\epsilon}(Z_{k})$
converges almost surely an interior mere VSS $x^{\star}$.

∎

## References

  * [1]
  * [2] P. Mertikopoulos and Z. Zhou, “Learning in games with continuous action sets and unknown payoff functions,” Mathematical Programming, vol. 173, no. 1-2, pp. 465-507, 2019.
  * [3] T. Migot and M.-G. Cojocaru, “On Minty-variational inequalities and evolutionary stable states of generalized monotone games,” Operations Research Letters, vol. 49, no. 1, pp. 96 - 100, 2021.
  * [4] A. S. Nemirovsky and D. B. Yudin, Problem Complexity and Method Efficiency in Optimization (Discrete Mathematics). Wiley, 1983.
  * [5] M. Raginsky and J. Bouvrie, “Continuous-time stochastic mirror descent on a network: Variance reduction, consensus, convergence”, in Proc. IEEE CDC, pp. 6793-6800, 2012.
  * [6] W. Krichene, A. Bayen, and P. Bartlett, “Accelerated Mirror Descent in Continuous and Discrete Time”, in NIPS 29, 2015.
  * [7] A. Beck, M. Teboulle, “Mirror descent and nonlinear projected subgradient methods for convex optimization,” Oper. Research Lett., vol. 31, no. 3, pp. 167-175, 2003.