LEARNING GRAPHON MEAN FIELD GAMES
Year: 2022

Abstract:
Recent advances at the intersection of dense large graph limits and mean field games have begun to enable the scalable analysis of a broad class of dynamical sequential games with large numbers of agents. So far, results have been largely limited to graphon mean field systems with continuous-time diffusive or jump dynamics, typically without control and with little focus on computational methods. We propose a novel discrete-time formulation for graphon mean field games as the limit of non-linear dense graph Markov games with weak interaction. On the theoretical side, we give extensive and rigorous existence and approximation properties of the graphon mean field solution in sufficiently large systems. On the practical side, we provide general learning schemes for graphon mean field equilibria by either introducing agent equivalence classes or reformulating the graphon mean field system as a classical mean field system. By repeatedly finding a regularized optimal control solution and its generated mean field, we successfully obtain plausible approximate Nash equilibria in otherwise infeasible large dense graph games with many agents. Empirically, we are able to demonstrate on a number of examples that the finite-agent behavior comes increasingly close to the mean field behavior for our computed equilibria as the graph or system size grows, verifying our theory. More generally, we successfully apply policy gradient reinforcement learning in conjunction with sequential Monte Carlo methods.

heading: 1 INTRODUCTION
text: Today, reinforcement learning (RL) finds application in various application areas such as robotics (Kober et al., 2013), autonomous driving (Kiran et al., 2021) or navigation of stratospheric balloons (Bellemare et al., 2020) as a method to realize effective sequential decision-making in complex problems. RL remains a very active research area, and there remain many challenges in multi-agent reinforcement learning (MARL) as a generalization of RL such as learning goals, non-stationarity and scalability of algorithms (Zhang et al., 2021). Nonetheless, potential applications for MARL are manifold and include e.g. teams of unmanned aerial vehicles (Tožička et al., 2018; Pham et al., 2018) or video games (Berner et al., 2019; Vinyals et al., 2017). While the domain of MARL is somewhat empirically successful, problems quickly become intractable as the number of agents grows, and methods in MARL typically miss a theoretical foundation. A recent tractable approach to handling the scalability problem in MARL are competitive mean field games and cooperative mean field control (Gu et al., 2021). Instead of considering generic multi agent Markov games, one considers many agents under the weak interaction principle, i.e. each agent alone has a negligible influence on all other agents. This class of models naturally contains a large number of real world scenarios and can find application e.g. in analysis of power network resilience (Bagagiolo & Bauso, 2014), smart heating (Kizilkale & Malhame, 2014), edge computing (Banez et al., 2019) or flocking (Perrin et al., 2021b). See also Djehiche et al. (2017) for a review of other engineering applications.
Mean field games. Mean field games (MFGs) were first popularized in the independent seminal works of Huang et al. (2006) and Lasry & Lions (2007) for the setting of differential games with diffusion-type dynamics given by stochastic differential equations. See also Guéant et al. (2011); Bensoussan et al. (2013) for a review. Since then, extensions have been manifold and include e.g. discrete-time (Saldi et al., 2018), partial observability (Saldi et al., 2019), major-minor formulations
(Nourian & Caines, 2013) and many more. In the learning community, there has been immense recent interest in finding and analyzing solution methods for mean field equilibria (Cardaliaguet & Hadikhanloo, 2017; Mguni et al., 2018; Guo et al., 2019; Subramanian & Mahajan, 2019; Elie et al., 2020; Cui & Koeppl, 2021; Pasztor et al., 2021; Perolat et al., 2021; Perrin et al., 2021a), solving the inverse reinforcement learning problem (Yang et al., 2018a) or applying related approximations directly to MARL (Yang et al., 2018b). In contrast to our work, Yang et al. (2018a) consider states instead of agents on a graph, while Yang et al. (2018b) requires restrictive assumptions and only considers averages instead of distributions of the neighbors. Recently, focus increased also on the cooperative case of mean field control (Carmona et al., 2019b; Mondal et al., 2021), for which dynamic programming holds on an enlarged state space, resulting in a high-dimensional Markov decision process (Pham & Wei, 2018; Motte & Pham, 2019; Gu et al., 2020; Cui et al., 2021).
Mean field systems on graphs. For mean field systems on dense graphs, prior work mostly considers mean field systems without control (Vizuete et al., 2020) or time-dynamics, i.e. the static case (Parise & Ozdaglar, 2019; Carmona et al., 2019a). To the best of our knowledge, Gao & Caines (2017) and Caines & Huang (2019) are the first to consider general continuous-time diffusion-type graphon mean field systems with control, the latter proposing many clusters of agents as well as proving an approximate Nash property as the number of clusters and agents grows. There have since been efforts to control cooperative graphon mean field systems with diffusive linear dynamics using spectral methods (Gao & Caines, 2019a;b). On the other hand, Bayraktar et al. (2020); Bet et al. (2020) consider large non-clustered systems in a continuous-time diffusion-type setting without control, while Aurell et al. (2021b) and Aurell et al. (2021a) consider continuous-time linear-quadratic systems and continuous-time jump processes respectively. To the best of our knowledge, only Vasal et al. (2021) have considered solving and formulating a graphon mean field game in discrete time, though requiring analytic computation of an infinite-dimensional value function defined over all mean fields and thus being inapplicable to arbitrary problems in a black-box, learning manner. In contrast, we give a general learning scheme and also provide extensive theoretical analysis of our algorithms and (slightly different) model. Finally, for sparse graphs there exist preliminary results (Gkogkas & Kuehn, 2020; Lacker & Soret, 2020), though the setting remains to be developed.
Our contribution. In this work, we propose a dense graph limit extension of MFGs in discrete time, combining graphon mean field systems with mean field games. More specifically, we consider limits of many-agent systems with discrete-time graph-based dynamics and weak neighbor interactions. In contrast to prior works, we consider one of the first general discrete-time formulations as well as its controlled case, which is a natural setting for many problems that are inherently discrete in time or to be controlled digitally at discrete decision times. Our contribution can be summarized as: (i) formulating one of the first general discrete-time graphon MFG frameworks for approximating otherwise intractable large dense graph games; (ii) providing an extensive theoretical analysis of existence and approximation properties in such systems; (iii) providing general learning schemes for finding graphon mean field equilibria, and (iv) empirically evaluating our proposed approach with verification of theoretical results in the finite N -agent graph system, finding plausible approximate Nash equilibria for otherwise infeasible large dense graph games with many agents.

