2.1 Entropic regularisation for spatially inhomogeneous evolutionary games

In the model (1.3), mixed strategies $\sigma_i$ are prone to converging exponentially to very singular distributions so that agents cannot react quickly to changes in the environment and exhibit strong inertia, contradicting the notion of rational agents. This behaviour can be illustrated with a simple example.

As a remedy, we propose an adaptation of the model (1.3) by adding an entropic term to the dynamics of the mixed strategies $\sigma_i$ . The update rule (1.3b) modifies $\sigma_i$ towards maximising the average benefit

\begin{equation*} \frac{1}{N} \sum_{j=1}^N \int_{U \times U} J(x_i,u,x_j,u^{\prime})\,\textrm{d} \sigma_i(u)\,\textrm{d} \sigma_j\!\left(u^{\prime}\right).\end{equation*}

To this objective, we will now add the entropy $\varepsilon\cdot \int_U \log\left(\frac{\textrm{d} {\sigma_i}}{\textrm{d} {\eta}}\right)\!\textrm{d} \sigma_i$ where $\eta \in \mathcal{M}_1(U)$ is a suitable ‘uniform’ reference measure and $\varepsilon>0$ is a weighting parameter. This will pull $\sigma_i$ towards $\eta$ and thus prevent exponential degeneration.

Entropic regularisation. The entropic regularisation implies that mixed strategies are no longer general probability measures over U, but they become densities with respect to a reference measure $\eta \in \mathcal{M}_1(U)$ , which we can assume without loss of generality to have full support, that is, $\textrm{spt}(\eta) = U$ . Thus, for a fixed $\varepsilon > 0$ , set

(2.1) \begin{align} \mathcal{S}_{a,b} \;:=\; \left\{ \sigma\colon U \to \mathbb{R}_+, \;\sigma \text{ measurable,} \;\sigma(u) \in [a,b]\text{ } \eta\text{-a.e.},\; \int_U \sigma\,\textrm{d} \eta = 1 \right\},\end{align}

where $0 < a < 1 < b < \infty$ (the bounds a and b are required for technical reasons, see below). For $\lambda > 0$ , the modified entropic dynamics is then given by

(2.2a) \begin{align} \partial_t x_i(t) & = \int_U e(x_i(t),u)\, \sigma_i(t)(u)\,\textrm{d} \eta(u),\qquad\qquad\qquad\qquad\qquad \end{align}

(2.2b) \begin{align} \partial_t \sigma_i(t) & = \lambda \cdot \left[ \frac{1}{N} \sum_{j=1}^N f^J\!\big(x_i(t),\sigma_i(t),x_j(t),\sigma_j(t)\big) +f^\varepsilon\big(\sigma_i(t)\big) \right],\end{align}

where the function $f^J$ now formally needs to be given in terms of densities, instead of measures:

(2.3) \begin{align} & f^J\!\big(x,\sigma,x^{\prime},\sigma^{\prime}\big) \;:=\;\nonumber\\[4pt]& \quad\left[\int_U J(x,\cdot,x^{\prime},u^{\prime})\,\sigma^{\prime}\left(u^{\prime}\right)\!\textrm{d} \eta\left(u^{\prime}\right) - \int_U \int_U J(x,v,x^{\prime},u^{\prime})\,\sigma^{\prime}\left(u^{\prime}\right)\,\sigma(v)\,\textrm{d} \eta\left(u^{\prime}\right)\!\textrm{d} \eta(v) \right] \cdot \sigma.\end{align}

The additional term $f^\varepsilon$ , corresponding to entropy regularisation, is given by:

(2.4) \begin{align} f^\varepsilon(\sigma) & \;:=\; \varepsilon \cdot \left[ -\log(\sigma(\cdot)) + \int_U \log(\sigma(v))\,\sigma(v)\,\textrm{d} \eta(v) \right] \cdot \sigma.\end{align}

In (2.2b), we have explicitly added the factor $\lambda \in (0,\infty)$ to control the relative time scale of the dynamics for mixed strategies and locations.

Remark 2.3 (Entropic gradient flow). Formally, the contribution of $f^\varepsilon$ to the dynamics (2.2b) corresponds to a gradient flow in the Hellinger–Kakutani metric over $\mathcal{S}_{a,b}$ of the (negative) entropy

(2.5) \begin{align} H(\sigma) \;:=\; \int_U \log(\sigma)\,\sigma\,\textrm{d} \eta\end{align}

of the density $\sigma$ with respect to the reference measure $\eta$ .

Well-posedness and mean-field limit of the regularised system. Following [Reference Ambrosio, Fornasier, Morandotti and Savaré6] (see (1.5)), we analyse (2.2) as an interaction system in the Banach space $Y = \mathbb{R}^d \times L^p_\eta(U)$ , $1 \leq p < \infty$ . For $y = (x,\sigma) \in Y$ we set $\|y\|_Y = \|x\| + \|\sigma\|_{L^p_\eta(U)}$ . Define

(2.6) \begin{align} \mathcal{C}_{a,b} \;:=\; \mathbb{R}^d \times \mathcal{S}_{a,b},\end{align}

and for $y, y^{\prime} \in \mathcal{C}_{a,b}$ let us set

(2.7) \begin{align}f \colon \mathcal{C}_{a,b} \times \mathcal{C}_{a,b} \to Y, \quad f\!\left(y,y^{\prime}\right) = \left[f^{e}(x,\sigma), \lambda f^{J,\varepsilon}\left(x,\sigma,x^{\prime},\sigma^{\prime}\right)\right],\end{align}

where

(2.8) \begin{align}f^{e}(x,\sigma) = \int_U e(x,u)\, \sigma(u)\,\textrm{d} \eta(u) \quad \text{and} \quad f^{J,\varepsilon}\!\!\left(x,\sigma,x^{\prime},\sigma^{\prime}\right) = f^J\!\big(x,\sigma,x^{\prime},\sigma^{\prime}\big) +f^\varepsilon(\sigma),\end{align}

so that (2.2) takes the equivalent form

(2.9) \begin{equation}\partial_t y_i(t) = \frac{1}{N} \sum_{j=1}^N f(y_i(t), y_j(t)) \quad \text{ for }i = 1,\dots, N.\end{equation}

Similar to [Reference Ambrosio, Fornasier, Morandotti and Savaré6], well-posedness of (2.9) relies on the Lipschitz continuity of f and on the following compatibility condition (for the corresponding proofs see Appendix A, Lemmas A.1, A.2 and A.4).

Lemma 2.5 (Compatibility condition). For $J \in \mathcal{X}$ and $\varepsilon > 0$ , let $f^J$ and $f^\varepsilon$ be defined as in (2.3) and (2.4). Then, there exist a,b with $0 < a < 1 < b < \infty$ such that for any $(x,\sigma),(x^{\prime},\sigma^{\prime}) \in \mathbb{R}^d \times \mathcal{S}_{a,b}$ there exists some $\theta>0$ such that

(2.10) \begin{align} \sigma + \theta \lambda \!\left[f^J\!\big(x,\sigma,x^{\prime},\sigma^{\prime}\big) +f^\varepsilon(\sigma) \right] \in \mathcal{S}_{a,b}. \end{align}

Intuitively, for specific choices of the bounds a and b, (2.10) states that moving from $\sigma \in \mathcal{S}_{a,b}$ into the direction generated by $f^J$ and $f^\varepsilon$ , we will always remain within $\mathcal{S}_{a,b}$ for some finite time. Eventually, similar to (1.6), a mean-field limit description of (2.9) is formally given by

(2.11) \begin{align} \partial_t \Sigma(t) + \textrm{Div}\big(b(\Sigma(t)) \cdot \Sigma(t)\big) = 0 \quad \text{with} \quad b(\Sigma(t))(y) \;:=\; \int_{\mathcal{C}_{a,b}}\!\! f\!\left(y,y^{\prime}\right)\!\textrm{d}\Sigma(t)(y^{\prime}).\end{align}

Like (1.6), this is a PDE whose domain is a Banach space and we refer to [Reference Ambrosio, Fornasier, Morandotti and Savaré6] for the technical details. We can then summarise the main result in the following theorem.

Proof. The theorem follows by invoking Theorem 4.1 from [Reference Ambrosio, Fornasier, Morandotti and Savaré6]. On a more technical level, Theorem 5.3 of [Reference Ambrosio, Fornasier, Morandotti and Savaré6] provides the uniqueness of the solution to (2.11) in a Eulerian sense. We now show that the respective requirements are met.

First, the set $\mathcal{C}_{a,b}$ , (2.6), is a closed convex subset of Y with respect to $\|\cdot\|_Y$ for any $0 < a < b < +\infty$ . Likewise, for any $0 < a < b < +\infty$ the map f driving the evolution (2.9) is Lipschitz continuous: indeed, one can prove that $f^J$ and $f^\varepsilon$ are Lipschitz continuous (see Lemmas A.1 and A.2 in Appendix A) and Lipschitz continuity of $f^e$ follows from the fact that e is Lipschitz continuous and bounded. Furthermore, as a consequence of Lemma 2.5, one can choose a and b so that the extended compatibility condition

\begin{align*} \forall\,R > 0 \, \exists \, \theta > 0\,:\, \quad y,y^{\prime} \in \mathcal{C}_{a,b} \cap B_R(0) \quad \Rightarrow \quad y + \theta\,f\!\left(y,y^{\prime}\right) \in \mathcal{C}_{a,b} \end{align*}

holds, cf. [Reference Ambrosio, Fornasier, Morandotti and Savaré6, Theorem B.1]. Therefore, we may invoke [Reference Ambrosio, Fornasier, Morandotti and Savaré6, Theorem 4.1]. Combining this with [Reference Ambrosio, Fornasier, Morandotti and Savaré6, Theorem 5.3], which is applicable since $L^p_\eta(U)$ is separable and thus so is Y, the result follows.

2.2 Undisclosed setting and fast reactions

In the model (2.2), the decision process of each agent potentially involves knowledge of the mixed strategies of the other agents. Often it is reasonable to assume that there is no such knowledge, which can be reflected in the model by assuming that the payoff J(x, u, x ^′, u ^′) does not actually depend on u ^′, the other agent’s strategy. Note that J(x, u, x ^′, u ^′) may still depend on x ^′, the other agent’s location. We call this the undisclosed setting.

Additionally, often it is plausible to assume that the (regularised) dynamic that governs the mixed strategies $(\sigma_i(t))_{i=1}^{N}$ of the agents runs at a much faster time scale compared to the physical motion of the spatial locations $(x_i(t))_{i=1}^N$ . This corresponds to a large value of $\lambda$ in (2.2). Therefore, for the undisclosed setting we study the fast-reaction limit $\lambda \to \infty$ .

The main results of this section are Theorems 2.13 and 2.14 which establish that the undisclosed fast-reaction limit is in itself a well-defined model (with a consistent mean-field limit as $N \to \infty$ ) and that the undisclosed model converges to this limit as $\lambda \to \infty$ .

Undisclosed setting. In the undisclosed setting, the general formulas (2.2) for the dynamics can be simplified as follows:

(2.13a) \begin{align} \partial_t x_i(t) & = \int_U e(x_i(t),u)\, \sigma_i(t)(u)\,\textrm{d} \eta(u),\qquad\qquad\ \ \qquad \end{align}

(2.13b) \begin{align} \partial_t \sigma_i(t) & = \lambda \cdot \left[ \frac{1}{N} \sum_{j=1}^N f^J\!\big(x_i(t),\sigma_i(t),x_j(t)\big) +f^\varepsilon\big(\sigma_i(t)\big) \right] \end{align}

where $f^\varepsilon$ is as given in (2.4) and $f^J$ simplifies to

(2.14) \begin{align} f^J\!\big(x,\sigma,x^{\prime}\big) \;:=\; \left[ J(x,\cdot,x^{\prime}) - \int_U J(x,v,x^{\prime})\,\sigma(v)\,\textrm{d} \eta(v) \right] \cdot \sigma.\end{align}

