This paper was converted on www.awesomepapers.org from LaTeX by an anonymous user.
Want to know more? Visit the Converter page.

Metastability in Stochastic Replicator Dynamics

Konstantin Avrachenkov111INRIA Sophia Antipolis, 2004 Route des Lucioles, B.P.93, 06902, Sophia Antipolis Cedex, France. Tel: +33 4 92.38.77.51 E-mail: k.avrachenkov@inria.fr  and  Vivek S. Borkar222Department of Electrical Engineering, Indian Institute of Technology Bombay, Powai, Mumbai 400076, India. E-mail: borkar.vs@gmail.com

Abstract: We consider a novel model of stochastic replicator dynamics for potential games that converts to a Langevin equation on a sphere after a change of variables. This is distinct from the models of stochastic replicator dynamics studied earlier. In particular, it is ill-posed due to non-uniqueness of solutions, but is amenable to a natural selection principle that picks a unique solution. The model allows us to make specific statements regarding metastable states such as small noise asymptotics for mean exit times from their domain of attraction, and quasi-stationary measures. We illustrate the general results by specializing them to replicator dynamics on graphs and demonstrate that the numerical experiments support theoretical predictions.

Key words: stochastic replicator dynamics; Langevin equation on sphere; metastable states; intermittency; small noise asymptotics; mean exit time; quasi-stationary distributions

MSC: Primary: 91A22 Secondary: 91A25; 60H10; 34F05

1 Introduction

This work is motivated by the studies on intermittent dynamics of complex large-scale systems [23], [54], [62]. We model the dynamics of a large-scale complex system by replicator dynamics with noise, a classical model in evolutionary biology. Replicator dynamics per se is an extensively studied dynamics [36], [55], [65]. Its noisy version is much less so. In their seminal work, Foster and Young [27] considered small noise asymptotics of the stationary regime of replicator dynamics perturbed by small noise and used it to formulate the notion of stochastically stable equilibria, i.e., Nash equilibria that arise as limit(s) when noise is reduced to nil in the stationary regime. This provided a selection principle for Nash equilibria (see also [67]). This is of interest because multiple Nash equilibria is a commonplace situation and one does need a selection principle of some sort to narrow down the relevant possibilities. While early work in this direction was based on static refinements of the notion of Nash equilibrium [34], [64], later developments sought dynamic models where some Nash equilibria would naturally emerge as stable equilibria of the dynamics [30], [55]. This, however, only partially solves the problem, because, as mentioned above, even stable equilibria can be many. The important observation of Foster and Young was that only some of them may be stable under stochastic perturbations in the stationary regime, a physically relevant concern. (See also [45], [47] for interesting takes on behavioral aspects of noise and [14], [25], [32] for connections with finite population models. The general philosophy of using small noise limits for selection principles in ill-posed dynamics has been attributed to Kolmogorov in [26], p. 626.)

The elegant machinery for small noise asymptotics of noise-perturbed dynamics due to Freidlin and Wentzell [28] gives us a handle to tag these equilibria in a precise manner, viz., as minimizers of the so-called Freidlin-Wentzell potential associated with the process. This is precisely what Foster and Young did. But this still leaves wide open the issue of metastability, a familiar phenomenon in statistical physics [18]. Well known examples are supercooled water, air supersaturated with vapor, etc., which on suitable perturbation result in resp. freezing (ice) or condensation (rain). There have been various analytic approaches to quantitative analysis of metastability, see the opening chapters of [18] for a succinct overview and historical account. Our approach follows that formulated by Kramers [43] who viewed them as equilibria of a gradient system perturbed by noise. The later seminal work of Freidlin and Wentzell [28] generalized this viewpoint considerably, encompassing non-gradient systems and attractors more general than equilibria. In this viewpoint, metastable states are stable equilibria for the deterministic dynamics that will not arise as a small noise limit in the stationary regime, but are still important because of the celebrated quasi-stationarity phenomenon [24]: the process may spend long time near such equilibria on time scales relevant to us before approaching true stationarity. An alternative though related viewpoint is that these are equilibria which the system gets out of when given a perturbation of a certain magnitude. In case of gradient systems, the size of the perturbation will be related to the barrier between the current metastable state and its neighboring equilibria as captured by the optimal trajectory (i.e., one minimizing the ‘rate function’ of the associated large deviations) needed to climb in order to get out of the current ‘valley’ or domain of attraction for the noiseless system. Analogous ideas exist for non-gradient systems in terms of the aforementioned Freidlin-Wentzell potential. See [50] for an extensive account of this approach and [8] for some recent developments.

A related phenomenon is that of intermittency, where one sees the process spend nontrivial time in the neighborhood of a succession of stable equilibria with abrupt noise-induced episodic transitions between them (see [23], [54], [62] for examples from physics and biology). This will be the case when the noise is small but not too small so that these transitions occur on an observable time scale. This should not be confused with the phenomenon of intermittency in chaotic systems wherein a deterministic dynamics exhibits long periods of slow variation punctuated by spurts of chaotic behavior [51]. The analysis of Fredlin-Wentzell shows that the process may be viewed as approximating a Markov chain over stable equilibria (more generally, stable attractors) with transition times and probabilities characterized in terms of a suitable ‘action’ functional [28]. The small noise asymptotics gives a feel for the noise levels for which one expects intermittency instead of quasi-stationarity.

In view of the foregoing one important numerical quantity associated with such a metastable equilibrium is the mean exit time from its domain of attraction, or rather, the scaling properties thereof in the small noise limit. Our aim here is to present one such result in the case of replicator dynamics with potential. While bounds on this quantity are available [39] for the specific model considered by [19], [29], we give a finer result using the Freidlin–Wentzell theory [28]. A highlight of our approach is that it is based on an interesting connection between replicator dynamics for potential games and gradient descent on a sphere, a compact manifold without boundary which allows the application of Freidlin-Wentzell theory without facing the complications encountered by Foster and Young in doing so - see, e.g., the correction note to [27]. It provides an alternative treatment to such issues distinct from that of [39] that is simpler in many ways. We also gain some important insights into quasi-stationary distributions.

In related works [57], [58], replicator dynamics is viewed as a large population ‘mean field’ limit of finite population stochastic models (see, e.g., [16]) and large deviations phenomena in the stationary regime are analyzed. Our emphasis, in contrast, is on quasi-stationary behavior. As noted in [52], the quasi-stationary behavior is closely associated with the ‘hardness’ of the underlying optimization problem of minimizing the potential. This also reflects in the slow convergence of the simulated annealing algorithm for non-convex potentials, which is a time-inhomogeneous stochastic gradient descent with slowly decreasing noise [2], [20], [21]. Other related works on asymptotic behavior of noisy replicator dynamics are [35], [48], [49]. Mathematically, what sets our model apart is a diffusion matrix involving in non-diagonal positions square-roots of state variables that puts it beyond available well-posedness results for stochastic differential equations. In fact it is ill-posed due to non-uniqueness of solutions, thus we need to resort to a well motivated selection principle to pick a unique choice.

To summarize, the main contributions of the present work are: (a) we introduce a new variant of the stochastic replicator dynamics, which in particular allows us to work around the problem of “corners” in [27], (b) we propose a transformation that converts our model to a Langevin equation on a sphere; (c) a selection principle to overcome the ill-posedness of the original model is then proposed; (c) we characterize for our model the exit time asymptotics from metastable states and also characterize its quasi-stationary distributions; (d) we apply these general results to study the intermittent behaviour of large complex systems described by graphs; and finally, (e) we illustrate main theoretical findings by numerical examples.

The article is organized as follows. The next section introduces our model of stochastic replicator dynamics and discusses in detail its relationship to other stochastic replicator dynamics in literature. Section 3 uses a transformation to map it into the Langevin equation on a sphere and exploits the structure therein to analyze the small noise asymptotics. Section 4 specializes the results to the case of quadratic potentials on graphs and makes some relevant observations for this special case which has important applications such as the study of the intermittent behaviour in complex large-scale systems [23], [54], [62] and clique detection problem [4], [10], [11], [12]. Section 5 presents some numerical experiments to support the theoretical findings. We conclude the paper with suggestions for future research in Section 6.

2 Stochastic replicator dynamics with potential

Replicator dynamics, introduced by Taylor and Jonker [61], is a model for behavior adjustment of population games, i.e., a game theoretic scenario in which we analyze collective behavior of a population of agents belonging to certain ‘types’ or ‘species’ in conflict situation, retaining only their species/type tag and erasing their individual identity. This is because their strategic behavior is dictated by the type, as also the environment they face is dictated by the collective sizes and strategic behavior of other types. The origin of this particular dynamics of population games (among many, see [55] for an extensive treatment) lies in models of Darwinian evolution in biology [36]. Issues of interest typically are convergence to an equilibrium behavior, persistence or non-extinction, etc. Among others, one very significant outcome of this theory is the link between the so-called evolutionarily stable states, a biologically relevant equilibrium notion introduced by John Maynard Smith, and asymptotically stable equilibrium points of replicator dynamics. See [36], [55], [65] for a detailed treatment. We further restrict to a special subclass of these, viz., potential games wherein the payoff function is restricted to be of a particular form, see again [55] for an extensive account of the deterministic case. This subclass, defined later, is much more analytically amenable and has important applications, e.g., in congestion games.

Our focus will be on stochastic replicator dynamics with potential. We begin with a motivating example of deterministic replicator dynamics with linear payoffs and then introduce the more general paradigm.

2.1 Deterministic replicator dynamics and potential games

Let

Pn:={x=[x1,,xn]:xi0i,i=1nxi=1}P^{n}:=\left\{x=[x_{1},\cdots,x_{n}]:x_{i}\geq 0\ \forall i,\ \sum_{i=1}^{n}x_{i}=1\right\}

denote the (n1)(n-1)-dimensional probability simplex. The deterministic replicator dynamics with linear payoffs is the following differential equation in PnP^{n}:

dxidt(t)=xi(t)([Mx(t)]ixT(t)Mx(t)),\frac{dx_{i}}{dt}(t)=x_{i}(t)\left([Mx(t)]_{i}-x^{T}(t)Mx(t)\right), (1)

where the n×nn\times n matrix MM is a payoff matrix, i.e., the iith element [Mx]i[Mx]_{i} of MxMx is the payoff to species ii when the population profile is xx. We assume that MM is symmetric.

In evolutionary biology, xi(t)x_{i}(t) is interpreted as the fraction of species ii in the population, x(t)x(t) thus being the overall population profile. The quantity [Mx(t)]i[Mx(t)]_{i} is the payoff for species ii at time tt when the population profile is x(t)x(t). The quantity xT(t)Mx(t)x^{T}(t)Mx(t) is then the average payoff for the population at time tt. Thus, the population share of species ii waxes or wanes at a rate proportional to its excess, resp., deficit of payoff over the population average. The standard equilibrium notion for population game is that of Nash equilibrium defined as follows: xx^{*} is a Nash equilibrium if xTMx(x)TMx,xPnx^{T}Mx^{*}\leq(x^{*})^{T}Mx^{*},\ \forall\ x\in P^{n}. An analogous definition is also used when MxMx is replaced by a nonlinear payoff function V(x)V(x), as we do a bit later. A population profile x^\hat{x} is evolutionarily stable if x^TMx>xTMx\hat{x}^{T}Mx>x^{T}Mx for xPnx\in P^{n} sufficiently close to x^\hat{x}. Then Nash equilibria are the equilibrium points of the deterministic replicator dynamics and the evolutionarily stable states are asymptotically stable equilibria, though neither converse holds (see e.g., [36], [55], [65].) While natural and appealing, the problem with stopping with these definitions is that of multiple equilibria as already mentioned, calling for a selection principle which stochastic replicator dynamics tries to provide.

