Phase transitions in non-linear urns with interacting types

In general terms, an urn model is a system containing a number of particles of different types (often regarded as balls of different colours, for ease of visualisation). At each time step, a set of particles is sampled from the system, whose contents are then altered depending on the sample which was drawn. Pemantle [10] surveys several ways to approach this model framework.

This paper is limited to models with a single urn from which a single ball is drawn, its colour is noted and it is then returned to the urn along with one new ball of that same colour. In addition, we introduce a graph-based interaction according to which the probability of choosing a ball of a given colour is reinforced not only by its own proportion in the urn, but also by the proportions of balls of other colours. Therefore, the interaction arises among balls of different colours, as opposed to the so-called interacting urn models consisting of systems of multiple urns (e.g. Benaïm et al [2], Launay and Limic [7]) in which different urns (each containing balls of different colours) interact with one another. Our model is also different from the graph-based competition described in van der Hofstad et al [11], where the colours correspond to edges of the graph, which compete with, as opposed to being reinforced by, other edges incident on the same vertices.

We now formally define our model. Consider an urn containing balls of $d$ colours. The vector $x(n)=(x_{1}(n),\ldots,x_{d}(n))\in\mathbb{N}^{d}$ denotes the number of balls of each colour at time $n=0,1,2,...$. The strength of the reinforcement is given by a positive real number $\beta>0$ and we denote by $x^{\beta}(n)$ the coordinate-wise $\beta$ power of the column vector $x(n)$. The interaction is defined as follows. Given a non-negative matrix $A=(a_{ij})_{i,j=1}^{d}$, define the column vector

$$u(n):=Ax^{\beta}(n)$$

(1)

Let $\mathcal{F}_{n}$ be the $\sigma$-algebra generated by the $x(m)$ for $0\leq m\leq n$ and write $u_{i}(n)$ for the $i$-th component of $u(n)$. The transition probabilities are then

$$\mathbb{P}(x(n+1)-x(n)=\mathbf{e}_{i}\>|\>\mathcal{F}_{n})=\frac{u_{i}(n)}{\sum_{j=1}^{d}u_{j}(n)},\quad i=1,\ldots,d,$$

(2)

where $\mathbf{e}_{i}$ is the unit vector in direction $i$. That is, one ball is added to the urn at each time step, and the right hand side of (2) gives the probability that it is of colour $i$.

Now, let $n_{0}>0$ be the initial number of balls so that at time $n$ the urn contains $n+n_{0}$ balls. Then, the proportion of balls of each colour is a process in the $(d-1)$-dimensional simplex $\Delta^{d-1}$ given by the vector

$$\bar{x}(n)=x(n)/(n+n_{0}).$$

(3)

When $A$ is a multiple of the identity matrix (therefore, $\bar{x}(n)$ having no interaction) it is well-known (see Oliveira [8]) that the process $\bar{x}$ undergoes a phase transition as follows. For $\beta<1$, the process converges almost surely to the ‘centre’ of the simplex, that is, the asymptotic proportion of balls of each colour are all the same. For $\beta=1$, commonly referred to as the Pólya urn model, the process converges almost surely to a non-trivial random variable supported in the interior of the simplex. For $\beta>1$, the process converges almost surely to one of the corners of the simplex. In this case that a single type dominates was proved by Khanin and Khanin [5] following on from the two-type case which can be covered using Rubin’s Theorem in Davis [4].

For a two-colour urn model $d=2$ and symmetric interaction,

$$A=\left(\begin{matrix}1&a\\ a&1\end{matrix}\right),\quad a>0,$$

it was proved by the first author in Theorem 2.2.1, [3], that there was a phase transition as follows.

	$\displaystyle(i)\quad\text{ if }\quad\left(\tfrac{1-a}{1+a}\right)\beta\leq 1,\quad\text{then}\quad\bar{x}(n)\rightarrow(\tfrac{1}{2},\tfrac{1}{2})\quad a.s.$
	$\displaystyle(ii)\quad\text{if }\quad\left(\tfrac{1-a}{1+a}\right)\beta>1,\quad\text{then}\quad\bar{x}(n)\rightarrow\Psi\quad a.s.,$