For finite $\lambda < +\infty,$ the dynamics of this model are still covered by Theorem 2.6.

Fast reactions of the agents. Intuitively, as $\lambda \to \infty$ in (2.13), at any given time t the mixed strategies $(\sigma_i(t))_{i=1}^N$ will be in the unique steady state of the dynamics (2.13b) for fixed spatial locations $(x_i(t))_{i=1}^N$ . For given locations $\mathbf{x} = (x_1,\ldots,x_N) \in [\mathbb{R}^d]^N$ this steady state is given by

(2.15a) \begin{align} \sigma_{i}^J(\mathbf{x}) &\equiv \sigma_{i}^J(x_1,\ldots,x_N) \;:=\; \frac{\exp\left(\dfrac{1}{\varepsilon N}\sum_{j=1}^N J\!\left(x_i,\cdot,x_j\right)\right)}{ \int_U \exp\left(\dfrac{1}{\varepsilon N}\sum_{j=1}^N J\!\left(x_i,v,x_j\right)\right)\!\textrm{d} \eta(v) }. \end{align}

(This computation is explicitly shown in the proof of Theorem 2.14.) The spatial agent velocities associated with this steady state are given by

(2.15b) \begin{align} v_{i}^J(\mathbf{x}) &\equiv v_i^J\!\left(x_1,\dots,x_N\right) \;:=\; \int_U e(x_i,u)\, \sigma_{i}^J(x_1,\ldots,x_N)(u)\,\textrm{d} \eta(u)\end{align}

and system (2.13) turns into a purely spatial ODE in the form

(2.15c) \begin{align} \partial_t x_i(t) &= v_i^J(x_1(t),\dots,x_N(t)) \quad \text{for } i = 1,\dots,N.\end{align}

Unlike in Newtonian models over $\mathbb{R}^d$ , (1.1), here the driving velocity field $v_i^J(x_1(t),\ldots,x_N(t))$ is not a linear superposition of the contributions by each $x_j(t)$ . This nonlinearity is a consequence of the fast-reaction limit and allows for additional descriptive power of the model that cannot be captured by the Newton-type model (1.1). This is illustrated in the two subsequent Examples 2.7 and 2.8.

Example 2.7 (Describing Newtonian models as undisclosed fast-reaction models). Newtonian models can be approximated by undisclosed fast-reaction entropic game models. We give a sketch for this approximation procedure. For a model as in (1.1), choose $U \subset \mathbb{R}^d$ such that it contains the range of f (for simplicity we assume that it is compact), let $e(x,u) \;:=\; u$ and then set

\begin{equation*}J\!\left(x,u,x^{\prime}\right) \;:=\; -\|u-f\!\left(x,x^{\prime}\right)\!\|^2.\end{equation*}

Accordingly, the stationary mixed strategy of agent i is given by

\begin{align*}\sigma_{i}^J(\mathbf{x})(u) & = \mathcal{N}\left({\exp\left(-\frac{1}{\varepsilon N}\sum_{j=1}^N \|u-f\!\left(x_i,x_j\right)\!\|^2\right)}\right),\end{align*}

where $\mathcal{N}\left({\cdot}\right)$ denotes the normalisation operator for a density with respect to $\eta$ . Observe now that

\begin{align*}\frac{1}{N} \sum_{j=1}^N \|u-f\!\left(x_i,x_j\right)\!\|^2 &= \|u\|^2 - \frac2N \sum_{j=1}^N u \cdot f\!\left(x_i,x_j\right) + \frac1N \sum_{j=1}^N \|{\kern1pt}f\!\left(x_i,x_j\right)\!\|^2\end{align*}

while

\begin{align*}\left\|u-\frac{1}{N}\sum_{j=1}^N f\!\left(x_i,x_j\right)\right\|^2 &= \|u\|^2 - \frac2N \sum_{j=1}^N u \cdot f\!\left(x_i,x_j\right) + \frac{1}{N^2}\!\left\| \sum_{j=1}^N f\!\left(x_i,x_j\right)\right\|^2.\end{align*}

Hence, since the terms not depending on u are cancelled by the normalisation, we also have

\begin{align*}\sigma_{i}^J(\mathbf{x})(u) & = \mathcal{N}\left({\exp\left(-\frac{1}{\varepsilon}\left\|u-\frac{1}{N}\sum_{j=1}^N f(x_i,x_j)\right\|^2\right)}\right).\end{align*}

So the maximum of the density $\sigma_{i}^J(\mathbf{x})$ is at $\tfrac{1}{N}\sum_{j=1}^N f(x_i,x_j)$ , which is the velocity imposed by the Newtonian model. For small $\varepsilon$ , $\sigma_{i}^J(\mathbf{x})$ will be increasingly concentrated around this point, and hence, the game model will impose a similar dynamic. Numerically, this is demonstrated in Figure 2 for a 1-dimensional Newtonian model as in (1.1) driven by

(2.16) \begin{align}f(x,x^{\prime}) = - x - \frac{\tanh\left(5\left(x^{\prime}-x\right)\right)}{\left(1+\|x^{\prime}-x\|\right)^2}.\end{align}

Figure 2. Approximation of a Newtonian model by an undisclosed fast-reaction entropic game model. Solid lines: original model, dashed lines: approximation. The Newtonian model is driven by (2.16), the approximation procedure is described in Example 2.7. The approximation becomes more accurate as $\varepsilon$ decreases.

Example 2.8 (Undisclosed fast-reaction models are strictly richer than Newtonian models). In the previous example, each agent j tried to persuade agent i to move with velocity $f(x_i,x_j)$ . Deviations from this velocity were penalised in the payoff function with a quadratic function. By minimising the sum of these functions, the compromise of the linear average of all $f(x_i,x_j)$ then yields the best payoff. By picking different penalties we may obtain other interesting behaviour that cannot be described by Newtonian systems.

As above, let U be a (sufficiently large) compact subset of $\mathbb{R}^d$ , $e(x,u) \;:=\; u$ . Let $g \,:\, \mathbb{R}^d \to \mathbb{R}$ be a Lipschitz non-negative ‘bump function’ with compact support, that is, g is maximal at 0 with $g(0)=1$ , and $g(u)=0$ for $\|u\|\geq \delta$ for some small $\delta >0$ . Then we set

\begin{equation*}J\!\left(x,u,x^{\prime}\right) \;:=\; \left(C-\|x-x^{\prime}\|\right) \cdot g\left(u-\frac{x^{\prime}-x}{\|x^{\prime}-x\|} \right)-C \cdot g(u),\end{equation*}

where $C>0$ is a suitable, sufficiently large constant. Then approximately, the function $u \mapsto \frac{1}{N}\sum_j J(x_i,u,x_j)$ is maximal at $u=\frac{x_j-x_i}{\|x_j-x_i\|}$ where j is the agent that is closest to i (if it is closer than C and not closer than $\delta$ , in particular the last term in J is added to avoid that $u=0$ is the maximiser). Therefore, approximately, we have created a model, where agents are attracted with unit speed by their nearest neighbours, but not by other agents. Such a nonlinear trade-off between the influences of the other agents is not obtainable by a Newtonian model.

Well-posedness and mean-field limit of the undisclosed fast-reaction system. Similar to the full model (2.2), we are also interested in a mean-field limit of the undisclosed fast-reaction limit as the number of agents tends to infinity, $N \to \infty$ . A formal limiting procedure leads to the equation

(2.17) \begin{equation} \partial_t \mu(t) +\textrm{div} \left( v^J(\mu(t))\cdot \mu(t) \right) = 0\end{equation}

where for $x \in \mathbb{R}^d$ and $\nu \in \mathcal{M}_1(\mathbb{R}^d)$ we set

(2.18) \begin{align} v^J(\nu)(x) & \;:=\; \int_U e(x,u)\sigma^J(\nu)(x, u)\,d\eta(u) \qquad\qquad\,\,\end{align}

(2.19) \begin{align} \sigma^J(\nu)(x,\cdot) & \;:=\; \dfrac{\exp\left(\dfrac{1}{\varepsilon}\int_{\mathbb{R}^d} J(x,\cdot,x^{\prime})\,d\nu(x^{\prime})\right)}{ \int_U \exp\left(\dfrac{1}{\varepsilon}\int_{\mathbb{R}^d} J(x,v,x^{\prime})\,d\nu(x^{\prime})\right)\!\textrm{d} \eta(v)}.\end{align}

Given $\mathbf{x} = \left(x_1,\dots,x_N\right) \in [\mathbb{R}^d]^N$ and setting $\mu^N = \frac1N \sum_{j=1}^N \delta_{x_j}$ , (2.15), (2.18) and (2.19) are related through

\begin{align*}\sigma_i^J(\mathbf{x}) = \sigma_i^J\!\left(x_1,\dots,x_N\right) = \sigma^J\!\left(\mu^N\right)(x_i,\cdot) \quad \text{and} \quad v_i^J(\mathbf{x}) = v_i^J\!\left(x_1,\dots,x_N\right) = v^J\!\left(\mu^N\right)(x_i).\end{align*}

The key ingredient for the study of (2.15) and its limiting behaviour is to establish Lipschitz continuity of the fast-reaction equilibrium strategies and velocities with respect to payoff function and particle locations.

Proof. Let us define two continuous functions $g, g^{\prime} \colon U \to \mathbb{R}$ as

\begin{equation*} \begin{aligned} g(u) = \int_{\mathbb{R}^d} J(x, u, y)\,d\mu(y) \quad \text{and} \quad g^{\prime}(u) = \int_{\mathbb{R}^d} J^{\prime}(x^{\prime}, u, y)\,d\mu^{\prime}(y) \end{aligned} \end{equation*}

For every $u \in U$ , using $M \geq \|J\|_{\infty}$ and $M \geq \|J^{\prime}\|_{\infty}$ , we immediately obtain the global bounds $-M \leq g(x) \leq M$ and $-M \leq g^{\prime}(x) \leq M$ . Using the triangle inequality and the dual definition for $W_1$ in (1.9), we also estimate

(2.22) \begin{equation} |g(u)-g^{\prime}(u)| \leq C\left( W_1(\mu,\mu^{\prime}) + \|J - J^{\prime}\|_\infty + \|x-x^{\prime}\|\right) \end{equation}

for $C = C(M)$ . Recall now that

\begin{equation*} \sigma^{J}(\mu)(x,u) = \frac{\exp( g(u)/\varepsilon )}{\int_U \exp\left( g(v)/\varepsilon \right)\,d\eta(v)} \quad \text{and} \quad \sigma^{J^{\prime}}(\mu^{\prime})(x^{\prime},u) = \frac{\exp( g^{\prime}(u)/\varepsilon )}{\int_U \exp\left( g^{\prime}(v)/\varepsilon \right)\,d\eta(v)}. \end{equation*}

Thus, (2.20) follows combining global boundedness of g, g ^′ with (2.22), the uniform bounds $1/C < \sigma^{J}(\mu)(x,u) < C$ are obtained in the same way, while (2.21) follows from (2.20) combined with $e \in \textrm{Lip}_b(\mathbb{R}^d \times U; \mathbb{R}^d)$ .

For the discrete fast-reaction system, the following Lemma adapts the estimate (2.21) as one in terms of discrete particle locations and their velocities.

Lemma 2.11 (Map to undisclosed fast-reaction agent velocity is Lipschitz continuous). For any $N > 0$ , consider the map

\begin{equation*} \mathbf{v}^J \colon \left[\mathbb{R}^d\right]^N \to \left[\mathbb{R}^d\right]^N, \quad \left(x_1,\dots,x_N\right) \mapsto \left(v_1^J\!\left(x_1,\dots,x_N\right), \dots, v_N^J\!\left(x_1,\dots,x_N\right)\right) \end{equation*}