heading: 2 DENSE GRAPH MEAN FIELD GAMES
text: In the following, we will give a dense graph N -agent model as well as its corresponding mean field system, where agents are affected only by the overall state distribution of all neighbors, as visualized in Figure 1. As a result of the law of large numbers, this distribution will become deterministic – the mean field – as N →∞. We begin with graph-theoretical preliminaries, see also Lovász (2012) for a review. The study of dense large graph limits deals with the limiting representation of adjacency matrices called graphons. Define I := [0, 1] and W0 as the space of all bounded, symmetric and measurable functions (graphons) W ∈ W0, W : I × I → R bounded by 0 ≤ W ≤ 1. For any simple graph G = ({1, . . . , N}, E), we define its step-graphon a.e. uniquely by
WG(x, y) = ∑
i,j∈{1,...,N}
1(i,j)∈E · 1x∈( i−1N , iN ] · 1y∈( j−1N , jN ], (1)
see e.g. Figure 1. We equipW0 with the cut (semi-)norm ‖·‖ and cut (pseudo-)metric δ
‖W‖ := sup S,T
∣∣∣∣ ∫
S×T W (x, y) dxdy
∣∣∣∣ , δ (W,W ′) := infϕ ‖W −W ′ ϕ‖ , (2)
for graphons W,W ′ ∈ W0 and W ′ϕ(x, y) := W ′(ϕ(x), ϕ(y)), where the supremum is over all measurable subsets S, T ⊆ I and the infimum is over measure-preserving bijections ϕ : I → I. To provide motivation, note that convergence in δ is equivalent to e.g. convergence of probabilities of locally encountering any fixed subgraph by randomly sampling a subset of nodes. Many such properties of graph sequences (GN )N∈N converging to some graphon W ∈ W0 can then be described by W , and we point to Lovász (2012) for details. In this work, we will primarily use the analytical fact that for converging graphon sequences ‖WGN −W‖ → 0, we equivalently have
‖WGN −W‖L∞→L1 = sup ‖g‖∞≤1
∫
I
∣∣∣∣ ∫
I (WGN (α, β)−W (α, β))g(β) dβ
∣∣∣∣ dα→ 0 (3)
under the operator norm of operators L∞ → L1, see e.g. Lovász (2012), Lemma 8.11. By Lovász (2012), Theorem 11.59, the above is equivalent to convergence in the cut metric δ (WGN ,W )→ 0 up to relabeling. In the following, we will therefore assume sequences of simple graphs GN = (VN , EN ) with vertices VN = {1, . . . , N}, edge sets EN , edge indicator variables ξNi,j := 1(i,j)∈EN for all nodes i, j ∈ VN , and associated step graphons WN converging in cut norm. Assumption 1. The sequence of step-graphons (WN )N∈N converges in cut norm ‖·‖ or equivalently in operator norm ‖·‖L∞→L1 as N →∞ to some graphon W ∈ W0, i.e.
‖WN −W‖ → 0, ‖WN −W‖L∞→L1 → 0 . (4)
Next, we define W -random graphs to consist of vertices VN := {1, . . . , N} with adjacency matrices ξN generated by sampling graphon indices αi uniformly from I and edges ξNi,j ∼ Bernoulli(W (αi, αj)) for all vertices i, j ∈ VN . For experiments, by Lovász (2012), Lemma 10.16, we can thereby generate a.s. converging graph sequences by sampling W -random graphs for any fixed graphon W ∈ W0. In principle, one could also consider arbitrary graph generating processes whenever a valid relabeling function ϕ is known.
In our work, the usage of graphons enables us to find mean field systems on dense graphs and to extend the expressiveness of classical MFGs. As examples, we will use the limiting graphons of uniform attachment, ranked attachment and p-Erdős–Rényi (ER) random graphs given by Wunif(x, y) = 1−max(x, y), Wrank(x, y) = 1− xy and Wer(x, y) = p respectively (Borgs et al., 2011; Lovász, 2012), each of which exhibits different node connectivities as shown in Figure 2.
Finite agent graph game. For simplicity of analysis, we consider finite state and action spaces X ,U as well as times T := {0, 1, . . . , T − 1}. On a metric space A, define the spaces of all Borel probability measures P(A) and all Borel measures B1(A) bounded by 1, equipped with the L1 norm. For simplified notation, we denote both a measure ν and its probability mass function by ν(·). Define the space of policies Π := P(U)T ×X , i.e. agents apply Markovian feedback policies πi = (πit)t∈T ∈ Π that act on local state information. This allows for the definition of weakly interacting agent state and action random variables
Xi0 ∼ µ0, U it ∼ πit(· | Xit), Xit+1 ∼ P (· | Xit , U it ,Git), ∀t ∈ T ,∀i ∈ VN (5) under some transition kernel P : X ×U×B1(X )→ P(X ), where the empirical neighborhood mean field Git of agent i is defined as the B1(X )-valued (unnormalized) neighborhood state distribution
Git := 1
N
∑
j∈VN
ξNi,jδXjt , (6)
where δ is the Dirac measure, i.e. each agent affects each other at most negligibly with factor 1/N . Finally, for each agent i we define separate, competitive objectives
JNi (π 1, . . . , πN ) := E
[ T−1∑
t=0
r(Xit , U i t ,Git)
] (7)
to be maximized over πi, where r : X × U × B1(X )→ R is an arbitrary reward function. Remark 1. We can also consider infinite horizons, J̃Ni (π1, . . . , πN ) ≡ E [∑∞ t=0 γ tr(Xit , U i t ,Git)
] with all results but Theorem 4 holding. One may also extend to state-action distributions, heterogeneous starting conditions and time-dependent r, P , though we avoid this for expositional simplicity. Remark 2. Note that it is straightforward to extend to heterogeneous agents by modelling agent types as part of the agent state, see also e.g. Mondal et al. (2021). It is only required to model agent states in a unified manner, which does not imply that there can be no heterogeneity.
With this, we can give a typical notion of Nash equilibria as found e.g. in Saldi et al. (2018). However, under graph convergence in Assumption 1, it is always possible for a finite number of nodes to have an arbitrary neighborhood differing from the graphon as N → ∞. Thus, it is impossible to show approximate optimality for all nodes and only possible to show for an increasingly large fraction 1 − p → 1 of nodes. For this reason, we slightly weaken the notion of Nash equilibria by restricting to a fraction 1− p of agents, as e.g. considered in (Carmona, 2004; Elie et al., 2020). Definition 1. An (ε, p)-Markov-Nash equilibrium (almost Markov-Nash equilibrium) for ε, p > 0 is defined as a tuple of policies (π1, . . . , πN ) ∈ ΠN such that for any i ∈ WN , we have
JNi (π 1, . . . , πN ) ≥ sup π∈Π JNi (π 1, . . . , πi−1, π, πi+1, . . . , πN )− ε, (8)
for some setWN ⊆ VN containing at least b(1− p)Nc agents, i.e. |WN | ≥ b(1− p)Nc.
The minimal such ε > 0 for any fixed policy tuple (and typically p = 0) is also called its exploitability. Whilst we ordain ε-optimality only for a fraction 1− p of agents, if the fraction p is negligible, it will have negligible impact on other agents as a result of the weak interaction property. Thus, the solution will remain approximately optimal for almost all agents for sufficiently small p regardless of the behavior of that fraction p of agents. In the following, we will give a limiting system that shall provide (ε, p)-Markov-Nash equilibria with ε, p→ 0 as N →∞.
Graphon mean field game. The formal N → ∞ limit of the N -agent game constitutes its graphon mean field game (GMFG), which shall be rigorously justified in Section 3. We define the space of measurable state marginal ensemblesMt := P(X )I and measurable mean field ensembles M := P(X )T ×I , in the sense that α 7→ µαt (x) is measurable for any µ ∈ M, t ∈ T , x ∈ X . Similarly, we define the space of measurable policy ensembles Π ⊆ ΠI , i.e. with measurable α 7→ παt (u | x) for any π ∈ Π, t ∈ T , x ∈ X , u ∈ U . In the GMFG, we will consider infinitely many agents α ∈ I instead of the finitely many i ∈ VN . As a result, we will have infinitely many policies πα ∈ Π – one for each agent α – through some measurable policy ensemble π ∈ Π. We again define state and action random variables
Xα0 ∼ µ0, Uαt ∼ παt (· | Xαt ), Xαt+1 ∼ P (· | Xαt , Uαt ,Gαt ), ∀(α, t) ∈ I × T (9)
where we introduce the (now deterministic) B1(X )-valued neighborhood mean field of agents α as
Gαt := ∫
I W (α, β)µβt dβ (10)
for some deterministic µ ∈M. Under fixed π ∈ Π, µαt should be understood as the law of Xαt , µαt ≡ L(Xαt ). Finally, define the maximization objective of agent α over πα for fixed µ ∈M as
Jµα (π α) ≡ E
[ T−1∑
t=0
r(Xαt , U α t ,Gαt )
] . (11)
To formulate the limiting version of Nash equilibria, we define a map Ψ: Π→M mapping from a policy ensemble π ∈ Π to the corresponding generated mean field ensemble µ = Ψ(π) ∈M by
µα0 ≡ µ0, µαt+1(x′) ≡ ∑
x∈X µαt (x)
∑ u∈U παt (u | x)P (x′ | x, u,Gαt ), ∀α ∈ I (12)
where integrability in (10) holds by induction, and note how then µαt = L(Xαt ). Similarly, let Φ: M → 2Π map from a mean field ensemble µ to the set of optimal policy ensembles π characterized by πα ∈ arg maxπ∈Π Jµα (π) for all α ∈ I, which is particularly fulfilled if παt (u | x) > 0 =⇒ u ∈ arg maxu′∈U Qµα(t, x, u′) for all α ∈ I, t ∈ T , x ∈ X , u ∈ U , where Qµα is the optimal action value function under fixed µ ∈M following the Bellman equation
Qµα(t, x, u) = r(x, u,Gαt ) + ∑
x′∈X P (x′ | x, u,Gαt ) arg max u′∈U Qµα(t+ 1, x ′, u′) (13)
with Qµα(T, x, u) ≡ 0 and generally time-dependent, see Puterman (2014) for a review. We can now define the GMFG version of Nash equilibria as policy ensembles π generating mean field ensembles µ under which they are optimal, as µαt = L(Xαt ) if all agents α ∈ I follow πα. Definition 2. A Graphon Mean Field Equilibrium (GMFE) is a pair (π,µ) ∈ Π ×M such that π ∈ Φ(µ) and µ = Ψ(π).

heading: 3 THEORETICAL ANALYSIS
text: To obtain meaningful optimality results beyond empirical mean field convergence, we will need a Lipschitz assumption as in the uncontrolled, continuous-time case (cf. Bayraktar et al. (2020), Condition 2.3) and typical in mean field theory (Huang et al., 2006). Assumption 2. Let r, P , W be Lipschitz continuous with Lipschitz constants Lr, LP , LW > 0.
Note that all proofs but Theorem 1 also hold for only block-wise Lipschitz continuous W , see Appendix A.1. Since X × U ×B1(X ) is compact, r is bounded by the extreme value theorem. Proposition 1. Under Assumption 2, r will be bounded by |r| ≤Mr for some constant Mr > 0.
We then obtain existence of a GMFE by reformulating the GMFG as a classical MFG and applying existing results from Saldi et al. (2018). More precisely, we consider the equivalent MFG with extended state space X × I, action space U , policy π̃ ∈ P(U)T ×X×I , mean field µ̃ ∈ P(X × I)T , reward function r̃((x, α), u, µ̃) := r(x, u, ∫ IW (αt, β)µ̃t(·, β) dβ) and transition dynamics such that the states (X̃t, αt) follow (X̃0, α0) ∼ µ̃0 := µ0 ⊗Unif([0, 1]) and
Ũt ∼ π̃t(· | X̃t, αt), X̃t+1 ∼ P (· | X̃t, Ũt, ∫
I W (αt, β)µ̃t(·, β) dβ), αt+1 = αt . (14)
Theorem 1. Under Assumption 2, there exists a GMFE (π,µ) ∈ Π×M.
Meanwhile, in finite games, even showing the existence of Nash equilibria in local feedback policies is problematic (Saldi et al., 2018). Note however, that while this reformulation will be useful for learning and existence, it does not allow us to conclude that the finite graph game is well approximated, as classical MFG approximation theorems e.g. in Saldi et al. (2018) do not consider the graph structure and directly use the limiting graphon W in the dynamics (14).
As our next main result, we shall therefore show rigorously that the GMFE can provide increasingly good approximations of the N -agent finite graph game as N → ∞. As mentioned, the following also holds for only block-wise Lipschitz continuous W instead of fully Lipschitz continuous W . Complete mathematical proofs together with additional theoretical supplements can be found in Appendix A.1 and A.2. To obtain joint N -agent policies as approximate Nash equilibria from a GMFE (π,µ), we define the map ΓN (π) := (π1, π2, . . . , πN ) ∈ ΠN , where
πit(u | x) := παit (u | x), ∀(α, t, x, u) ∈ I × T × X × U (15) with αi = iN , as by Assumption 1 the agents are correctly labeled such that they match up with their limiting graphon indices αi ∈ I. In our experiments, we use the αi generated during the generation process of the W -random graphs, though for arbitrary finite systems one would have to first identify the graphon as well as an appropriate assignment of agents to graphon indices αi ∈ I, which is a separate, non-trivial problem requiring at least graphon estimation, e.g. Xu (2018).
For theoretical analysis, we propose to lift the empirical distributions and policy tuples to the continuous domain I, i.e. under an N -agent policy tuple (π1, . . . , πN ) ∈ ΠN , we define the step policy ensemble πN ∈ Π and the random empirical step measure ensemble µN ∈M by
πN,αt := ∑
i∈VN
1α∈( i−1N , i N ] · πit, µN,αt :=
∑
i∈VN
1α∈( i−1N , i N ] · δXjt , ∀(α, t) ∈ I × T . (16)
In the following, we consider deviations of the i-th agent from (π1, π2, . . . , πN ) = ΓN (π) ∈ ΠN to (π1, . . . , πi−1, π̂, πi+1, . . . , πN ) ∈ ΠN , i.e. the i-th agent deviates by instead applying π̂ ∈ Π. Note that this includes the special case of no agent deviations. For any f : X × I → R and state marginal ensemble µt ∈ Mt, define µt(f) := ∫ I ∑ x∈X f(x, α)µ α t (x) dα. We are now ready to state our first result of convergence of empirical state distributions to the mean field, potentially at the classical rate O(1/ √ N) and consistent with results in uncontrolled, continuous-time diffusive graphon mean field systems (cf. Bayraktar et al. (2020), Theorem 3.2). Theorem 2. Consider Lipschitz continuous π ∈ Π up to a finite number of discontinuities Dπ , with associated mean field ensemble µ = Ψ(π). Under Assumption 1 and the N -agent policy (π1, . . . , πi−1, π̂, πi+1, . . . , πN ) ∈ ΠN with (π1, π2, . . . , πN ) = ΓN (π) ∈ ΠN , π̂ ∈ Π, t ∈ T , we have for all measurable functions f : X × I → R uniformly bounded by some Mf > 0, that
E [∣∣µNt (f)− µt(f) ∣∣]→ 0 (17) uniformly over all possible deviations π̂ ∈ Π, i ∈ VN . Furthermore, if the graphon convergence in Assumption 1 is at rate O(1/ √ N), then this rate of convergence is also O(1/ √ N).
In particular, the technical Lipschitz requirement of π typically holds for neural-network-based policies (Mondal et al., 2021; Pasztor et al., 2021) and includes also the case of finitely many optimality regimes over all graphon indices α ∈ I, which is sufficient to achieve arbitrarily good approximate Nash equilibria through our algorithms as shown in Section 4. We would like to remark that the above result generalizes convergence of state histograms to the mean field solution, since the state marginals of agents are additionally close to each of their graphon mean field equivalents. The above will be necessary to show convergence of the dynamics of a deviating agent to
X̂ i N 0 ∼ µ0, Û i N t ∼ π̂t(· | X̂ i N t ), X̂ i N t+1 ∼ P (· | X̂ i N t , Û i N t ,G i N t ), ∀t ∈ T (18)
for almost all agents i, i.e. the dynamics are approximated by using the limiting deterministic neighborhood mean field G iN , see Appendix A.1. This will imply the approximate Nash property: Theorem 3. Consider a GMFE (π,µ) with Lipschitz continuous π up to a finite number of discontinuities Dπ . Under Assumptions 1 and 2, for any ε, p > 0 there exists N ′ such that for all N > N ′, the policy (π1, . . . , πN ) = ΓN (π) ∈ ΠN is an (ε, p)-Markov Nash equilibrium, i.e.
JNi (π 1, . . . , πN ) ≥ max π∈Π JNi (π 1, . . . , πi−1, π, πi+1, . . . , πN )− ε (19)
for all i ∈ WN and someWN ⊆ VN : |WN | ≥ b(1− p)Nc.
In general, Nash equilibria are highly intractable (Daskalakis et al., 2009). Therefore, solving the GMFG allows obtaining approximate Nash equilibria in theN -agent system for sufficiently largeN , since ε, p→ 0 as N →∞. As a side result, we also obtain first results for the uncontrolled discretetime case by considering trivial action spaces with |U| = 1, see Corollary A.2 in the Appendix.