Equation (1) can be written in the following matrix-vector form

x˙(t)=diag(x(t))(Mx(t)x(t)TMx(t)1),\dot{x}(t)=\mbox{diag}(x(t))(Mx(t)-x(t)^{T}Mx(t)\textbf{1}), (2)

where we introduce the notation:

  • diag(z)(z) for z=[z1,,zn]Tnz=[z_{1},\cdots,z_{n}]^{T}\in\mathcal{R}^{n} is the diagonal matrix whose iith diagonal entry is zii;z_{i}\ \forall i; and,

  • 1 is the constant vector of all 11’s.

Then, we observe that (2) is a particular case of

x˙(t)=diag(x(t))(F(x(t))x(t)TF(x(t))1),\dot{x}(t)=\mbox{diag}(x(t))(\nabla F(x(t))-x(t)^{T}\nabla F(x(t))\textbf{1}), (3)

with F(x):=12xTMxF(x):=\frac{1}{2}x^{T}Mx. Given x(t)Pnx(t)\in P^{n}, the right hand side summed over its components adds to zero, implying that ixi(t)\sum_{i}x_{i}(t)\equiv a constant on PnP^{n}. Since the iith component of this dynamics is of the form

x˙i(t)=xi(t)gi(t)\dot{x}_{i}(t)=x_{i}(t)g_{i}(t)

for a suitably defined gig_{i}, its solution is of the form xi(0)exp(0tgi(s)𝑑s)0x_{i}(0)exp(\int_{0}^{t}g_{i}(s)ds)\geq 0, which together with the preceding observation implies that the simplex of probability vectors is invariant under (3). In fact, it is also clear that if a component is zero for some tt, it is always so - this is immediate from the uniqueness of solutions. So the faces of PnP^{n} are also invariant under (3) and the trajectory in any face cannot enter a lower dimensional face in finite time.

For a replicator dynamics of form (3), FF is said to be its associated potential function, assumed continuously differentiable (which it indeed is in our case). Then it is easy to check (see [65], p. 110, or [55], Section 7.1, for a more general statement and an extensive account of the deterministic dynamics (3)) that

ddtF(x(t))=ixi(t)(Fxi(x(t)))2(ixi(t)(Fxi(x(t))))20,\frac{d}{dt}F(x(t))=\sum_{i}x_{i}(t)\left(\frac{\partial F}{\partial x_{i}}(x(t))\right)^{2}-\left(\sum_{i}x_{i}(t)\left(\frac{\partial F}{\partial x_{i}}(x(t))\right)\right)^{2}\geq 0,

with equality on the set

Aeq:={x:Fxi=Fxj,ij,i,jsupport(x)}.A_{eq}:=\left\{x:\frac{\partial F}{\partial x_{i}}=\frac{\partial F}{\partial x_{j}},\ \forall\ i\neq j,\ i,j\in\mbox{support}(x)\right\}.

This includes the critical points of FF. Then the above inequality implies that FF is non-decreasing along any trajectory of (3) and strictly increasing outside AeqA_{eq}. This ensures convergence of all trajectories to AeqA_{eq}. Note that AeqA_{eq} is specified by a set of nonlinear equations. We shall assume that it is a union of lower dimensional manifolds, which can usually be verified by using the implicit function theorem. In that case, it is clear that for any points in AeqA_{eq} other than the local maxima of FF, there will be a direction transversal to AeqA_{eq} along which FF must increase (otherwise it would have been a local maximum). (See [55], p. 273, for a related discussion.) Thus F-F serves as a Liapunov function for (3), implying convergence to the local maxima of FF for generic (i.e., belonging to an open dense set) initial conditions. Then (2) is a special case with F(x)=12xTMxF(x)=\frac{1}{2}x^{T}Mx. The entire development to follow goes through for general continuously differentiable F:PnF:P^{n}\mapsto\mathcal{R}, so we shall consider this general case henceforth.

2.2 Stochastic replicator dynamics

The stochastic or noise-perturbed replicator dynamics is the following stochastic differential equation in PnP^{n}:

dxi(t)=xi(t)(Fxi(x(t))x(t)TF(x(t))1)dt+εxi(t)Γi(x(t))dW(t),dx_{i}(t)=x_{i}(t)\left(\frac{\partial F}{\partial x_{i}}(x(t))-x(t)^{T}\nabla F(x(t))\textbf{1}\right)dt+\varepsilon x_{i}(t)\Gamma_{i}(x(t))dW(t), (4)

where

  • W(t)=[W1(t),,Wn(t)]T,t0,W(t)=[W_{1}(t),...,W_{n}(t)]^{T},\ t\geq 0, is the standard nn-dimensional Wiener process;

  • ϵ>0\epsilon>0 is a small parameter;

  • F():PnF(\cdot):P^{n}\mapsto\mathcal{R} is a continuously differentiable potential function and Fxi(x)\frac{\partial F}{\partial x_{i}}(x) is a payoff to species ii when the population profile is xx;

  • Γ()=[[Γij()]]1i,jn:Pnn×n\Gamma(\cdot)=[[\Gamma_{ij}(\cdot)]]_{1\leq i,j\leq n}:P^{n}\mapsto\mathcal{R}^{n\times n} with Γi():=\Gamma_{i}(\cdot):= the iith row thereof.

In order that

x(t)=[x1(t),,xn(t)]Tx(t)=[x_{1}(t),\cdots,x_{n}(t)]^{T}

truly evolves in a probability simplex, we require that

xTΓ(x)0,xPn.x^{T}\Gamma(x)\equiv 0,\quad\forall x\in P^{n}. (5)

In [27], there is an additional requirement of continuity on Γ()\Gamma(\cdot) which ensures well-posedness of (4). It is worth noting that in the model we shall propose, we do away with this: our choice is continuous in the interior of PnP^{n}, but singular at its boundaries. The product xiΓi(x)x_{i}\Gamma_{i}(x) is, however, continuous, but not Lipschitz at the boundary. Not surprisingly, well-posedness fails - the solutions are not unique. Note also that our state space is the simplex PnP^{n} and hence the above applies only to PnP^{n}, the rest of n\mathcal{R}^{n} is irrelevant. Except for the aforementioned regularity requirement, our model does fit the Foster-Young framework [27] in its overall structure.

In [27], the noise component of (4) is attributed to mutation and migration. The latter is also used to justify a reflected diffusion formulation. The continuity assumption on Γ\Gamma suffices to ensure the existence of a weak solution, but since (5) implies in particular that the ensuing diffusion is not uniformly nondegenerate, some additional care is needed to ensure uniqueness. For example, Lipschitz condition would suffice. That said, the variant we introduce later does away with the uniqueness requirement and opts for a selection principle instead – more on that later.

The standard framework for analyzing diffusions described by stochastic differential equations that allows one to handle a sufficiently rich class of problems is the Stroock-Varadhan martingale formulation, which is a ‘weak’ formulation as opposed to the original Ito formulation in terms of ‘strong’ solutions. In this formulation, one is given a densely defined local333Extensions to non-local operators for, e.g., discontinuous processes, exist. operator 𝒜\mathcal{A} on the space of bounded continuous functions on the state space, with suitable technical conditions. By a solution for a given initial condition xx, one means a process X()X(\cdot) defined on some probability space so that

f(X(t))0t𝒜f(X(s))𝑑s,t0,f(X(t))-\int_{0}^{t}\mathcal{A}f(X(s))ds,\ t\geq 0,

is a martingale with respect to the natural filtration of σ\sigma-fields generated by X(t),t0X(t),t\geq 0, for bounded continuous ff on the state space in the domain of 𝒜\mathcal{A}. The problem is well-posed if this process is unique in law for each xx and the resulting laws, tagged by xx, satisfy the Chapman-Kolmogorov equation. In case of non-uniqueness, under mild continuity requirements one has a nonempty compact family of solutions for each xx. A selection principle then entails picking one solution per xx so that the Chapman-Kolmogorov equation is satisfied. The classic work in this direction is that of [44], which, however, has limitations that we shall point out later. We shall employ a different selection principle suited for our purposes, which can be developed in a more principled manner. In particular, it is in the spirit of the Kolmogorov philosophy alluded to in the introduction and thus well-motivated for game theoretic application, somewhat in the spirit of ‘trembling hand perfect equilibria’, a refinement of Nash equilibria in static framework defined in terms of stability under random perturbations (see, e.g., [31], [67]).

On purely mathematical grounds, an advantage of our formulation is that it is the natural extension of deterministic replicator dynamics for which the celebrated coordinate transformation of [3], [59] that maps deterministic replicator dynamics for potential games to a gradient flow on sphere, extends to the stochastic case. In fact we exploit this fact in a major way in our analysis. Furthermore, it avoids artificial restrictions on the state space in the passage from deterministic to stochastic and when converted to gradient dynamics on sphere, it allows for seamless conversion from Ito to Stratonovich form and back due to absence of a correction term, allowing either formulation to be used according to convenience. All these issues are explained later in this work, as they require some preparatory material.

Note that PnP^{n} is an (n1)(n-1)-dimensional manifold with boundary and corners. This presented some problem in the analysis of small noise limit for the noisy replicator dynamics in [27], which required the authors to use the extension of the Freidlin-Wentzell theory due to Anderson and Orey [5] that handles such cases (see the correction note to [27]). We propose here a simpler route by confining ourselves to an alternative of Γ\Gamma, which we motivate next.

Let Vi(x)=Fxi(x)V_{i}(x)=\frac{\partial F}{\partial x_{i}}(x) denote the payoff accumulation rate for species ii when the population profile is xPnx\in P^{n}. Consider first the unnormalized population dynamics for the iith species given by

dri(t)=ri(t)Vi(x(t))dt+2ϵri(t)r(t)dWi(t),t0,dr_{i}(t)=r_{i}(t)V_{i}(x(t))dt+2\epsilon\sqrt{r_{i}(t)r(t)}dW_{i}(t),\ t\geq 0, (6)

where ri(t)r_{i}(t) is the population of species ii at time tt and

r(t):=i=1nri(t),x~i(t):=ri(t)r(t),r(t):=\sum_{i=1}^{n}r_{i}(t),\ \ \tilde{x}_{i}(t):=\frac{r_{i}(t)}{r(t)}, (7)

is resp. the total population at time tt and the population share (as a fraction of total population) of species ii at time tt (00=0\frac{0}{0}=0 by convention). We assume Vi:PnnV_{i}:P^{n}\mapsto\mathcal{R}^{n} to be Lipschitz so that well-posedness of the corresponding deterministic dynamics

r˙(t)=ri(t)Vi(x(t))\dot{r}(t)=r_{i}(t)V_{i}(x(t))