associated with the fast-reaction ODE (2.15). Then, $\mathbf{v}^J$ is Lipschitz continuous under the distance induced by $\|\cdot\|_N$ , with Lipschitz constant $L = L(J, e, \varepsilon)$ .

Proof. Fix any $\mathbf{x},\mathbf{x}^{\prime} \in \left[\mathbb{R}^d\right]^N$ and define $\mu = \frac{1}{N} \sum_{j=1}^N \delta_{x_i}$ and $\mu^{\prime} = \frac{1}{N} \sum_{j=1}^N \delta_{x^{\prime}_i}$ . Then, by Lemma 2.10, with $J^{\prime} = J$ , $x = x_i$ and $x^{\prime} = x^{\prime}_i$ , we have

\begin{equation*} \left\| v^{J}(\mu)(x_i) - v^{J}(\mu^{\prime})\left(x^{\prime}_i\right) \right\| \leq C\left( W_1(\mu,\mu^{\prime}) + \|x_i-x^{\prime}_i\|\right) \end{equation*}

for every $i = 1,\dots, N$ , with $C = C(J,e,\varepsilon)$ . Using the definition of $W_1$ in (1.8), we observe

\begin{equation*} W_1(\mu,\mu^{\prime}) \leq \frac{1}{N}\sum_{j=1}^N \|x_i - x^{\prime}_i\| = \|x - x^{\prime}\|_N \end{equation*}

so that

\begin{align} &\|\mathbf{v}^J(\mathbf{x})-\mathbf{v}^J\!\left(\mathbf{x}^{\prime}\right)\!\|_N =\frac{1}{N}\sum_{i=1}^N \| v_i^J\!\left(x_1,\dots,x_N\right) - v_i^J\!\left(x^{\prime}_1,\dots,x^{\prime}_N\right) \| \nonumber \\ &\quad = \frac{1}{N}\sum_{i=1}^N\left\| v^{J}(\mu)(x_i) - v^{J}(\mu^{\prime})\left(x^{\prime}_i\right) \right\|\leq \frac{C}{N} \sum_{i=1}^N \left( \|x-x^{\prime}\|_N + \|x_i-x^{\prime}_i\| \right) = 2C \|\mathbf{x}-\mathbf{x}^{\prime}\|_N \nonumber \end{align}

which is the sought-after estimate.

After clarifying what we mean by a solution of (2.17), we summarise the mean-field result in Theorem 2.13, whose proof then builds on Lipschitz continuity of the velocity field in the discrete system and follows by fairly standard arguments (see, e.g. [Reference Cañizo, Carrillo and Rosado12, Reference Carrillo, Choi and Hauray14]).

Definition 2.12. We say a curve $\mu \in C([0,T]; (\mathcal{M}_1(\mathbb{R}^d), W_1) )$ solves (2.17) if $\mu(t)$ has uniformly compact support for $t \in [0,T]$ and

\begin{align*} \frac{d}{dt} \int_{\mathbb{R}^d} \phi(x)\textrm{d}\mu(t)(x) = \int_{\mathbb{R}^d} \nabla\phi(x) \cdot v^J(\mu(t))(x)\,\textrm{d}\mu(t)(x) \quad \text{for every $\phi \in C^\infty_c(\mathbb{R}^d)$}.\end{align*}

Proof outline. The proof is rather standard and builds on the Lipschitz continuity of the finite agents system (2.15). We highlight the main steps without in-depth details and refer to Proposition 2.1 and Theorem 2.4 of [Reference Bongini, Fornasier, Hansen and Maggioni10], and references therein, for further information.

Part 1: finite agents setting. For given $\bar{\mathbf{x}}^N = (\bar{x}_{1}^N,\dots, \bar{x}_{N}^N) \in [B^d({\bar{R}})]^N$ , using the Lipschitz continuity provided in Lemma 2.11, there exists a unique curve $\mathbf{x}^N \colon [0,T] \to [\mathbb{R}^d]^N$ of class $C^1$ solving (2.15) with $\mathbf{x}^N(0) = \bar{\mathbf{x}}^N$ . A direct estimate provides $\|\partial x_i^N(t)\| \leq \|e\|_\infty$ , so that $\|x_i^N(t)\| \leq \bar{R} + T \cdot \|e\|_\infty$ for every $t \in [0,T]$ , $i = 1,\dots,N$ . Furthermore, setting $\mu^N(t) \;:=\; \tfrac{1}{N} \sum_{i=1}^N \delta_{x^N_i(t)}$ one can explicitly verify that $\mu^N$ solves (2.17) with initial condition $\bar{\mu}^N \;:=\; \tfrac{1}{N} \sum_{i=1}^N \delta_{\bar{x}^N_i}$ . This establishes point .

Part 2: mean-field solution as limit of finite agents solutions. Fix $\bar{\mu} \in \mathcal{M}_1(B^d({\bar{R}}))$ . For $N > 0$ , let $(\bar{\mu}^N)_N \subset \mathcal{M}_1(B^d({\bar{R}}))$ be a sequence of empirical measures such that $W_1(\bar{\mu}^N,\bar{\mu}) \to 0$ as $N \to \infty$ and let $(\mu^N)_N$ be the respective sequence of solutions of (2.17) with initial conditions $\bar{\mu}^N$ . This sequence of curves $(\mu^N)_{N \in \mathbb{N}} \subset C([0,T]; (\mathcal{M}_1(B^d({R})), W_1))$ is equicontinuous and equibounded. Therefore, an application of the Ascoli–Arzelà theorem provides a cluster point $\mu \in C([0,T]; (\mathcal{M}_1(B^d({R})), W_1))$ such that, up to a subsequence, we have

\begin{equation*} \lim_{N\to\infty} W_1\!\left(\mu^N(t), \mu(t)\right) = 0 \quad \text{ uniformly for } t \in [0,T]. \end{equation*}

Invoking Lemma 2.10 for $x=x^{\prime}$ and $J=J^{\prime}$ this implies $\| v^J\!\left(\mu^N(t)\right)(x) - v^J(\mu(t))(x)\| \to 0$ uniformly in $t \in [0,T]$ , $x \in B^d({R})$ as $N \to \infty$ . Consequently, the cluster point $\mu$ is a solution to (2.17) with $\bar{\mu}$ as initial condition. This establishes the existence part of point .

Part 3: stability estimates. For fixed N, a stability estimate of the form $\|\mathbf{x}^N_1(t)-\mathbf{x} ^N_2(t)\|_N \leq \exp(C\,(t-s)) \cdot \|\mathbf{x}^N_1(s)-\mathbf{x} ^N_2(s)\|_N$ for solutions to the discrete system (2.15) follows quickly from Grönwall’s lemma. Since the $\|\cdot\|_N$ -norm between two point clouds provides an upper bound for the $W_1$ distance between the respective empirical measures, this would provide point for empirical measures. The extension to arbitrary measures can be done as in the proof of Theorem 2.4 in [Reference Bongini, Fornasier, Hansen and Maggioni10]. This then provides uniqueness of the solution $\mu$ to (2.17) for initial condition $\bar{\mu} \in \mathcal{M}_1(B^d({\bar{R}}))$ , completing point .

Finally, we establish convergence of (2.13) to (2.15) (and their respective mean-field versions) as $\lambda \to \infty$ . The proof is given in Section A.2 of the Appendix, a simple numerical example is shown in Figure 3.

Figure 3. Convergence of the undisclosed model to the fast-reaction limit as $\lambda \to \infty$ . The payoff function encourages agents to move to the origin, but penalises small pairwise distances (see Example 4.1 for a detailed description). With $\lambda$ small, agents cannot adjust their strategies fast enough. They overshoot the origin and completely fail to avoid each other. The situation improves as $\lambda$ increases. For $\lambda=100,$ the model closely resembles the fast-reaction limit, in accordance with (2.24).

3.1 Discussion

Our motivation is as follows: we observe a set of rational agents, for example, pedestrians in a confined space, and want to learn about their interaction rules. As argued above (e.g. in Examples 2.7 and 2.8), undisclosed fast-reaction entropic games can be used to parametrise a rich class of interaction behaviours, in particular subsuming Newtonian models. We may therefore hope that, when feeding our observations into an inference functional for J, that by analysing and interpreting the resulting payoff function we can learn something about the interaction rules of the agents. We focus on the inference of J and assume that U, e and $\varepsilon$ are known. We will demonstrate during the numerical examples (Section 4) that it is usually possible to make plausible choices for U, e and that qualitatively equivalent choices yield qualitatively equivalent results for the inferred J. One may also assume that $\varepsilon=1$ , since re-scaling $\varepsilon$ can be compensated by re-scaling J in the same way, see (2.15).

In Section 3.2, we discuss a differential inference functional, proposed in analogy to (1.7), for the full entropic game model, that is, without the undisclosed assumption and not in the fast-reaction limit. For this we assume that we are able to observe the full dynamics, that is, physical locations and velocities of the agents, as well as their mixed strategies and temporal derivatives, $(x_i,\partial_t x_i,\sigma_i,\partial_t \sigma_i)_{i=1}^N$ and propose a functional that seeks to infer J by comparing its predictions with the observed $\partial_t \sigma_i$ .

In many cases, it may not be possible to observe the mixed strategies $\sigma_i$ , let alone their temporal derivatives $\partial_t \sigma_i$ , since these may correspond to internal states of the agents that we, as external observers, cannot perceive. In Section 3.3, we turn to the undisclosed fast-reaction setting, where physical velocities and mixed strategies are functions of current physical locations, and show that there one can still perform inference with limited observations. In Section 3.3.1, we propose an inference functional, again in analogy to (1.7), that only requires the observation of the physical locations $x_i$ and velocities $\partial_t x_i$ . Although it is nonconvex, we observe in our numerical examples that it provides plausible payoff functions J, indicating that nonconvexity does not seem to be a practical issue.

In Section 3.3.2, we propose a convex inference functional that requires knowledge of the mixed strategies $\sigma_i$ , but not of their temporal derivatives. This may be practical in two cases: In Section 4.1, we propose a scheme to provide such data from observed physical velocities, intuitively by ‘inverting’ the approximation scheme of Example 2.7. Also, mixed strategies could represent a stochastic movement of agents and the ‘smooth’ physical velocity $\partial_t x_i$ is merely the result of an averaging over a smaller time scale. Should we indeed be able to observe the fluctuating agent dynamics at small time scales, this could be used to approximate the mixed strategies. This was one of our motivations for providing Example 4.3.

For both functionals, there are natural candidates for the mean-field limit. We analyse these two functionals in more detail in Section 3.3.3, providing an estimate of the approximation quality (Theorem 3.3), as well as existence of minimisers and consistency in the mean-field limit (Theorem 3.6).

Now, we formalise our notion of admissible observations for inference and a notion of consistency of observations in the mean-field limit.

Assumption 3.1 (Admissible observations).

1. 1. Discrete observation: For a fixed number of agents N, a time horizon $T \in (0,\infty)$ , and some radius $R \in (0,\infty)$ we observe the agents’ physical paths $\textbf{x}^N \;:=\; \left(x^N_1,\ldots,x^N_N\right) \in C^1([0,T],B^d(R))^N$ with velocities $\textbf{v}^N \;:=\; \left(v^N_1,\ldots,v^N_N\right)$ , $v^N_i = \partial_t x^N_i$ .

   Optionally, in addition we may also observe the agents’ mixed strategies $\boldsymbol{\sigma}^N \;:=\; \left(\sigma^N_1,\ldots,\sigma^N_N\right) \in C([0,T],\mathcal{S}_{a,b})^N$ for some bounds $0<a\leq b<\infty$ in the definition of $\mathcal{S}_{a,b}$ , see (2.1). The mixed strategies are consistent with the observed velocities, that is,

   (3.1) \begin{equation}v_i^N(t) = \int_U e\left(x_i^N(t),u\right)\sigma_i^N(t)(u)\,\textrm{d} \eta(u).\end{equation}
   Note: The assumption that the agent’s velocity is exactly consistent with their mixed strategy may seem unrealistic, as the velocity may be subject to other external influences and noise. However, we will subsequently often assume that the mixed strategies are not observed directly, but are only inferred from the velocities. In this case, satisfying assumption (3.1) is quite natural. In other cases, assumption (3.1) will be used in Theorem 3.3 to bound the error on trajectories based on the error on strategies. If (3.1) only holds approximately, an additional corresponding error term will appear there.