heading: 4 LEARNING GRAPHON MEAN FIELD EQUILIBRIA
text: By learning GMFE, one may potentially solve otherwise intractable largeN -agent games. For learning, we can apply any existing techniques for classical MFGs (e.g. Mguni et al. (2018); Subramanian & Mahajan (2019); Guo et al. (2019)), since by (14) we have reformulated the GMFG as a classical MFG with extended state space. Nonetheless, it may make sense to treat the graphon index α ∈ I separately, e.g. when treating special cases such as block graphons, or by grouping graphically similar agents. We repeatedly apply two functions Φ̂, Ψ̂ by beginning with the mean field µ0 = Ψ̂(π0) generated by the uniformly random policy π0, and computing πn+1 = Φ̂(µn), µn+1 = Ψ̂(πn+1) for n = 0, 1, . . . until convergence using one of the following two approaches:
• Equivalence classes method. We introduce agent equivalence classes, or discretization, of I for the otherwise uncountably many agents α ∈ I by partitioning I into M subsets. For example, in the special case of block graphons (block-wise constant W ), one can solve separately for each block equivalence class (type) of agents, since all agents in the class share the same dynamics. Note that in contrast to multi-class MFGs (Huang et al., 2006), GMFGs are rigorously connected to finite graph games and can handle an uncountable number of classes α. To deal with general graphons, we choose equidistant representatives αi ∈ I, i = 1, . . . ,M covering the whole interval I, and approximate each agent α ∈ Ĩi by the nearest αi for the intervals Ĩi ⊆ I of points closest to that αi to obtain M approximate equivalence classes. Formally, we approximate mean fields Ψ̂(π) =
∑M i=1 1α∈Ĩi µ̂ αi
recursively computed over all times for any fixed policy ensemble π, and similarly policies Φ̂(µ) = ∑M i=1 1α∈Ĩiπ
αi where παi is the optimal policy of αi for fixed µ. We solve the optimal control problem for each equivalence class using backwards induction (alternatively, one may use reinforcement learning), and solve the evolution equation for the representatives αi of the equivalence classes recursively. For space reasons, the details are found in Appendix A.3. Note that this does not mean that we consider theN -agent problem with N = M , but instead we approximate the limiting problem with the limiting graphon W , and the solution will be near-optimal for all sufficiently large finite systems at once.
• Direct reinforcement learning. We directly apply RL as Φ̂. The central idea is to consider the GMFG as a classical MFG with extended state space X × I, i.e. for fixed mean fields, we solve the MDP defined by (14). Agents condition their policy not only on their own state, but also their node index α ∈ I and the current time t ∈ T , since the mean fields are non-stationary in general and require time-dependent policies for optimality. Here, we assume that we can sample from a simulator of (9) for a given fixed mean field as commonly assumed in MFG learning literature (Guo et al., 2019; Subramanian & Mahajan, 2019). For application to arbitrary finite systems, one could apply a model-based RL approach coupled with graphon estimation, though this remains outside the scope of this work. For solving the mean field evolution equation (12), we can again use any applicable numerical method and choose a conventional sequential Monte Carlo method for Ψ̂. While it is possible to exactly solve optimal control problems for each agent equivalence class with finite stateaction spaces, this is generally not the case for e.g. continuous state-action spaces. Here, a general reinforcement learning solution can solve otherwise intractable problems in an elegant manner, since the graphon index α simply becomes part of a continuous state space.
For convergence, we begin by stating the classical feedback regularity condition (Huang et al., 2006; Guo et al., 2019) after equipping Π, M e.g. with the supremum metric. Proposition 2. Assume that the maps Ψ̂, Φ̂ are Lipschitz with constants c1, c2 and c1 · c2 < 1. Then the fixed point iteration µn+1 = Ψ̂(Φ̂(µn)) converges.
Feedback regularity is not assured, and thus there is no general convergence guarantee. Whilst one could attempt to apply fictitious play (Mguni et al., 2018), additional assumptions will be needed for convergence. Instead, whenever necessary for convergence, we regularize by introducing Boltzmann policies παt (u | x) ∝ exp( 1ηQµα(t, x, u)) with temperature η, provably converging to an approximation for sufficiently high temperatures (Cui & Koeppl, 2021). Theorem 4. Under Assumptions 1 and 2, the equivalence classes algorithm with Boltzmann policies Φ̂(µ)αt (u | x) ∝ exp( 1ηQµα(t, x, u)) converges for sufficiently high temperatures η > 0.
Importantly, even an exact solution of the GMFG only constitutes an approximate Nash equilibrium in the finite-graph system. Furthermore, even the existence of exact finite-system Nash equilibria in local feedback policies is not guaranteed, see the discussion in Saldi et al. (2018) and references therein. Therefore, little is lost by introducing slight additional approximations for the sake of a tractable solution, if at all needed (e.g. the Investment-Graphon problem in the following converges without introducing Boltzmann policies), since near optimality holds for small temperatures (Cui & Koeppl, 2021). Indeed, we find that we can show optimality of the equivalence classes approach for sufficiently fine partitions of I, giving us a theoretical foundation for our proposed algorithms. Theorem 5. Under Assumptions 1 and 2, for a solution (π,µ) ∈ Π×M, π ∈ Φ̂(µ), µ = Ψ̂(π) of the M equivalence classes method and for any ε, p > 0 there exists N ′,M ∈ N such that for all N > N ′, the policy (π1, . . . , πN ) = ΓN (π) ∈ ΠN is an (ε, p)-Markov Nash equilibrium.
A theoretically rigorous analysis of the elegant direct reinforcement learning approach is beyond our scope and deferred to future works, though we empirically find that both methods agree.

heading: 5 EXPERIMENTS
text: In this section, we will give an empirical verification of our theoretical results. As we are unaware of any prior discrete-time GMFGs (except for the example in Vasal et al. (2021), which is similar to the first problem in the following), we propose two problems adapted from existing non-graphbased works on the three graphons in Figure 2. For space reasons, we defer detailed descriptions of problems and algorithms, plots as well as further analysis, including exploitability and a verification of stability of our solution with respect to the number of equivalence classes – to Appendix A.3.
The SIS-Graphon problem was considered in Cui & Koeppl (2021) as a classical discrete-time MFG. We impose an epidemics scenario where people (agents) are infected with probability proportional to the number of infected neighbors and recover with fixed probability. People may choose to take precautions (e.g. social distancing), avoiding potential costly infection periods at a fixed cost.
In the Investment-Graphon problem – an adaptation of a problem studied by Chen et al. (2021), where it was in turn adapted from Weintraub et al. (2010) – we consider many firms maximizing profits, where profits are proportional to product quality and decrease with total neighborhood product quality, i.e. the graph models overlap in e.g. product audience or functionality. Firms can invest to improve quality, though it becomes more unlikely to improve quality as their quality rises.
Learned equilibrium behavior. For the SIS-Graphon problem, we apply softmax policies for each approximate equivalence class to achieve convergence, see Appendix A.3 for details on temperature choice and influence. In Figure 3, the learned behavior can be observed for various α. As expected, in the ER graphon case, behavior is identical over all α. Otherwise, we find that agents take more precautions with many connections (low α) than with few connections (high α). For the uniform attachment graphon, we observe no precautions in case of negligible connectivity (α→ 1), while for the ranked attachment graphon there is no such α ∈ I (cf. Figure 2). Further, the fraction of infected agents at stationarity rises as α falls. A similar analysis holds for Investment-Graphon without need for regularization, see Appendix A.3.
Note that the specific method of solution is not of central importance here, as in general any RL and filtering method can be substituted to handle 1. otherwise intractable or 2. inherently samplebased settings. Indeed, we achieve similar results using PPO (Schulman et al., 2017) in InvestmentGraphon, enabling a general RL-based methodology for GMFGs. In Appendix A.3, we find that PPO achieves qualitatively and quantitatively similar behavior to the equivalence classes method, with slight deviations due to the approximations from PPO and Monte Carlo. In particular, the PPO exploitability ε ≈ 2 remains low compared to ε > 30 for the uniform random policy, see Appendix A.3. In Appendix A.3, we also show how, due to the non-stationarity of the environment, a naive application of MARL (Yu et al., 2021) fails to converge, while existing mean field MARL techniques (Yang et al., 2018b) remain incomparable as agents must observe the average actions of all neighbors. On SIS-Graphon, we require softmax policies to achieve convergence, which is not possible with PPO as no Q-function is learned. In general, one could use entropy regularized policies, e.g. SAC (Haarnoja et al., 2018), or alternatively use any value-based reinforcement learning method, though an investigation of the best approach is outside of our scope.
Quantitative verification of the mean field approximation. To verify the rigorously established accuracy of our mean field system empirically, we will generate W -random graphs. Note that there are considerable difficulties associated with an empirical verification of (19), since 1. for any N one must check the Nash property for (almost) all N agents, 2. finding optimal π̂ is intractable, as no dynamic programming principle holds on the non-Markovian local agent state, while acting on the full state fails by the curse of dimensionality, and 3. the inaccuracy from estimating all JNi , i = 1, . . . , N at once increases with N due to variance, i.e. cost scales fast with N for fixed variance. Instead, we verify (26) in Appendix A.1 using the GMFE policy on systems of up to N = 100 agents, i.e. π̂ = παi for the closest αi and comparing for all agents at once (p = 0). Shown in Figure 4, forW -random graph sequences, at eachN we performed 10000 runs to estimate maxi |JNi − Jαi |. We find that the maximum deviation between achieved returns and mean field return decreases as N → ∞, verifying that we obtain an increasingly good approximation of the finite N -agent graph system. The oscillations in Figure 4 stem from the randomly sampled graphs.

heading: 6 CONCLUSION
text: In this work, we have formulated a new framework for dense graph-based dynamical games with the weak interaction property. On the theoretical side, we have given one of the first general discretetime GMFG formulations with existence conditions and approximate Nash property of the finite graph system, thus extending classical MFGs and allowing for a tractable, theoretically well-founded solution of competitive large-scale graph-based games on large dense graphs. On the practical side, we have proposed a number of computational methods to tractably compute GMFE and experimentally verified the plausibility of our methodology on a number of examples. Venues for further extensions are manifold and include extensions of theory to e.g. continuous spaces, partial observability or common noise. So far, graphons assume dense graphs and cannot properly describe sparse graphs (W = 0), which remain an active frontier of research. Finally, real-world application scenarios may be of interest, where estimation of agent graphon indices becomes important for model-based MARL. We hope that our work inspires further applications and research into scalable MARL using graphical dynamical systems based on graph limit theory and mean field theory.