in {r=[r1,,rn]T:ri0i;iri>0}\{r=[r_{1},\cdots,r_{n}]^{T}:r_{i}\geq 0\ \forall i;\ \sum_{i}r_{i}>0\} is not an issue. We do not consider the well-posedness of (6) which is a degenerate diffusion with a non-Lipschitz diffusion coefficient. But we do so later in this article for the stochastic replicator dynamics we get after normalization. Our point of departure from earlier models of stochastic replicator dynamics is the stochastic forcing term, which is to be interpreted as follows. The noise seen by the iith species as a whole has incremental variance proportional to its own population and the population of each of other species (inclusive of itself, because the ‘infinitesimal’ agents of a species also interact with other such agents in the same species), summed over the latter. In other words, it is j=1nri(t)rj(t)Δt\propto\sum_{j=1}^{n}r_{i}(t)r_{j}(t)\Delta t where Δt\Delta t is the time increment. This is motivated by the fact that the variance of a sum of i.i.d. random variables scales linearly with the number thereof, therefore the net variance in the stochastic interaction of two populations should be proportional to the product of their respective sizes, assuming that each ‘infinitesimal’ interaction between two agents picks independent noise. This leads to the noise term 2ϵri(t)r(t)dWi(t)2\epsilon\sqrt{r_{i}(t)r(t)}dW_{i}(t) above where ϵ>0\epsilon>0 is a scale parameter and Wi()W_{i}(\cdot) are independent standard Brownian motions. This is in contrast to the earlier models that add noise directly to the normalized (replicator) dynamics or to the payoffs in the unnormalized population dynamics - the noise is now attached to each (infinitesimal) interaction episode and not to any overall gross feature such as the net payoff. This makes ours a game theoretically distinct model.

A direct application of the Ito formula to (6) and (7) leads to:

Lemma 2.1 The process x~()=[x~(1),,x~(n)]T\tilde{x}(\cdot)=[\tilde{x}(1),\cdots,\tilde{x}(n)]^{T} satisfies

dx~i(t)\displaystyle d\tilde{x}_{i}(t) =\displaystyle= x~i(t)(Vi(x~(t))j=1nx~j(t)Vj(x~(t)))dt+\displaystyle\tilde{x}_{i}(t)\Big{(}V_{i}(\tilde{x}(t))-\sum_{j=1}^{n}\tilde{x}_{j}(t)V_{j}(\tilde{x}(t))\Big{)}dt\ + (9)
2ϵ(x~i(t)dWi(t)x~i(t)j=1nx~j(t)dWj(t))\displaystyle 2\epsilon\Big{(}\sqrt{\tilde{x}_{i}(t)}dW_{i}(t)-\tilde{x}_{i}(t)\sum_{j=1}^{n}\sqrt{\tilde{x}_{j}(t)}dW_{j}(t)\Big{)}
=\displaystyle= x~i(t)[(Fx~i(x~(t))jx~j(t)Fx~j(x~(t)))dt\displaystyle\tilde{x}_{i}(t)\Big{[}\Big{(}\frac{\partial F}{\partial\tilde{x}_{i}}(\tilde{x}(t))-\sum_{j}\tilde{x}_{j}(t)\frac{\partial F}{\partial\tilde{x}_{j}}(\tilde{x}(t))\Big{)}dt
+ 2ϵx~i(t)(dWi(t)x~i(t)jx~j(t)dWj(t))].\displaystyle+\ 2\frac{\epsilon}{\sqrt{\tilde{x}_{i}(t)}}\Big{(}dW_{i}(t)-\sqrt{\tilde{x}_{i}(t)}\sum_{j}\sqrt{\tilde{x}_{j}(t)}dW_{j}(t)\Big{)}\Big{]}.

This will be the stochastic replicator dynamics we analyze. Note that the diffusion coefficients in

x~i(t)dWi(t)x~i(t)jx~j(t)dWj(t)\sqrt{\tilde{x}_{i}(t)}dW_{i}(t)-\tilde{x}_{i}(t)\sum_{j}\sqrt{\tilde{x}_{j}(t)}dW_{j}(t)

are not Lipschitz near the boundary Pn\partial P^{n} of PnP^{n}, so the standard proof of well-posedness of stochastic differential equations with Lipschitz coefficients does not apply. Even the celebrated work of Watanabe and Yamada [66] that establishes well-posedness under weaker (e.g., certain Hölder continuity) conditions does not apply because of the presence of non-diagonal terms. In fact (9) is not well-posed - its solution is not unique as we shall see later. We shall therefore propose a selection principle to choose a specific solution as the ‘natural’ or ‘physical’ solution. For the time being, note that in the notation of (4),

Γi(x~(t))dW(t)=1x~i(t)(dWi(t)x~i(t)jx~j(t)dWj(t)),\Gamma_{i}(\tilde{x}(t))dW(t)=\frac{1}{\sqrt{\tilde{x}_{i}(t)}}\left(dW_{i}(t)-\sqrt{\tilde{x}_{i}(t)}\sum_{j}\sqrt{\tilde{x}_{j}(t)}dW_{j}(t)\right),

or,

Γij(x):=1xi(δijxixj),\Gamma_{ij}(x):=\frac{1}{\sqrt{x_{i}}}(\delta_{ij}-\sqrt{x_{i}x_{j}}),

where δij\delta_{ij} is the Kronecker delta. Letting x:=[x1,,xn]T\sqrt{x}:=[\sqrt{x_{1}},\cdots,\sqrt{x_{n}}]^{T}, we have xTΓ(x)=xT(IxxT)=0x^{T}\Gamma(x)=\sqrt{x}^{T}(I-\sqrt{x}\sqrt{x}^{T})=0 for xPnx\in P^{n} as desired, because IxxTI-\sqrt{x}\sqrt{x}^{T} is the projection to the space (x)(\sqrt{x})^{\perp}.

We take up the ill-posedness issues and subsequent analysis of this dynamics in the next section.

As was already mentioned, our model is a particular instance of a general formulation in [27], but with different regularity conditions. We now compare it with other models somewhat similar in spirit. An important predecessor is the work of Fudenberg and Harris [29] mentioned earlier, followed by a series of works [19, 39, 35, 42]. In that model, they begin with a ‘raw’ (i.e., unnormalized) growth model for the populations with the payoffs perturbed by scaled Brownian noise, in the spirit of the classical derivation of the deterministic replicator dynamics [65]. This leads to the equation (6) modified by using geometric Brownian motion as follows:

dri(t)=ri(t)(Vi(x(t))dt+ϵidWi(t)),t0.dr_{i}(t)=r_{i}(t)(V_{i}(x(t))dt+\epsilon_{i}dW_{i}(t)),\quad t\geq 0.

When the Ito transformation is applied to obtain the s.d.e. for the population fractions xi(t):=ri(t)(jrj(t))1x_{i}(t):=r_{i}(t)(\sum_{j}r_{j}(t))^{-1}, this leads to

dxi(t)\displaystyle dx_{i}(t) =\displaystyle= xi(t)(Vi(x(t))j=1nxj(t)Vj(x(t)))dt\displaystyle x_{i}(t)\Big{(}V_{i}(x(t))-\sum_{j=1}^{n}x_{j}(t)V_{j}(x(t))\Big{)}dt
xi(t)(ϵi2xi(t)j=1nϵj2xj2(t))\displaystyle-x_{i}(t)\Big{(}\epsilon_{i}^{2}x_{i}(t)-\sum_{j=1}^{n}\epsilon_{j}^{2}x_{j}^{2}(t)\Big{)}
+xi(t)(ϵidWi(t)j=1nxj(t)ϵjdWj(t)).\displaystyle+x_{i}(t)\Big{(}\epsilon_{i}dW_{i}(t)-\sum_{j=1}^{n}x_{j}(t)\epsilon_{j}dW_{j}(t)\Big{)}.

Note that new terms appear in the deterministic part.

Recently, Mertikopoulos and Viossat [49] proposed to introduce directly the geometric-type Brownian motion in the s.d.e. for the population fractions, which resulted in the following version of the stochastic replicator dynamics:

dxi(t)\displaystyle dx_{i}(t) =\displaystyle= xi(t)(Vi(x(t))j=1nxj(t)Vj(x(t)))dt\displaystyle x_{i}(t)\Big{(}V_{i}(x(t))-\sum_{j=1}^{n}x_{j}(t)V_{j}(x(t))\Big{)}dt
+xi(t)(σi(x(t))dWi(t)j=1nxjσj(x(t))dWj(t))\displaystyle+x_{i}(t)\Big{(}\sigma_{i}(x(t))dW_{i}(t)-\sum_{j=1}^{n}x_{j}\sigma_{j}(x(t))dW_{j}(t)\Big{)}

with the diffusion coefficients σi(x(t)),i=1,,n\sigma_{i}(x(t)),i=1,...,n being Lipschitz. This is a well-posed stochastic differential equation unlike ours and while also a particular case of [27], it is qualitatively different in various aspects as will become apparent as we proceed. Most importantly, their model shares with classical models the feature that no species can become extinct (i.e., the corresponding component hit zero) in finite time, whereas ours can, and get reflected, implying fresh arrivals of the species into the system. Also, in [49] the authors address issues different from the ones we do.

3 A transformation of the state space and analysis of metastability

Our first objective is to overcome the difficulty in treating corners in [27] by using dynamics (9).

This is achieved through a change of variable y(t)=x(t)y(t)=\sqrt{x(t)}. For ySn:=y\in S^{n}:= the (n1)(n-1)-dimensional sphere, define F~(y)=F([y12,,yn2])\tilde{F}(y)=F([y_{1}^{2},\cdots,y_{n}^{2}]). For the deterministic o.d.e. (3), this change of variable leads to

y˙(t)=14(F~(y(t))F~(y(t)),y(t)y(t))=14F~(y(t)),\dot{y}(t)=\frac{1}{4}\left(\nabla\tilde{F}(y(t))-\langle\nabla\tilde{F}(y(t)),\ y(t)\rangle y(t)\right)=\frac{1}{4}\nabla^{*}\tilde{F}(y(t)), (10)

where \nabla^{*} denotes the projected gradient operator that gives gradient in the tangent space of SnS^{n}. The new dynamics (10) is then nothing but gradient ascent on SnS^{n}, a compact manifold without boundary.

This transformation goes back to [3], [59], see [56] for a recent treatment which also discusses other related projection dynamics. The dynamics (10) as defined above is well defined on the entire SnS^{n}, not just on its positive orthant. Setting xi=yi2,i,x_{i}=y_{i}^{2},\ \forall i, establishes a one-to-one correspondence of each orthant of SnS^{n} with PnP^{n}. Note that all local maxima, resp. minima of FF in PnP^{n} get mapped to local maxima, resp. minima of F~\tilde{F} in a one to many map. At the boundary of the orthants, the vector fields on either side will exhibit a mirror symmetry with respect to this boundary because of the nature of our transformation, in particular, the component of the driving vector field 14F~(y(t))\frac{1}{4}\nabla^{*}\tilde{F}(y(t)) normal to the boundary will be symmetric across the boundary and being continuous, must vanish at the boundary. As observed earlier, the replicator dynamics in a face of PnP^{n} does not reach a lower dimensional face in finite time. The same will be true for the gradient ascent above vis-a-vis boundaries of orthants, because the normal component of the above vector field vanishes at the boundary, slowing down its approach to the boundary. This, however, is not the case for the stochastic version as we shall see later.

Consider now the process y~()\tilde{y}(\cdot) satisfying

dy~(t)\displaystyle d\tilde{y}(t) =\displaystyle= 14(F~(y~(t))F~(y~(t)),y~(t)y~(t))dt\displaystyle\frac{1}{4}\left(\nabla\tilde{F}(\tilde{y}(t))-\langle\nabla\tilde{F}(\tilde{y}(t)),\ \tilde{y}(t)\rangle\tilde{y}(t)\right)dt (11)
+ϵ(Iy~(t)y~(t)T)dW(t)\displaystyle+\ \epsilon(I-\tilde{y}(t)\tilde{y}(t)^{T})\circ dW(t)
=\displaystyle= 14F~(y~(t))dt+ϵdW~(t),\displaystyle\frac{1}{4}\nabla^{*}\tilde{F}(\tilde{y}(t))dt+\epsilon d\widetilde{W}(t), (12)

where ‘\circ’ denotes the Stratonovich integral and