2. 2. Consistent mean-field behaviour: For fixed T, R, for an increasing sequence of N we make observations as specified in part , and there is a limit observation $t \mapsto \mu^\infty(t) \in \mathcal{M}_1(B^d(R))$ with a velocity field $t \mapsto v^\infty(t)$ such that

   (3.2) \begin{align} W_1\!\left(\mu^N(t),\mu^\infty(t)\right) \rightarrow 0 \quad \text{uniformly in $t \in [0,T]$ where}\quad \mu^N(t) \;:=\; \frac1N \sum_{i=1}^N \delta_{x^N_i(t)},\end{align}
   (3.3) \begin{align} \int_0^T \frac{1}{N} \sum_{i=1}^N \left \langle \phi\left(t,x^N_i(t)\right), v^N_i(t) \right \rangle \,\textrm{d} t \rightarrow \int_0^T \int_{\mathbb{R}^d} \left \langle \phi(t,x) , v^\infty(t)(x) \right \rangle \,\textrm{d} \mu^\infty(t)(x)\,\textrm{d} t \end{align}
   for all $\phi \in \left[C\!\left([0,T] \times \mathbb{R}^d\right)\right]^d$ , and
   (3.4) \begin{align} \int_0^T \frac{1}{N} \sum_{i=1}^N \|v^N_i(t)\|^2 \,\textrm{d} t \rightarrow \int_0^T \int_{\mathbb{R}^d} \|v^\infty(t)(x)\|^2 \,\textrm{d} \mu^\infty(t)(x)\,\textrm{d} t<\infty.\end{align}
   (3.2) implies that physical locations are consistent, (3.3) implies that observed velocities converge in a weak sense, and (3.4) implies that they are consistent. Intuitively, for (3.4) to hold, mass particles that converge to the same limit location x as $N \to \infty$ need to converge to the same velocity, otherwise Jensen’s strict inequality contradicts (3.4).

   If we also observe mixed strategies in part , then the bounds $0<a \leq b<\infty$ are uniform in N and there also is a mixed strategy mean-field $(t,x) \mapsto \sigma^\infty(t)(x) \in \mathcal{S}_{a,b}$ such that

(3.5) \begin{multline} \int_0^T \frac{1}{N} \sum_{i=1}^N \int_U \phi\left(t,x^N_i(t),u\right) \cdot \sigma^N_i(t)(u)\, \textrm{d} \eta(u) \,\textrm{d} t\\ \rightarrow \int_0^T \int_{\mathbb{R}^d} \int_U \phi(t,x,u) \cdot \sigma^\infty(t)(x)(u) \textrm{d} \eta(u)\,\textrm{d} \mu^\infty(t)(x)\,\textrm{d} t \end{multline}

for all $\phi \in \left[C\!\left([0,T] \times \mathbb{R}^d \times U\right)\right]^d$ and

(3.6) \begin{align} & \int_0^T \frac{1}{N} \sum_{i=1}^N \int_U \sigma^N_i(t)(u)\log\!\left(\sigma^N_i(t)(u)\right)\, \textrm{d} \eta(u) \,\textrm{d} t\nonumber\\[4pt] & \qquad\qquad \rightarrow \int_0^T \int_{\mathbb{R}^d} \int_U \sigma^\infty(t)(x)(u)\log(\sigma^\infty(t)(x)(u)) \textrm{d} \eta(u)\,\textrm{d} \mu^\infty(t)(x)\,\textrm{d} t.\end{align}

These are direct equivalents of (3.3) and (3.4).

In particular, observations are admissible when they were generated by an entropic game model with some ground truth payoff J.

The proof is quite straightforward and provided in Section A.3.

In the following, we assume that observations are admissible in the sense of Assumption 3.1.

3.2 Differential inference functional

In this section, we discuss an inference functional for payoff function J, in close analogy to (1.7). For now, we assume that in addition to $\left(\textbf{x}^N, \textbf{v}^N, \boldsymbol{\sigma}^N\right)$ we can even observe $\partial_t \sigma_i^N$ for all $i \in \{1,\ldots,N\}$ . Our differential inference functional is therefore aimed at recovering J, by comparing the observed $\partial_t \sigma_i^N$ with the predictions by the model (2.2b). In (1.7), the discrepancy between observed $v_i^N$ and predicted velocities is penalised by the squared Euclidean distance. As metric to compare the observed $\partial_t \sigma_i^N$ and the prediction by (2.2b), we choose the (weak) Riemannian tensor of the Hellinger–Kakutani and Fisher–Rao metrics. That is, we set for a base point $\sigma \in \mathcal{S}_{a,b}$ , and two tangent directions $\delta \mu, \delta \nu \in L^p_{\eta}(U)$ attached to it,

(3.7) \begin{align} d_\sigma(\delta \mu,\delta \nu)^2 \;:=\; \frac{1}{4} \int_U \frac{\left(\delta \mu- \delta \nu\right)^2}{\sigma}\,\textrm{d} \eta.\end{align}

A potential inference functional for J could then be:

(3.8) \begin{align} & \mathcal{E}^{N}_{\dot{\sigma}}({J}) \;:=\; \nonumber\\ & \frac{1}{T} \int_0^T \left[ \frac{1}{N} \sum_{i=1}^N d_{\sigma_i^N(t)}\left( \partial_t \sigma_i^N(t), \frac{1}{N} \sum_{j=1}^N f^{{J}}\left(x_i^N(t),\sigma_i^N(t),x_j^N(t),\sigma_j^N(t)\right) +f^\varepsilon\left(\sigma_i^N(t)\right) \right)^2 \right] \textrm{d} t.\end{align}

Minimisation should be restricted to a sufficiently regular and compact class of J. Since $d_\sigma(\cdot,\cdot)^2$ is quadratic in its arguments and $f^{{J}}$ is linear in J, the mismatch term is again quadratic in J and thus restricted to a suitable class, this is a convex optimisation problem.

In Section 4.3, we provide a simple example where we simulate trajectories $(\textbf{x}^N,\boldsymbol{\sigma}^N)$ according to the dynamics (2.2), based on a given J, and then obtain a corresponding minimiser $\kern1pt\hat{{\kern-1pt}J}$ of (3.8). We then show that trajectories simulated with J and $\kern1pt\hat{{\kern-1pt}J}$ are close even when the initial configurations are not drawn from training data used to infer $\kern1pt\hat{{\kern-1pt}J}$ . Indeed, upon selecting a suitable functional framework, one can in principle extend the analysis in [Reference Bongini, Fornasier, Hansen and Maggioni10] (and the analysis presented here below) to prove existence of minimisers, limiting behaviour as $N \to \infty$ and ability of the inferred model to generalise.

However, application to real data seems challenging, as one usually only observes $\left(\textbf{x}^N, \textbf{v}^N\right)$ but not $\boldsymbol{\sigma}^N$ , let alone variations $\partial_t \sigma_i^N$ . To address this, we propose in the next section two new inference functionals for the undisclosed fast-reaction setting.

3.3 Undisclosed fast-reaction inference functionals

3.3.1 Penalty on velocities

The inference situation is simplified in the undisclosed fast-reaction setting since now, for fixed payoff function J, mixed strategies and physical velocities are a direct function of physical locations. We can therefore attempt to minimise the discrepancy between observed physical velocities and those predicted by J.

For a discrete admissible observation $\left(\textbf{x}^N,\textbf{v}^N\right)$ , as in Assumption 3.1, we define the inference functional on hypothetical payoff functions ${J} \in X$

(3.9a) \begin{align} \mathcal{E}^{N}_{v}({J}) & \;:=\; \frac{1}{T} \int_0^T \left[ \frac{1}{N} \sum_{i=1}^N \left\|v_i^N(t)-v_i^{{J}}\left(x_1^N(t),\ldots,x_N^N(t)\right) \right\|^2 \right] \textrm{d} t. \end{align}

The natural candidate for the limit functional is

(3.9b) \begin{align} \mathcal{E}_v({J}) & \;:=\; \frac{1}{T} \int_0^T \int_{\mathbb{R}^d} \left\| v^\infty(t)(x) - v^{{J}}\left(\mu^\infty(t)\right)(x) \right\|^2 \textrm{d} \mu^\infty(t)(x) \textrm{d} t, \end{align}

where the couple $\left(\mu^\infty, v^\infty\right)$ is the corresponding limit of $\left(\textbf{x}^N,\textbf{v}^N\right)_N$ as introduced in Assumption 3.1, point .

This functional is intuitive and requires only the observations of $\left(\textbf{x}^N,\textbf{v}^N\right)$ and no observations about the mixed strategies. However, it is not convex. This could potentially lead to poor local minima, but we did not observe such problems in our numerical examples.

The motivation for studying mean-field inference functionals is twofold. First, it establishes asymptotic consistency of inference in the limit of many agents. Second, the mean-field equation yields an approximation for the expected behaviour of the finite-agent system under many repetitions with random initial conditions. While we do not prove this approximation property, the existence of a consistent mean-field inference functional is still an encouraging indicator that inference over the collection of many finite-agent samples will yield a reasonable result.

3.3.2 Penalty on mixed strategies

In the case where information about the mixed strategies of the agents is available (see Section 3.1 for a discussion), the discrepancy between observed mixed strategies and those predicted by J can be used for inference. We choose to measure this discrepancy with the Kullback–Leibler (KL) divergence.

Similar to above, for a discrete admissible observation $\left(\mathbf{x}^N,\mathbf{v}^N,\boldsymbol{\sigma}^N\right)$ , we define the inference functional for hypothetical payoff functions ${J} \in X$

(3.10a) \begin{align} \mathcal{E}^{N}_\sigma({J}) & \;:=\; \frac{1}{T} \int_0^T \left[ \frac{1}{N} \sum_{i=1}^N \textrm{KL}\left(\sigma_i^N(t)|\sigma_i^{J}\left(x_1^N(t),\ldots,x_N^N(t)\right)\right) \right] \textrm{d} t. \end{align}

The natural candidate for the limit functional is

(3.10b) \begin{align} \mathcal{E}_\sigma({J}) & \;:=\; \frac{1}{T} \int_0^T \int_{\mathbb{R}^d} \textrm{KL}\big(\sigma^\infty(t)(x,\cdot)|\sigma^{{J}}(\mu^\infty(t))(x,\cdot)\big) \textrm{d}\mu^\infty(t)(x)\, \textrm{d} t \nonumber \\[5pt] & = \frac1T \int_0^T \int_{\mathbb{R}^d} \int_U \log\left( \frac{\sigma^\infty(t)(x,u)}{\sigma^{{J}}(\mu^\infty(t))(x,u)} \right) \sigma^\infty(t)(x,u) \textrm{d} \eta(u) \textrm{d} \mu^\infty(t)(x)\,\textrm{d} t, \end{align}

where $(\mu^\infty, \sigma^\infty)$ are as introduced in Assumption 3.1, point . Due to the particular structure of stationary mixed strategies, (2.15), and the KL divergence, functionals (3.10) are convex on X (Proposition A.5), which is a potential advantage over (3.9).

3.3.3 Analysis of the inference functionals

We now provide some theoretical results about the inference functionals of Sections 3.3.1 and 3.3.2. The first result establishes that the obtained minimal inference functional value provides an upper bound on the accuracy of the trajectories that are simulated with the inferred $\kern1pt\hat{{\kern-1.5pt}J}$ (the proof builds upon the proof of [Reference Bongini, Fornasier, Hansen and Maggioni10, Proposition 1.1] and is given in Section A.3).