heading: ACKNOWLEDGMENTS
text: This work has been funded by the LOEWE initiative (Hesse, Germany) within the emergenCITY center. The authors also acknowledge support by the German Research Foundation (DFG) via the Collaborative Research Center (CRC) 1053 – MAKI.

heading: ETHICS STATEMENT
text: Existing mean field methodologies, including ours, currently require manual modeling and have not been applied in a model-based reinforcement learning setting for given finite agent systems. As a result, we do not foresee any immediate ethical issues stemming from this work.

heading: REPRODUCIBILITY STATEMENT
text: For reproducibility, in the supplement we provide all code required to reproduce all results in this work. This includes but is not limited to our models and problems, algorithms as well as all plotting scripts for all of the figures found in this work.

heading: A APPENDIX

heading: A.1 THEORETICAL DETAILS
text: In this section, we will give all intermediate results required to prove the results in the main text, as well as additional results, e.g. for the uncontrolled case. For convenience, we first state all obtained theoretical result. Proofs for each of the theorems and corollaries can be found in their own sections further below.
Note that except for Theorem 1, as mentioned in the main text we can also slightly weaken Assumption 2 to block-wise Lipschitz continuous W , i.e. there exist LW > 0 and disjoint intervals {I1, . . . , IQ}, ∪iIi = I s.t. ∀i, j ∈ {1, . . . , Q},
|W (x, y)−W (x̃, ỹ)| ≤ LW (|x− x̃|+ |y − ỹ|), ∀(x, y), (x̃, ỹ) ∈ Ii × Ij (20) which is fulfilled e.g. for block-wise Lipschitz-continuous or block-wise constant graphons.
For α ∈ I, define the α-neighborhood maps Gα : Mt → P(X ) and empirical α-neighborhood maps GαN : Mt → P(X ) as
Gα(µt) := ∫
I W (α, β)µβt dβ, GαN (µt) :=
∫
I WN (α, β)µ
β t dβ (21)
and note how we naturally have Gαt = Gα(µt) in the mean field system and Git = G i N N (µ N t ) in the finite system. Finally, for ν,ν′ ∈ Mt, π ∈ Π and graphon W , define the ensemble transition kernel operator Pπt,ν,W : Mt →Mt via
( νPπt,ν′,W )α ≡ ∑
x∈X να(x)
∑ u∈U παt (u | x)P
( · | x, u, ∫
I W (α, β)ν′β dβ
) (22)
and note how we have µt+1 = µtPπt,µt,W in the mean field system.
After showing Theorem 2, we continue by showing convergence of the law of deviating agent state Xit to the law of the corresponding auxiliary mean field systems given by
X̂ i N 0 ∼ µ0, Û i N t ∼ π̂t(· | X̂ i N t ), X̂ i N t+1 ∼ P (· | X̂ i N t , Û i N t ,G i N t ), ∀t ∈ T (23)
for almost all agents i as N →∞. Lemma A.1. Consider Lipschitz continuous π ∈ Π up to a finite number of discontinuities Dπ , with associated mean field ensemble µ = Ψ(π). Under Assumptions 1 and 2 and the N -agent policy (π1, . . . , πi−1, π̂, πi+1, . . . , πN ) ∈ ΠN where (π1, π2, . . . , πN ) = ΓN (π) ∈ ΠN , π̂ ∈ Π arbitrary, for any uniformly bounded family of functions G from X to R and any ε, p > 0, t ∈ T , there exists N ′ ∈ N such that for all N > N ′ we have
sup g∈G
∣∣∣E [ g(Xit) ] − E [ g(X̂ i N t ) ]∣∣∣ < ε (24)
uniformly over π̂ ∈ Π, i ∈ WN for someWN ⊆ VN with |WN | ≥ b(1− p)Nc. Similarly, for any uniformly Lipschitz, uniformly bounded family of measurable functions H from X × B1(X ) to R and any ε, p > 0, t ∈ T , there exists N ′ ∈ N such that for all N > N ′ we have
sup h∈H
∣∣∣E [ h(Xit ,G i N N (µ N t )) ] − E [ h(X̂ i N t ,G i N (µt)) ]∣∣∣ < ε (25)
uniformly over π̂ ∈ Π, i ∈ WN for someWN ⊆ VN with |WN | ≥ b(1− p)Nc.
As a direct implication of the above results, the objective functions of almost all agents converge uniformly to the mean field objectives. Corollary A.1. Consider Lipschitz continuous π ∈ Π up to a finite number of discontinuities Dπ , with associated mean field ensemble µ = Ψ(π). Under Assumptions 1 and 2 and the N -agent policy (π1, . . . , πi−1, π̂, πi+1, . . . , πN ) ∈ ΠN where (π1, π2, . . . , πN ) = ΓN (π) ∈ ΠN , π̂ ∈ Π arbitrary, for any ε, p > 0, there exists N ′ ∈ N such that for all N > N ′ we have∣∣∣JNi (π1, . . . , πi−1, π̂, πi+1, . . . , πN )− Jµi
N
(π̂) ∣∣∣ < ε (26)
uniformly over π̂ ∈ Π, i ∈ WN for someWN ⊆ VN with |WN | ≥ b(1− p)Nc.
The approximate Nash property (Theorem 3) of a GMFE (π,µ) then follows immediately from the definition of a GMFE, since π is by definition optimal under µ.
As a corollary, we also obtain results for the uncontrolled case without actions, which is equivalent to the case where |U| = 1, i.e. there being only one trivial policy that is always optimal. Corollary A.2. Under Assumption 1 and |U| = 1, we have for all measurable functions f : X×I → R uniformly bounded by |f | ≤Mf and all t ∈ T that
E [∣∣µNt (f)− µt(f) ∣∣]→ 0 . (27)
Furthermore, if the convergence in Assumption 1 is at rate O(1/ √ N), the rate of convergence is also at O(1/ √ N).
If further Assumption 2 holds, then for any uniformly bounded family of functions G from X to R and any ε, p > 0, t ∈ T , there exists N ′ ∈ N such that for all N > N ′ we have
sup g∈G
∣∣∣E [ g(Xit) ] − E [ g(X̂ i N t ) ]∣∣∣ < ε (28)
uniformly over i ∈ WN for some WN ⊆ VN with |WN | ≥ b(1− p)Nc, and similarly for any uniformly Lipschitz, uniformly bounded family of measurable functions H from X × B1(X ) to R and any ε, p > 0, t ∈ T , there exists N ′ ∈ N such that for all N > N ′ we have
sup h∈H
∣∣∣E [ h(Xit ,G i N N (µ N t )) ] − E [ h(X̂ i N t ,G i N (µt)) ]∣∣∣ < ε (29)
uniformly over i ∈ WN for someWN ⊆ VN with |WN | ≥ b(1− p)Nc.

heading: A.2 PROOFS
text: In this section, we will give full proofs to the statements made in the main text and in Appendix A.1.

heading: A.2.1 PROOF OF THEOREM 1
text: Proof. First, we will verify Saldi et al. (2018), Assumption 1 for the MFG with dynamics given by (14). For this purpose, as in Saldi et al. (2018) let us metrize the product space with the sup-metric, and equip the space P(X × I) with the weak topology. Note that the results hold for both the finite and infinite horizon setting, see Saldi et al. (2018), Remark 6.
(a) The reward function r̃((x, α), u, µ̃) := r(x, u, ∫ IW (αt, β)µ̃t(·, β) dβ) is continuous,
since for ((xn, αn), un, µ̃n)→ ((x, α), u, µ̃) we have ∫
I W (αn, β)µ̃n(·, β) dβ →
∫
I W (α, β)µ̃(·, β) dβ
by Lipschitz continuity of W and weak convergence of µ̃n, and therefore
r(xn, un,
∫
I W (αn, β)µ̃n(·, β) dβ)→ r(x, u,
∫
I W (α, β)µ̃(·, β) dβ)
by Assumption 2.
(b) The action space is compact and the state space is locally compact.
(c) Consider the moment function w(x, α) ≡ 2. In this case, we can choose ζ = 1 (we use ζ instead of α in Saldi et al. (2018)).
(d) The stochastic kernel P̃ that fulfills (14) such that (X̃t+1, αt+1) ∼ P̃ (· | (X̃t, αt), Ũt, µ̃t) is weakly continuous, since for ((xn, αn), un, µ̃n)→ ((x, α), u, µ̃) we again have
∫
I W (αn, β)µ̃n(·, β) dβ →
∫
I W (α, β)µ̃(·, β) dβ
and therefore for any continuous bounded f : X × I → R,∫
X×I f dP̃ (· | (xn, αn), un, µ̃n) =
∫
I
∫
X f dPα(· | (xn, α), un, µ̃n) δαn(dα)
=
∫
X f dP (· | xn, un,
∫
I W (αn, β)µ̃n(·, β) dβ)
→ ∫
X f dP (· | x, u,
∫
I W (α, β)µ̃(·, β) dβ)
=
∫
X×I f dP̃ (· | (x, α), u, µ̃)
by disintegration of P̃ and (14).
(e) By boundedness of r, we trivially have v(x) ≡ 1 ≤ ∞. (f) By boundedness of r, we can trivially choose β = 1 (we flip the usage of β and γ).
(g) As a result of the above choices, ζγβ = γ < 1 trivially.
By Saldi et al. (2018), Theorem 3.3 we have the existence of a mean field equilibrium (π̃, µ̃) with some Markovian feedback policy π̃ acting on the state (X̃t, αt). By defining the mean field and policy ensembles π, µ via παt (u | x) = π̃t(u | x, α), µαt = µ̃t(·, α), we obtain existence of the α-a.e. optimal policy ensemble π, since at any time t ∈ T , the joint state-action distribution µ̃t⊗ π̃t puts mass 1 on optimal state-action pairs (see Saldi et al. (2018), Theorem 3.6), implying that for a.e. α the policy must be optimal, as otherwise there exists a non-null set Ĩ0 ⊆ I such that for all α ∈ Ĩ0, there is some suboptimality ε > 0, which directly contradicts the prequel. For the remaining suboptimal α ∈ I0 in the null set I0 ⊆ I, we redefine π optimally for those α (always possible in our case, see e.g. Puterman (2014)). This policy ensemble generates µ = Ψ(π) α-a.e. uniquely, and we need only consider its α-a.e. unique equivalence class for optimality, implying π ∈ Φ(µ). Furthermore, µ is always measurable by definition, whereas π is measurable because π̃t is by definition a Markov kernel, and thus π̃t(u | ·, ·) = π̃t({u} | ·, ·) for Borel set {u} is a measurable function, which implies measurability of π̃t(u | x, ·) (see e.g. Yeh (2014), Appendix E). Therefore, we have proven existence of the GMFE (π,µ).