where $\Psi$ is a random vector supported on $\left\{\left(\frac{1}{1+r},\frac{r}{1+r}\right),\left(\frac{r}{1+r},\frac{1}{1+r}\right)\right\}$ and $r:=r(a,\beta)$ is the unique root in $(0,1)$ of $\mathcal{P}_{a,\beta}(z)=az^{\beta+1}-z^{\beta}+z-a=0.$ In case (ii), $\mathbb{P}[\bar{x}(n)\rightarrow(\frac{1}{2},\frac{1}{2})]=0.$ Note that for $\beta=1$, the process $(u(n))_{n\geq 0}$ is a Friedman’s urn model and statement (i) yields $u(n)/(u_{1}(n)+u_{2}(n))\rightarrow(\frac{1}{2},\frac{1}{2})$ $a.s.$ as expected.

Some similar results appear in Laruelle and Pagès [6]. The definition of the model there allows for a more general skewing function than $x^{\beta}$ as above, but in the $d=2$ case gives phase transitions similar to those in [3], and relates them to the eigenvalues of the generating matrix $H$, which is related to our matrix $A$. Furthermore they show that, for any $d$, for a concave skewing function (corresponding to $\beta<1$ in our model) and a bi-stochastic generating matrix there is almost sure convergence to the centre of the simplex, and they give conditions, related to the eigenvalues of the generating matrix, under which for a convex skewing function (our $\beta>1$) there is probability zero of converging to the centre of the simplex.

In this paper, we follow up on the results for $d=2$ in [3] and those in [6], with the aim to generalise from $d=2$ to larger values of $d$ and to see whether more types of behaviour emerge when this is done. We show that this is indeed the case when $d=3$, where we consider two particular choices of $A$. Our results can be seen as extensions of those of [6] in certain specific cases.

First of all, we consider a choice of $A$ with a symmetric interaction of the same strength $a$ for each pair of colours. The following theorem shows that in this system there are three phases, as opposed to two when $d=2$; there are phases where there is almost sure convergence to a symmetric limit and where there is almost sure convergence to one of a number of asymmetric limits, which are analogues of the phases when $d=2$, but there is also an intermediate phase where both symmetric and asymmetric limits are possible.

Theorem 1.1 presents the results in terms of phase transitions in $\beta$ with $a$ fixed. However, because both $\beta_{1}(a)$ and $\frac{1+2a}{1-a}$ are increasing functions of $a$ which converge to $1$ as $a\to 0$ and to $\infty$ as $a\to 1$, it is also possible to see them as phase transitions in $a$ with $\beta>1$ fixed: if $a<\frac{\beta-1}{2+\beta}$ then we will be in case (c), if $\frac{\beta-1}{2+\beta}<a<\beta_{1}^{-1}(\beta)$ then we will be in case (b), and if $a>\beta_{1}^{-1}(\beta)$ we will be in case (a).

Theorem 1.1 can be seen as an extension of the results of Proposition 2.15 of Laruelle and Pagès [6] in the case where the matrix $H=\frac{1}{1+2a}\begin{pmatrix}1&a&a\\ a&1&a\\ a&a&1\end{pmatrix}$ in that their result shows the non-convergence to $(\frac{1}{3},\frac{1}{3},\frac{1}{3})$ in the case (i)(c), as the second largest eigenvalue of $H$ is $\frac{1-a}{1+2a}$. We can also state the condition $\beta>\frac{1+2a}{1-a}$ in terms of the eigenvalues of $A$: the right hand side can be seen as the ratio of the two largest eigenvalues. It might be reasonable to conjecture that an extension to $d>3$ of Theorem 1.1 might involve a similar condition on the eigenvalues; however we note that the other phase transition does not appear to be related to the eigenvalues of $A$. We discuss the question of what happens with $d>3$ further in Section 5.4.

We also consider a system where each colour is reinforced by itself and by one other, in a cyclic way. For this system, the following theorem shows the existence of a phase transition between a phase with convergence with positive probability to a symmetric limit and a phase where there is no convergence to a limit and there is cycling behaviour.

In Section 2 we discuss the stochastic approximation methods we use in the proofs, while the proofs themselves are in Section 3 for Theorem 1.1 and Section 4 for Theorem 1.2. In the final Section 5, we illustrate the results with some examples and simulations, including some examples beyond those covered by Theorems 1.1 and 1.2.

In this section we introduce some of the stochastic approximation ideas which appear in our proofs.