Theorem 3.3. Let $\kern1pt\hat{{\kern-1pt}J} \in X$ and $\left(\mathbf{x}^N,\mathbf{v}^N,\boldsymbol{\sigma}^N\right)$ be an admissible discrete observation for N agents (cf. Assumption 3.1, point ). Let $\hat{\mathbf{x}}^N$ be the solution of (2.15) induced by $\kern1pt\hat{{\kern-1pt}J}$ for the initial condition $\hat{\textbf{x}}^N(0) = \textbf{x}^N(0) = \left(\bar{x}_{1}^N,\dots,\bar{x}_{N}^N\right) \in [\mathbb{R}^d]^N$ . Then,

\begin{equation*} \|\mathbf{x}^N(t)-\hat{\mathbf{x}}^N(t)\|_N \leq C \sqrt{\mathcal{E}^{N}_v(\kern1pt\hat{{\kern-1pt}J})} \quad \text{and} \quad \|\mathbf{x}^N(t)-\hat{\mathbf{x}}^N(t)\|_N \leq C \sqrt{\mathcal{E}^{N}_\sigma(\kern1pt\hat{{\kern-1pt}J})} \end{equation*}

for all $t \in [0,T]$ , with $C = C(T, \kern1pt\hat{{\kern-1pt}J}, e, \varepsilon)$ . Analogously, let $(\mu^\infty, \sigma^\infty)$ as introduced in Assumption 3.1, point and let $\hat{\mu}$ be the solution of (2.17) induced by $\kern1pt\hat{{\kern-1pt}J}$ for the initial condition $\bar{\mu} = \mu^\infty(0)$ . Then,

\begin{equation*} W_1(\mu^\infty(t),\hat{\mu}(t)) \leq C \sqrt{\mathcal{E}_v(\kern1pt\hat{{\kern-1pt}J})} \quad \text{and} \quad W_1(\mu^\infty(t),\hat{\mu}(t)) \leq C \sqrt{\mathcal{E}_\sigma(\kern1pt\hat{{\kern-1pt}J})} \end{equation*}

for all $t \in [0,T]$ and the same constant C as above.

Next, we address the existence of minimising payoff functions, both in theory and numerical approximation. At the theoretical level, we need to ensure compactness of minimising sequences, which we obtain here by restriction to a suitable compact space. At the numerical level, we are interested in finite-dimensional approximations of this space that asymptotically are dense as the discretisation is refined.

Remark 3.4 (Compactness and finite-dimensional approximation). In order to obtain compactness, we restrict the variational problems to suitable subspaces of X. In particular, for $R, M > 0$ , let us define

\begin{equation*}X_{R,M} = \left\{ J \in \textrm{Lip}_b\!\left(B^d({R}) \times U \times B^d({R})\right) \mid \|J\|_{\infty} + \textrm{Lip}(J) \leq M \right\}.\end{equation*}

The parameter R bounds the learning domain in space: inference where we have no data available is simply meaningless. The parameter M will ensure compactness with respect to uniform convergence.

For each $R,M > 0,$ we consider a family of closed convex subsets $\left(X_{R,M}^N\right)_{N \in \mathbb{N}} \subset X_{R,M}$ satisfying the uniform approximation property: for every $J \in X_{R,M}$ there exists a sequence $(J^N)_{N \in \mathbb{N}} $ uniformly converging to J with $J^N \in X_{R,M}^N$ for each $N \in \mathbb{N}$ . These approximating spaces can be selected to be finite dimensional.

An alternative way to obtain compactness would be via the introduction of additional regularising functionals. This is discussed in Section 4.1. Conceptually, both approaches serve the same purpose and we find that the former is more convenient for analysis while the latter is more attractive numerically.

The following key lemma ensures convergence of the error functionals under uniform convergence of the payoff functions and/or $W_1$ convergence of measures.

Lemma 3.5. Let $\left(\mathbf{x}^N,\mathbf{v}^N,\boldsymbol{\sigma}^N\right)_N$ be a family of admissible observations and $\mathcal{E}^{N}_v$ , $\mathcal{E}_v$ , $\mathcal{E}^{N}_\sigma$ and $\mathcal{E}_\sigma$ be the corresponding discrete and continuous inference functionals (3.9) and (3.10). For $M>0$ , let $\left({J}^N\right)_N$ be a sequence of payoff functions in $X_{R,M}$ , converging uniformly to some ${J} \in X_{R,M}$ as $N \to \infty$ . Then,

(3.11) \begin{align} \lim_{N\to\infty} \mathcal{E}^{N}_v\!\!\left({J}^N\right) = \mathcal{E}_v({J}) \quad \text{ and } \quad \lim_{N\to\infty} \mathcal{E}^{N}_\sigma\!\!\left({J}^N\right) = \mathcal{E}_\sigma({J}). \end{align}

Proof. Thanks to Assumption 3.1, point , we have $W_1\left(\mu^N(t),\mu^\infty(t)\right) \to 0$ uniformly in $t \in [0,T]$ as $N \to \infty$ . Lemma 2.10 provides

\begin{align*} \left\| v^{{J}^N}\left(\mu^N(t)\right)(x) - v^{{J}}(\mu^\infty(t))(x) \right\| & \leq C \cdot \left( W_1\left(\mu^N(t),\mu^\infty(t)\right) + \|\kern1pt\hat{{\kern-1pt}J}^N-\kern1pt\hat{{\kern-1pt}J}\| \right) \end{align*}

uniformly in $x \in B^d({R}), t \in [0,T]$ and N. Hence, $v^{{J}^N}\!\left(\mu^N(t)\right)(x)$ converges uniformly to $v^J(\mu^\infty(t))(x)$ in $C([0,T]\times B^d({R}); \mathbb{R}^d)$ . Since $\|v^J(\mu^\infty(t))(x)\| \leq \|e\|_\infty$ for every $x \in B^d({R})$ and $t \in [0,T]$ , we also have $\|v^{J^N}\left(\mu^N(t)\right)(x) \|^2 \to \|v^J(\mu^\infty(t))(x)\|^2$ uniformly in $C([0,T]\times B^d({R}))$ . Recalling that $v_i^{{J}^N}\left(x_1^N(t),\ldots,x_N^N(t)\right)(x) = v^{{J}^N}\left(\mu^N(t)\right)\left(x_i^N(t)\right)$ by definition and using the convergence of velocities provided by Assumption 3.1, point , the result follows by passing to the limit in the definition:

\begin{align*} \lim_{N \to \infty} \mathcal{E}^{N}_v\!\left({J}^N\right) &= \lim_{N \to \infty} \frac{1}{T} \int_0^T \left[ \frac{1}{N} \sum_{i=1}^N \left\|v_i^N(t)-v_i^{{J}^N}\left(x_1^N(t),\ldots,x_N^N(t)\right) \right\|^2 \right] \textrm{d} t \\[4pt] &= \frac{1}{T} \int_0^T \int_{\mathbb{R}^d} \left\|v^\infty\left(\mu^\infty(t)\right)(x) - v^{{J}}\left(\mu^\infty(t)\right)(x)\right\|^2 \textrm{d} \mu^\infty(t)(x) \textrm{d} t = \mathcal{E}_v({J}). \end{align*}

An analogous derivation applies to the $\sigma$ -energies (3.10) where we use the Lipschitz estimate (2.20) and the $\sigma$ -part of Assumption 3.1, point .

We are now ready to state the main result concerning the inference functionals and convergence of minimisers for $\mathcal{E}_v^{N}$ and $\mathcal{E}_\sigma^{N}$ .

Theorem 3.6. Let $\left(\mathbf{x}^N,\mathbf{v}^N,\boldsymbol{\sigma}^N\right)_N$ be a family of admissible observations and fix $M > 0$ . Then:

1. 1. The functionals (3.9a) and (3.10a) have minimisers over $X_{R,M}^N$ and (3.9b) and (3.10b) have minimisers over $X_{R,M}$ .

2. 2. For $N > 0$ , let

   \begin{equation*} {J}^N_v \in arg min_{{J} \in X_{R,M}^N} \mathcal{E}^{N}_v({J}) \quad {and} \quad {J}^N_\sigma \in arg min_{{J} \in X_{R,M}^N} \mathcal{E}^{N}_\sigma({J}). \end{equation*}
   The sequence $\left({J}^N_v\right)_N$ , resp. $\left({J}^N_\sigma\right)_N$ , has a subsequence converging uniformly to some continuous function ${J}_v \in X_{R,M}$ , resp. ${J}_\sigma \in X_{R,M}$ , such that
   (3.12) \begin{align} \lim_{N \to \infty} \mathcal{E}^{N}_v\!\left({J}^N_v\right) = \mathcal{E}_v({J}_v) \quad {and} \quad \lim_{N \to \infty} \mathcal{E}^{N}_\sigma\left({J}^N_\sigma\right) = \mathcal{E}_\sigma({J}_\sigma). \end{align}
   In particular, ${J}_v$ (resp. ${J}_\sigma$ ) is a minimiser of $\mathcal{E}_v^J$ (resp. $\mathcal{E}_\sigma^J$ ) over $X_{R,M}$ .

Proof. Part 1: existence of minimisers. All functionals (3.9) and (3.10) are continuous in J with respect to uniform convergence which follows from Lemma 3.5 setting the initial measure $\bar{\mu}$ to be the empirical one (i.e. just looking at converging J for fixed trajectories). The sets $X_{R,M}$ and $X_{R,M}^N$ are compact with respect to uniform convergence. Therefore, minimisers exist.

Part 2: convergence of minimisers. The proof is inspired by [Reference Bongini, Fornasier, Hansen and Maggioni10, Section 4.2]. Since $X_{R,M}^N \subset X_{R,M}$ and the latter is compact, the sequence of minimisers $\left({J}_v^N\right)_N$ with ${J}^N_v \in X_{R,M}^N$ must have a cluster point ${J}_v$ in $X_{R,M}$ . For simplicity, denote the converging subsequence again by $\left({J}^N_v\right)_N$ . Let $\tilde{J} \in X_{M,R}$ . By the uniform approximation property of $X_{R,M}^N$ let $\left(\tilde{J}^N\right)_N$ be a sequence converging uniformly to $\tilde{J}$ and $\tilde{J}^N \in X_{R,M}^N$ . By Lemma 3.5 and optimality of each ${J}^N_v$ ,

\begin{align*} \mathcal{E}_v({J}_v) = \lim_{N\to\infty} \mathcal{E}^{N}_v\!\left({J}^N_v\right) \leq \lim_{N\to\infty} \mathcal{E}^{N}_v\!\left(\tilde{J}^N\right) = \mathcal{E}_v(\tilde{J}). \end{align*}

Therefore, ${J}_v$ minimises $\mathcal{E}_v$ and the optimal values converge. The same argument covers the $\sigma$ -functionals (3.10) and the corresponding sequence of minimisers $\left({J}_\sigma^N\right)_N$ .

We conclude this section with a couple of comments and remarks on the interpretation of results in Theorem 3.6.

Remark 3.7 (Realisability of data). In general, observed data may or may not be realisable by the model. This distinguishes two regimes.

• • Assume there is a true regular J that generated the observed data (as for example in Lemma 3.2). If we pick M large enough in Theorem 3.6, so that $J \in X_{R,M}$ , one has $\mathcal{E}_v({J}_v) \leq \mathcal{E}_v({J}) = 0$ and analogously $\mathcal{E}_v({J}_\sigma) \leq \mathcal{E}_\sigma({J}) = 0$ . Hence, minimisers of $\mathcal{E}_v$ and $\mathcal{E}_\sigma$ over $X_{R,M}$ reproduce exactly the trajectories of the system thanks to Theorem 3.3.