heading: A.2.2 PROOF OF THEOREM 2
text: Proof. The proof is by induction as follows.
Initial case. For t = 0, we trivially have for all measurable functions f : X × I → R uniformly bounded by |f | ≤Mf a law of large numbers result
E [∣∣µN0 (f)− µ0(f) ∣∣]
= E [∣∣∣∣∣ ∫
I
∑ x∈X µN,α0 (x) f(x, α)− ∑ x∈X µα0 (x) f(x, α) dα ∣∣∣∣∣ ]
= E [∣∣∣∣∣ 1 N ∑
i∈VN
(∫
( i−1N , i N ]
f(Xi0, α) dα− E [∫
( i−1N , i N ]
f(Xi0, α) dα ])∣∣∣∣∣ ]
≤  E   ( 1
N
∑
i∈VN
(∫
( i−1N , i N ]
f(Xi0, α) dα− E [∫
( i−1N , i N ]
f(Xi0, α) dα
]))2    1 2
=   1 N2 ∑
i∈VN
E
  (∫
( i−1N , i N ]
f(Xi0, α) dα− E [∫
( i−1N , i N ]
f(Xi0, α) dα
])2    1 2
≤ 2Mf√ N
by definition of µNt , independence of {Xi0}i∈VN and Xi0 ∼ µ0 = µα0 for all i ∈ VN , α ∈ I, where the second equality follows from Fubini’s theorem.
Induction step. Assume that the induction assumption holds at t. Then by definition of µNt , for all bounded function f : X × I → R with |f | ≤Mf ,
E [∣∣µNt+1(f)− µt+1(f) ∣∣] ≤ E [∣∣∣µNt+1(f)− µNt Pπ N t,µNt ,WN (f) ∣∣∣ ]
+ E [∣∣∣µNt Pπ N
t,µNt ,WN (f)− µNt Pπ
N t,µNt ,W (f) ∣∣∣ ]
+ E [∣∣∣µNt Pπ N
t,µNt ,W (f)− µNt Pπt,µNt ,W (f)
∣∣∣ ]
+ E [∣∣∣µNt Pπt,µNt ,W (f)− µ N t P π t,µt,W (f) ∣∣∣ ] + E [∣∣µNt Pπt,µt,W (f)− µt+1(f) ∣∣] .
First term. We have by definition of µNt E [∣∣∣µNt+1(f)− µNt Pπ N t,µNt ,WN (f) ∣∣∣ ]
= E [∣∣∣∣∣ ∫
I
∑ x∈X µN,αt+1 (x) f(x, α) dα
− ∫
I
∑ x∈X µN,αt (x) ∑ u∈U πN,αt (u | x) ∑ x′∈X P
( x′ | x, u, ∫
I WN (α, β)µ
N,β t dβ
) f(x′, α) dα ∣∣∣∣∣ ]
= E [∣∣∣∣∣ 1 N ∑
i∈VN
(∫
( i−1N , i N ]
f(Xit+1, α) dα− E [∫
( i−1N , i N ]
f(Xit+1, α) dα ∣∣∣∣∣ Xt ])∣∣∣∣∣ ]
≤  E   ( 1
N
∑
i∈VN
(∫
( i−1N , i N ]
f(Xit+1, α) dα− E [∫
( i−1N , i N ]
f(Xit+1, α) dα ∣∣∣∣∣ Xt ]))2    1 2
=   1 N2 ∑
i∈VN
E
  (∫
( i−1N , i N ]
f(Xit+1, α) dα− E [∫
( i−1N , i N ]
f(Xit+1, α) dα ∣∣∣∣∣ Xt ])2    1 2
≤ 2Mf√ N
where the last equality follows from conditional independence of {Xit+1}i∈VN given Xt ≡ {Xit}i∈VN and the law of total expectation.
Second term. We have E [∣∣∣µNt Pπ N
t,µNt ,WN (f)− µNt Pπ
N t,µNt ,W (f) ∣∣∣ ]
= E [∣∣∣∣∣ ∫
I
∑ x∈X µN,αt (x) ∑ u∈U πN,αt (u | x) ∑ x′∈X P
( x′ | x, u, ∫
I WN (α, β)µ
N,β t dβ
) f(x′, α) dα
− ∫
I
∑ x∈X µN,αt (x) ∑ u∈U πN,αt (u | x) ∑ x′∈X P
( x′ | x, u, ∫
I W (α, β)µN,βt dβ
) f(x′, α) dα ∣∣∣∣∣ ]
≤ |X |MfLP E [∫
I
∥∥∥∥ ∫
I WN (α, β)µ
N,β t dβ −
∫
I W (α, β)µN,βt dβ
∥∥∥∥ dα ]
≤ |X |2MfLP sup x∈X
E [∫
I
∣∣∣∣ ∫
I WN (α, β)µ
N,β t (x)−W (α, β)µN,βt (x) dβ ∣∣∣∣ dα ] → 0
by Assumption 1 and µN,βt (x) trivially being bounded by 1. If the convergence in Assumption 1 is at rate O(1/ √ N), then this convergence is also at rate O(1/ √ N).
Third term. We have
E [∣∣∣µNt Pπ N
t,µNt ,W (f)− µNt Pπt,µNt ,W (f)
∣∣∣ ]
= E [∣∣∣∣∣ ∫
I
∑ x∈X µN,αt (x) ∑ u∈U πN,αt (u | x) ∑ x′∈X P
( x′ | x, u, ∫
I W (α, β)µN,βt dβ
) f(x′, α) dα
− ∫
I
∑ x∈X µN,αt (x) ∑ u∈U παt (u | x) ∑ x′∈X P
( x′ | x, u, ∫
I W (α, β)µN,βt dβ
) f(x′, α) dα ∣∣∣∣∣ ]
≤ |X ||U|Mf E [∫
I
∣∣∣πN,αt (u | x)− παt (u | x) ∣∣∣ dα
]
= |X ||U|Mf E
  ∑
j∈VN\{i}
∫
( j−1N , j N ]
∣∣∣∣π dNαe N t (u | x)− παt (u | x) ∣∣∣∣ dα
 
+ |X ||U|Mf E [∫
( i−1N , i N ]
|π̂t(u | x)− παt (u | x)| dα ]
≤ |X ||U|Mf · Lπ N + |X ||U|Mf · 2|Dπ| N + |X ||U|Mf · 2 N
by assumption of Lipschitz continuous π up to a finite number of discontinuities Dπ as well as the deviating agent i’s error term, for which the integrands are bounded by 2.
Fourth term. We have
E [∣∣∣µNt Pπt,µNt ,W (f)− µ N t P π t,µt,W (f) ∣∣∣ ]
= E [∣∣∣∣∣ ∫
I
∑ x∈X µN,αt (x) ∑ u∈U παt (u | x) ∑ x′∈X P
( x′ | x, u, ∫
I W (α, β)µN,βt dβ
) f(x′, α) dα
− ∫
I
∑ x∈X µN,αt (x) ∑ u∈U παt (u | x) ∑ x′∈X P
( x′ | x, u, ∫
I W (α, β)µβt dβ
) f(x′, α) dα ∣∣∣∣∣ ]
≤Mf |X |E [ sup x,u ∫
I
∣∣∣∣P ( x′ | x, u, ∫
I W (α, β)µN,βt dβ
)
− P ( x′ | x, u, ∫
I W (α, β)µβt dβ
)∣∣∣∣ dα ]
≤Mf |X |LP ∑ x′∈X E [∫ I ∣∣∣∣ ∫ I W (α, β)µN,βt (x ′) dβ − ∫ I W (α, β)µβt (x ′) dβ ∣∣∣∣ dα ] ≤Mf |X |2LP · C ′(1)√ N
in the case of rate O(1/ √ N), or uniformly to zero otherwise, from Lipschitz P by defining the functions f ′x′,α(x, β) = W (α, β)·1x=x′ for any (x′, α) ∈ X×I and using the induction assumption on f ′x′,α to obtain
E [∫
I
∣∣∣∣ ∫
I W (α, β)µN,βt (x
′) dβ − ∫
I W (α, β)µβt (x ′) dβ
∣∣∣∣ dα ]
=
∫ I E [∣∣∣∣ ∫ I W (α, β)µN,βt (x ′) dβ − ∫ I W (α, β)µβt (x ′) dβ ∣∣∣∣ ] dα
=
∫
I E [∣∣µNt (f ′x′,α)− µt(f ′x′,α) ∣∣] dα ≤ C ′(1)√ N
for some C ′(1) > 0 uniformly over all f ′ bounded by 1 if the convergence in Assumption 1 is at rate O(1/ √ N), or uniformly to zero otherwise.
Fifth term. We have
E [∣∣µNt Pπt,µt,W (f)− µtPπt,µt,W (f) ∣∣]
= E [∣∣∣∣∣ ∫
I
∑ x∈X µN,αt (x) ∑ u∈U παt (u | x) ∑ x′∈X P
( x′ | x, u, ∫
I W (α, β)µβt dβ
) f(x′, α) dα
− ∫
I
∑ x∈X µαt (x) ∑ u∈U παt (u | x) ∑ x′∈X P
( x′ | x, u, ∫
I W (α, β)µβt dβ
) f(x′, α) dα ∣∣∣∣∣ ]
= E [∣∣∣∣∣ ∫
I
∑ x∈X µN,αt (x)f ′(x, α) dα− ∫ I ∑ x∈X µαt (x)f ′(x, α) dα ∣∣∣∣∣ ]
= E [∣∣µNt (f ′)− µt(f ′) ∣∣] ≤ C ′(Mf )√ N .
in the case of rateO(1/ √ N), or uniformly to zero otherwise, again by induction assumption applied to the function
f ′(x, α) = ∑
u∈U παt (u | x)
∑
x′∈X P
( x′ | x, u, ∫
I W (α, β)µβt dβ
) f(x′, α)
bounded by Mf . This completes the proof by induction.