W~(t):=0t(Iy~(s)y~(s)T)𝑑W(s),t0,\widetilde{W}(t):=\int_{0}^{t}(I-\tilde{y}(s)\tilde{y}(s)^{T})\circ dW(s),\ t\geq 0,

is a Brownian motion in SnS^{n}. (This is a consequence of the fact that IyyTI-yy^{T} for ySny\in S^{n} is a projection to the tangent space of SnS^{n} at yy. Note that this is an extrinsic description of W~()\widetilde{W}(\cdot) which requires that SnS^{n} be viewed as a manifold embedded in n\mathcal{R}^{n}.) The state space of y~()\tilde{y}(\cdot) is the entire SnS^{n}, with full support because as a diffusion in SnS^{n}, y~()\tilde{y}(\cdot) is non-degenerate. In particular, unlike the deterministic case, the trajectories do cross over from one orthant to another because of the Brownian component: at a point on the boundary, the normal component of the gradient drift is zero, so the normal component of the diffusion is a pure Brownian motion, implying that it will be on either side of the boundary with probability one in an arbitrarily small time interval.

Lemma 3.1 We can equivalently write (11) as

dy~(t)\displaystyle d\tilde{y}(t) =\displaystyle= 14(F~(y~(t))F~(y~(t)),y~(t)y~(t))dt\displaystyle\frac{1}{4}\left(\nabla\tilde{F}(\tilde{y}(t))-\langle\nabla\tilde{F}(\tilde{y}(t)),\ \tilde{y}(t)\rangle\tilde{y}(t)\right)dt (13)
+ϵ(Iy~(t)y~(t)T)dW(t).\displaystyle+\ \epsilon(I-\tilde{y}(t)\tilde{y}(t)^{T})dW(t).

Proof: The Ito to Stratonovich conversion of a vector s.d.e.

dX(t)=m(X(t))dt+σ(X(t))dW(t),dX(t)=m(X(t))dt+\sigma(X(t))dW(t),

written componentwise as

dXi(t)=mi(X(t))dt+jσij(X(t))dWj(t),dX_{i}(t)=m_{i}(X(t))dt+\sum_{j}\sigma_{ij}(X(t))dW_{j}(t),

goes as follows. For any twice continuously differentiable f:nf:\mathcal{R}^{n}\mapsto\mathcal{R}, we have by Ito formula,

df(X(t))\displaystyle df(X(t)) =\displaystyle= kfxk(X(t))(mk(X(t))dt+σk(X(t))dW(t))+\displaystyle\sum_{k}\frac{\partial f}{\partial x_{k}}(X(t))\Big{(}m_{k}(X(t))dt+\sum_{\ell}\sigma_{k\ell}(X(t))dW_{\ell}(t)\Big{)}+
12k,,mσkm(X(t))σm(X(t))2fxkx(X(t))dt\displaystyle\frac{1}{2}\sum_{k,\ell,m}\sigma_{km}(X(t))\sigma_{\ell m}(X(t))\frac{\partial^{2}f}{\partial x_{k}\partial x_{\ell}}(X(t))dt

and

df(X()),W(t)=kfxk(X(t))σk(X(t))dt.d\langle f(X(\cdot)),W_{\ell}\rangle(t)\ =\ \sum_{k}\frac{\partial f}{\partial x_{k}}(X(t))\sigma_{k\ell}(X(t))dt.

Hence

σij(X(t))dWj(t)\displaystyle\sigma_{ij}(X(t))\circ dW_{j}(t) =\displaystyle= σij(X(t))dWj(t)+12dσij(X()),Wj(t)\displaystyle\sigma_{ij}(X(t))dW_{j}(t)+\frac{1}{2}d\langle\sigma_{ij}(X(\cdot)),W_{j}\rangle(t)
=\displaystyle= σij(X(t))dWj(t)+12kσijxk(X(t))σkj(X(t))dt.\displaystyle\sigma_{ij}(X(t))dW_{j}(t)+\frac{1}{2}\sum_{k}\frac{\partial\sigma_{ij}}{\partial x_{k}}(X(t))\sigma_{kj}(X(t))dt.

Applying this to (11), we have

(Iy~(t)y~(t)T)dW~(t)=(Iy~(t)y~(t)T)dW~(t),(I-\tilde{y}(t)\tilde{y}(t)^{T})\circ d\widetilde{W}(t)=(I-\tilde{y}(t)\tilde{y}(t)^{T})d\widetilde{W}(t), (14)

which proves the claim. \Box

Remark: Ito and Stratonovich444The work [42] suggests the use of Stratonovich formulation for a different reason, viz., that it correctly captures the limiting behavior of ordinary differential equations with rapidly varying random right hand side. An older work [63], however, is critical of Stratonovich formulation because it misses out on some biologically relevant or realistic aspects, but as our main motivation is not from biology, we do not address these issues. See [37] for general background on diffusions on manifolds. integrals are certainly not always identical and there can be an O(ϵ2)O(\epsilon^{2}) error term. The lemma states that this term is zero for us. We emphasize this need not always be the case. The reason it works out for us is that our stochastic gradient dynamics on sphere is driven by a pure Brownian motion for which the diffusion coefficient is identity and therefore there is no Stratonovich correction. Our aim in stating (11) in Stratonovich form is purely to avoid technicalities in handling the Ito version (13) in the manifold framework. The advantage of Stratonovich formulation is that it is more amenable to smooth coordinate transformations: even though we eventually work on the probability simplex, a compact manifold with boundary that admits a single local chart for the entire manifold, our passage to it is via a stochastic differential equation on a sphere, a compact manifold without boundary that does not admit a single local chart.


We may identify a vector [y1,,yn]T[y_{1},\cdots,y_{n}]^{T} in SnS^{n} with [|y1|,,|yn|]T[|y_{1}|,\cdots,|y_{n}|]^{T} in its positive orthant. Under this equivalence relation, the quotient space is the positive orthant S+nS^{n}_{+} of SnS^{n} and y~()\tilde{y}(\cdot) maps to a reflected diffusion yˇ()\check{y}(\cdot) in S+nS^{n}_{+}. This requires that we modify (9)-(9) to include a boundary term to account for the reflection, specifically a bounded variation increasing process supported on the boundary that confines the original process to the prescribed domain. We also need to specify a direction of reflection. This is because the Brownian and hence the diffusion trajectory is non-differentiable everywhere with probability one, therefore its direction of reflection on the boundary cannot be defined in a conventional manner. As in the case of simple reflected Brownian motion, pure reflection across a codimension one smooth surface entails normal reflection, as can be justified also by symmetry considerations: yi-y_{i}^{\prime} is the mirror image of yiy_{i}^{\prime} across the plane yi=0y_{i}=0 and the line joining the two is normal to this plane. The other specification needed is whether the reflection is instantaneous or ‘sticky’. Our definition of yˇ()\check{y}(\cdot) implies instantaneous reflection. Both these specifications are achieved by adding to the dynamics a bounded variation process which confines the process to the desired domain while respecting the prescribed reflection direction. We describe this in detail below. This material is standard and a nice exposition can be found in [38]. Thus, if yˇ()\check{y}(\cdot) is a reflected diffusion with normal reflection in S+nS^{n}_{+}, it satisfies

dyˇ(t)=14F~(yˇ(t))dt+ϵdW~(t)+dZ(t),d\check{y}(t)=\frac{1}{4}\nabla^{*}\tilde{F}(\check{y}(t))dt+\epsilon d\widetilde{W}(t)+dZ(t), (15)

where Z()=[Z1(),,Zn()]TZ(\cdot)=[Z_{1}(\cdot),\cdots,Z_{n}(\cdot)]^{T} is a process (called the ‘local time’ at the boundary) satisfying:

  1. 1.

    Z(0)=Z(0)= the zero vector;

  2. 2.

    If |Zi|(T)|Z_{i}|(T) denotes the total variation of Zi(t)Z_{i}(t) on [0,T][0,T], then |Zi|(t)<t0|Z_{i}|(t)<\infty\ \forall t\geq 0 (i.e., it is a bounded variation process) and

    |Zi|(t)=0tI{yˇi(s)=0}d|Z|(s),t0|Z_{i}|(t)=\int_{0}^{t}I\{\check{y}_{i}(s)=0\}d|Z|(s),\ t\geq 0

    (i.e., it increases only when yˇi(t)=0\check{y}_{i}(t)=0, implying instantaneous reflection);

  3. 3.

    For yS+ny\in\partial S_{+}^{n}, let C(y):={C(y):=\{ the inward normal }\} if yy belongs to one of the primary faces of S+nS^{n}_{+}, and :=:= the cone formed by the inward normals to the adjoining faces if yy belongs to the intersection of two or more primary faces of S+nS^{n}_{+}. Then there exists a measurable γ:[0,)n\gamma:[0,\infty)\mapsto\mathcal{R}^{n} such that γ(t)C(yˇ(t)),t0,\gamma(t)\in C(\check{y}(t)),\ \forall\ t\geq 0, and

    Z(t)=0tγ(yˇ(s))d|Z|(s).Z(t)=\int_{0}^{t}\gamma(\check{y}(s))d|Z|(s).

    (This specifies the direction of reflection.)

Theorem 3.1 The stochastic differential equation (15) has a unique strong solution.

Proof: Consider the stochastic differential equation in the convex set

𝒞:={y=[y1,,yn]n:yi0i,y1η}\mathcal{C}:=\{y=[y_{1},\cdots,y_{n}]\in\mathcal{R}^{n}:\ y_{i}\geq 0\ \forall i,\ \|y\|_{1}\geq\eta\}

for some η(0,1)\eta\in(0,1), given by

dy˘(t)\displaystyle d\breve{y}(t) =\displaystyle= 14(F~(y˘(t))F~(y˘(t)),y˘(t)y˘(t)y˘(t)y˘(t))dt\displaystyle\frac{1}{4}\left(\nabla\tilde{F}(\breve{y}(t))-\left\langle\nabla\tilde{F}(\breve{y}(t)),\ \frac{\breve{y}(t)}{\|\breve{y}(t)\|}\right\rangle\frac{\breve{y}(t)}{\|\breve{y}(t)\|}\right)dt (16)
+ϵ(Iy˘(t)y˘(t)Ty˘(t)2)dW(t)+dZ(t),\displaystyle+\ \epsilon\left(I-\frac{\breve{y}(t)\breve{y}(t)^{T}}{\|\breve{y}(t)\|^{2}}\right)\circ dW(t)+dZ(t),

with Z()Z(\cdot) defined analogously to the above. This is a degenerate reflected diffusion in 𝒞\mathcal{C} with normal reflection, with the property that it stays a.s. on the spherical shell 𝒞{y:y=y˘(0)}\mathcal{C}\cap\{y:\|y\|=\|\breve{y}(0)\|\} whenever y˘(0)>η\|\breve{y}(0)\|>\eta. In particular, if y˘(0)S+n\breve{y}(0)\in S^{n}_{+}, it remains in S+nS^{n}_{+} and reduces to yˇ()\check{y}(\cdot). By Theorem 4.1, p. 175, of [60], the claim follows for (16) and ipso facto for (15), since η>0\eta>0 is arbitrary. \Box

Recall that for an ill-posed diffusion, a selection principle involves picking one solution law per initial condition so that collectively they satisfy the Chapman-Kolmogorov conditions, see, e.g., the classical work of Krylov [44].


Theorem 3.2 Setting

x~(t):=yˇ(t)2\tilde{x}(t):=\check{y}(t)^{2} (17)

gives a unique selection for a strong Markov solution to the ill-posed s.d.e. (9)-(9).