• • Assume there is no true regular J that generated the observed data. This means that, no matter how large M is taken in Theorem 3.6, the minimal limit energies $\mathcal{E}_\sigma({J}_\sigma)$ and $\mathcal{E}_v({J}_v)$ may not be equal to 0. However, the bound on the trajectories generated by the inferred model, Theorem 3.3, still holds and so the remaining residual values $\mathcal{E}_\sigma({J}_\sigma)$ and $\mathcal{E}_v({J}_v)$ are then an indication for how well our parametrised family of interaction models was able to capture the structure in the observed data.

Remark 3.8 (Non-uniqueness of minimisers). We point out that the parametrisation of an entropic game model by a payoff function J is not unique. Replacing J by $\kern1pt\hat{{\kern-1pt}J}\left(x,u,x^{\prime},u^{\prime}\right) \;:=\; J\!\left(x,u,x^{\prime},u^{\prime}\right) + g\left(x,x^{\prime},u^{\prime}\right)$ yields the same agent trajectories for arbitrary $g \,:\, \Omega \times \Omega \times U$ . This can be seen from the fact that $f^{J+g}=f^J$ for $f^J$ and $f^{J+g}$ as given in (2.3). Intuitively, given some x, x^′ and u^′, which strategy u is most attractive does not change, if we add a constant benefit for all potential choices. The analogous observation holds for the undisclosed case. This implies that minimisers of (3.8), (3.9) and (3.10) are not unique. This has to be taken into account when we try to interpret an inferred J. In our numerical examples with the undisclosed fast-reaction model, we usually made additional structural assumptions on J (see Section 4.1) and as a consequence did not observe issues with non-unique minimisers. The phenomenon is encountered in a proof-of-concept example on the differential inference functional, Section 4.3.

4.1 Numerical setup

Preliminaries and discretisation. We assume for simplicity that the strategy space U is finite, that is, $U=\{u_1,\ldots,u_K\}$ for some $K>0$ and fix the reference measure $\eta = \sum_{k=1}^K \frac{1}{K}\delta_{u_k}$ as the normalised uniform measure over U.

For a system of N agents starting at positions $x_1(0), \dots, x_N(0)$ , we observe S snapshots of the evolution at times $t_s = s\cdot \frac{T}{S}$ for $s = 1,\dots,S$ , so that we are given locations $\{x_i(t_s)\}_{s,i}$ , velocities $\{\partial_t x_i(t_s)\}_{s,i}$ and (sometimes) mixed strategies $\{\sigma_i(t_s)\}_{s,i}$ where indices s and i run from 1 to S and 1 to N respectively. Here, $\sigma_{i}(t_s) = (\sigma_{i,1}(t_s),\dots,\sigma_{i,K}(t_s))$ is a discrete probability density with respect to $\eta$ .

Let $\Omega = \prod_{i=1}^d [\ell_i, r_i]$ , $\ell_i,r_i \in \mathbb{R}$ , be the smallest hypercube containing all observed locations. We discretise $\Omega$ by a regular Cartesian grid and describe the unknown payoff J by its values $J(\cdot,u_k,\cdot)$ at grid points for $u_k \in U$ (or $J(\cdot,u_k,\cdot,u_{k^{\prime}})$ in the fully general setting). Between grid points, J is extended by coordinate-wise linear interpolation, that is, we consider piecewise d-linear finite elements over hypercubes.

Within this setting, the inference functional (3.10) reduces to the discrete error functional

(4.1) \begin{align}\mathcal{E}_\sigma^N(J) = \frac{1}{SNK} \sum_{s=0}^S \sum_{i=1}^N \sum_{k=1}^K\log\left( \frac{\sigma_{i,k}(t_s)}{\sigma^{J}_{i,k}(t_s)} \right)\sigma_{i,k}(t_s),\end{align}

where for $s \in \{0,\ldots,S\}$ , $i \in \{1,\dots,N\}$ and $k \in \{1,\dots,K\}$ , the optimal mixed strategy at time $t_s$ is given by (2.15) as

\begin{equation*}\sigma^{J}_{i,k}(t_s) = \sigma^{J}_{i}(x_1(t_s),\dots,x_N(t_s))(u_k) = \frac{\exp\left(\dfrac{1}{\varepsilon N}\sum_{j=1}^N J(x_i(t_s),u_k,x_j(t_s))\right)}{\frac{1}{K}\sum_{k^{\prime}=1}^K \exp\left(\dfrac{1}{\varepsilon N}\sum_{j=1}^N J(x_i(t_s),u_{k^{\prime}},x_j(t_s))\right)}.\end{equation*}

Similarly, the inference functional (3.9) reduces to the discrete functional

(4.2) \begin{align}\mathcal{E}_v^N(J) = \frac{1}{SN} \sum_{s=0}^S \sum_{i=1}^N \left\| v_i^J(t_s) - \partial_t x_i(t_s)\right\|^2,\end{align}

where

\begin{equation*}v_i^J(t_s) = v^{J}_{i}\left(x_1(t_s),\dots,x_N(t_s)\right) = \frac{1}{K}\sum_{k=1}^K e\left(x_i(t_s), u_k\right) \sigma_{i,k}^J(t_s).\end{equation*}

Reducing dimensionality. When $d \geq 2$ , describing a general function $J \,:\, \Omega \times U \times \Omega \to \mathbb{R}$ where $\Omega \subset \mathbb{R}^d$ (or even $\Omega \times U \times \Omega \times U \to \mathbb{R}$ in the fully general setting) requires many degrees of freedom (one value for each node on the grid) and inferring them from data in a reliable way would require a large set of observations. After inference, interpreting such a general function to understand the interactions between agents can be a daunting task. For these two reasons, we may wish to incorporate some prior knowledge on the structure of J via a suitable ansatz. A prototypical example is

(4.3) \begin{align}J\!\left(x,u,x^{\prime}\right)=J_1(x,u) + J_2(x^{\prime}-x,u),\end{align}

where the first term models how a single agent prefers to move by itself depending on its absolute location, and the second term models a translation invariant pairwise interaction with other agents. This reduces the dimensionality from one function over $\mathbb{R}^d \times U \times \mathbb{R}^d$ to two functions over $\mathbb{R}^d \times U$ . Further simplifications can be made, for instance, by assuming that the interaction is also rotation invariant. Such parametrisations are very common in Newtonian models, see for instance [Reference Bongini, Fornasier, Hansen and Maggioni10], where radially symmetric interaction functions are assumed. Unless noted otherwise, we will in the following use the ansatz (4.3) in our experiments.

Regularisation and dropping the Lipschitz constraints. Even after the reduction (4.3), inferring the coefficients of J is a high-dimensional nonlinear inverse problem. We may observe no (or very few) agents near some grid points and thus have no reliable data to set the corresponding coefficients of J directly. A common approach to avoiding ill-posedness is to add regularisation to the minimisation problem. Intuitively, this will diffuse information about observations over the grid to some extent. For regularisation in $\Omega,$ we use the squared $L^2$ -norm of the gradient of J, which boils down to weighted sums of the finite differences between grid points. For the strategy space U, one can in principle do the same, although this was not necessary in our examples. When J is split as in (4.3), a regulariser can be applied separately to each term. The regularised version of (4.1) then becomes

(4.4) \begin{align}\mathcal{E}_\sigma^N(J) + \lambda_1 \cdot \mathcal{R}_1(J_1) + \lambda_2 \cdot \mathcal{R}_2(J_2).\end{align}

Of course, over-regularisation will lead to loss of contrast and high-frequency features of J. The effect of under-regularisation is illustrated in Example 4.5.

While in Theorem 3.6 we consider increasingly finer discretisations $X_{R,M}^N$ , $N \to \infty$ , of the payoff function space $X_{R,M}$ , in our numerical examples a fixed resolution is entirely sufficient (approximately 30 grid points along any spatial dimension). On a fixed grid, the Lipschitz and $L^2$ -gradient semi-norms are equivalent. Consequently, for simplicity, in our numerical experiments we may drop the explicit constraint on the Lipschitz constant of J and only use the regularisation term to impose regularity. This has the advantage that our finite-dimensional discrete problem is then unconstrained, and we may approach it with quasi-Newton methods such as L-BFGS, for which we use the Julia package Optim [Reference Mogensen and Riseth52].

Inferring ${\boldsymbol\sigma}_{\kern-1pt\textbf{i}}(\textbf{t}_{\textbf{s}})$ from observations. In many situations, we are not able to observe the mixed strategies of agents, but merely their locations and velocities. We can then decide to use the nonconvex velocity-based error functional (4.2). Alternatively, in case we want to optimise the convex inference functional $\mathcal{E}_\sigma^N$ (4.1) instead, we can try to infer some meaningful $\sigma_i(t_s)$ from the observed velocities. In both cases, we need to make suitable choices for U and e. Unless stated otherwise, we pick U as a discrete subset of $\mathbb{R}^d$ such that observed velocities lie in the (interior of the) convex hull of U and pick $e(x,u) \;:=\; u$ . This means that every observed velocity can be reproduced by some mixed strategy. As a heuristic method for inferring a viable mixed strategy $\sigma$ from a velocity v, we propose to solve the following minimisation problem:

(4.5) \begin{align}\arg \min\left\{\int_U \left[ \varepsilon \cdot \log(\sigma(u)) + \|u-v\|^2 \right]\,\sigma(u)\,\textrm{d} \eta(u)\middle|\sigma \in \mathcal{S}_{a,b}, \,\int_U u \cdot \sigma(u)\,\textrm{d} \eta(u)=v\right\}\end{align}

The constraint enforces that $\sigma$ reproduces the observed velocity v. The term $\|u-v\|^2\,\sigma(u)$ encourages $\sigma$ to be concentrated in the vicinity of v in U, whereas the entropic term $\varepsilon \cdot \log(\sigma(u))\,\sigma(u)$ keeps $\sigma$ from collapsing onto a single point. It is easy to see that minimisers of (4.5) are of the form $\sigma(u)=A \cdot \exp(-\|u-\tilde{v}\|^2/\varepsilon)$ where A is a normalisation factor and $\tilde{v}$ is a velocity close to v, potentially slightly shifted to make sure that the constraint $\int_U u \cdot \sigma(u)\,\textrm{d} \eta(u)=v$ is satisfied. (4.5) is readily solved by standard methods. We consider this approach a heuristic ‘inverse’ of the approximation of Newtonian models by game models, as outlined in Example 2.7.

Observations from forward models: Subsampling and multiple realisations. In all the examples here below, observations are obtained by means of simulating a forward model with explicit Euler-stepping with time-step size $\Delta t=0.02$ . Since the numerical complexity of our inference functionals grows linearly with the number of observations, we usually choose to sub-sample the trajectories in time, keeping only every $\delta$ -th step, since intuitively the information contained in subsequent steps is largely redundant as agents are confronted with a very similar configuration. Unless noted otherwise, we set $\delta=2$ .

The quality of the reconstructed model is highly related to the data quality: more observations in a given region imply a better reconstruction of the payoff function in that specific region. Generally, if we consider as inputs only data coming from a single realisation of the system we have highly correlated observations, which tend to concentrate only on a small subset of the reconstruction domain. Vaster regions of the domain are then explored taking $N \to \infty$ . However, another option to enhance the quality of the reconstructed model is to combine data from multiple realisations with different initial conditions, that is, we fix the number of agents N and simulate the system for multiple different initial conditions. This provides a good coverage of the domain. Extending the functionals in Section 3 to multiple observations is straightforward.

4.2 Learning from undisclosed fast-reaction game models

We now begin with relatively simple numerical examples: an explicit game model is fixed (including U, e, J), simulations according to this model are performed and subsequently the payoff function is inferred from the observations. By possibly re-scaling J, we may w.l.o.g. set $\varepsilon=1$ for all numerical experiments, see (2.15). In this section we focus on the undisclosed fast-reaction setting. A proof of concept of inference for the full entropic game models via the differential functional of Section 3.2 is discussed in Section 4.3.