heading: A.2.3 PROOF OF LEMMA A.1
text: Proof. First, we will show that (24) implies (25).
Proof of (24) =⇒ (25). We consider a uniformly Lipschitz, uniformly bounded family of measurable functionsH from X ×B1(X ) to R. Let Mh be the uniform bound of functions inH and Lh be the uniform Lipschitz constant. Then, for arbitrary h ∈ H we have
∣∣∣E [ h(Xit ,G i N N (µ N t )) ] − E [ h(X̂ i N t ,G i N (µt)) ]∣∣∣
= ∣∣∣E [ h(Xit ,G i N N (µ N t )) ] − E [ h(Xit ,G i N N (µt)) ]∣∣∣
+ ∣∣∣E [ h(Xit ,G i N N (µt)) ] − E [ h(Xit ,G i N (µt)) ]∣∣∣ + ∣∣∣E [ h(Xit ,G i N (µt)) ] − E [ h(X̂ i N t ,G i N (µt)) ]∣∣∣
which we will analyze in the following.
First term. We have ∣∣∣E [ h(Xit ,G i N N (µ N t )) ] − E [ h(Xit ,G i N N (µt)) ]∣∣∣
≤ E [ E [∣∣∣h(Xit ,G i N N (µ N t ))− h(Xit ,G i N N (µt)) ∣∣∣ ∣∣∣ Xit ]] ≤ Lh E [∥∥∥G i N N (µ N t )−G i N N (µt) ∥∥∥ ] = Lh ∑
x∈X E [∣∣∣∣ ∫ I WN ( i N , β)µN,βt (x) dβ − ∫ I WN ( i N , β)µβt (x) dβ ∣∣∣∣ ] ≤ C(1)√ N
by Theorem 2 applied to the functions f ′N,i,x(x ′, β) = WN ( i N , β) · 1x=x′ uniformly bounded by 1.
Second term. Similarly, we have ∣∣∣E [ h(Xit ,G i N N (µt)) ] − E [ h(Xit ,G i N (µt)) ]∣∣∣
≤ Lh‖G i N N (µt)−G i N (µt)‖1
≤ Lh ∑
x∈X
∣∣∣∣ ∫
I
( WN ( i
N , β)−W ( i N , β)
) µβt (x) dβ ∣∣∣∣
≤ Lh ∑
x∈X
∣∣∣∣∣ ∫
I
( WN ( i
N , β)−N
∫
( i−1N , i N ]
W (α, β) dα ) µβt (x) dβ ∣∣∣∣∣
+ Lh ∑
x∈X
∣∣∣∣∣ ∫
I
( N ∫
( i−1N , i N ]
W (α, β) dα−W ( i N , β)
) µβt (x) dβ ∣∣∣∣∣
where the latter term can be bounded as
Lh ∑
x∈X
∣∣∣∣∣ ∫
I
( N ∫
( i−1N , i N ]
W (α, β) dα−W ( i N , β)
) µβt (x) dβ ∣∣∣∣∣
≤ Lh ∑
x∈X
∣∣∣∣∣ ∫ I N ∫
( i−1N , i N ]
( W (α, β)−W (dNαe
N , β)
) µβt (x) dα dβ ∣∣∣∣∣
≤ Lh|X |N · 1 N · LW N = LWLh|X | N
by Assumption 2.
Alternatively, if we assumed the weaker block-wise Lipschitz condition onW in (20), we can obtain the same result for almost all i ∈ VN , i.e. for any p0 > 0 there exists N ′ ∈ N such that for any N > N ′, there exists a setW0N , |W0N | ≥ b(1− p0)Nc such that for all i ∈ W0N the above is true: Since by (20) there exist only a finite number Q of intervals and therefore jumps, there can be only Q many i for which the above fails, while for all other i we again have
∣∣∣∣∣ ∫ I N ∫
( i−1N , i N ]
( W (α, β)−W (dNαe
N , β)
) µβt (x) dα dβ ∣∣∣∣∣
≤ ∑
j∈{1,...,Q}
∣∣∣∣∣ ∫
Ij N
∫
( i−1N , i N ]
( W (α, β)−W (dNαe
N , β)
) µβt (x) dα dβ ∣∣∣∣∣
≤ N · 1 N · LW N = LWLh|X | N
by (20), as ( i−1N , i N ]× Ij ⊆ Ik × Ij for some k ∈ {1, . . . , Q}.
For the former term we observe that
Lh ∑
x∈X
∣∣∣∣∣ ∫
I
( WN ( i
N , β)−N
∫
( i−1N , i N ]
W (α, β) dα ) µβt (x) dβ ∣∣∣∣∣
≤ Lh ∑
x∈X N
∫
( i−1N , i N ]
∣∣∣∣ ∫
I (WN (α, β)−W (α, β))µβt (x) dβ
∣∣∣∣dα
and by defining for any x ∈ X the terms INi (x) via
INi (x) := N
∫
( i−1N , i N ]
∣∣∣∣ ∫
I (WN (α, β)−W (α, β))µβt (x) dβ
∣∣∣∣ dα
and noticing that we have
1
N
N∑
i=1
INi (x) =
∫
I
∣∣∣∣ ∫
I (WN (α, β)−W (α, β))µβt (x) dβ
∣∣∣∣dα→ 0
by Assumption 1, we can conclude that for any ε1, p1 > 0 there exists N ′ ∈ N such that for any N > N ′, there exists a setW1N , |W1N | ≥ b(1− p1)Nc such that for all i ∈ W1N we have
INi (x) < ε1,
since by the above we can choose N ′ ∈ N such that for any N > N ′ we have 1N ∑N i=1 I N i (x) <
ε1p1, and from INi (x) ≥ 0 it would otherwise follow that 1N ∑N i=1 I N i (x) ≥ 1N · dp1Ne ε1 ≥ ε1p1 which would be a direct contradiction. Therefore, for all i ∈ W1N , we have uniformly
Lh ∑
x∈X N
∫
( i−1N , i N ]
∣∣∣∣ ∫
I (WN (α, β)−W (α, β))µβt (x) dβ
∣∣∣∣dα = Lh ∑
x∈X INi (x)→ 0 .
Third term. By (24), for any ε2, p2 > 0 there exists a setW2N , |W2N | ≥ b(1− p2)Nc such that for all i ∈ W2N we have∣∣∣E [ h(Xit ,G i N (µt)) ] − E [ h(X̂ i N t ,G i N (µt))
]∣∣∣ < ε2 independent of π̂ ∈ Π. The intersection ofW0N ,W1N ,W2N has at least N −dp0Ne−dp1Ne−dp2Ne agents fulfilling (25), which completes the proof of (24) =⇒ (25) for almost all agents by choosing ε1, ε2 sufficiently small and p0, p1, p2 < p3 such thatN−dp0Ne−dp1Ne−dp2Ne ≥ b(1− p)Nc, which is equivalent to 1 − dp0NeN − dp1Ne N − dp2Ne N ≥ b(1−p)Nc N and is true for sufficiently large N , since in the limit, 1− dp0NeN − dp1Ne N − dp2Ne N → 1− p0 − p1 − p2 and b(1−p)Nc N → 1− p as N →∞.
Proof of (24). All that remains is to show (24), which will automatically imply (25) at all times t ∈ T by the prequel. We will show (24) by induction.
Initial case. At t = 0, L(Xit) = µ0 = L(X̂ i N t ) by definition. Thus, trivially∣∣∣E [ g(Xi0) ] − E [ g(X̂ i N 0 ) ]∣∣∣ = 0 < ε .
Induction step. For any uniformly bounded family of functions G from X to R with bound Mg , we will show that for any ε, p > 0, there exists N ′ ∈ N such that for all N > N ′ we have
∣∣∣E [ g(Xit+1) ] − E [ g(X̂ i N t+1) ]∣∣∣ < ε
uniformly over π̂ ∈ Π, i ∈ WN for someWN ⊆ VN with |WN | ≥ b(1− p)Nc. Observe that ∣∣∣E [ g(Xit+1) ] − E [ g(X̂ i N t+1) ]∣∣∣ = ∣∣∣E [ lN,t(X i t ,G i N N (µ N t )) ] − E [ lN,t(X̂ i N t ,G i N (µt)) ]∣∣∣
where we defined the uniformly bounded, uniformly Lipschitz functions
lN,t(x, ν) ≡ ∑
u∈U π̂t(u | x)
∑
x′∈X P (x′ | x, u, ν)g(x′)
with Lipschitz constant |X |MgLP and uniform bound Mg . By the induction assumption and (24) =⇒ (25) from the prequel, there exists N ′ ∈ N such that for all N > N ′ we have
∣∣∣E [ lN,t(X i t ,G i N N (µ N t )) ] − E [ lN,t(X̂ i N t ,G i N (µt)) ]∣∣∣ < ε
uniformly over π̂ ∈ Π, i ∈ WN for someWN ⊆ VN with |WN | ≥ b(1− p)Nc, which completes the proof by induction.

heading: A.2.4 PROOF OF COROLLARY A.1
text: Proof. Define the uniformly bounded, uniformly Lipschitz functions
rπ̂(x, ν) ≡ ∑
u∈U r(x, u, ν)π̂t(u | s)
with Lipschitz constant |U |Lr and uniform bound Mr such that by Lemma A.1 and Fubini’s theorem, there exists N ′ ∈ N such that for all N > N ′ we have
∣∣∣JNi (π1, . . . , πi−1, π̂, πi+1, . . . , π̂)− Jµi N (π̂) ∣∣∣
≤ T−1∑
t=0
∣∣∣E [ rπ̂t(X i t ,G i N (µt)) ] − E [ rπ̂t(X̂ i N t ,G i N (µt)) ]∣∣∣ < ε .
uniformly over π̂ ∈ Π, i ∈ WN for some WN ⊆ VN with |WN | ≥ b(1− p)Nc by choosing the maximum over all N ′ at each finite time step from Lemma A.1.
In case of the infinite horizon discounted objective, we instead first cut off at a time T > log
ε(1−γ) 4Mr
log γ
such that trivially
T−1∑
t=0
γt ∣∣∣E [ rπ̂t(X i t ,G i N (µt)) ] − E [ rπ̂t(X̂ i N t ,G i N (µt)) ]∣∣∣
+ γT ∞∑
t=T
γt−T ∣∣∣E [ rπ̂t(X i t ,G i N (µt)) ] − E [ rπ̂t(X̂ i N t ,G i N (µt)) ]∣∣∣
<
T−1∑
t=0
γt ∣∣∣E [ rπ̂t(X i t ,G i N (µt)) ] − E [ rπ̂t(X̂ i N t ,G i N (µt)) ]∣∣∣+ ε 2
and then handle the remaining term analogously to the finite horizon case.

heading: A.2.5 PROOF OF THEOREM 3
text: Proof. By Corollary A.1, for any ε > 0 there exists N ′ ∈ N such that for all N > N ′ we have max π∈Π ( JNi (π 1, . . . , πi−1, π, πi+1, . . . , πN )− JNi (π1, . . . , πN ) )
≤ max π∈Π
( JNi (π
1, . . . , πi−1, π, πi+1, . . . , πN )− Jµi N
(π) )
+ max π∈Π ( Jµi N (π)− Jµi N (π i N ) ) + ( Jµi N (π i N )− JNi (π1, . . . , πN ) )
< ε
2 + 0 +
ε 2 = ε
uniformly over i ∈ WN for some WN ⊆ VN with |WN | ≥ b(1− p)Nc, since π i N ∈ arg maxπ J µ i N (π) by definition of a GMFE. Reordering completes the proof.

heading: A.2.6 PROOF OF COROLLARY A.2
text: Proof. The proof follows immediately from Theorem 2 and Lemma A.1 by considering the trivial policy π that always chooses the only action available together with its generated mean field µ = Ψ(π).

heading: A.2.7 PROOF OF PROPOSITION 2
text: Proof. The set Π is a complete metric space, since existence of limits follows from completeness of R, pointwise limits of measurable functions are measurable, and policies will remain normalized. Banach’s fixed point theorem applied to Φ̂ ◦ Ψ̂ gives us the desired result.

heading: A.2.8 PROOF OF THEOREM 4
text: Proof. Formally, we approximate mean fields by Ψ̂(π) = ∑M i=1 1α∈Ĩi µ̂
αi for any fixed policy ensemble π, and similarly policies Φ̂(µ) =
∑M i=1 1α∈Ĩiπ
αi where παi is the softmax policy of αi for fixed µ, i.e.
παt (u | x) = exp
( Qµα(t,x,u)
η
)
∑ u∈U exp
( Qµα(t,x,u)
η
) . (30)
By the Bellman equation (13), Qµα(t, x, u) is Lipschitz in µ for all (t, x, u) ∈ T × X × U under Assumption 2. Since the Lipschitz constants are shared over all α, by Cui & Koeppl (2021), Lemma B.7.5, (30) is therefore Lipschitz with Lipschitz constant proportional to 1/η, which immediately implies that Φ̂ is also Lipschitz with Lipschitz constant c1/η. By its recursive definition as compositions of Lipschitz functions, Ψ̂ is Lipschitz as well with some constant c2. Therefore, the composition of both functions Ψ̂ ◦ Φ̂ is Lipschitz with constants c1c2/η, which will be less than 1 for sufficiently large η. By Proposition 2, the equivalence classes algorithm Ψ̂ ◦ Φ̂ converges to a fixed point.