Proof: By Ito formula,

dx~i(t)\displaystyle d\tilde{x}_{i}(t) =\displaystyle= x~i(t)(Vi(x~(t))j=1nx~j(t)Vj(x~(t)))dt\displaystyle\tilde{x}_{i}(t)\Big{(}V_{i}(\tilde{x}(t))-\sum_{j=1}^{n}\tilde{x}_{j}(t)V_{j}(\tilde{x}(t))\Big{)}dt
+ 2ϵ(x~i(t)dWi(t)x~i(t)j=1nx~j(t))dWj(t)\displaystyle+\ 2\epsilon\Big{(}\sqrt{\tilde{x}_{i}(t)}dW_{i}(t)-\tilde{x}_{i}(t)\sum_{j=1}^{n}\sqrt{\tilde{x}_{j}(t)}\Big{)}dW_{j}(t)
+ 2x~i(t)dZi(t),\displaystyle+\ 2\sqrt{\tilde{x}_{i}(t)}dZ_{i}(t),
=\displaystyle= x~i(t)(Vi(x~(t))j=1nx~j(t)Vj(x~(t)))dt\displaystyle\tilde{x}_{i}(t)\Big{(}V_{i}(\tilde{x}(t))-\sum_{j=1}^{n}\tilde{x}_{j}(t)V_{j}(\tilde{x}(t))\Big{)}dt
+ 2ϵ(x~i(t)dWi(t)x~i(t)j=1nx~j(t))dWj(t)\displaystyle+\ 2\epsilon\Big{(}\sqrt{\tilde{x}_{i}(t)}dW_{i}(t)-\tilde{x}_{i}(t)\sum_{j=1}^{n}\sqrt{\tilde{x}_{j}(t)}\Big{)}dW_{j}(t)
=\displaystyle= x~i(t)[(Fx~i(x~(t))jx~j(t)Fx~j(x~(t)))dt\displaystyle\tilde{x}_{i}(t)\Big{[}\Big{(}\frac{\partial F}{\partial\tilde{x}_{i}}(\tilde{x}(t))-\sum_{j}\tilde{x}_{j}(t)\frac{\partial F}{\partial\tilde{x}_{j}}(\tilde{x}(t))\Big{)}dt
+ 2ϵx~i(t)(dWi(t)x~i(t)jx~j(t)dWj(t))],\displaystyle+\ 2\frac{\epsilon}{\sqrt{\tilde{x}_{i}(t)}}\Big{(}dW_{i}(t)-\sqrt{\tilde{x}_{i}(t)}\sum_{j}\sqrt{\tilde{x}_{j}(t)}dW_{j}(t)\Big{)}\Big{]},

for t0t\geq 0, where we have used the fact that 2x~i(t)dZi(t)02\sqrt{\tilde{x}_{i}(t)}dZ_{i}(t)\equiv 0 because Zi(t)Z_{i}(t) increases only on the set {t:x~i(t)=0}\{t:\tilde{x}_{i}(t)=0\}. Thus x~()\tilde{x}(\cdot) defined by (17) satisfies (9) and (9). Furthermore, by Theorem 3.1, (17) uniquely specifies x~()\tilde{x}(\cdot). On the other hand, y~()\tilde{y}(\cdot) hits the boundary of Sn+S^{+}_{n} and therefore x~()\tilde{x}(\cdot) must hit the boundary of PnP^{n}. Then any solution that remains confined to a face of PnP^{n} after hitting the same with some prescribed probability >0>0 is also a solution. This shows that the above x~()\tilde{x}(\cdot) is not a unique solution. In other words, (9) and (9) are ill-posed. Hence (17) implies a genuine selection process. The strong Markov property follows from that of the well-posed reflected diffusion yˇ()\check{y}(\cdot) because the two are interconvertible by a continuous invertible transformation S+nPnS^{n}_{+}\leftrightarrow P_{n}. \Box

Remarks 1. This interconvertibility also implies that x~()\tilde{x}(\cdot) is adapted to the increasing filtration of σ\sigma-fields generated by the driving Brownian motion (because y~()\tilde{y}(\cdot) is) and since it has continuous paths, is also predictable with respect to it. Thus it remains amenable to the use of standard tools of Ito calculus such as the Ito formula, a fact we shall use implicitly throughout what follows.

2. Observe that by the well-posedness of (15) in the strong and therefore weak sense, yˇ()\check{y}(\cdot) is also the unique limit in law as 0<η00<\eta\downarrow 0 of the well-posed non-degenerate reflected diffusion yˇη()\check{y}^{\eta}(\cdot) in (n)+(\mathcal{R}^{n})^{+} given by

dyˇη(t)=14F~(yˇη(t))dt+ϵdW~(t)+ηdW^(t)+dZ(t),d\check{y}^{\eta}(t)=\frac{1}{4}\nabla^{*}\tilde{F}(\check{y}^{\eta}(t))dt+\epsilon d\widetilde{W}(t)+\eta d\widehat{W}(t)+dZ(t),

where W^()\widehat{W}(\cdot) is an independent Brownian motion in n\mathcal{R}^{n} and η>0\eta>0. Setting x~iη(t)=yˇiη(t)2i\tilde{x}_{i}^{\eta}(t)=\check{y}_{i}^{\eta}(t)^{2}\ \forall i, x~η()\tilde{x}^{\eta}(\cdot) satisfies the stochastic differential equation in (n)+(\mathcal{R}^{n})^{+} given by

dx~iη(t)\displaystyle d\tilde{x}^{\eta}_{i}(t) =\displaystyle= x~iη(t)(Vi(x~η(t))j=1nx~jη(t)Vj(x~η(t)))dt+ 2(ϵx~iη(t)dWi(t)+\displaystyle\tilde{x}^{\eta}_{i}(t)\Big{(}V_{i}(\tilde{x}^{\eta}(t))-\sum_{j=1}^{n}\tilde{x}^{\eta}_{j}(t)V_{j}(\tilde{x}^{\eta}(t))\Big{)}dt+\ 2\Big{(}\epsilon\sqrt{\tilde{x}^{\eta}_{i}(t)}dW_{i}(t)+
ηx~iη(t)dW^i(t)ϵx~iη(t)j=1nx~jη(t)dWj(t))+η2dt.\displaystyle\eta\sqrt{\tilde{x}^{\eta}_{i}(t)}d\widehat{W}_{i}(t)-\epsilon\tilde{x}^{\eta}_{i}(t)\sum_{j=1}^{n}\sqrt{\tilde{x}^{\eta}_{j}(t)}dW_{j}(t)\Big{)}+\eta^{2}dt.

This is indeed well posed as an (n)+(\mathcal{R}^{n})^{+}-valued diffusion in the sense of Stroock-Varadhan martingale problem, in particular it has a unique weak solution [7]. Furthermore, we recover x~()\tilde{x}(\cdot) as the unique limit in law thereof as η0\eta\downarrow 0 (because the corresponding limit in law of yˇη()\check{y}^{\eta}(\cdot) is unique). This selection is along the lines of the selection principle enunciated in [17], but cannot be directly deduced from the results of [17] because (9) and (9) do not satisfy the conditions stipulated therein. Note also that the selection procedure employed above is in the spirit of the Kolmogorov philosophy of ‘selection through small noise limit’ alluded to in the introduction and leads to a unique choice, unlike the procedure of [44] where the result depends on the choice and ordering of the ‘test functionals’ employed. It is also worth noting that a key condition, condition (1.4) of [7], fails here for η=0\eta=0 and sure enough, there is no uniqueness as already observed.


We also lose the phenomenon of extinction in this framework. In fact, our model describes a community of fixed size where different sub-communities wax or wane in their relative share depending on the selection mechanism, but do not become extinct - they can bounce back from zero population.


We next consider the asymptotic behavior of the above mentioned processes.


Theorem 3.3 The process {y~(t),t[0,)}\{\tilde{y}(t),t\in[0,\infty)\} is ergodic with stationary distribution

ξϵ(dy)=φϵ(y)κ(dy):=Z1eF~(y)8ϵ2κ(dy)\xi^{\epsilon}(dy)=\varphi^{\epsilon}(y)\kappa(dy):=Z^{-1}e^{\frac{\tilde{F}(y)}{8\epsilon^{2}}}\kappa(dy) (18)

for κ:=\kappa:= the uniform surface measure on SnS^{n}, where ZZ is the normalizing factor

Z:=eF~(y)8ϵ2κ(dy).Z:=\int e^{\frac{\tilde{F}(y^{\prime})}{8\epsilon^{2}}}\kappa(dy^{\prime}).

Proof: Note that ξϵ(dy)\xi^{\epsilon}(dy) is the Gibbs distribution with potential F~-\tilde{F}. Its invariance is easily established by verifying that the stationary Fokker-Planck equation φ0\mathcal{L}^{*}\varphi\equiv 0, where

:=14F~+ϵ22ΔS\mathcal{L}:=\frac{1}{4}\nabla\widetilde{F}\cdot\nabla+\frac{\epsilon^{2}}{2}\Delta_{S}

is the generator of the diffusion y~()\tilde{y}(\cdot) and \mathcal{L}^{*} is its formal adjoint, given by

(f):=14(fF~1)+ϵ22ΔSf\mathcal{L}^{*}(f):=-\frac{1}{4}\nabla(f\nabla\widetilde{F}\cdot\textbf{1})+\frac{\epsilon^{2}}{2}\Delta_{S}f

for fC(Sn)f\in C_{\infty}(S^{n}), 1:=[1,1,,1]T\textbf{1}:=[1,1,\cdots,1]^{T}. (ΔS\Delta_{S} denotes the Laplace-Beltrami operator on SnS^{n}.) Ergodicity follows by classical arguments: As a non-degenerate diffusion on SnS^{n}, y~()\tilde{y}(\cdot) has for each t>0t>0 a transition probability that is absolutely continuous with respect to κ(dy)\kappa(dy) with strictly positive density (say) p(y|x,t)p(y|x,t). Then every invariant distribution η(dy)\eta(dy) must satisfy

p(y|z,t)η(dz)κ(dy)=η(dy)\int p(y|z,t)\eta(dz)\kappa(dy)=\eta(dy)

by invariance. It follows that any two candidate η()\eta(\cdot) are mutually absolutely continuous w.r.t. κ\kappa and therefore w.r.t.  each other. From the Doeblin decomposition from ergodic theory of Markov processes, we know that extremal invariant measures are singular with respect to each other. The foregoing then implies that there is only one invariant measure. \Box

For an extensive treatment of invariant measures of reflected diffusions, see [41]. In particular, we have

Corollary 3.1 The process x~()\tilde{x}(\cdot) has a stationary distribution supported on the interior of PnP^{n} given by

Z^1(m=1n1xm)eF(x)8ϵ2ζ(dx),\hat{Z}^{-1}\left(\prod_{m=1}^{n}\frac{1}{\sqrt{x_{m}}}\right)e^{\frac{F(x)}{8\epsilon^{2}}}\zeta(dx), (19)

with ζ:=\zeta:= the normalized Lebesgue measure on PnP^{n}, Z^\hat{Z} being the normalizing factor. In particular, as ϵ0\epsilon\downarrow 0, the stationary distribution concentrates on the global maxima of FF, rendering them stochastically stable.

Proof: The first claim follows by an application of the change of variables formula to (18). The second claim then follows by standard arguments. \Box

Note that we proved ergodicity for the process y~()\tilde{y}(\cdot). That of x~()\tilde{x}(\cdot) follows. The latter, however, is a statement regarding the specific selection. Other solutions remain, e.g., those restricted to the boundary, but we do not consider them. An important advantage of the explicit Gibbs distribution above is that it lays bare the metastable states as local maxima of FF and flags the global maxima as stochastically stable states.