Example 4.1 (Game model, $d=1$ ). Let $d=1$ , $U=\{-1,+1\}$ and $e(x,u)=u$ . We set

(4.6) \begin{align} J\!\left(x,u,x^{\prime}\right) & \;:=\; J_1(x,u) + J_2(x^{\prime}-x,u), \qquad J_1(x,u) \;:=\; -u \cdot x, \nonumber \\ J_2(\Delta x,u) & \;:=\; -u \cdot \tanh(5\Delta x) \cdot \left(\max\left\{1-|\Delta x|^2,0\right\}\right)^2. \end{align}

$J_1$ encourages agents to move towards the origin $x=0$ ; for example, for $x>0$ one has $J_1(x,-1)>J_1(x,+1)$ , thus preferring strategy (and velocity) $u=-1$ over $u=+1$ . $J_2$ implements repulsion between nearby agents: when $\Delta x \in (0,1)$ , that is, the ‘other agent’ is to the right of ‘me’ and relatively close, then $J_2(\Delta x,-1)>J_2(\Delta x,+1)$ , and I am encouraged to move left.

We set $N=8$ , simulated 100 instances over a time span of $0.2$ , collecting 5 data points per instance, that is, a total of 500 observed configurations (consisting of locations and mixed strategies). In each instance, initial locations are sampled uniformly from $[-1,1]$ . Observed relative locations between two agents are thus distributed over $[-2,2]$ . Describing the discrete J required 178 coefficients. Since observed mixed strategies are available in this example, we use the energy (4.1) for inference augmented with a regulariser as discussed in Section 4.1, with $\lambda_{1}=\lambda_2=10^{-6}$ , see (4.4). The results are illustrated in Figure 4. We find that, as intended, the model describes agents moving towards the origin while avoiding getting too close to each other. The functions $J_1$ and $J_2$ are accurately recovered from the data. For newly generated initial conditions, that is, not taken from the training data, the trajectories simulated with the inferred J are in excellent agreement with the underlying true model, demonstrating that the inferred model is able to generalise.

Example 4.2 (Towards the mean-field regime). We can also explore the mean-field regime numerically. Figure 5 shows trajectories for the same model as in the previous example for an increasing number of particles where the initial locations are sampled from some underlying distribution. As anticipated, the behaviour approaches a consistent limit.

Inference can also be performed for larger numbers of particles. We applied the same approach as in the previous example with $N=100$ agents. Since we are dealing with a large number of agents, already a small number of configurations carries enough information for the inference. Thus, we simulate 10 instances over a time span of $0.02$ with time-step size $\Delta t = 0.002$ , collecting 2 data points per instance; that is, a total of 20 observed configurations (for locations and mixed strategies). In each instance, initial locations are sampled uniformly from $[-1,1]$ . The inferred payoff function is essentially identical to the one obtained in Figure 4.

Figure 4. Inference from 1d observations of an undisclosed fast-reaction system. (a,b): reconstruction of the payoff function. (c): distribution of the data in the training set. (d): comparison between true trajectories (solid lines) and trajectories generated by the inferred model (dashed lines) for two new realisations (not part of the training data).

Figure 5. Numerical approximation of the mean-field limit. The forward model of Example 4.1 is simulated for an increasing number of agents. Each panel shows a discrete histogram of the particle distribution over time. As N increases, the histograms approach a consistent limit.

Example 4.3 (Game model, $d=1$ , noisy data). We repeat the previous example with some added noise. We consider the same set of observations, but the observed mixed strategies are corrupted by re-sampling: For a given ‘true’ simulated $\sigma$ , we draw 20 times from U, according to $\sigma$ and now use the resulting empirical distribution for inference. These new corrupted $\sigma$ provide corrupted velocities predictions: in our example, the overall standard deviation between exact and corrupted velocities is $\approx 0.2$ (with the range of admissible velocities being $[-1,1]$ ). This could model the situation where observation of the mixed strategies is imperfect or by observing very noisy (stochastic) agent movement and using empirical distributions of agent velocities over short time intervals as substitute for mixed strategies. The results are shown in Figure 6, for parameters $\lambda_{1}=\lambda_2=10^{-5}$ to enforce more regularity due to noisy inputs. The inferred payoff functions are similar to the previous example with just a few spurious fluctuations. The trajectories simulated with the inferred J are close to the (unperturbed) true trajectories, the error is consistently smaller than the one caused by the re-sampling noise, indicating that inference also works well in noisy settings.

Figure 6. Inference from noisy 1d observations of an undisclosed fast-reaction system. (a): reconstruction of the payoff function for $u=-1$ . (b): comparison between true trajectories (solid lines) and trajectories generated by the inferred model (dashed lines) for multiple new realisations (not part of the training data).

Figure 7. Inference from 2d observations of an undisclosed fast-reaction system. (a): inferred payoff function and ground truth for pure strategies (1,1) and $(-1,1)$ . (b): distribution of the data in the training set (agents locations and relative locations), dark blue : low, yellow: high. (c): comparison between exact trajectories (solid lines) and trajectories generated by the inferred model (dashed lines) for two new realisations (not part of the training data).

Example 4.5 (Under-regularisation). As is common in inverse problems, if the regularisation parameters $\lambda_1, \lambda_2$ in (4.4) are chosen too high, one obtains an over-regularised result, where the large penalty for spatial variation leads to a loss of contrast and high-frequency features of J. Conversely, we may experience under-regularisation when the parameters are small and the observed data does not cover the learning domain sufficiently, for example when we only observe few realisations of the system. Such a case is illustrated in Figure 8(c) where the observations are concentrated on a few one-dimensional trajectories. For weak regularisation, the reconstructed J then tends to be concentrated on these trajectories as well, see Figure 8(a). And while the inferred J may predict the agent behaviour with high accuracy on these trajectories, it will generalise very poorly to different data. With stronger regularisation, J is extended beyond the trajectories in a meaningful way, see Figure 8(b).

Figure 8. Influence of the regularisation parameter. Analogous to Figure 7 but with fewer observations, as displayed in (c). (a): reconstruction for low regularisation parameters, $\lambda_1=\lambda_2=10^{-15}$ . (a): reconstruction for higher regularisation parameters, $\lambda_1=\lambda_2=10^{-5}$ . See Example 4.5 for more details.

4.3 An example for differential inference

The ‘differential inference’ functional $\mathcal{E}_{\dot{\sigma}}^{J,N}$ introduced in Section 3.2 was constructed in close analogy to (1.7). We have argued that it may be challenging to apply to real data, as it would require the observation of mixed strategies $\{\sigma_{i}(t_s)\}_{s,i}$ and their temporal derivatives $\{\partial_t \sigma_{i}(t_s)\}_{s,i}$ . As a proof of concept, we now demonstrate that if these observations are indeed available, the numerical setup described up to now can be easily adapted to $\mathcal{E}_{\dot{\sigma}}^{J,N}$ . After discretisation, $\mathcal{E}_{\dot{\sigma}}^{J,N}$ reduces to a finite-dimensional quadratic objective. For a quadratic regulariser, such as the squared $L^2$ -norm of the gradient of J, minimisation then corresponds to solving a linear system. Alternatively, we can implement the Lipschitz constraint on J which amounts to adding linear constraints to the quadratic problem. We solve the resulting finite-dimensional optimisation problem using the operator splitting approach provided by the OSQP solver [Reference Stellato, Banjac, Goulart, Bemporad and Boyd69].

As an example, we simulate the entropic regularised model in (2.2) for $U = \{-1,1\}$ and $J \colon \mathbb{R} \times \{-1,1\} \times \mathbb{R} \times \{-1,1\}$ defined by

\begin{align*}J\!\left(x,u,x^{\prime},u^{\prime}\right) &= -\frac{1}{2}\left(u+u^{\prime}\right)\left((u+1)x^5+(u-1)(x+x^{\prime})^3\right)\end{align*}

Our choice of J is not motivated by providing a meaningful model for anything, but simply to illustrate that inference works in principle. We collect data using the same sub-sampling approach presented above: 8 agents, 100 realisations for a time span of 0.2, initial conditions uniformly sampled from $[-1,1]$ , only every other observation is kept for a total of 500 input configurations. Observed positions (x, x ^′) with $x \neq x^{\prime}$ are displayed in Figure 9(a). No particular structure is assumed for J.

Figure 9. Inference from 1d observations generated with (2.2). (a): distribution of observed pairs (x, x ^′), $x \neq x^{\prime}$ . (b): comparison between exact trajectories (solid lines) and inferred model (dashed lines) on new realisations (not part of the training data).

As discussed in Remark 3.8, the optimal inferred J is not unique. As a remedy, we normalise the reconstructed $J^r$ such that $J^r(x,1,x^{\prime},-1) = J^r(x,-1,x^{\prime},1) = 0$ . Comparisons for trajectories are reported in Figure 9(b) for new randomly selected initial conditions with $T = 8$ : trajectories with the reconstructed $J^r$ follow closely the trajectories generated with the true J.

4.4 Learning for Newtonian models

In Example 2.7, we have shown that entropic game models with fast reaction can also capture Newtonian interaction models. In this section, we perform inference with an entropic game model ansatz on data generated by a Newtonian model. As an example, we pick

(4.7) \begin{align}\partial_t x_i(t) = \frac{1}{N} \sum_{j=1}^N f\!\left(x_i(t),x_j(t)\right)\qquad \text{with} \qquad f(x,x^{\prime}) = - x - \frac{\tanh\left(5\left(x^{\prime}-x\right)\right)}{\left(1+\|x^{\prime}-x\|\right)^2},\end{align}

where the function $\tanh$ is assumed to be applied componentwise. The first term in f encourages agents to move towards the origin, the second term implements a pairwise repulsion at close range. Trajectory data are collected for a two dimensional system: for $N=8$ agents, we simulate 100 instances of (4.7) over a time span of $0.2$ , collecting 5 data points per instance to total 500 observed configurations. In each instance, initial locations are sampled uniformly from $[-0.75,0.75]^2$ . Observed relative locations between two agents are thus distributed over $[-1.5,1.5]^2$ , as displayed in Figure 10(b).

Figure 10. Inference from 2d observations of a Newtonian system. Top rows: inferred payoff function for pure strategies (1,1) and $(-1,1)$ , bottom line: distribution of the data in the training set (agents locations and pairwise distances) and comparison between exact Newtonian trajectories (solid lines) and trajectories generated by the inferred game model (dashed lines).

We use the parametrisation (4.3) for the payoff function. Unlike in the previous examples, the original model does not involve a strategy set U or a map e. Thus, prior to inference, we have to make informed choices for these. As in earlier examples, we choose $U \subset \mathbb{R}^2$ and $e(x,u)=u$ (i.e. pure strategies correspond directly to velocities) and in particular we pick $U = \{u_1,\dots,u_K\}$ such that the convex hull of U in $\mathbb{R}^2$ contains (almost) all observed velocities, meaning that they can be reproduced by suitable mixed strategies. So for (almost) every i and s, we can write $\partial_t x_i(t_s) = \sum_k u_k\sigma_{k}$ for some mixed strategy $\sigma$ over U.

Example 4.7. In most cases, there will not be a unique preferred choice for U and various options could be reasonable. We now demonstrate robustness of the inference results between different choices that are essentially equivalent. Thus, we repeat Example 4.6 changing U into $U = \{(1,0),(0,1),(-1,0),(0,-1)\}$ . The results are illustrated in Figure 11. Essentially, a rotation of U caused a corresponding rotation in the reconstructed functions, without changing the underlying behaviour that is being described.

Figure 11. Inferred payoff function via $\mathcal{E}_v^N$ -based energy for pure strategies (1,0) and (0,1) (same Newtonian input data as in Figure 10(c)).