The type of stochastic approximation process we will be interested in is a Robbins-Monro algorithm in $\mathbb{R}^{d}$, following section 4.2 of Benaïm [1]. This is a stochastic process $(y(n))_{n\in\mathbb{N}}$, with natural filtration $(\mathcal{F}_{n})$, taking values in $\mathbb{R}^{d}$ which satisfies

$$y(n+1)-y(n)=\gamma_{n}(F(y(n))+\xi_{n+1}),$$

(5)

where $F:\mathbb{R}^{d}\to\mathbb{R}^{d}$ is a deterministic vector field, $\gamma_{n}$ is a step size satisfying certain conditions, and $\mathbb{E}(\xi_{n+1}|\mathcal{F}_{n})=0$ with $\xi_{n}$ being $\mathcal{F}_{n}$-measurable. For most results it is required that $\gamma_{n}\to 0$ as $n\to\infty$ but that $\sum_{n=1}^{\infty}\gamma_{n}=\infty$, and some conditions are also needed on the noise term $\xi_{n+1}$, typically that it is bounded in $L^{q}$ for some $q$ which depends on the $\gamma_{n}$: see for example Proposition 4.2 of Benaïm [1].

The above sequence can be thought as a numerical approximating method with varying step size $\gamma_{n}$ for solving the ODE $dx/dt=F(x)$. Under the conditions discussed above, the asymptotic behavior of $(\bar{x}_{n})_{n\in\mathbb{N}}$ and the underlying ODE are closely connected: as described in [1], define an interpolated version $(X(t))_{t\geq 0}$ by defining $\tau_{0}=0$ and for $n\in\mathbb{N}$ $\tau_{n}=\sum_{i=1}^{n}\gamma_{i}$, then defining $X(\tau_{n}+s)=\bar{x}_{n}+\frac{\bar{x}_{n+1}-\bar{x_{n}}}{\gamma_{n+1}}$ for $0\leq s<\gamma_{n+1}$. Then the interpolated process $(X(t))_{t\geq 0}$ is an asymptotic pseudotrajectory for the ODE. This is called the ODE method or the dynamical system approach, which alongside some probabilistic techniques, is applied to examine almost sure dynamics of stochastic approximation processes.

Let $\Phi:\mathbb{R}_{+}\times\mathbb{R}^{d}\to\mathbb{R}^{d}$ be the semiflow induced by the vector field $F$, so that $(\Phi(t,x))_{t\geq 0})$ is the trajectory of $F$ started at $x$. A useful situation for analysis of stochastic approximation processes is where there is a Lyapunov function for the vector field $F$: a function $V:\mathbb{R}^{d}\to\mathbb{R}^{d}$ such that on a trajectory $(\Phi(t,x))_{t\geq 0}$ of the vector field, $V(\Phi(t,x))$ is strictly decreasing in $t$ except where $\Phi(t,x)$ is a stationary point of $F$. If this holds, then under mild conditions then the Robbins-Monro algorithm will converge almost surely to a (possibly random) limit, which will be a stationary point of $F$. See section 6.2 of Benaïm [1], and in particular Corollary 6.6 therein.

We now show that our process $(\bar{x}_{n})$ can be put into Robbins-Monro form, and that the assumptions of Proposition 4.2 of Benaïm [1] are satisfied, allowing the theory of asymptotic pseudotrajectories to be used. For a general matrix $A$ and a given configuration of balls $x(n)$ at time $n$, let $i_{n+1}\in\{1,\ldots,d\}$ be the random colour of the ball to be added in the urn at time $n+1$. Then, note that

	$\displaystyle\bar{x}(n+1)$	$\displaystyle=\frac{x(n)+{\bf e}_{i_{n+1}}}{n_{0}+n+1}=\frac{(n_{0}+n)\bar{x}(n)+{\bf e}_{i_{n+1}}}{n_{0}+n+1}$
		$\displaystyle=\left(1-\frac{1}{n_{0}+n+1}\right)\bar{x}(n)+\frac{{\bf e}_{i_{n+1}}}{n_{0}+n+1},$		(6)

implying

$$\bar{x}(n+1)-\bar{x}(n)=\frac{1}{n_{0}+n+1}({\bf e}_{i_{n+1}}-\bar{x}(n)).$$

(7)

To put (7) into Robbins-Monro form, we rearrange the right-hand side into a deterministic part and a zero mean “noise”. More specifically, let

$$F(\bar{x}(n)):=\mathbb{E}[{\bf e}_{i_{n+1}}\,|\,\mathcal{F}_{n}]-\bar{x}(n),$$