We analyze next the important phenomenon of metastability. The advantage of taking the present route is twofold. First, working with a compact manifold without boundary, we are spared the technicalities due to a non-smooth boundary which occurred in [27]. Secondly, as in [27], we shall be seeking small noise asymptotics for concentration on the minima of the associated Freidlin-Wentzell potential ([28], Chapter 6), which for the above gradient ascent, is simply proportional to F-F. Hence for this special case, one can get a good handle on exit times from the domain of attraction of a stable equilibrium xx^{*} of (3):


Lemma 3.2 Let G~\tilde{G} be the domain of attraction of a stable equilibrium ySny^{*}\in S^{n} of (10) with boundary G~\partial\tilde{G} and let

τ~ϵ:=inf{t0:y~(t)G~},\tilde{\tau}_{\epsilon}:=\inf\{t\geq 0:\tilde{y}(t)\notin\tilde{G}\},

where y~(t)\tilde{y}(t) is given by o.d.e. (13). Then for yG~y\in\tilde{G},

limϵ0ϵ2logEy[τ~ϵ]=12minzG~(F~(y)F~(z)).\lim_{\epsilon\downarrow 0}\epsilon^{2}\log E_{y}\left[\tilde{\tau}_{\epsilon}\right]=\frac{1}{2}\min_{z\in\partial\tilde{G}}(\tilde{F}(y^{*})-\tilde{F}(z)). (20)

Proof: Let G~0:={yG~:F~(y)F~0:=maxyG~F~(y)}\tilde{G}_{0}:=\{y\in\tilde{G}:\tilde{F}(y)\geq\tilde{F}_{0}:=\max_{y^{\prime}\in\partial\tilde{G}}\tilde{F}(y^{\prime})\}. Then it is easy to see that supG~F~=supG~0F~\sup_{\tilde{G}}\tilde{F}=\sup_{\tilde{G}_{0}}\tilde{F}. By combining Theorem 3.1, p. 100, and Theorem 4.1, p. 106, of [28], we then have

limϵ0ϵ2logEx[τ~ϵ]=12minyG~(F~(y)F~(y)).\lim_{\epsilon\downarrow 0}\epsilon^{2}\log E_{x}\left[\tilde{\tau}_{\epsilon}\right]=\frac{1}{2}\min_{y\in\partial\tilde{G}}(\tilde{F}(y^{*})-\tilde{F}(y)).

The claim follows. \Box

Theorem 3.4 If GG is the domain of attraction of a stable equilibrium xx^{*} (possibly on the boundary) of (3) and let

τϵ:=inf{t0:x~(t)G}.\tau_{\epsilon}:=\inf\{t\geq 0:\tilde{x}(t)\notin G\}.

where x~(t)\tilde{x}(t) is given by (17). Then for xGx\in G,

limϵ0ϵ2logEx[τϵ]=12minzG(F(x)F(z)).\lim_{\epsilon\downarrow 0}\epsilon^{2}\log E_{x}\left[\tau_{\epsilon}\right]=\frac{1}{2}\min_{z\in\partial G}(F(x^{*})-F(z)). (21)

Proof: In Lemma 3.2, let y=xy^{*}=\sqrt{x^{*}}. The claim then follows from the relationship between FF and F~\tilde{F}. \Box


It should be emphasized that this captures the dominant part of the small noise asymptotics for mean exit times. A finer analysis is possible by combining Theorem 11.2, p. 267; Proposition 11.7, p. 275, and Theorem 11.12, p. 280, of [18].

Note also that we are able to make a statement about boundary equilibria for which a general theory appears unavailable. We are able to do so because we first derive our results for a process on a manifold without boundary and then deduce the result for the desired process by a change of variables.

Following [28], Chapter 6, we can infer considerable additional information regarding the behavior of the process for small ϵ\epsilon. For example, we know that with probability approaching 11 as ϵ0\epsilon\downarrow 0, the exit from the domain of attraction of xx^{*} occurs near a minimizer of the r.h.s. in (21). This allows us to consider possible paths from xx^{*} to a ‘ground state’ or global minimum of the potential via a succession of stable equilibria, subject to the transition between two consecutive ones being compatible with the foregoing. By Theorem 4.2, Chapter 4 of [28], we also know that the exit time distribution for each successive jump is approximately exponential. Furthermore, the process behaves like a continuous time Markov chain on the set of metastable states, i.e., the stable equilibria. Then each path is traversed with time approximately distributed as sums of independent (by the strong Markov property) approximately exponential random variables, of which the smallest exponential rate will dominate in the ϵ0\epsilon\downarrow 0 limit. Thus the transit time from xx^{*} to a ground state will be approximately exponential with exponential rate equal to the largest, over the aforementioned paths, of the smallest exponential rate along the path. This is akin to the considerations that characterize the optimal constant for logarithmic cooling in simulated annealing [33].

We conclude this section with a remark on quasi-stationary distributions. In the small noise regime, the process may spend a long time in the domain of attraction GG of a stable equilibrium, and thereby be well approximated by another process confined to G¯\bar{G}. The stationary distribution of the latter then well approximates the behavior of the original process in the time scales of interest. In discrete Markov chain case, several notions of quasi-stationary distributions exist, but they agree in the small noise limit [6]. In case of diffusions, one of these notions is the most convenient to work with: the limit in law as TT\uparrow\infty of the conditional law of the process conditioned on never exiting the region of interest till TT. This is the so called ‘Q-process’ of [24], Chapter 7. As shown in [53], this turns out to be the stationary distribution of another diffusion which lives in G¯\bar{G} and whose extended generator can be exactly characterized, viz., it is given as the unique minimizer of the (convex, lower semi-continuous) Donsker-Varadhan rate function for large deviations of empirical measures of the original diffusion from its stationary distribution, when the minimization is performed over probability measures supported in G¯\bar{G}. The results of [53] are for diffusions in d\mathcal{R}^{d}, but similar results will hold for diffusions on a sphere. A recent work [22] provides an alternative route to the Q-process for diffusions which is more explicit. Translated into our framework, their Theorems 3.1-3.2 state that the extended generator of the Q-process is given by

˘f=(ϕf)ϕ+λ0f,\breve{\mathcal{L}}f=\frac{\mathcal{L}(\phi f)}{\phi}+\lambda_{0}f,

where \mathcal{L} is the extended generator of the original diffusion (in our case, (12)) and ϕ,λ0\phi,\lambda_{0} are respectively the principal eigenfunction and eigenvalue associated with -\mathcal{L} with Dirichlet boundary condition on G\partial G. In particular, a simple application of Girsanov transformation shows that this amounts to changing the drift of (12) by an additional term

ϵ2logϕ.\epsilon^{2}\nabla^{*}\log\phi.

We summarize this below as a theorem.

Theorem 3.5 The law of the Q-process corresponds to the diffusion in GG given by

dy~(t)=14F~(y~(t))dt+ϵ2logϕ(y~(t))+ϵdW~(t),d\tilde{y}(t)=\frac{1}{4}\nabla^{*}\tilde{F}(\tilde{y}(t))dt+\epsilon^{2}\nabla^{*}\log\phi(\tilde{y}(t))+\epsilon d\widetilde{W}(t),

where ϕ\phi is the principal eigenfunction of the operator -\mathcal{L} above with Dirichlet boundary condition on G\partial G.

4 Intermittent dynamics on large graphs

We next describe our general results for the case of quadratic potential F(x):=12xTMxF(x):=\frac{1}{2}x^{T}Mx with xx belonging to a high dimensional space and MM representing a network structure. We know at least two important applications motivating such a setting. The first motivation comes from the Tangled Nature model [23] or its high level abstraction based on replicator dynamics [54, 62]. These models have been proposed to explain the intermittent behaviour of complex ecological systems as well as sudden, fundamental changes in economic systems. The authors of [54, 62] considered a discretized version of the standard replicator dynamics with matrix MM representing a large weighted random graph. They add to the standard replicator dynamics the processes of extinction and mutation. The extinction happens when the population density of a species drops below a certain threshold. Also, at each time step, a fraction of a species can mutate to some other species with a certain probability distribution, which can be tuned to encourage or discourage mutation to similar species. We feel that our continuous time model share important features with the models of [23, 54, 62]. As in [23, 54, 62], we too observe in our model an intermittent behaviour when a state stable for relatively long time changes suddenly to a new state, which can be quite different in terms of its set of strategies that have a significant population share relative to the noise level. As in [23, 54, 62], our model is not confined to the interior of the simplex. Also, as in the above cited models, the metastable states in our model have only a small number of strategies with significant population share. The reader will be able to see these features in the numerical examples described in the ensuing section.

Our second motivating application comes from the series of works [10, 11, 12] on the application of the replicator dynamics to Maximum Clique Detection Problem (MCP). The goal of MCP is to find in a graph a clique, i.e., a complete subgraph, with maximum size. Let AA be the adjacency matrix of a graph G=(V,E)G=(V,E), with VV as the vertex set, |V|=n|V|=n, and EE as the edge set, |E|=m|E|=m. In other words, AA is an n×nn\times n symmetric matrix whose (i,j)(i,j)th element =1=1 if (i,j)E(i,j)\in E and =0=0, otherwise. We assume that there are no self-loops, i.e., Aii=0,iA_{ii}=0,\ \forall i. If we take M=A+12IM=A+\frac{1}{2}I, the strict local maxima of F(x)F(x) over the simplex and therefore Evolutionarily Stable Strategies (ESS), will correspond to maximal cliques, see e.g. [10, 11, 12]. Specifically, if CC is a maximal clique (i.e., not a proper subset of another clique) and xCx_{C} is an associated characteristic vector with mass 1/|C|1/|C| at each vertex belonging to CC and with zero mass at all other vertices, then FF defined over the simplex achieves a local maximum at xCx_{C}. Clearly, FF achieves the global maximum over the simplex at the characteristic vector corresponding to the maximum clique (the largest clique in the graph.)

Note that maximal cliques will typically correspond to boundary equilibria, for which our model is particularly suited, because they become interior equilibria when the dynamics is lifted to the sphere. In general we expect the analysis of such equilibria to be much harder in models where the boundary is not reached in finite time.

Now let us provide an upper bound on the expected intermittence times by specializing further the results of Theorem 3.4. Let x=xCx^{*}=x_{C} in Theorem 3.4 with GG its domain of attraction and τϵ\tau_{\epsilon} defined correspondingly.

Theorem 4.1 For any graph we have the following bound on the asymptotics of the expected exit times from metastable states

limϵ0ϵ2logEx[τϵ]12[(112|C|)14n].\lim_{\epsilon\downarrow 0}\epsilon^{2}\log E_{x}\left[\tau_{\epsilon}\right]\leq\frac{1}{2}\left[\left(1-\frac{1}{2|C|}\right)-\frac{1}{4n}\right].

Proof: Recall from [13] the following lower bound for the quadratic function F(x)F(x) over the simplex:

F(x)12(i=1nmii1)1.F(x)\geq\frac{1}{2}\left(\sum_{i=1}^{n}m_{ii}^{-1}\right)^{-1}.

Since in our case mii=1/2m_{ii}=1/2, we obtain

F(x)14n.F(x)\geq\frac{1}{4n}.

Using the above bound, we can write