Example 4.8. In the previous two examples, we used the velocity-based functional $\mathcal{E}_v^N$ (4.2) for inference. As discussed in Section 4.1, we can also optimise the error functional $\mathcal{E}_\sigma^N$ (4.1) upon providing some suitable mixed strategies via (4.5). We repeat Example 4.6 with this approach. The results are illustrated in Figure 12: the reconstructed functions closely resemble the previous reconstruction provided in Figure 10(a).

Figure 12. Inferred payoff function via $\mathcal{E}_\sigma^N$ -based energy for pure strategies (1,1) and $(-1,1)$ (same Newtonian input data as in Figure 10(c)). To compare with $\mathcal{E}_v^N$ -based reconstruction in Figure 10(a).

4.5 Pedestrian dynamics

The examples presented in Sections 4.2 and 4.3 were based on observations generated with the underlying entropic game model with known U and e. The examples in Section 4.4 were based on observations generated by Newtonian models, which by virtue of Example 2.7 are known to be covered by the entropic game models, if suitable U and e are chosen. Here we consider a more challenging scenario, trying to learn interaction rules between pedestrians in a model adapted from [Reference Bailo, Carrillo and Degond7, Reference Degond, Appert-Rolland, Moussaïd, Pettré and Theraulaz32]. This model attempts to describe ‘rational’ agents (pedestrians) that avoid collisions with other agents not merely via a repulsion term involving the other agents’ current locations, but by explicit anticipation of their trajectories. Related datasets and models in the computer vision community can be found, for instance, in [Reference Pellegrini, Ess, Schindler and van Gool59, Reference Rahmati, Talebpour, Mittal and Fishelson62, Reference Yang, Li, Redmill and Özgüner76]. It will be interesting future work to see if our model can be trained on such data. Note that this differs from the setting where pedestrians in high-density and high-stress situations are modelled as a non-sentient fluid that is solely driven by external gradients and internal pressure, see for instance [Reference Bellomo and Dogbe9, 41] for reviews on the literature of micro- and macroscopic crowd motion models.

Forward model. We now describe the model for pedestrian movement that we use to generate observations. It is slightly adapted from [Reference Bailo, Carrillo and Degond7, Reference Degond, Appert-Rolland, Moussaïd, Pettré and Theraulaz32], to which we refer for more details and motivation of the quantities we re-introduce here. We consider N pedestrians in $\mathbb{R}^2$ , each of which with desired target velocity $\bar{v}_i \in \mathbb{R}^2$ , $i=1,\ldots,N$ . The dynamic is described by the coupled system of second-order equations

(4.8) \begin{align}\left\{\!\begin{aligned}\partial_t x_i(t) &= v_i(t), \\[4pt]\partial_t v_i(t) &= -\nabla_{v_i} \Phi_i(x_1(t),\ldots,x_N(t),v_1(t),\ldots,v_N(t)), \\[4pt]\|v_i(t)\| & = 1,\end{aligned}\right.\end{align}

where the third line means that we confine the gradient movement of the velocities in the second line to the unit sphere. This reduces the dimension of the strategy space and thus simplifies inference (see below). The potential $\Phi_i$ is designed to steer the velocity towards the target velocity $\bar{v}_i$ while at the same time trying to avoid as much as possible close encounters with other pedestrians. The potential is computed from pairwise interaction terms that involve anticipated (near-)collision events. A key assumption is that agents are able to evaluate at each time the current locations and velocities of the other agents in their cone of vision; that is, they fully know the relevant status of the system. Thus, for locations $x_1, \dots, x_N$ and velocities $v_1, \dots, v_N$ , each agent i assumes the others’ velocities to be fixed and can evaluate ‘optimality’ of a new velocity v through the following quantities:

(4.9) \begin{align} D_{i,j}(v) &= -\frac{\left(x_j-x_i\right)\cdot\left(v_j-v\right)}{\|v_j-v\|^2}\|v\|\qquad\qquad\qquad\qquad\qquad\qquad \end{align}

(4.10) \begin{align} \ \ \ \ \ C_{i,j}(v) &= \left( \|x_j-x_i\|^2 - \frac{(\left(x_j-x_i\right)\cdot\left(v_j-v\right))^2}{\|v_j-v\|^2}\right)^{1/2}\qquad\qquad\qquad\quad \end{align}

(4.11) \begin{align} S_{i,j}(v) &=\left\{\begin{aligned}1 \quad&\text{ if } \frac{\left(x_j-x_i\right)\cdot v}{\|x_j-x_i\|\|v\|} > \cos\left(\frac{7}{12}\pi\right) \text{ and } D_{i,j}(v) > 0 \\[5pt]0 \quad&\text{ else}.\end{aligned}\right.\end{align}

These are illustrated in Figure 13. $D_{i,j}$ is the travelling distance of i to the closest encounter with j, based on velocities v and $v_j$ . It is negative if the closest encounter would have happened in the past. $C_{i,j}$ is their distance realised at that closest encounter and $S_{i,j}$ counts whether agent j is visible to i and getting closer. Note that $D_{i,j}$ , $C_{i,j}$ and $S_j$ are functions of the other agent’s location and velocity. From the perspective of agent i, these are considered parameters, and thus, we drop them in the explicit notation of arguments. For a parameter $a > 0$ define the quadratic penalty

\begin{align*}\phi_a(x) =\begin{cases}\left(\dfrac{x}{a}-1\right)^2 & \text{if } x\leq a, \\[5pt]0 & \text{else,}\end{cases}\end{align*}

and then set the i-th potential to be

(4.12) \begin{align}\Phi_i(v) = k_1 \|v-\bar{v}_i\|^2 + \frac{k_2}{\sum_{j \neq i} S_{i,j}(v)}\sum_{j \neq i} S_{i,j}(v) \cdot \phi_R\left(C_{i,j}(v)\right) \cdot \phi_L\left(D_{i,j}(v)\right)\end{align}

for some parameters $k_1, k_2, L, R > 0$ . Agent i will then evolve its velocity based on the gradient of $\Phi_i$ , which balances the desire to stay close to the target velocity $\bar{v}_i$ and the avoidance of close encounters with other agents. The second term penalises anticipated minimal distances $C_{i,j}$ below R, that happen on a distance horizon of L, with a stronger penalty if the anticipated encounter is closer.

Figure 13. Adapted from [Reference Bailo, Carrillo and Degond7].

State space and strategy space. System (4.8) is simulated with potentials $\Phi_i$ defined as in (4.12). By the third line of (4.8) we can write $v_i(t) = (\cos(\theta_i(t)), \sin(\theta_i(t)))$ at each t for a suitable time depended angle $\theta_i$ . Observations consist of locations $\{x_i(t_s)\}_{s,i}$ and velocities $\{v_i(t_s)\}_{s,i}$ for $s \in \{0,\dots,S\}$ and $i \in \{1,\dots,N\}$ . In particular, observations on velocities can be recast into observations of the angles $\{\theta_i(t_s)\}_{s,i}$ , $\theta_i(t_s) \in S^1$ , so that $v_i(t_s) = (\cos(\theta_i(t_s)),\sin(\theta_i(t_s)))$ . Likewise, the desired velocity $\bar{v}_i$ can be encoded by an orientation $\bar{\theta}_i \in S^1$ .

We want to model the system as an undisclosed fast-reaction entropic game. Consequently, we must choose a suitable ‘location’ space, strategy space and map e for the agents. In the undisclosed setting, the choice of strategy of an agent does not depend on the choices of strategies of the other agents. However, for pedestrians, the movement choices of one agent clearly depend on the current velocities of other agents. This is explicit in the above model, and of course also true in reality, where we can learn about other persons’ current movement by their orientation (people walk forward). Consequently, the orientations $\theta_i$ cannot be interpreted as strategies and should be considered as part of the agents’ locations. Instead, the second line of (4.8) suggests that strategies describe the agents’ choice to change their orientations. Finally, different agents may have different desired orientations $\bar{\theta}_i$ ; that is, we may have different ‘species’ of agents. As discussed in Remark 2.4, this can be formally incorporated into our model by enriching the physical space.

Consequently, as physical space we use $\mathbb{R}^2 \times S^1 \times S^1$ where the first component describes the physical location of the agent, the second component their current orientation and the third component their desired orientation.

As strategy space, we pick $U \subset \{-C,+C\}$ for some $C>0$ , where a pure strategy $u \in U$ corresponds to the desire to change the orientation with rate u. C should be picked sufficiently large to capture the maximally observed changes of orientation of the agents.

In summary, in this particular instance of (2.15) we choose

\begin{align*}e\left((x,\theta,\bar{\theta}),u\right) =\begin{pmatrix}\cos(\theta),\sin(\theta),u,0\end{pmatrix},\end{align*}

and consequently obtain the full system

(4.13) \begin{align}\left\{\!\begin{aligned}\partial_t x_i(t) &= (\cos(\theta_i(t)),\sin(\theta_i(t))) \\[3pt]\partial_t \theta_i(t) &= \vartheta_i^J\big((x_1(t), \theta_1(t),\bar{\theta}_1(t)),\dots,(x_N(t), \theta_N(t), \bar{\theta}_N(t))\big)\\[3pt]&\;:=\; \frac1K \sum_{k=1}^K u_k \cdot \sigma_{i}^J\big((x_1(t), \theta_1(t),\bar{\theta}_1(t)),\dots,(x_N(t), \theta_N(t), \bar{\theta}_N(t))\big)(u_k) \\[3pt]\partial_t \bar{\theta}_i(t) &= 0\end{aligned}\right.\end{align}

Ansatz for payoff function. For the payoff function, in analogy to (4.3), we make the following ansatz:

(4.14) \begin{align}J\!\left(\left(x,\theta, \bar{\theta}\right)\!,u,\left(x^{\prime},\theta^{\prime},\bar{\theta}^{\prime}\right)\right) = J_1\left(\theta, u, \bar{\theta}\right) + J_2\left(R_{-\theta} \cdot \left(x^{\prime}-x\right), \theta^{\prime}-\theta, u\right),\end{align}

where $R_{-\theta}$ denotes the rotation matrix by angle $-\theta$ . $J_1$ is intended to recover the drive towards the desired orientation $\bar{\theta}_i$ and $J_2$ the pairwise interaction between agents. We have added the natural assumption that $J_2$ is invariant under translations and rotations. For each agent, we may choose a coordinate system such that they are currently in the origin ( $x=0$ ), heading right ( $\theta=0$ ), whereas the other agent is located at $y=R_{-\theta} \left(x^{\prime}-x\right)$ with orientation $\Delta \theta=\theta^{\prime}-\theta$ . Further dimensionality reduction can be obtained by setting $J_2$ to zero, when the other agent is not within the agent’s field of vision. Alternatively, we can attempt to infer the field of vision from the observations during learning.

Discrete inference functional. We approximate the missing velocities $\{\partial_t \theta_i(t_s)\}_{s,i}$ by finite differences,

\begin{equation*}\partial_t \theta_i(t_s) \approx \frac{\theta_i(t_s)-\theta_i(t_{s-1})}{t_s-t_{s-1}} \quad \text{ for all }s=1,\dots,S, \text{ and } i = 1,\dots,N,\end{equation*}

and the corresponding velocity-based inference functional (3.9) is given by

(4.15) \begin{align}\mathcal{E}_v^N(J) = \frac{1}{SN} \sum_{s=1}^S \sum_{i=1}^N \left| \vartheta_i^J(t_s) - \partial_t \theta_i(t_s)\right|^2,\end{align}

where $\vartheta_i^J$ is introduced in (4.13).

Figure 14. Inference from 2d observations of a pedestrian model. (a): distribution of the data in the training set. (b,c): reconstruction of the payoff function (section for $J_2$ at $\Delta \theta = \pi$ ). (d): comparison between exact trajectories (solid lines) and trajectories generated by the inferred model (dashed lines) for four new realisations (not part of the training data).