heading: A.2.9 PROOF OF THEOREM 5
text: Proof. First, note that under the equivalence classes method, the distance between any α and its representant αi uniformly shrinks to zero as M →∞, i.e. maxi=1,...,M supα∈Ĩi |α− αi| → 0. We begin by showing that a solution of the M equivalence classes method (π,µ) ∈ Π ×M, π ∈ Φ̂(µ), µ = Ψ̂(π) following (31), (32) fulfills approximate optimality, i.e. for any ε > 0 there exists M ′ s.t. for all M > M ′
sup α∈I max π∈Π
( J µ̄α (π)− J µ̄α (πα) ) < ε,
where we introduced the true, exact mean field ensemble µ̄ = Ψ(π) following (12) generated by the block-wise solution policy
∑M i=1 1α∈Ĩiπ
αi of the M equivalence classes method, as well as the true mean field system under µ̄ and any policy π ∈ Π
X̄α0 ∼ µ0, Ūαt ∼ πt(· | X̄αt ), X̄αt+1 ∼ P (· | X̄αt , Ūαt , Ḡαt ), ∀(α, t) ∈ I × T
with B1(X )-valued Ḡαt := ∫ IW (α, β)µ̄ β t dβ and J µ̄ α (π) ≡ E [∑T−1 t=0 r(X̄ α t , Ū α t , Ḡαt ) ] , while system (9) is to be understood as the system under the approximate mean field ensemble µ.
To see this, we will analyze
sup α∈I max π∈Π
( J µ̄α (π)− J µ̄α (πα) ) ≤ max i=1,...,M
sup α∈Ĩi max π∈Π
( J µ̄α (π)− Jµα (π) )
+ max i=1,...,M sup α∈Ĩi max π∈Π
( Jµα (π)− Jµαi(π) )
+ max i=1,...,M sup α∈Ĩi max π∈Π
( Jµαi(π)− Jµαi(παi) )
+ max i=1,...,M sup α∈Ĩi
( Jµαi(π αi)− Jµα (παi) )
+ max i=1,...,M sup α∈Ĩi
( Jµα (π α)− J µ̄α (πα) ) .
First term. For any π ∈ Π, define the uniformly bounded, uniformly Lipschitz functions
rπ(x, ν) ≡ ∑
u∈U r(x, u, ν)πt(u | s)
with Lipschitz constant |U |Lr and uniform bound Mr such that for the first term, we have ( J µ̄α (π)− Jµα (π) ) ≤ ∣∣J µ̄α (π)− Jµα (π) ∣∣
≤ T−1∑
t=0
∣∣E [ rπ(X̄ α t , Ḡαt ) ] − E [rπ(Xαt ,Gαt ] ∣∣
≤ T−1∑
t=0
∣∣E [ rπ(X̄ α t , Ḡαt )− rπ(X̄αt ,Gαt ]∣∣+ T−1∑
t=0
∣∣E [ rπ(X̄ α t ,Gαt ] − E [rπ(Xαt ,Gαt ] ∣∣
≤ T−1∑
t=0
|U |Lr ∥∥∥∥ ∫
I W (α, β)(µ̄βt − µβt ) dβ)
∥∥∥∥+ T−1∑
t=0
∣∣E [ rπ(X̄ α t ,Gαt ] − E [rπ(Xαt ,Gαt ] ∣∣ .
For the former term, note that ∥∥∥∥ ∫
I W (α, β)(µ̄βt − µβt ) dβ) ∥∥∥∥ = ∥∥∥∥∥ ∑
i
∫
Ĩi W (α, β)(µ̄βt − µαit ) dβ) ∥∥∥∥∥ ≤ max i=1,...,M
sup α∈Ĩi |X | ‖µ̄αt − µαit ‖
and we will show by induction over t = 0, 1, . . . , T that supα∈Ĩi ‖µ̄αt − µ αi t ‖ → 0 over all α uniformly over all equivalence classes Ĩi. At t = 0, we have trivially µ̄0 = µ0. Assume that supα∈Ĩi ‖µ̄αt − µ αi t ‖ → 0. Then for t+ 1, we have
sup α∈Ĩi
∥∥µ̄αt+1 − µαit+1 ∥∥
= sup α∈Ĩi
∥∥∥∥∥ ∑
x∈X µ̄αt (x)
∑ u∈U παt (u | x)P (· | x, u, Ḡαt )− ∑ x∈X µαit (x) ∑ u∈U παit (u | x)P (· | x, u,Gαit ) ∥∥∥∥∥
≤ sup α∈Ĩi
∥∥∥∥∥ ∑
x∈X µ̄αt (x)
∑ u∈U παt (u | x)P (· | x, u, Ḡαt )− ∑ x∈X µ̄αit (x) ∑ u∈U παit (u | x)P (· | x, u, Ḡαit ) ∥∥∥∥∥
+ ∥∥∥∥∥ ∑
x∈X µ̄αit (x)
∑ u∈U παit (u | x)P (· | x, u, Ḡαit )− ∑ x∈X µαit (x) ∑ u∈U παit (u | x)P (· | x, u, Ḡαit ) ∥∥∥∥∥
+ ∥∥∥∥∥ ∑
x∈X µ̄αit (x)
∑ u∈U παit (u | x)P (· | x, u, Ḡαit )− ∑ x∈X µαit (x) ∑ u∈U παit (u | x)P (· | x, u,Gαit ) ∥∥∥∥∥
≤ sup α∈Ĩi
∥∥∥∥∥ ∑
x∈X µ̄αt (x)
∑ u∈U παt (u | x)P (· | x, u, Ḡαt )− ∑ x∈X µ̄αit (x) ∑ u∈U παit (u | x)P (· | x, u, Ḡαit ) ∥∥∥∥∥
+ |X |2 ‖µ̄αit − µαit ‖+ |X |2|U|LP ‖µ̄αit − µαit ‖ → 0 asM →∞, since the first term is uniformly Lipschitz in α by (12) as a recursive composition, finite multiplication and addition of Lipschitz functions, whereas the other terms tend to zero by induction hypothesis. Since the Lipschitz constants do not depend on Ĩi, the convergence is uniform. To bound the latter term, we first note that rπ(·,Gαt ) is always bounded by Mr regardless of t, α, π, i.e. it again suffices to show that for any family of functions G from X to R uniformly bounded by Mr, we have
sup g∈G
∣∣E [ g(X̄αt ) ] − E [g(Xαt )] ∣∣→ 0 .
The proof is by induction. At t = 0, we trivially have L(X̄αt ) = µ0 = L(Xαt ). Assuming that the induction hypothesis holds at t, then at t+ 1 we have
sup g∈G
∣∣E [ g(X̄αt+1) ] − E [ g(Xαt+1) ]∣∣ = sup g∈G ∣∣E [ lt(X̄ α t ,Gαt )) ] − E [lt(Xαt ,Gαt )] ∣∣→ 0
by the induction hypothesis, where we defined the uniformly bounded functions
lt(x, ν) ≡ ∑
u∈U πt(u | x)
∑
x′∈X P (x′ | x, u, ν)g(x′)
with uniform bound Mr. Therefore, |J µ̄α (π)− Jµα (π)| → 0 uniformly over all α, π.
Second term. For the second term, we analogously have ( Jµα (π)− Jµαi(π) ) ≤ ∣∣Jµα (π)− Jµαi(π) ∣∣
≤ T−1∑
t=0
|U |Lr ∥∥∥∥ ∫
I (W (α, β)−W (αi, β))µβt dβ)
∥∥∥∥+ T−1∑
t=0
|E [rπ(Xαt ,Gαit ]− E [rπ(Xαit ,Gαit ]|
where the former term uniformly tends to zero as M → ∞ over all α by Lipschitz W from Assumption 2 and increasingly fine partition intervals Ĩi, while for the latter term we again show that for any family of functions G from X to R uniformly bounded by Mr, we have
sup g∈G |E [g(Xαt )]− E [g(Xαit )]| → 0 .
The proof is by induction. At t = 0, we trivially have L(Xαt ) = µ0 = L(Xαit ). Assuming that the induction hypothesis holds at t, then at t+ 1 we have
sup g∈G
∣∣E [ g(Xαt+1) ] − E [ g(Xαit+1) ]∣∣ = sup g∈G |E [lt(Xαt ,Gαit ))]− E [lt(Xαit ,Gαit )]| → 0
by the induction hypothesis, where we defined the uniformly bounded functions
lt(x, ν) ≡ ∑
u∈U πt(u | x)
∑
x′∈X P (x′ | x, u, ν)g(x′)
with uniform bound Mr. Therefore, ∣∣Jµα (π)− Jµαi(π) ∣∣→ 0 uniformly over all α, π.
Third term. By definition, we have optimality of π ∈ Φ̂(µ) under the approximate mean field µ at each representative αi. Therefore, the term maxπ∈Π ( Jµαi(π)− Jµαi(παi) ) is upper bounded by 0, as there is no policy π that improves over παi .
Fourth and fifth term. The results follow from the first and second term by inserting πα for π.
Variations on the setting. The infinite horizon discounted case is handled as in the proof of Corollary A.1, i.e. repeating the above up to some chosen time horizon T and trivially bounding all terms with t ≥ T . The block-wise Lipschitz graphon case (20) is handled by choosing the equivalence classes Īi ⊆ Ij such that they are part of at most one block Ij of the graphon.
Proof of Theorem 5. Now fix any ε, p > 0. As a result of the prequel, we have that there exists M ′ s.t. for all M > M ′
sup α∈I max π∈Π |Jµα (πα)− Jµα (π)| <
ε 3 .
Pick any such M > M ′. By Corollary A.1 (for the first and third term, since π is constant with at most M discontinuities) and the prequel (for the second term), there exists N ′ ∈ N such that for all N > N ′ we have
max π∈Π
( JNi (π 1, . . . , πi−1, π, πi+1, . . . , πN )− JNi (π1, . . . , πN ) )
≤ max π∈Π
( JNi (π
1, . . . , πi−1, π, πi+1, . . . , πN )− J µ̄i N
(π) )
+ max π∈Π ( J µ̄i N (π)− J µ̄i N (π i N ) ) + ( J µ̄i N (π i N )− JNi (π1, . . . , πN ) )
< ε 3 + ε 3 + ε 3 = ε
which holds uniformly over i ∈ WN for someWN ⊆ VN with |WN | ≥ b(1− p)Nc, since π i N ∈ arg maxπ J µ i N (π) by definition of a GMFE. Reordering completes the proof.

heading: A.3 EXPERIMENTAL DETAILS
text: In this section, we will give a full description of all the algorithms and hyperparameters we used during our experiments. For reinforcement learning, we use PPO (Schulman et al., 2017).
For the approximate equivalence classes, we shall consider grids (αm ∈ [0, 1])m=1,...,M with associated policies (παm ∈ Π)m=1,...,M and mean fields (µαm ∈ P(X )T )m=1,...,M . For the grid, we
Algorithm 1 Fixed point iteration 1: Initialize µ0 as the mean field induced by the uniformly random policy q. 2: for k = 0, 1, . . . do 3: Compute πk ∈ Π either directly by PPO on (14), or by computing Qµ via Algorithm 2 and
using (38) to obtain a softmax policy. 4: Compute µk+1 induced by πk using Algorithm 3, or for RL the neighborhood mean fields Gαt directly using Algorithm 4. 5: end for
Algorithm 2 Backwards induction 1: Input: Grid (αm ∈ [0, 1])m=1,...,M , mean field µ ∈M. 2: for m = 1, . . . ,M do 3: Initialize terminal condition Qµα(T, x, u) ≡ 0 for all (x, u) ∈ X × U . 4: for t = T − 1, . . . , 0 do 5: for (x, u) ∈ X × U do 6: Qµαm(t, x, u)← r(x, u,G αm t )+ ∑ x′∈X P (x
′ | x, u,Gαmt ) maxu′∈U Qµαm(t+1, x′, u′). 7: end for 8: end for 9: end for
10: Return (Qµαm)m=1,...,M
choose the points αm = m100 with m = 0, . . . , 100. Here, an agent α shall use the policy π αm with the closest αm.
To be precise, for the approximate mean field µ = Ψ̂(π) we define µα ≡ µ̂αm for the αm closest to α, i.e. formally, we thus have
Ψ̂(π) =
M∑
m=1
1α∈Ĩm µ̂ αm (31)
for any fixed policy ensemble π, with µ̂ defined through the recursive equation
µ̂αm0 ≡ µ0, µ̂αmt+1(x′) ≡ ∑
x∈X µ̂αmt (x)
∑ u∈U παmt (u | x)P (x′ | x, u, Ĝαmt ), m = 1, . . . ,M (32)
where under the assumption of equivalence classes Ĩm ≡ [am, bm] of size (bm − am), we obtain neighborhood mean fields via
Ĝαt = M∑
m=1
(bm − am)W (α, αm)µ̂αmt . (33)
Note that in our algorithms, we shall assume equisized partitions and use (bm−am) = 1M . Similarly, the policy ensemble is approximated by
Φ̂(µ) =
M∑
i=1
1α∈Ĩiπ αi (34)
where παi is the optimal policy of αi for any fixed µ, i.e. the optimal policy and mean field of each α is approximated by the optimal solution and mean field of the closest αi, which is an increasingly good approximation for sufficiently fine grids under the standing Lipschitz assumptions. In the case of block-wise Lipschitz continuous graphons via (20), a similar justification holds as long as each equivalence class remains constrained to one of the blocks of the graphon, see Theorem 5.
In Algorithm 1, the learning scheme is described on a high level. In our experiments, we either use approximate equivalence classes via Algorithms 2 and 3, or reinforcement learning in the form of PPO together with sequential Monte Carlo in Algorithm 4, though in principle one can mix arbitrary methods.
Algorithm 3 Forward simulation 1: Input: Grid (αm ∈ [0, 1])m=1,...,M , policy π ∈ Π. 2: Initialize starting condition µαm0 ≡ µ0 for all m = 1, . . . ,M . 3: for t = 0, . . . , T − 2 do 4: for m = 1, . . . ,M do 5: µαmt+1 ← ∑ x∈X µ αm t (x) ∑ u∈U π αm t (u | x)P (· | x, u, 1M ∑M n=1W (αm, αn)µ αn t ).
6: end for 7: end for 8: Return (µαm)m=1,...,M
Algorithm 4 Sequential Monte Carlo 1: Input: Number of trajectories K = 5, number of particles L = 200, policy π ∈ Π. 2: for k = 1, . . . ,K do 3: Initialize particles αm ∼ Unif([0, 1]), xm,k0 ∼ µ0 for all m = 1, . . . , L. 4: for t = 1, . . . , T − 1 do 5: for m = 1, . . . , L do 6: Sample action u ∼ παmt (· | xm,kt ). 7: Sample new particle state xm,kt+1 ∼ P (· | xm,kt , u, 1L ∑L n=1W (αm, αn)δxn,kt
). 8: end for 9: end for
10: end for 11: return neighborhood mean fields Gαt ≈ 1K ∑K k=1 1 L ∑L m=1W (α, αm)δxm,kt .
We ran each trial of our experiments on a single conventional CPU core, with typical wall-clock times reaching up to at most a few days. We estimate the required compute to approximately 6500 core hours. We did not use any GPUs or TPUs. More specifically, the training of our approximate equivalence class approach took on average approximately 24 hours for 250 iterations in SIS and 50 iterations in Investment. As a result, Figure 5 for the selection of appropriate temperatures took around 2500 core hours. The PPO experiments took approximately 3 days for each configuration, resulting in approximately 200 core hours for Figure 7. Finally, for the N -agent evaluations in Figure 4, each run up to 100 agents takes up to 4 core hours. Adding on top of that around 250 core hours for the rest of the experiments results in a total of approximately 4000 core hours.
For PPO, we used the RLlib implementation by Liang et al. (2018) (version 1.2.0, Apache-2.0 license). To allow for time-dependent policies, we append the current time to the network inputs. Further, discrete-valued observations are one-hot encoded. Any other parameter configurations are given in Algorithms 1, 2, 3 and 4, as well as in Table 2.
As for the specific configurations used in the PPO experiments, we give the hyperparameters in Table 1 and used with a feedforward neural network policy consisting of two hidden layers with 256 nodes and tanh activations, outputting a softmax policy over all actions.

heading: A.3.1 PROBLEM DEFINITIONS
text: For each possible problem setting, we list the applied temperature setting in Table 2. In the following, let G ∈ P(X ).
SIS-Graphon. In the SIS-Graphon game as described in the main text, we have X = {S, I}, U = {U,D}, µ0(I) = 0.5, r(x, u,G) = −2 · 1{I}(x) − 0.5 · 1{D}(u) and T = {0, . . . , 49}. Similar parameters produce similar results, and we set the transition probabilities as
P(S | I, ·, ·) = 0.2, P(I | S,U,G) = 0.8 ·G(I), P(I | S,D, ·) = 0 .
Investment-Graphon. Similarly, in the Investment-Graphon game we have X = {0, 1, . . . 9}, U = {I,O}, µ0(0) = 1, r(x, u,G) = 0.3x1+∑x′∈X x′G(x′) − 2 · 1{I}(u) and T = {0, . . . , 49}. We set
the transition probabilities for x = 0, 1, . . . , 8 as
P(x+ 1 | x, I, ·) = 9− x 10 ,
P(x | x, I, ·) = 1 + x 10 ,
P(x | x,O, ·) = 1, while for x = 9 the next state is always x = 9.

heading: A.3.2 EXPLOITABILITY AND TEMPERATURE CHOICE
text: In the following, we will explain our choice of temperatures in Table 2 by approximately evaluating the average exploitability of GMFE candidates (π,µ) – as it is intractable to approximately evaluate the maximum exploitability over all α ∈ I – defined by
∆J(π,µ) =
∫
I sup π∗∈Π
Jµα (π ∗)− Jµα (πα) dα . (35)
More specifically, when using approximate equivalence classes, we compute the exploitability of some policy π by computing the optimal policy π∗ obtained via Algorithm 2, under the fixed mean field µ generated by π via Algorithm 3, inserting π∗,α into (35) and then approximating by
∫
I Jµα (π) dα ≈
1
M
∑
m=1,...,M
∑ x∈X µ0(x) ∑ u∈U π0(u | x)Qµ,παm (0, x, u) . (36)
Here, we defined for any policy π ∈ Π and α ∈ I the policy evaluation functions Qµ,πα as usual via
Qµ,πα (t, x, u) = r(x, u,Gαt ) + ∑
x′∈X P (x′ | x, u,Gαt )
∑
u′∈U π0(u
′ | x)Qµ,πα (t+ 1, x′, u′) (37)
with terminal condition Qµ,πα (T, x, u) ≡ 0, which can be computed as in Algorithm 2, see also Puterman (2014) for a review.
To achieve convergence of fixed point iterations to approximate equilibria, for previous µn ∈M we compute the action value function Qµnα via Algorithm 2 using approximate equivalence classes
and then define the next policy πn+1 = Φ̂(µn) for every α ∈ I via the softmax function
πn+1,αit (u | x) = exp
( Qµ n
αi (t,x,u)
η
)
∑ u∈U exp
( Qµ
n αi (t,x,u)
η
) (38)
for the closest αi with some temperature η > 0 chosen minimally for convergence.
For choosing the temperature, we evaluate the approximate final exploitability at various temperatures. The results can be seen in Figure 5, where we plot the average, minimum and maximum exploitability over the last 10 iterations of the fixed point learning scheme. The reasoning behind choosing our temperatures as in Table 2 is that we can see no fluctuations (indicating convergence of our learning scheme) together with a low approximate exploitability at the indicated temperatures.

heading: A.3.3 ADDITIONAL EXPERIMENTS
text: In Figure 6, we plot investment behavior at quality x = 0 as well as expected quality for each α of the approximate equivalence class solution, and similarly in Figure 7 for the PPO solution with sequential Monte Carlo. Here, for each α we averaged quality over all particles within a distance of 0.05 to α. We can see that PPO achieves qualitatively and quantitatively similar behavior, deviating slightly due to the approximate optimality of the PPO algorithm. To be precise, when evaluating exploitability via either solution, we find that the learned policy exploitability remains around ε ≈ 2, compared to ε > 30 for the uniform random policy.
In Figure 8 the equilibrium behavior is shown for the Investment-Graphon problem without softmax policy regularization (except for the ER graphon case), as we find that the problem already converges to a very good equilibrium with low approximate exploitability, see Figure 5. In this problem, we find that the resulting (deterministic without regularization) policy will let agents invest up to a certain quality, after which any further investment is avoided. The agents with higher connectivity will invest up to a lower quality, as they are in competition with more products.
In Figure 9 and 10, we have performed ablations over the number of equivalence classes for the SISGraphon problem. As can be observed, the solution obtained by approximate equivalence classes remains stable regardless of the particular number of equivalence classes, showing the stability of discretization approach and supporting Theorem 5.
Finally, in Figure 11 we exemplarily show training results of applying state-of-the-art multi-agent reinforcement learning methods such as multi-agent PPO (MAPPO, Yu et al. (2021)) on the finiteagent system with observed, randomized-per-episode graphon indices and W -random graphs. Here, we use the same hyperparameters as shown in Table 1. As can be seen, due to the non-stationarity of the other agents, a naive application of MARL techniques fails to converge at all.

Review summaries:
This paper proposes a discrete-time graphon mean-field game (GMFG) framework as a limiting model to approximate large population dense graph games. Existence of graphon mean-field equilibrium (GMFE) is proved by reformulating GMFG as a classical MFG, and a graphon mean-field approximation result is also obtained. The authors then propose two learning algorithms for computing GMFE, one based on discretizing the graphon index, while the other based on the aforementioned classical MFG reformulation. Some preliminary numerical experiments are also provided to demonstrate the theoretical findings and to validate the proposed algorithms. 

This paper studies a class of games with a continuum of agents that connected by a graphon. This corresponds to the limit of games with a finite number of players connected by a graph. In the recent literature, such games have been considered in continuous time, and here the authors focus on a discrete time version. They first analyze the game by showing existence of a Nash equilibrium (defined in a suitable sense), and they prove that using the equilibrium policy from the game with a continuum of players in a game with a finite number of players provides an approximate Nash equilibrium. Then, two numerical approaches are proposed, both based on fixed point iterations that alternate between updating the distribution and updating the optimal control. The first method relies on backward induction to compute the optimal control of representative agents in a set of equivalence classes. The second method relies on reinforcement learning to compute the optimal control. Last, numerical examples are provided. 

This paper studies discrete time dense graph mean field games (GMFG). It discusses the settings of both finite-agent graph game and GMFG and shows that under mild Lipschitz continuous assumptions, the Nash equilibrium of the GMFG exists. It establishes that the NE of the GMFG is a good approximation of the NE of an N-agent game and empirically verifies this result. It then proposes two algorithms for solving the NE of GMFG, one with the idea of discretizing the graphon index, and the other utilizing the equivalence of GMFG to a classical MFG with an extended state space. The proposed algorithms are tested on an SIS-Graphon problem and an Investment-Graphon problem to validate their performances. 

This paper studies Markov games with large number of agents. The authors start by proposing a novel discrete-time formulation for graphon mean field games. They show that the mean field game exists, and that finite graph games converge to the mean field game when the number of agents goes to infinity. They also provide an algorithm that learns an approximate Nash equilibria.