limϵ0ϵ2logEx[τϵ]=12minyG(F(x)F(y))12[(112|C|)14n],\lim_{\epsilon\downarrow 0}\epsilon^{2}\log E_{x}\left[\tau_{\epsilon}\right]=\frac{1}{2}\min_{y\in\partial G}(F(x^{*})-F(y))\leq\frac{1}{2}\left[\left(1-\frac{1}{2|C|}\right)-\frac{1}{4n}\right], (22)

where CC is the clique associated with a local maximum. The claim follows. \Box

This suggests that the maximum mean exit time from the domain of attraction of the local maxima grows very slowly with nn, implying that while the number of local maxima increases, the sizes of the barrier heights for their domains of attraction do not increase drastically.

We can make the result even more precise if we consider a particular model of large graphs, e.g., the Erdős-Rényi G(n,p)G(n,p) random graph model [9]. Recall that in G(n,p)G(n,p) random graph of size nn, the edges between any two vertices are drawn at random with probability pp.

Theorem 4.2 For sufficiently large nn, with high probability we have

limϵ0ϵ2logEx[τϵ]=12minyG(F(x)F(y))12[(1122log1/p(n))14n].\lim_{\epsilon\downarrow 0}\epsilon^{2}\log E_{x}\left[\tau_{\epsilon}\right]=\frac{1}{2}\min_{y\in\partial G}(F(x^{*})-F(y))\leq\frac{1}{2}\left[\left(1-\frac{1}{2\lceil 2\log_{1/p}(n)\rceil}\right)-\frac{1}{4n}\right].

Proof: It is known [9, Chapter 11] that in G(n,p)G(n,p) the maximum clique size is either 2log1/p(n)\lfloor 2\log_{1/p}(n)\rfloor or 2log1/p(n)\lceil 2\log_{1/p}(n)\rceil, with high probability. By high probability we mean that when nn goes to infinity, the probability of an event goes to one. The claim then is a simple rephrasing of (22) in view of preceding remarks. \Box

5 Numerical examples

Before proceeding to the description of two particular numerical examples, let us describe the overall discretization scheme that is used for numerical experiments:

Y^((n+1)Δ)\displaystyle\widehat{Y}((n+1)\Delta) =\displaystyle= Y(nΔ)\displaystyle Y(n\Delta)
Δ4{F~(Y(nΔ))Y(nΔ),F~(Y(nΔ))Y(nΔ)}\displaystyle-\ \frac{\Delta}{4}\Big{\{}\nabla\tilde{F}(Y(n\Delta))-\langle Y(n\Delta),\nabla\tilde{F}(Y(n\Delta))\rangle Y(n\Delta)\Big{\}}
+ϵΔ(IY(nΔ)Y(nΔ)T)(W((n+1)Δ)W(nΔ)),\displaystyle+\ \epsilon\sqrt{\Delta}(I-Y(n\Delta)Y(n\Delta)^{T})(W((n+1)\Delta)-W(n\Delta)),
Yi((n+1)Δ)\displaystyle Y_{i}((n+1)\Delta) =\displaystyle= Y^i((n+1)Δ)Y^((n+1)Δ),\displaystyle\frac{\widehat{Y}_{i}((n+1)\Delta)}{\|\widehat{Y}((n+1)\Delta)\|},
Xi((n+1)Δ)\displaystyle X_{i}((n+1)\Delta) =\displaystyle= Yi2((n+1)Δ).\displaystyle Y_{i}^{2}((n+1)\Delta).

Here the components of W((n+1)Δ)W(nΔ)W((n+1)\Delta)-W(n\Delta) are i.i.d. N(0,1)N(0,1) for all n0n\geq 0. The second step above is a renormalization step to ensure that the iterates remain on the unit sphere despite numerical approximations. This is an example of a retraction, a common operation for, e.g., optimization algorithms on an embedded manifold (see, e.g., [1], section 4.1).

Let us first consider the simplest non-trivial example with two edges connected in one node. Equivalently, it is described by the following interaction matrix:

M=[1/21011/21011/2].M=\left[\begin{array}[]{ccc}1/2&1&0\\ 1&1/2&1\\ 0&1&1/2\end{array}\right].

An example of system trajectory is presented in Figure 1 with a zoom of one of transitions between the metastable states shown in Figure 2.

Refer to caption
(a) The fraction of the first subpopulation, X1(t)X_{1}(t).
Refer to caption
(b) The fraction of the second subpopulation, X2(t)X_{2}(t).
Refer to caption
(c) The fraction of the third subpopulation, X3(t)X_{3}(t).
Figure 1: Sample of trajectory in two edge example. The x-axis corresponds to the time and the y-axis to the subpopulation fraction.
Refer to caption
Figure 2: Sample of trajectory in two edge example (zoom). The x-axis corresponds to the time and the y-axis to the subpopulation fractions.

The associated potential function F(x)=12xTMxF(x)=\frac{1}{2}x^{T}Mx has two global maxima at [1/2 1/2 0]T[1/2\ 1/2\ 0]^{T} and [0 1/2 1/2]T[0\ 1/2\ 1/2]^{T}, corresponding to the two maximum cliques. In game theoretic terminology, these are two evolutionarily stable strategies. In this simple example, two maximal, and also maximum, cliques are just the two edges.

At the maxima, the potential achieves the value of 3/8. It is easy to see from the symmetry of the problem and one-dimensional optimization that the easiest path to reach one maximum from the other goes through the point [1/5 3/5 1/5][1/5\ 3/5\ 1/5] with the value of the potential 7/20. Thus, in formula (21) we have that

minyG(F~(x)F~(y))=3/87/20.\min_{y\in\partial G}(\tilde{F}(\sqrt{x^{*}})-\tilde{F}(y))=3/8-7/20.

In Figure 3 we plot ϵ2logτ¯ϵexper\epsilon^{2}\log\bar{\tau}_{\epsilon}^{exper} vs (3/87/20)/2(3/8-7/20)/2, where τ¯ϵexper\bar{\tau}_{\epsilon}^{exper} is the average exit time obtained from experiments and ϵ\epsilon changes from 0.1 down to 0.05. We note that for values of ϵ\epsilon smaller than 0.05, the exit times become excessively large and we cannot collect enough reliable statistics. Figure 3 suggests that formula (21) describes correctly the asymptotics of the system for small values of ϵ\epsilon.

Refer to caption
Figure 3: Comparison of experiments with asymptotics. The x-axis corresponds to the value of ϵ\epsilon and the y-axis to the scaled logarithm of the average exit time.

We also plot in Figure 4 the empirical complementary cumulative distribution function of the exit times for ϵ=0.1\epsilon=0.1 in semi-log scale. An approximately linear shape of the curve is in agreement with our expectation that the exit times are exponentially distributed. In all numerical experiments of this section we took the step size Δ=0.05\Delta=0.05 and the number of steps between one and ten million depending on the value of ϵ\epsilon.

Refer to caption
Figure 4: Empirical complementary cumulative distribution function of the exit times. The x-axis corresponds to the value of the exit time.

Our second numerical example is based on a sample of Erdős-Rényi random graph with n=100n=100 vertices, with the edge probability p=0.25p=0.25, and with the noise level ϵ=0.02\epsilon=0.02 in the stochastic replicator dynamics. According to [9, Chapter 11], with high probability, the size of the maximum clique in Erdős-Rényi graph is either 2log1/p(n)\lfloor 2\log_{1/p}(n)\rfloor or 2log1/p(n)\lceil 2\log_{1/p}(n)\rceil. The substitution of n=100n=100, p=0.25p=0.25 into the latter formulas indicates that in our example the maximum clique size is either 6 or 7. Furthermore, using [9, Eq. (11.6)], we calculate that the expected number of maximal cliques is 1273 in G(100,0.25)G(100,0.25). The cliques of sizes 3 and 4 contribute most to this number. Thus, it is no surprise that most often we observe the cliques of sizes 3 and 4 and rarely 5. In our experiments, we could never observe a clique of size 6 or 7. The system moves from one metastable state to another and it can take a very long time for the system to come across the metastable state corresponding to the maximum clique. A sample of the intermittent evolution of one subpopulation fraction is shown in Figure 5.

Next, we planted a clique of size 10 in the same realization of G(100,0.25)G(100,0.25). We have chosen the size 10, because n\sqrt{n} is believed to be the critical scaling for detecting the planted clique [4, 40]. It is interesting to observe that in most simulation runs the system state is concentrated on the subpopulations corresponding to the planted clique (see as example Figures 6 and 7).

The planted clique experiments, in agreement with the theoretical results from [4, 40], indicate that there is a critical threshold such that, if cardinality of a tightly interacting subpopulation (clique) exceeds this threshold, the system will fairly quickly converge to the metastable state corresponding to the planted clique and will remain in that metastable state over long time intervals.

Refer to caption
Figure 5: A sample of trajectory of one subpopulation in the G(100,0.25)G(100,0.25) example. The x-axis corresponds to the time and the y-axis to the subpopulation fraction.
Refer to caption
Figure 6: Trajectories of subpopulations in the G(100,0.25)G(100,0.25) example with the planted clique of size 10. The x-axis corresponds to the time and the y-axis to the subpopulation fractions.
Refer to caption
Figure 7: The subpopulation fractions at discrete time step 7000 in G(100,0.25)G(100,0.25) with the planted clique of size 10. The x-axis corresponds to the subpopulation index and the y-axis to the subpopulation fraction. The first 10 indices correspond to the planted clique.

6 Conclusions and future research

We studied the phenomena of metastability and quasi-invariance as well as intermittent dynamics of a class of stochastic replicator equation in the small noise regime by invoking the Freidlin-Wentzell theory for small noise asymptotics of diffusions and some recent works on quasi-stationarity. We further specialized these results to a noisy replicator dynamics on graphs, with numerical experiments that confirmed our theoretical observations. We see a number of interesting future research directions worth note.

  1. 1.

    The exit time from the domain of attraction of a stable equilibrium is approximately exponentially distributed. Our example of a stochastic replicator dynamics on graphs shows that an exponential distribution still remains a very good approximation even for the empirical distribution of exit times when the identity of the metastable state is erased, so that each exit episode could be from a different metastable state. In some highly irregular potential functions, however, the exit distribution is observed to be closer to Pareto [23]. This could be because of clusters of closely spaced metastable states at many length scales, which leads to aggregate effects from a large number of individually exponential random variables. Clearly a more refined analysis is needed to explain this observation.

  2. 2.

    In the context of graph problems, estimating the depth of the deepest well of the potential and finding the precise scaling behavior of the mixing time with the size of the graph (see, e.g., [52]) are interesting problems, as is that of characterizing the mixing time in terms of other related stopping times such as hitting times [46]. Specializing further the graph dynamics to the planted clique problem [4, 40], it would be interesting to explore for what scales and after what time the clique becomes detectable by the replicator dynamics. These questions have been investigated for a different Metropolis dynamics in [40]. It will be interesting to carry out a parallel development for our model emphasizing the role of metastability.

Acknowledgements

This work was partly supported by EU Project Congas FP7-ICT-2011-8-317672, CEFIPRA Grants 5100-IT and IFC/DST-Inria-2016-01/448 and ARC Discovery Grant DP 160101236. We thank the referees and the editor for their careful reading and valuable suggestions which have immensely improved this article. We would like to acknowledge the fruitful discussions with Henrik Jensen and Jelena Grujić. We also thank Professors Siva Athreya, K. Suresh Kumar, Laurent Miclo, Kavita Ramanan, K. S. Mallikarjuna Rao and Alexander Veretennikov for valuable pointers to the literature.

This is the author version of the paper accepted to Dynamic Games and Applications journal.