(8)

and

$$\xi_{n+1}:={\bf e}_{i_{n+1}}-\mathbb{E}[{\bf e}_{i_{n+1}}\,|\,\mathcal{F}_{n}].$$

(9)

By setting $\gamma_{n}=1/(n_{0}+n+1)$, we obtain

$$\bar{x}(n+1)-\bar{x}(n)=\gamma_{n}(F(\bar{x}(n))+\xi_{n+1}).$$

(10)

We note that the components of $\xi_{n+1}$ are uniformly bounded in modulus by $2$, meaning that the conditions of Proposition 4.2 of Benaím [1] are satisfied, meaning that the interpolated process $(X(t))_{t\geq 0}$ is indeed an asymptotic pseudotrajectory of $F$.

We can write the components of $F$ as

$$F_{i}(x)=\frac{u_{i}}{\sum_{j=1}^{d}u_{j}}-x_{i},\quad i=1,\ldots,d,$$

where $u_{i}$ has the same relationship to $x$ as $u_{i}(n)$ to $x(n)$.

Typically convergence happens with positive probability to a stable stationary point but with probability zero to a linearly unstable one. For the former, Theorem 7.3 of Benaïm [1], which holds under the assumptions on the noise of Proposition 4.2, shows that a stable stationary point has positive probability of being a limit as long as its basin of attraction $W$ contains points which are attainable in the sense that the process has positive probability of being indefinitely close to them at indefinitely large times, a condition which is usually satisfied for all points in the simplex for urn type processes under mild conditions; see Example 7.2 of [1].

To show that linearly unstable stationary points have zero probability of being limits it is usually necessary to check that there is expectation bounded away from zero of the positive part of the component of the noise in any given direction; see Pemantle [9], Theorem 1. This condition ensures that the process has zero probability of being trapped on a stable manifold and converging to the unstable point. Theorem 1 of [9] also requires the noise term $\xi_{n+1}$ (in our notation) to be uniformly bounded; as noted above this is satisfied in our process.

In this section we prove Theorem 1.1.

Throughout this section we let $A$ be the matrix

$$\begin{pmatrix}1&a&a\\ a&1&a\\ a&a&1\end{pmatrix},\quad a>0.$$

(11)

In this case the vector field $F$ is given by

	$\displaystyle F_{1}(x_{1},x_{2},x_{3})$	$\displaystyle=$	$\displaystyle\frac{x_{1}^{\beta}+ax_{2}^{\beta}+ax_{3}^{\beta}}{(1+2a)(x_{1}^{\beta}+x_{2}^{\beta}+x_{3}^{\beta})}-x_{1}$		(12)
	$\displaystyle F_{2}(x_{1},x_{2},x_{3})$	$\displaystyle=$	$\displaystyle\frac{ax_{1}^{\beta}+x_{2}^{\beta}+ax_{3}^{\beta}}{(1+2a)(x_{1}^{\beta}+x_{2}^{\beta}+x_{3}^{\beta})}-x_{2}$		(13)
	$\displaystyle F_{3}(x_{1},x_{2},x_{3})$	$\displaystyle=$	$\displaystyle\frac{ax_{1}^{\beta}+ax_{2}^{\beta}+x_{3}^{\beta}}{(1+2a)(x_{1}^{\beta}+x_{2}^{\beta}+x_{3}^{\beta})}-x_{3}$		(14)

The stochastic approximation approach indicates that the possible limits of our process will be stationary points of $F$, so we start by identifying these. Noting that the lines $x_{1}=x_{2}$, $x_{1}=x_{3}$ and $x_{2}=x_{3}$ are each invariant under $F$, define the function

$$\mathcal{P}_{a,\beta}(z)=az^{\beta+1}-z^{\beta}+(1+a)z-2a,$$

(15)

which we will see is related to the dynamics restricted to one of these lines. The following result shows that all stationary points of $F$ are located on at least one of these lines and expresses them in terms of solutions to $\mathcal{P}_{a,\beta}(z)=0$.

We now show that the process will indeed, almost surely, converge to one of the stationary points identified in Proposition 3.1.

Note that the Lyapunov function $L$ generalises in an obvious way to more than three types, as long as all off-diagonal entries are equal.

We now investigate the stability of these roots, for which recall the terminology in Definition 1.