References

  • [1] Absil, P.-A., Mahoney, R. and Sepulchre, R.  (2008) Optimization Algorithms on Matrix Manifolds, Princeton Uni. Press, Princeton, NJ.
  • [2] Aiba, K. (2015) “Waiting times in evolutionary dynamics with time-decreasing noise”, International J. Game Theory 44(2), 499-514.
  • [3] Akin, E. (1979) The Geometry of Population Genetics, Springer Verlag, Berlin.
  • [4] Alon, N., Krivelevich, M. and Sudakov, B. (1998) “Finding a large hidden clique in a random graph”, Random Structures and Algorithms 13(3-4), 457-466.
  • [5] Anderson, R. F. and Orey, S. (1976) “Small random perturbation of dynamical systems with reflecting boundary”, Nagoya Math. J. 60, 189-216.
  • [6] Avrachenkov, K., Borkar, V. S. and Nemirovski, D. (2010) “Quasi-stationary distributions as centrality measures for the giant strongly connected component of a reducible graph”, J. Computational and Appl. Math. 234(11), 3075-3090.
  • [7] Bass, R. F. and Perkins, E. A. (2002) “Degenerate stochastic differential equations with Hölder continuous coefficients and super-Markov chains”, Trans. of the American Math. Society 355(1), 373-405.
  • [8] Bianchi, A. and Gaudilliere, A. (2016) “Metastable states, quasi-stationary and soft measures”, Stochastic Processes and Their Applications 126(6), 1622-1680.
  • [9] Bollobas B. (2001) Random Graphs, 2nd ed., Cambridge University Press.
  • [10] Bomze, I. M. (1997) “Evolution towards the maximum clique”, J. Global Optimization 10(2), 143-164.
  • [11] Bomze, I. M., Budinich, M., Pardalos, P. M. and Pelillo, M. (1999) “The maximum clique problem”. In Handbook of Combinatorial Optimization (pp. 1-74). Springer.
  • [12] Bomze, I. M., Budinich, M., Pelillo, M. and Rossi, C. (2002) “Annealed replication: a new heuristic for the maximum clique problem”, Discrete Applied Mathematics, 121(1), 27-49.
  • [13] Bomze, I. M., Locatelli, M. and Tardella, F. (2008) “New and old bounds for standard quadratic optimization: dominance, equivalence and incomparability”, Mathematical Programming (Ser. A), 115(1), 31-64.
  • [14] Börgers, T. and Sarin, R. (1997) “Learning through reinforcement and replicator dynamics”, Journal of Economic Theory 77(1), 1-14.
  • [15] Borkar, V. S. (2008) Stochastic Approximation: A Dynamical Systems Viewpoint, Hindustan Publ. Agency, New Delhi, and Cambridge Uni. Press, Cambridge, UK.
  • [16] Borkar, V.S. and Sundaresan R. (2012) “Asymptotics of the invariant measure in mean field model with jumps”, Stochastic Systems, 2, 322-380.
  • [17] Borkar, V. S. and Kumar, K. S. (2010) “A new Markov selection procedure for degenerate diffusions”, Journal of Theoretical Probability 23(3), 729-747.
  • [18] Bovier, A. and Den Hollander, F. (2015) Metastability: A Potential Theoretic Approach, Springer Verlag, Berlin-Heidelberg.
  • [19] Cabrales, A. (2000). Stochastic replicator dynamics. International Economic Review, 41(2), 451-481.
  • [20] Catoni, O. (1991) “Sharp large deviations estimates for simulated annealing algorithms.” Annales de l’IHP Probabilités et statistiques. Vol. 27. No. 3.
  • [21] Catoni, O. (1992) “Rough large deviation estimates for simulated annealing: Application to exponential schedules.” Annals of Probability, 20(3), 1109-1146.
  • [22] Champagnat, N. and Villemonais, D. (2016) “Exponential convergence to quasi-stationary distribution and Q-process”, Prob. Theory and Related Fields 164(1), 243-283.
  • [23] Christensen, K., Di Collobiano, S. A., Hall, M. and Jensen, H. J. (2002) “Tangled nature: A model of evolutionary ecology”, Journal of Theoretical Biology 216(1), 73-84.
  • [24] Collet, P., Martinez, S. and San Martin J. (2014) Quasi-stationary Distributions: Markov Chains, Diffusions and Dynamical Systems, Springer Verlag, Berlin-Heidelberg.
  • [25] Corradi, V. and Sarin, R. (2000) “Continuous approximations of stochastic evolutionary game dynamics”, J. Economic Theory 94(2), 163-191. (Corrigendum in ibid., 140(1) (2008), e2-e4.)
  • [26] Eckmann, J. P. and Ruelle, D. (1985) “Ergodic theory of chaos and strange attractors”, Reviews of Modern Physics 57(3), 617-656.
  • [27] Foster, D. and Young H. P. (1990) “Stochastic evolutionary game dynamics”, Theoretical Population Biology 38(2), 219-232. (Correction note in ibid. (1997) 51(1), 77-78.)
  • [28] Freidlin, M. I. and Wentzell, A. D. (2012) Random Perturbations of Dynamical Systems (3rd ed.). Springer Verlag, New York.
  • [29] Fudenberg, D. and Harris, C. (1992) “Evolutionary dynamics with aggregate shocks”, J. Economic Theory 57(2), 420-441.
  • [30] Fudenberg, D. and Levine, D. K. (1998) Theory of Learning in Games, MIT Press, Cambridge, Mass.
  • [31] Fudenberg, D. and Tirole, J. (1991) Game Theory, MIT Press, Cambridge, Mass.
  • [32] Gao, X., Zhong, W. and Mei, S. (2013) “Stochastic evolutionary game dynamics and their selection mechanisms”, Computational Economics 41(2), 233-247.
  • [33] Hajek, B. (1988) “Cooling schedules for optimal annealing”, Mathematics of Operations Research 13(2), 311-329.
  • [34] Harsanyi, J. C. and Selten, R. (1988) A General Theory of Equilibrium Selection in Games, MIT Press, Cambridge, Mass.
  • [35] Hofbauer, J. and Imhof, L. A. (2009) “Time averages, recurrence and transience in stochastic replicator dynamics”, Annals Appl. Prob. 19(4), 1347-1368.
  • [36] Hofbauer, J. and Sigmund, K. (1998) Evolutionary Games and Population Dynamics, Cambridge Uni. Press, Cambridge, UK.
  • [37] Hsu, E. P. (2002) Stochastic Analysis on Manifolds, Graduate Studies in Mathematics Vol. 38, American Math. Society.
  • [38] Ikeda, N. and Watanabe, S. (1989) Stochastic Differential Equations and Diffusion Processes (2nd edition), North Holland–Kodansha, Amsterdam/Tokyo.
  • [39] Imhof, L. A. (2005) “The long-run behavior of the stochastic replicator dynamics”, Annals Appl. Prob., 15(1B), 1019-1045.
  • [40] Jerrum, M. (1992). “Large cliques elude the Metropolis process”, Random Structures and Algorithms 3(4), 347-359.
  • [41] Kang, W.  and Ramanan, K. (2014) “Characterization of stationary distributions of reflected diffusions”, Ann.  Applied Prob. 24(4), 1329-1374.
  • [42] Khasminskii, R. and Potsepun, N. (2006) “On the replicator dynamics behavior under Stratonovich type random perturbations”, Stochastics and Dynamics 6(2), 197-211.
  • [43] Kramers, H. A. (1940) “Brownian motion in a field of force and the difffusion model of chemical reactions”, Physica 7, 284-304.
  • [44] Krylov, N. V. (1973) “The selection of a Markov proccess from a Markov system of processes” (in Russian) Izv. Akad. Nauk SSSR, Ser. Mat. 37, 691-708.
  • [45] Manapat, M. L., Rand, D. G., Pawlowitsch, C. and Nowak, M. A. (2012) “Stochastic evolutionary dynamics resolve the traveler’s dilemma.”, Journal of Theoretical Biology 303, 119-127.
  • [46] Manzo, F. and Scoppola, E. (2016) “Strong times and first hitting”, arXiv preprint arXiv:1606.07244.
  • [47] Mäs, M. and Nax, H. H. (2016) “A behavioral study of “noise” in coordination games”, Journal of Economic Theory 162, 195-208.
  • [48] Mertikopoulos, P. and Moustakas, A. L. (2010) “The emergence of rational behavior in the presence of stochastic perturbations”, Ann. Appl. Prob. 20(4), 1359-1388.
  • [49] Mertikopoulos, P. and Viossat, Y. (2016) “Imitation dynamics with payoff shocks”, Intern. J. Game Theory 45, 291-320.
  • [50] Olivieri, E. and Vares, M. E. (2004) Large Deviations and Metastability, Encyclopedia of Mathematics and its Applications Vol. 100, Cambridge Uni. Press, Cambridge, UK.
  • [51] Ott, E. (2002) Chaos in Dynamical Systems, Cambridge University Press, Cambridge, UK.
  • [52] Panageas I. and Vishnoi N.K. (2016) “Mixing time of Markov chains, dynamical systems and evolution”, in LIPIcs-Leibniz International Proceedings in Informatics (Vol. 55). Schloss Dagstuhl-Leibniz-Zentrum fuer Informatik.
  • [53] Pinsky, R. (1985) “On the convergence of diffusion processes conditioned to remain in a bounded region for large time to limiting positive recurrent diffusion processes”, Annals of Probability 13(2), 363-378.
  • [54] Piovani, D., Grujić, J. and Jensen, H. J. (2016) “Linear stability theory as an early warning sign for transitions in high dimensional complex systems”, Journal of Physics A: Mathematical and Theoretical 49(29), 295102.
  • [55] Sandholm, W. H. (2011) Population Games and Evolutionary Dynamics, MIT Press, Cambridge, Mass.
  • [56] Sandholm, W. H., Documaci, E. and Lahkar, R. (2008) “The projection dynamic and the replicator dynamic”, Games and Economic Behavior 64(2), 666-683.
  • [57] Sandholm, W. H. and Staudigl, M. (2016) “Large deviations and stochastic stability in the small noise double limit”, Theoretical Economics 11.1, 279-355.
  • [58] Sandholm, W. H. and Staudigl, M., “Sample path large deviations for stochastic evolutionary game dynamics”, ArXiv preprint arXiv:1511.07897.
  • [59] Shahshahani, S. (1979) A new mathematical framework for the study of linkage and selection, Mem. American Math. Society 211.
  • [60] Tanaka, H. (1979) “Stochastic differential equations with reflecting boundary condition in convex regions”, Hiroshima Math. J. 9, 163-177.
  • [61] Taylor, P. D. and Jonker, L. B. (1978) “Evolutionary stable strategies and game dynamics”, Mathematical Biosciences 40(1+2), 145-156.
  • [62] Tokita, K. and Yasutomi, A. (2003) “Emergence of a complex and stable network in a model ecosystem with extinction and mutation”, Theoretical Population Biology 63(2), 131-146.
  • [63] Turelli, M. (1977) “Random environments and stochastic calculus”, Theoretical Population Biology 12(2), 140-178.
  • [64] van Damme, E. (2002) Stability and Perfection of Nash Equilibria (2nd ed.), Springer Verlag, Berlin-Heidelberg.
  • [65] Weibull, J. (1995) Evolutionary Game Theory, MIT Press, Cambridge, Mass.
  • [66] Watanabe, S. and Yamada, T. (1971) “On the uniqueness of solutions of stochastic differential equations II”, Journal of Mathematics of Kyoto University 11(3), 553-563.
  • [67] Young, H. P. (1998) Individual Strategy and Social Structure, Princeton Uni. Press, Princeton, NJ.