Paper deep dive
The Dynamics of Policy Gradient in Social Dilemmas with Partner Selection
Benedict Russell, Chin-wing Leung, Paolo Turrini
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 93%
Last extracted: 7/8/2026, 5:14:11 PM
Summary
This paper provides an analytical solution for policy gradient dynamics in multi-agent social dilemmas with partner selection. It demonstrates that partner selection alters the opponent distribution and reward landscape to promote cooperation, proving that population variance is a necessary condition for cooperation to emerge. Using mean-field theory and a two-dimensional Wiener process, the authors derive sufficient conditions for cooperation and prove the existence of a stationary distribution, with simulations confirming the model's accuracy and the learning rate's impact.
Entities (10)
Relation Signals (8)
Population Variance → necessaryfor → Emergence of Cooperation
confidence 96% · we find that population variance is a necessary condition for cooperation to emerge
Partner Selection → modifies → Opponent Distribution
confidence 95% · We show how partner selection changes the opponent distribution and hence the reward landscape
Opponent Distribution → shapes → Reward Landscape
confidence 94% · changes the opponent distribution and hence the reward landscape
Out-for-Tat (OFT) → promotes → Cooperation
confidence 93% · prove this promotes cooperation under simple rules known from the literature.
Reverse Out-for-Tat (ROFT) → promotes → Cooperation
confidence 93% · prove this promotes cooperation under simple rules known from the literature.
Reward Landscape → influences → Policy Gradient Dynamics
confidence 92% · maps policy gradient (PG) updates onto the non-stationary reward landscapes induced by shifting opponent distributions
Learning Rate → affects → Emergence of Cooperation
confidence 91% · clarifies how the learning rate affects the emergence of cooperation
Two-Dimensional Wiener Process → captures → Stochastic Effects
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:In social dilemmas self-interested learning agents face the choice between the societal benefit of cooperation and the immediate reward of defection. Significant evidence exists on the benefits of assortment mechanisms such as partner selection for the emergence of cooperation, but this is largely available through agent-based simulations. In this paper, we provide an analytical solution to the problem, studying the policy-gradient dynamics in a multi-agent environment with partner selection. We show how partner selection changes the opponent distribution and hence the reward landscape, and prove this promotes cooperation under simple rules known from the literature. In particular, we find that population variance is a necessary condition for cooperation to emerge. Using a two-dimensional Wiener process, we extend the dynamics to capture the stochastic effects of partner selection and the resulting opponent distribution. We derive a sufficient condition for the population to be cooperation-promoting and prove the existence of a stationary distribution. Simulations confirm that the stochastic model accurately captures the policy-gradient dynamics and clarifies how the learning rate affects the emergence of cooperation.
Tags
Links
- Source: https://arxiv.org/abs/2605.18185v2
- Canonical: https://arxiv.org/abs/2605.18185v2
Trouble viewing inline? Open PDF directly →
Full Text
96,665 characters extracted from source content.
Expand or collapse full text
The Dynamics of Policy Gradient in Social Dilemmas with Partner Selection Benedict Russell Mathematics Institute Warwick University benedict.i.russell@warwick.ac.uk &Chin-wing Leung Department of Computer Science Warwick University chin-wing.leung@warwick.ac.uk &Paolo Turrini Department of Computer Science Warwick University p.turrini@warwick.ac.uk Abstract In social dilemmas self-interested learning agents face the choice between the societal benefit of cooperation and the immediate reward of defection. Significant evidence exists on the benefits of assortment mechanisms such as partner selection for the emergence of cooperation, but this is largely available through agent-based simulations. In this paper, we provide an analytical solution to the problem, studying the policy-gradient dynamics in a multi-agent environment with partner selection. We show how partner selection changes the opponent distribution and hence the reward landscape, and prove this promotes cooperation under simple rules known from the literature. In particular, we find that population variance is a necessary condition for cooperation to emerge. Using a two-dimensional Wiener process, we extend the dynamics to capture the stochastic effects of partner selection and the resulting opponent distribution. We derive a sufficient condition for the population to be cooperation-promoting and prove the existence of a stationary distribution. Simulations confirm that the stochastic model accurately captures the policy-gradient dynamics and clarifies how the learning rate affects the emergence of cooperation. 1 Introduction Social dilemmas capture the tension faced by self-interested agents when given the option to contribute to a common goal, or exploit the contributions of others. Ensuring that autonomous systems can work together for the greater good is currently one of the biggest challenges in the field of Artificial Intelligence. Recent research has explored the cooperation of LLMs [willis_will_2025, nguyen_navigating_nodate], multi-agent reinforcement learning (MARL) [leung_learning_2024], and evolutionary game theory [nowak_five_2006]; each providing unique insight into how cooperation can be sustained. In MARL, this problem is especially challenging, as agents experience a non-stationary environment induced by other learning agents [anastassacos_partner_2020]. Consequently, defection is a common outcome unless the agents are equipped with an additional mechanism such as behavioural signals of the opponent [priklopil_optional_2017, sabater-mir_reputation_2002] or dynamic interactions [bara_enabling_2022]. In his seminal work, nowak_five_2006 identified partner selection as a key mechanism for the emergence of cooperation. This has led to extensive research on cooperation on dynamic networks [bara_enabling_2022, wang_cooperation_2012, fehl_co-evolution_2011], as well as pairwise partner selection [anastassacos_partner_2020, leung_learning_2024]. Researchers in this area commonly find that partner selection rules which increase interaction frequency or duration between cooperative agents are key to the emergence of cooperation. A popular mechanism, coined Out-for-Tat (OFT), enforces cooperative pairs to stay together; whilst re-wiring pairs with a defector [izquierdo_leave_2014, zhang_opting_2016]. Recent contributions have found allowing agents to learn partner selection is sufficient for cooperation to emerge [defection_russell2026, fan_colearning]; the agents adopt OFT as a mechanism to avoid exploitation. When random matching is enforced, defection takes hold of the population. Despite the empirical success, the current understanding of partner selection remains heavily reliant on agent-based simulations [leung2026learning, leung_curiosity_2025]. While these simulations provide robust evidence that cooperation can emerge, they often fail to provide a formal theoretical justification for how these selection rules reshape the underlying learning dynamics. Specifically, there is a lack of analytical work that maps policy gradient (PG) updates onto the non-stationary reward landscapes induced by shifting opponent distributions. Without this theoretical foundation, it remains unclear which population conditions are strictly necessary for cooperation to take hold, or how stochasticity in action and partner selection influences long-term stability. Contribution In this paper, we provide a formal theoretical analysis of policy gradient dynamics in social dilemmas with partner selection. Using mean-field theory, we derive the conditional partner distribution to capture how partner selection rules alter the reward structure for learning agents. We prove that partner selection promotes cooperation under well-known rules (OFT and ROFT) and identify that population variance is a necessary condition for cooperation to emerge from an initially non-cooperative state. We extend the mean dynamics to a stochastic model using a two-dimensional Wiener process to capture the randomness inherent in action selection and partner selection. We derive a sufficient condition for the population to be cooperation-promoting and prove the existence of a stationary distribution under the derived model. Simulations corroborate our theoretical findings, demonstrating that our stochastic model accurately captures the evolution of strategies and clarifies the pivotal role of the learning rate in fostering cooperative clusters. Related Literature Partner selection. In evolutionary and network models, enabling cooperators to avoid defectors alters the payoff landscape, generating more assortative interaction patterns than a well-mixed population [pacheco_active_2006, santos_cooperation_2006, pacheco_coevolution_2006]. This can take many forms, including active linking [pacheco_active_2006, pacheco_repeated_2008, bara_enabling_2022], partner switching [fu_partner_2009, li_aspiration-based_2014], and conditional dissociation [izquierdo_leave_2014, izquierdo_option_2010]. Rewiring acts as a weak form of ostracism; defectors cannot repeatedly exploit cooperators, but cooperators have continued access to each other. These results carry over to human populations, where experimental work shows that dynamic partner updates can increase cooperation [rand_dynamic_2011, fehl_co-evolution_2011, zhang_opting_2016]. Recent research in MARL has implemented these ideas into learning environments. Anastassacos et al. [anastassacos_partner_2020] show that partner selection can induce cooperation amongst selfish reinforcement learning agents. This has been extended to enable agents to learn partner selection rules themselves, with behavioural signals [leung_learning_2024, leung_curiosity_2025] and without any prior information on opponents [defection_russell2026]. Learning Dynamics. The complexity associated with non-stationary environments has led researchers to model learning dynamics using tools from dynamical systems and evolutionary game theory. Prior work has linked learning to coupled replicator equations [sato_coupled_2003, galstyan_continuous_2013], selection-mutation models [tuyls_selection-mutation_2003, bloembergen_evolutionary_2015], and mean field approximations [Hu2019ModellingTD]. Learning dynamics have been considered in population games, using mean-field theory to derive the dynamics under the concurrent learning protocol [hu_dynamics_2022, Hu2019ModellingTD]. This has been extended to agents acting on graphs [chu_formal_2022], and considering stochastic interactions under the social learning protocol [leung_modelling_2022]. In a closely related work, zheng_simple_2017 uses a pair approximation in an evolutionary setting to show the emergence of cooperation under a partner selection rule; our work differs by analysing a learning population. Policy Gradient Dynamics. In 2-player symmetric games, srinivasan_actor-critic_nodate provides the mean dynamics of policy gradient and compares against the replicator dynamics. bernasconi_evolutionary_2025 formally establishes the connection between the replicator dynamics and the soft-max policy gradient, and proves asymptotic stability for evolutionarily stable strategies (ESS) in symmetric games. This setting has been extended to analyse the stochastic effects of exploration [leung_stochastic_evolution_2024], with stochastic stability shown under self-play. Stochastic stability has also been analysed when the payoff faces aggregate shocks [fudenberg_evolutionary_1992] and random perturbations [mertikopoulos_emergence_2010]. Whilst there has been extensive work on symmetric games, including the stochastic effects of shocks and action selection, there has not been theoretical research which shows the role of partner selection on the behaviour of reinforcement learning (RL) agents. This adds an additional layer of complexity, where the opponent distribution will influence both the mean and stochastic dynamics. 2 Background In this section, we provide the necessary background on social dilemmas with partner selection and the policy gradient algorithm in repeated Prisoner’s Dilemma. 2.1 Repeated Prisoner’s Dilemma In the Prisoner’s Dilemma (PD), two players simultaneously choose whether to cooperate (C) or defect (D), receiving payoffs based on their combined actions. Whilst mutual cooperation yields the greatest collective payoff, each player is incentivised to defect and exploit their opponent. Consequently, mutual defection remains the only Nash equilibrium [nash_non-cooperative_1951] and ESS in the one-shot game [smith_evolution_1982]. The general payoff matrix is shown in Table 1; the payoff is such that b>c>0b>c>0. Table 1: Prisoner’s Dilemma payoff matrix. C D C b 0 D b+cb+c c This paper examines a repeated Prisoner’s Dilemma, where for each episode, an agent plays H rounds of PD with their opponents. The opponent changes based on the actions of the pair, according to a predefined partner selection rule. We will consider 4 commonly studied rules [leung_learning_2024, defection_russell2026]: – Out-for-Tat (OFT): stay if and only if both cooperate – Reverse Out-for-Tat (ROFT): stay if and only if both defect – Always Stay (Stay): stay independent of the actions – Always Switch (Switch): switch independent of the actions These provide a baseline to analyse how partner selection can be used to promote cooperation in a population of self-interested learning agents. 2.2 Policy Gradient We study independent reinforcement learners where each agent aims at optimising their expected return. The policy gradient (PG) algorithm searches over the parameter space ∈ℝm ψ ^m and updates the policy through stochastic gradient ascent to find parameters that optimise the agent’s decision-making. We consider softmax parameterisation for each agent’s policy, which enables balancing exploration and exploitation. The policy for action aka_k is given by π(ak|(t)):=exp(ψk)∑l=1mexp(ψl).π(a_k| ψ(t)):= ( _k) _l=1^m ( _l). (1) The policy update follows the REINFORCE algorithm [williams_simple_1992]. If action aha_h is selected in round h of an episode and the resulting accumulated rewards from step h onwards are RahhR^h_a_h, where we assume the discount factor is 11, the approximated gradient of the score-function is (Rahh−β)(ah−)(R^h_a_h-β)( e_a_h- π), where β∈ℝβ is a baseline [Sutton2018ReinforcementL] and we set it to 0. The parameter update across an episode of length H is ∂t:=(t+1)−(t)=α∑h=0H−1Rahh(ah−). _t ψ:= ψ(t+1)- ψ(t)=α _h=0^H-1R^h_a_h( e_a_h- π). (2) where aha_h is the action sampled at round h, and ah=(0,…1,…0) e_a_h=(0,...1,...0) is the basis vector. In the PD, the actions are limited to a∈C,Da∈\C,D\. Denoting the probability of cooperating x:=π(C|(t))x:=π(C| ψ(t)) and defecting 1−x:=π(D|(t))1-x:=π(D| ψ(t)), the parameter mean dynamics are given by ∂tψC=αx(1−x)∑h=0H−1(GCh−GDh),∂tψD=−∂tψC, _t _C=α x(1-x) _h=0^H-1(G^h_C-G^h_D), _t _D=- _t _C, where GCh:=[Rh|ah=C]G^h_C:=E[R^h|a_h=C] and GDh:=[Rh|ah=D]G^h_D:=E[R^h|a_h=D]. Applying the chain rule to Equation (1), the policy mean dynamics are dxdt dxdt =2αx2(1−x)2∑h=0H−1(GCh−GDh), =2α x^2(1-x)^2 _h=0^H-1(G^h_C-G^h_D), (3) which is a one-dimensional non-linear ODE, where the sign of the episodic reward difference between cooperating and defecting will determine the evolution. Defining Ga=∑h=0H−1GahG_a= _h=0^H-1G_a^h, we denote the overall difference by ΔG=GC−GD G=G_C-G_D. 3 Model Formulation In this section, we introduce the agent-based model and derive how partner selection alters the opponent distribution and resulting reward structure. 3.1 Setup We consider a population of agents, where the distribution of strategies is denoted by ρ(x)ρ(x). Each strategy, x, corresponds to the probability of cooperating by that agent, and as such is bounded between 0 and 1. Therefore, if x=0x=0 the agent always defects, whilst if x=1x=1, the agent will always cooperate. For each episode, an agent i is randomly selected to be the focal agent. At round h=0h=0, sample an opponent j from the population; the pair interact by playing the PD game, and the payoff is given according to Table 1. The opponent j is then redrawn according to the partner selection rule. This repeats for a fixed number of rounds, H. At the end of the episode, the focal agent i updates their strategy according to the REINFORCE algorithm. The partner selection rule chosen will play a significant role in the evolutionary dynamics. For example, under the Out-For-Tat mechanism, an agent will continue to play with their partner if and only if they both cooperate. Otherwise, a new opponent is sampled (with replacement) from the population. The distribution of these opponents is inherently different: if the opponent stayed, they are more likely to have a higher cooperation rate x than a new draw from the population. 3.2 Conditional Opponent Distribution & Reward Structure Consider an arbitrary agent i in the population, where the strategy of this agent is denoted by x. At step h of the episode, the distribution of the opponent given agent i’s strategy is ρh(y|x)ρ^h(y|x). For any partner selection rule, the distribution at the next step is ρh+1(y|x) ρ^h+1(y|x) =∫01ρh+1(y|x,z)ρh(z|x)z = _0^1ρ^h+1(y|x,z)ρ^h(z|x)dz The one-step change between the distributions is subject to the chosen Markovian partner selection mechanism. For example, under the OFT rule, the pair will stay if and only if both agents cooperate. Hence, ρh+1(y|x) ρ^h+1(y|x) =∫01(δz(y)xz+ρh+1(y|x,switch)(1−xz))ρh(z|x)z = _0^1 ( _z(y)xz+ρ^h+1(y|x,switch)(1-xz) )ρ^h(z|x)dz =x∫01δz(y)zρh(z|x)z+ρh+1(y|x,switch)∫01(1−xz)ρh(z|x)z =x _0^1 _z(y)zρ^h(z|x)dz+ρ^h+1(y|x,switch) _0^1(1-xz)ρ^h(z|x)dz =xyρh(y|x)+ρ(y)(1−xmh(x)) =xyρ^h(y|x)+ρ(y)(1-xm^h(x)) where, upon breaking a connection, a new opponent is sampled from the population. Given the focal agent’s strategy x, the mean cooperation frequency of the opponents at step h is denoted by mh(x)=h[Y|x] m^h(x)=E^h[Y|x] =∫01yρh(y|x)y. = _0^1yρ^h(y|x)dy. For the 4 partner selection rules, the recursion is given in Table 2. The opponent distribution provides an explicit way to capture the reward at each step, h. We can utilise this structure to derive the conditional rewards across the episode. Table 2: Conditional partner distribution update step and per round reward difference for episode length H=2H=2. In each case, the initial distribution is ρ(y)ρ(y). Rule ρh+1(y|x)ρ^h+1(y|x) Δrh,h(x) r^h,h(x) Δr0,1(x) r^0,1(x) OFT xyρh(y|x)+ρ(y)(1−xmh(x))xyρ^h(y|x)+ρ(y)(1-xm^h(x)) −c-c b(μ2−μ12)b( _2- _1^2) ROFT (1−x)(1−y)ρh(y|x)+ρ(y)(1−(1−x)(1−mh(x)))(1-x)(1-y)ρ^h(y|x)+ρ(y)(1-(1-x)(1-m^h(x))) −c-c b(μ2−μ12)b( _2- _1^2) Stay ρ(y)ρ(y) −c-c 0 Switch ρ(y)ρ(y) −c-c 0 Following Equation (2), the episodic reward is the sum of all accumulated rewards conditioned on an action adopted at step h=kh=k, followed by the original policy x in future interactions. Consequently, we denote the agent’s strategy, conditioned on its action in round k, as x,Ckx,C^k and x,Dkx,D^k, for cooperating and defecting at k respectively. Likewise, the opponent distribution, expected cooperation, and reward at time h are denoted ρh(y|x,Ck)ρ^h(y|x,C^k), mCk,h(x)m_C^k,h(x), and rCk,hr_C^k,h. At step h=kh=k, the expected difference between cooperation and defection is always −c-c. The reward difference for h≥k+1h≥ k+1 is given by Δrk,h(x)=rCk,h−rDk,h r^k,h(x)=r_C^k,h-r_D^k,h =∫[by+(1−x)c](ρh(y|x,Ck)−ρh(y|x,Dk))y = [by+(1-x)c](ρ^h(y|x,C^k)-ρ^h(y|x,D^k))dy =b∫y(ρh(y|x,Ck)−ρh(y|x,Dk))y =b y(ρ^h(y|x,C^k)-ρ^h(y|x,D^k))dy =b(mCk,h(x)−mDk,h(x)) =b(m^k,h_C(x)-m^k,h_D(x)) =bΔmk,h(x) =b m^k,h(x) where given an agent has policy x and cooperated in the kkth round, mCk,h(x)m^k,h_C(x) is explicitly given by mCk,h(x)=h[Y|x,Ck] m^k,h_C(x)=E^h[Y|x,C^k] =∫01yρh(y|x,Ck)y. = _0^1yρ^h(y|x,C^k)dy. (4) Let μl=∫ylρ(y)y _l= y^lρ(y)dy denote the llth moment of the underlying population distribution ρ(y)ρ(y); the reward differences for the first two steps in the episode are shown in Table 2 (a full derivation is provided in Appendix A). Regardless of the partner selection rule, the expected reward difference between cooperation and defection in the kkth round is −c-c. Therefore, for cooperation to emerge the future rewards must compensate by pairing the focal agent with more cooperative agents which can build sustained partnerships. Interestingly, both OFT and ROFT provided the same expected reward across the first two rounds: keeping defective pairs together is as effective as keeping cooperative ones. We will now provide a more general result for when H>2H>2. The conditional reward difference for Always stay and Always switch is zero for all steps h≥k+1h≥ k+1. Here, we aim to show that the reward difference for all time steps h≥k+1h≥ k+1 is non-negative for the OFT and ROFT mechanisms. Proposition 3.1. For the OFT and ROFT partner selection rules, the reward difference Δrk,h(x)≥0 r^k,h(x)≥ 0 for all x∈[0,1]x∈[0,1] and h≥k+1h≥ k+1. Proof. All proofs are provided in Appendix B. ∎ Corollary 3.1. Under the OFT and ROFT partner selection rules, if ΔG[ρ]>0 G[ρ]>0 for H=2H=2 then ΔG[ρ]>0 G[ρ]>0 for all H≥2H≥ 2. The above result allows a much simpler reward structure to be studied, associated with the partner selection mechanisms. In particular, increasing the episodic length can only increase the marginal accumulated rewards for cooperation. Therefore, the results in Section 3.3 provide sufficient conditions for the emergence of cooperation for all H≥2H≥ 2. 3.3 Evolution of the Mean Dynamics In this section, we show how the above reward structure can be embedded within the policy gradient update. We then prove that always stay and always switch converge to the pure defection equilibrium, regardless of the initial distribution. Finally, we show that both OFT and R-OFT necessarily increase cooperation when the initial population has sufficient variance. Under the random matching and always stay mechanisms, the reward difference simplifies to −Hc-Hc for any episode length H. As such, when H=2H=2 the policy update is given by dxdt dxdt =−4αcx2(1−x)2 =-4α cx^2(1-x)^2 which generates the mean-field continuity PDE ∂tρ+∂x(−4αx2(1−x)2ρ) _tρ+ _x(-4α x^2(1-x)^2ρ) =0. =0. Theorem 3.1. Let ρ0∈L1(0,1),ρ0≥0 _0∈ L^1(0,1), _0≥ 0, and ∫01ρ0=1 _0^1 _0=1. Under the Stay and Switch rules, the population converges to pure defection. Consequently, we know that always-switch and always-stay will, under the mean-dynamics, cause mass defection to take over the population. Therefore, we seek to show that partner selection rules such as OFT and R-OFT can induce a more cooperative population. Under the OFT mechanism, the continuity equation is ∂tρ+∂x(u[ρ](x)ρ) _tρ+ _x(u[ρ](x)ρ) =0, =0, u[ρ](x)=2αΔG[ρ]x2(1−x)2, u[ρ](x)=2α G[ρ]x^2(1-x)^2, where for H=2H=2, we have ΔG[ρ] G[ρ] =∑h=0H−1ΔGh[ρ] = _h=0^H-1 G^h[ρ] =Δr0,0(x)+Δr0,1(x)+Δr1,1(x) = r^0,0(x)+ r^0,1(x)+ r^1,1(x) =−c+b(μ2−μ12)−c =-c+b( _2- _1^2)-c =bVar(ρ(y))−2c. =bVar(ρ(y))-2c. By Corollary 3.1, for H>2H>2 this value is a lower bound for the sign of the drift. As a consequence, if cooperation emerges for H=2H=2, it provides a sufficient condition for longer episode lengths. Denote ([0,1])P([0,1]) as the set of all probability measures on the domain. Proposition 3.2. Let ρ0∈([0,1]) _0 ([0,1]), then there exists a unique solution to the continuity equation. Furthermore, the unique solution can be expressed as ρ(t)=(Xt)#ρ0 ρ(t)=(X_t)_\# _0 where (Xt)#ρ0(X_t)_\# _0 is the push-forward of ρ0 _0 by the characteristic flow. Define the time integral of the non-local velocity induced by the episodic reward, K(t)=∫0tΔG[ρ(s)]s K(t)= _0^t G[ρ(s)]\;ds Then the continuity PDE depends on time through K, and hence satisfies ρ(t,⋅)=(XK(t))#ρ0.ρ(t,·)=(X_K(t))_\# _0. (5) See Appendix B for details. This enables the flow of the mass to be uniquely determined by the behaviour of K. Theorem 3.2. Let ρ0∈L1(0,1),ρ0≥0 _0∈ L^1(0,1), _0≥ 0, and ∫01ρ0=1 _0^1 _0=1. If ΔG[ρ0]>0 G[ _0]>0, then K(t)↗K∗K(t) K^* where K∗∈ℝK^* and ρ(t)⇀(XK∗)#ρ0ρ(t) (X_K^*)_\# _0. That is, if the initial population variance is high enough, the limiting distribution is shifted towards higher cooperation. We highlight that a finite K∗K^* is attained, which means XK∗(x)<1X_K^*(x)<1. This enforces that the population does not converge to a single strategy, and hence does not converge to pure cooperation. The positivity of the reward difference ΔG[ρ0] G[ _0] constrains the payoff; the variance on the domain [0,1][0,1] is bounded above by 14 14, meaning b>4cb>4c is a necessary condition when H=2H=2. This constraint relaxes as H increases, as the long term benefits of OFT and ROFT increase the reward difference. 4 Stochastic Dynamics To capture the stochasticity induced by action exploration and the opponent distribution, we model the episodic REINFORCE update as a random increment in the parameter space, ψ. To do this we follow the approach in [kifer_random_1988, leung_stochastic_evolution_2024], characterising the random increment by the first two moments. In particular, we approximate the discrete stochastic updates by a continuous-time diffusion process: d=dt+Σdt, d ψ= μdt+ \;d W_t, (6) where μ is the expected update, t W_t is a two-dimensional standard Wiener process, and Σ:=Σ(,t) := ( ψ,t) denotes the covariance matrix. In 2-action repeated games, we can exploit the logit difference to simplify the stochastic dynamics. Applying Ito’s lemma [ito_stochastic_1944], we derive the policy evolution as dx=[2αx2(1−x)2ΔG[ρ]+2x(1−x)(1−2x)ΣCC]dt+2x(1−x)ΣCCdWt dx=[2α x^2(1-x)^2 G[ρ]+2x(1-x)(1-2x) _C]\;dt+2x(1-x) _C\;dW_t (7) where ΣCC _C denotes the variance of the parameter update. Denote the hhth update step of the REINFORCE algorithm as Uh:=(Rahh−β)(1ah=C−x)U^h:=(R_a_h^h-β)(1\a_h=C\-x). The covariance is ΣCC _C :=α2∑h=0H−1(x(1−x)2SCh+x2(1−x)SDh−x2(1−x)2(ΔGh[ρ])2)+2α2∑h<lCov(Uh,Ul), :=α^2 _h=0^H-1 (x(1-x)^2S^h_C+x^2(1-x)S^h_D-x^2(1-x)^2( G^h[ρ])^2 )+2α^2 _h<lCov(U^h,U^l), where SCh=[(Rh−β)2|ah=C]S^h_C=E[(R^h-β)^2|a_h=C] and SDh=[(Rh−β)2|ah=D]S^h_D=E[(R^h-β)^2|a_h=D] are the second moments of the reward conditioned on actions C and D at round h. For any action a, this is given by Sah S^h_a =∑l=hH−1Var(rl|ah=a)+2∑h≤l<kCov(rl,rk|ah=a)+(Gah−β)2. = _l=h^H-1Var(r_l|a_h=a)+2 _h≤ l<kCov(r_l,r_k|a_h=a)+(G^h_a-β)^2. Details of the variance and covariance terms are provided in Appendix A. 4.1 Population Evolution By propagating the stochastic updates across the population, we can study the strategy evolution at a population scale. The Fokker-Planck equation (FPE) describes the evolution of the probability density function as a diffusion process [risken_fokker-planck_1996]. For the strategy space x, the FPE corresponding to (7) is ∂tρ(t,x)=−∂x[Aρ(t,x)ρ(t,x)]+12∂xx[Bρ2(t,x)ρ(t,x)] _tρ(t,x)=- _x[A_ρ(t,x)ρ(t,x)]+ 12 _x[B^2_ρ(t,x)ρ(t,x)] (8) where Aρ(t,x) A_ρ(t,x) =2αx2(1−x)2ΔG[ρ]+2x(1−x)(1−2x)ΣCC,Bρ2(t,x)=4x2(1−x)2ΣCC. =2α x^2(1-x)^2 G[ρ]+2x(1-x)(1-2x) _C, B^2_ρ(t,x)=4x^2(1-x)^2 _C. Given an initial strategy distribution ρ(0,x)=ρ0(x)ρ(0,x)= _0(x), the time evolution can be solved numerically. These numerical solutions will enable us to analyse the long-term dynamics; for the short-term we can provide an equivalent sufficient condition to the mean dynamics for an increasing level of cooperation. In particular, the next result shows that given a sufficiently low learning rate, the cooperation levels in the population will initially increase for sufficiently high variance in the initial population. Proposition 4.1. Suppose ρ0 _0 has non-zero interior mass and ΔG[ρ0]>0 G[ _0]>0. Then there exists a T∈ℝ+T ^+ and α∗α^* such that for all α<α∗α<α^*, the mean level of cooperation is increasing on [0,T][0,T]. It does not, however, guarantee that a stable cooperative equilibrium will be reached. When H=2H=2, the population variance remaining sufficiently high would prove a sufficient condition, but the evolution of the variance depends upon higher moments. To analyse the equilibrium which is attained in the stochastic case, we prove the limiting distribution is well-posed and turn to numerical solutions to analyse the behaviour. 4.2 Stationary Distribution The stationary distribution occurs at the point such that the time derivative in the PDE is zero. Applying no-flux boundary conditions, this condition then reduces to finding a solution to the ODE given by ∂x[Bρ2(x)ρ(x)] _x[B_ρ^2(x)ρ(x)] =2Aρ(x)ρ(x). =2A_ρ(x)ρ(x). (9) A major obstacle in proving such an equation is well-defined is the behaviour at the boundary. At x=0x=0 and x=1x=1, the noise term B2(x)=0B^2(x)=0 which can prevent using standard techniques. We first consider an ε− -regularised version, then extend the existence as ε→0 → 0. Fix η∈L∞([0,1])η∈ L^∞([0,1]) as a trial density such that Aη(x)A_η(x) and Bη2(x)B^2_η(x) are a function of x and η only. The regularised versions are Aηε(x)=Aη(x)A_η (x)=A_η(x) and Bη2,ε(x)=Bη2(x)+εB^2, _η(x)=B^2_η(x)+ . Define the set Sε=η∈C([0,1]):η≥0,∫01η(x)x=1,‖η‖∞≤Mε S =\η∈ C([0,1]):η≥ 0, _0^1η(x)dx=1,\|η\|_∞≤ M \ For a given η∈Sη∈ S, the unique solution ℱε[η](x)F [η](x) to the linear ODE is given by ℱε[η](x) [η](x) =wηε(x)∫01wηε(x)x,wηε(x)=1Bη2,ε(x)exp(∫0x2Aηε(y)Bη2,ε(y)y) = w _η(x) _0^1w _η(x)dx, w _η(x)= 1B^2, _η(x) ( _0^x 2A _η(y)B_η^2, (y)dy) By updating the value of η using the operator ℱF, we aim to converge to a fixed point of the system, which will correspond to a stationary distribution of the PDE in (8). The required regularity of the operator and space is proven in Appendix B. Proposition 4.2. Fix ε>0 >0, then there exists at least one fixed point of the mapping ℱε[ρ]=ρF [ρ]=ρ. That is, there is at least one solution to the regularised steady-state equations. Remark 4.1. One can consider a PDE with additional noise not captured by the first two moments as an independent stochastic fluctuation of size σ. The corresponding stationary distribution would be equivalent to the regularisation above, and therefore existence would follow. Moreover, uniqueness would follow for sufficiently high σ by analysing the Lipschitz constant of the operator, ℱF. We can now extend this to find an unregularised solution by sending the parameter ε→0 → 0. Theorem 4.1. There exists at least one stationary probability measure ρ∈([0,1])ρ ([0,1]) which solves (9) in the weak sense. The result implies that the steady state solutions can contain Dirac measures at the boundary. Consequently, we would expect in the long run for mass to accumulate in clusters of pure defection and cooperation. Moreover, the steady state solution is non-unique, explaining how the dynamics are strongly influenced by the initial population and parameter regimes. 5 Experiments In this section, we show how the theoretical analysis and model capture the dynamics displayed in the agent-based simulations. To do so, we solve Equation (8) for the 4 different partner selection rules. 5.1 Game Setting We use a finite volume method over the domain [0,1][0,1] to capture the time and spatial evolution. The agent-based simulations act as a ground truth in the experiments; we conduct 30 simulations of a population of 1,000 agents training over E=E=5e6 episodes. The population distribution at time t is obtained by averaging the simulations. For comparison, the simulation time when plotting is scaled such that t=E/Nt=E/N. The initial strategy values are sampled from a specified distribution ρ0 _0, with the parameter (0) ψ(0) obtained through the inverse of the softmax function (1). The payoff is given by Table 1, where b,c=3,0.1b,c=3,0.1 respectively. Unless specified, the parameters are α=0.01,H=2α=0.01,H=2. Evolution under a wide range of initial conditions is shown in Appendix C. 5.2 Result We compare the distribution of strategies in the population over time. Figure 1 presents the results of both the theoretical solution (solid line) and the simulations (histogram). We clearly see that the FPE derived captures the distribution of strategies as it evolves through time. This close alignment is present for all 4 partner selection rules, where the behaviour of the two cooperation-inducing rules (OFT and ROFT) displays very similar dynamics which lead to the emergence of a defecting and cooperating cluster. On the other hand, the rules Stay and Switch induce defection to take over the population. Figure 1: Evolution of the strategy distribution where the population is initialised with X∼β(2,2)X β(2,2) for the 4 rules. The theoretical solution (solid line) matches the simulations (histogram). As highlighted in Section 3.3, the underlying distribution requires sufficient variance for cooperation to emerge under the mean dynamics. This begs the question of how and why cooperation can emerge when the population is initialised to a singular policy [leung_learning_2024, defection_russell2026]. The answer lies in how the learning rate can induce population variance. With low learning rates, the stochastic effects are diminished and therefore policies are updated close to the underlying expectation. A higher value of α enables sampled trajectories to influence the dynamics and therefore induce a faster-growing population variance. Figure 2 shows how for an initially unbiased population, a cooperation cluster can be supported through only increasing the learning rate. Whilst higher learning rate can induce greater variance and support cooperation, the FPE approximation loses fidelity. Figure 2: Evolution of the strategy distribution under OFT where the population is initialised with a Dirac at 0.5. The effect of the learning rate is clear: higher values induce additional variance in the population which in turn induces cooperation. Time has been scaled with the learning rate for comparison. 6 Conclusion In this paper, we study the learning dynamics of policy gradient in a multi-agent environment when partner selection rules are enforced. We show how the opponent distribution alters the reward landscape, and prove that the underlying dynamics promote cooperation for both the OFT and ROFT partner selection rules. In particular, we find that population variance is a fundamental requirement for the emergence of cooperation under partner selection. We extend the model to capture the stochastic dynamics using a 2-dimensional Wiener process which models the stochasticity induced by action and opponent selection. In this setting, we show how the same variance requirement is sufficient for an initial increase in cooperation with a sufficiently small learning rate, before proving the existence of a stationary distribution. This analysis revealed that the structure of the long-term dynamics constitutes mass at the pure defection and pure cooperation strategies. In the experiments, we demonstrate that the stochastic model accurately describes the policy gradient dynamics and can be used to understand the role of the learning rate in inducing cooperation. Future directions include extending the model to other learning algorithms and considering a synchronous population partner selection model. References Appendix A Derivations A.1 Conditional Partner Distribution and Rewards For additional clarity, we provide the explicit iterative form for the first two steps of the episode with the 4 partner selection rules. All steps are computed using the iterative formulas provided in Table 2. For each, the initial opponent distribution is set as ρ(y)ρ(y). Table 3: Conditional partner distribution and reward difference for OFT h ρh(y|x,C0)ρ^h(y|x,C^0) ρh(y|x,D0)ρ^h(y|x,D^0) mC0,h(x)m^0,h_C(x) mD0,h(x)m^0,h_D(x) Δr0,h(x) r^0,h(x) 0 ρ(y)ρ(y) ρ(y)ρ(y) μ1 _1 μ1 _1 −c-c 1 yρ(y)+(1−μ1)ρ(y)yρ(y)+(1- _1)ρ(y) ρ(y)ρ(y) μ2−μ12+μ1 _2- _1^2+ _1 μ1 _1 b(μ2−μ12)( _2- _1^2) Table 4: Conditional partner distribution and reward difference for ROFT h ρh(y|x,C0)ρ^h(y|x,C^0) ρh(y|x,D0)ρ^h(y|x,D^0) mC0,h(x)m^0,h_C(x) mD0,h(x)m^0,h_D(x) Δr0,h(x) r^0,h(x) 0 ρ(y)ρ(y) ρ(y)ρ(y) μ1 _1 μ1 _1 −c-c 1 ρ(y)ρ(y) (1+μ1−y)ρ(y)(1+ _1-y)ρ(y) μ1 _1 μ1+μ12−μ2 _1+ _1^2- _2 b(μ2−μ12)( _2- _1^2) Table 5: Conditional partner distribution and reward difference for Stay and Switch h ρh(y|x,C0)ρ^h(y|x,C^0) ρh(y|x,D0)ρ^h(y|x,D^0) mC0,h(x)m^0,h_C(x) mD0,h(x)m^0,h_D(x) Δr0,h(x) r^0,h(x) 0 ρ(y)ρ(y) ρ(y)ρ(y) μ1 _1 μ1 _1 −c-c ≥1≥ 1 ρ(y)ρ(y) ρ(y)ρ(y) μ1 _1 μ1 _1 0 When H=2H=2, the only remaining reward difference to calculate is Δr1,1(x)=−c r^1,1(x)=-c for any partner selection rule. For larger values of H, the rewards depend on higher-order moments and can be calculated using the iterative formulas. Taking the sum across these rewards gives ΔG G for each rule. A.2 Stochastic Dynamics in 2-action Symmetric Games Recall the update to the parameter is ΔψC=α∑h=0H−1Uh, _C=α _h=0^H-1U^h, where Uh=(Rahh−β)(1ah=C−x)U^h=(R^h_a_h-β)(1\a_h=C\-x). Define the mean vector μ(,t):=[Δ(t)], μ( ψ,t):=E[ ψ(t)], and the covariance matrix by Σ(,t) ( ψ,t) =[(Δ(t)−μ)(Δ(t)−μ)T]. =E[( ψ(t)-μ)( ψ(t)-μ)^T]. In the two-action case, since ΔψC=−ΔψD _C=- _D, we have μ=(μC,−μC),Σ=(ΣCC−ΣCC−ΣCCΣCC) μ=( _C,- _C), = pmatrix _C&- _C\\ - _C& _C pmatrix where μC=[ΔψC]=αx(1−x)ΔG[ρ] _C=E[ _C]=α x(1-x) G[ρ], and ΔG[ρ]=GC−GD G[ρ]=G_C-G_D. To compute the covariance, note that ΣCC _C =Var(ΔψC)=α2∑h=0H−1Var(Uh)+2α2∑0≤h<k≤H−1Cov(Uh,Uk). =Var( _C)=α^2 _h=0^H-1Var(U^h)+2α^2 _0≤ h<k≤ H-1Cov(U^h,U^k). For each round h, the variance term is Var(Uh) (U^h) =x(1−x)2SCh+x2(1−x)SDh−x2(1−x)2(ΔGh)2 =x(1-x)^2S^h_C+x^2(1-x)S^h_D-x^2(1-x)^2( G^h)^2 where Sah=[(Rh−β)2|ah=a]S^h_a=E[(R^h-β)^2|a_h=a] and ΔGh=GCh−GDh G^h=G^h_C-G^h_D. For h<kh<k, define the conditional moment Ma,bh,k M^h,k_a,b =[(Rh−β)(Rk−β)|ah=a,ak=b], =E[(R^h-β)(R^k-β)|a_h=a,a_k=b], then each covariance term is Cov(Uh,Uk) (U^h,U^k) =x2(1−x)2(MC,Ch,k−MD,Ch,k−MC,Dh,k+MD,Dh,k−ΔGhΔGk). =x^2(1-x)^2(M^h,k_C,C-M^h,k_D,C-M^h,k_C,D+M^h,k_D,D- G^h G^k). We then approximate the discrete stochastic updates with the continuous time diffusion process as defined in Equation (6). To simplify the equations, consider the logit difference given by z=ψC−ψDz= _C- _D. Then Δz=ΔψC−ΔψD=2ΔψC z= _C- _D=2 _C. The corresponding one-dimensional parameter dynamics are dz=μzdt+σzdWt dz= _zdt+ _zdW_t (10) where μz=2αx(1−x)ΔG[ρ],σz2=4ΣCC. _z=2α x(1-x) G[ρ], _z^2=4 _C. The cooperation probability of an agent is then computed with the softmax function x(z)=11+e−z. x(z)= 11+e^-z. Applying Ito’s lemma with x′(z)=x(1−x),x′(z)=x(1−x)(1−2x) x (z)=x(1-x), x (z)=x(1-x)(1-2x) we obtain dx dx =[x(1−x)μz+12x(1−x)(1−2x)σz2]dt+x(1−x)σzdWt. =[x(1-x) _z+ 12x(1-x)(1-2x)σ^2_z]\;dt+x(1-x) _z\;dW_t. Substituting in the formulas of μz _z and σz _z gives the provided equation. A.3 Conditional Moments For action a∈C,Da∈\C,D\, the second moment of the episodic reward is Sah S^h_a =[(Rh−β)2|ah=a] =E[(R^h-β)^2|a_h=a] =Var(Rh−β|a)+[Rh−β|ah=a]2 =Var(R^h-β|a)+E[R^h-β|a_h=a]^2 =Var(Rh|ah=a)+(Gah−β)2 =Var(R^h|a_h=a)+(G^h_a-β)^2 =∑l=hH−1Var(rl|ah=a)+2∑h≤l<kCov(rl,rk|ah=a)+(Gah−β)2. = _l=h^H-1Var(r_l|a_h=a)+2 _h≤ l<kCov(r_l,r_k|a_h=a)+(G^h_a-β)^2. Let the opponent type at timestep l be given by YlY_l. For the focal agent, the action at the conditioned round is fixed: zh=1ah=C, z_h=1\a_h=C\, and for all round l>hl>h, the policy is given by x. Then the actions sampled in each round l>hl>h are ζl∼Ber(Yl),zl∼Ber(x) _l (Y_l), z_l (x) for the opponent and focal agent, respectively. Note that the focal agent’s actions are i.i.d. across the episode. The reward at round l is given by rl=bζl+c(1−zl) r_l=b _l+c(1-z_l) At the initial round l=hl=h, the action of the focal agent is fixed, so [rh|ah=a] [r_h|a_h=a] =b[ζh|ah=a]+c⋅1ah=D =bE[ _h|a_h=a]+c· 1\a_h=D\ =bmh(x)+c⋅1ah=D. =bm^h(x)+c· 1\a_h=D\. For all l>hl>h, the conditional mean for the focal agent taking action a at time h is [ζl|ah=a]=[Yl|ah=a]=mah,l(x) [ _l|a_h=a]=E[Y_l|a_h=a]=m_a^h,l(x) which means the expected conditional reward is [rl|ah=a] [r_l|a_h=a] =b[ζl|ah=a]+c[1−zl] =bE[ _l|a_h=a]+cE[1-z_l] =bmah,l(x)+c(1−x) =bm_a^h,l(x)+c(1-x) Likewise, the conditional variance when l=hl=h is Var(rh|ah=a) (r_h|a_h=a) =b2Var(ζh|ah=a) =b^2Var( _h|a_h=a) =b2mh(x)(1−mh(x)), =b^2m^h(x)(1-m^h(x)), and for l>hl>h it is Var(rl|ah=a) (r_l|a_h=a) =b2Var(ζl|ah=a)+c2Var(1−zl) =b^2Var( _l|a_h=a)+c^2Var(1-z_l) =b2mah,l(x)(1−mah,l(x))+c2x(1−x). =b^2m_a^h,l(x)(1-m_a^h,l(x))+c^2x(1-x). (11) The conditional means can be computed using repeated substitution of the iterative partner distribution update. For the covariance terms, let h≤l<k≤H−1h≤ l<k≤ H-1. Then Cov(rl,rk|ah=a) (r_l,r_k|a_h=a) =Cov(bζl+(1−zl)c,bζk+c(1−zk)|ah=a) =Cov(b _l+(1-z_l)c,b _k+c(1-z_k)|a_h=a) =b2Cov(ζl,ζk|ah=a)+bcCov(ζl,1−zk|ah=a)+bcCov(ζk,1−zl|ah=a) =b^2Cov( _l, _k|a_h=a)+bcCov( _l,1-z_k|a_h=a)+bcCov( _k,1-z_l|a_h=a) +c2Cov(1−zl,1−zk|ah=a) +c^2Cov(1-z_l,1-z_k|a_h=a) =b2Cov(Yl,Yk|ah=a)+bcCov(Yl,1−zk|ah=a)+bcCov(Yk,1−zl|ah=a) =b^2Cov(Y_l,Y_k|a_h=a)+bcCov(Y_l,1-z_k|a_h=a)+bcCov(Y_k,1-z_l|a_h=a) where the final term in the second line is zero due to the independence of action selection by the focal agent. Note that for H=2H=2, the final two terms are zero as l=hl=h is the only relevant term and the focal action is fixed. In this case, this covariance is given by Cov(Yl,Yk|ah=a) (Y_l,Y_k|a_h=a) =[YlYk|ah=a]−mah,l(x)mah,k(x). =E[Y_lY_k|a_h=a]-m_a^h,l(x)m_a^h,k(x). Computing the Covariance for H=2H=2 Beginning with the variance terms, we have for h=0h=0: [r0|Y0,a0=C] [r_0|Y_0,a_0=C] =b[ζ0|Y0]=bY0 =bE[ _0|Y_0]=bY_0 Var(r0|Y0,a0=C) (r_0|Y_0,a_0=C) =b2Var(ζ0|Y0)=b2Y0(1−Y0) =b^2Var( _0|Y_0)=b^2Y_0(1-Y_0) ⇒Var(r0|a0=C) (r_0|a_0=C) =[b2Y0(1−Y0)]+Var(bY0) =E[b^2Y_0(1-Y_0)]+Var(bY_0) =b2(μ1−μ12). =b^2( _1- _1^2). Likewise for defection: [r0|Y0,a0=D] [r_0|Y_0,a_0=D] =b[ζ0|Y0]+c=bY0+c =bE[ _0|Y_0]+c=bY_0+c Var(r0|Y0,a0=D) (r_0|Y_0,a_0=D) =b2Var(ζ0|Y0)=b2Y0(1−Y0) =b^2Var( _0|Y_0)=b^2Y_0(1-Y_0) ⇒Var(r0|a0=D) (r_0|a_0=D) =[b2Y0(1−Y0)]+Var(bY0+c) =E[b^2Y_0(1-Y_0)]+Var(bY_0+c) =b2(μ1−μ12). =b^2( _1- _1^2). For the following step: h=1h=1, [r1|Y1,a0=C] [r_1|Y_1,a_0=C] =b[ζ1|Y1]+c[1−z1]=bY1+c(1−x) =bE[ _1|Y_1]+cE[1-z_1]=bY_1+c(1-x) Var(r1|Y1,a0=C) (r_1|Y_1,a_0=C) =b2Var(ζ1|Y1)+c2Var(1−z1)=b2Y1(1−Y1)+c2x(1−x) =b^2Var( _1|Y_1)+c^2Var(1-z_1)=b^2Y_1(1-Y_1)+c^2x(1-x) ⇒Var(r1|a0=C) (r_1|a_0=C) =[b2Y1(1−Y1)+c2x(1−x)|a0=C]+Var(bY1+c(1−x)|a0=C) =E[b^2Y_1(1-Y_1)+c^2x(1-x)|a_0=C]+Var(bY_1+c(1-x)|a_0=C) =b2([Y1|a0=C]−[Y1|a0=C]2)+c2x(1−x). =b^2 (E[Y_1|a_0=C]-E[Y_1|a_0=C]^2 )+c^2x(1-x). Note this is the same form as Equation (11). Using the one-step Markov operator on the opponent distribution, we can derive the reward variances for each of the partner selection rules. The updated covariances can be computed using the conditional moment Ma,b0,1 M^0,1_a,b =[(R0−β)(R1−β)|a0=a,a1=b] =E[(R^0-β)(R^1-β)|a_0=a,a_1=b] =[(r0+r1−β)(r1−β)|a0=a,a1=b] =E[(r_0+r_1-β)(r_1-β)|a_0=a,a_1=b] =[r0r1|a0=a,a1=b]+[r12|a0=a,a1=b]−β[r0|a0=a,a1=b]−2β[r1|a0=a,a1=b]+β2. =E[r_0r_1|a_0=a,a_1=b]+E[r_1^2|a_0=a,a_1=b]- [r_0|a_0=a,a_1=b]-2 [r_1|a_0=a,a_1=b]+β^2. OFT The expected cooperation rate of the opponent in the next step given the current opponent and the focal agent’s action is [Y1|Y0,a0=C]=Y0⋅Y0+(1−Y0)μ1. [Y_1|Y_0,a_0=C]=Y_0· Y_0+(1-Y_0) _1. Then, computing the conditional mean we have [Y1|a0=C] [Y_1|a_0=C] =[[Y1|Y0,a0=C]] =E[E[Y_1|Y_0,a_0=C]] =[Y02+(1−Y0)μ1] =E[Y_0^2+(1-Y_0) _1] =μ2+μ1−μ12 = _2+ _1- _1^2 [Y1|a0=D] [Y_1|a_0=D] =μ1 = _1 For the covariance terms, we have [Y0Y1|a0=C] [Y_0Y_1|a_0=C] =[Y0[Y1|Y0,a0=C]] =E[Y_0E[Y_1|Y_0,a_0=C]] =[Y0(Y02+(1−Y0)μ1)] =E[Y_0(Y_0^2+(1-Y_0) _1)] =μ3+μ12−μ2μ1 = _3+ _1^2- _2 _1 [Y0Y1|a0=D] [Y_0Y_1|a_0=D] =[Y0[Y1|Y0,a0=D]] =E[Y_0E[Y_1|Y_0,a_0=D]] =[Y0μ1] =E[Y_0 _1] =μ12 = _1^2 The variance terms are Var(r1|a0=C)=b2(μ2+μ1−μ12)(1−μ2−μ1+μ12)+c2x(1−x),Var(r1|a0=D)=b2μ1(1−μ1)+c2x(1−x). (r_1|a_0=C)=b^2( _2+ _1- _1^2)(1- _2- _1+ _1^2)+c^2x(1-x),\;\;Var(r_1|a_0=D)=b^2 _1(1- _1)+c^2x(1-x). The covariance terms are Cov(r0,r1|a0=C) (r_0,r_1|a_0=C) =b2Cov(Y0,Y1|a0=C) =b^2Cov(Y_0,Y_1|a_0=C) =b2([Y0Y1|a0=C]−[Y0|a0=C][Y1|a0=C]) =b^2(E[Y_0Y_1|a_0=C]-E[Y_0|a_0=C]E[Y_1|a_0=C]) =b2(μ3+μ12−μ2μ1−μ1(μ2+μ1−μ12)) =b^2( _3+ _1^2- _2 _1- _1( _2+ _1- _1^2)) =b2(μ3−2μ2μ1+μ13) =b^2( _3-2 _2 _1+ _1^3) Cov(r0,r1|a0=D) (r_0,r_1|a_0=D) =b2Cov(Y0,Y1|a0=D) =b^2Cov(Y_0,Y_1|a_0=D) =b2([Y0Y1|a0=D]−[Y0|a0=D][Y1|a0=D]) =b^2(E[Y_0Y_1|a_0=D]-E[Y_0|a_0=D]E[Y_1|a_0=D]) =b2(μ12−μ12)=0 =b^2( _1^2- _1^2)=0 It remains to condition on the second step, h=1h=1. Explicitly, the terms are given by [r1|a1=C] [r_1|a_1=C] =bm1(x) =bm^1(x) =b∫01y[xyρ(y)+ρ(y)(1−xμ1)]y =b _0^1y[xyρ(y)+ρ(y)(1-x _1)]dy =b(μ1+x(μ2−μ12)) =b( _1+x( _2- _1^2)) [r1|a1=D] [r_1|a_1=D] =bm1(x)+c =bm^1(x)+c =b(μ1+x(μ2−μ12))+c =b( _1+x( _2- _1^2))+c Var(r1|a1=a) (r_1|a_1=a) =b2m1(x)(1−m1(x)) =b^2m^1(x)(1-m^1(x)) =b2(μ1+x(μ2−μ12))(1−(μ1+x(μ2−μ12))) =b^2( _1+x( _2- _1^2))(1-( _1+x( _2- _1^2))) Collating these terms, we can summarise the elements of the second moment in Table 6. Table 6: Second moments (SahS^h_a) of episodic reward when H=2H=2 under OFT. h a Var(Rh|ah=a)Var(R^h|a_h=a) (Gah−β)2(G^h_a-β)^2 0 C b2[(μ2+μ1−μ12)(1−μ2−μ1+μ12)+(μ1−μ12)+2(μ3−2μ2μ1+μ13)]+c2x(1−x)b^2[( _2+ _1- _1^2)(1- _2- _1+ _1^2)+( _1- _1^2) +2( _3-2 _2 _1+ _1^3)]+c^2x(1-x) [b(μ2+2μ1−μ12)+c(1−x)−β]2[b( _2+2 _1- _1^2)+c(1-x)-β]^2 D 2b2μ1(1−μ1)+c2x(1−x)2b^2 _1(1- _1)+c^2x(1-x) [2bμ1+c(2−x)−β]2[2b _1+c(2-x)-β]^2 11 C b2(μ1+x(μ2−μ12))(1−μ1−x(μ2−μ12))b^2( _1+x( _2- _1^2))(1- _1-x( _2- _1^2)) [b(μ1+x(μ2−μ12))−β]2[b( _1+x( _2- _1^2))-β]^2 D b2(μ1+x(μ2−μ12))(1−μ1−x(μ2−μ12))b^2( _1+x( _2- _1^2))(1- _1-x( _2- _1^2)) [b(μ1+x(μ2−μ12))+c−β]2[b( _1+x( _2- _1^2))+c-β]^2 ROFT The expected cooperation rate of the opponent in the next step given the current opponent and the focal agent’s action is [Y1|Y0,a0=C] [Y_1|Y_0,a_0=C] =μ1 = _1 [Y1|Y0,a0=D] [Y_1|Y_0,a_0=D] =Y0⋅(1−Y0)+Y0μ1 =Y_0·(1-Y_0)+Y_0 _1 Then computing the conditional mean, [Y1|a0=C] [Y_1|a_0=C] =μ1 = _1 [Y1|a0=D] [Y_1|a_0=D] =[[Y1|Y0,a0=D]] =E[E[Y_1|Y_0,a_0=D]] =[Y0(1−Y0)+Y0μ1] =E[Y_0(1-Y_0)+Y_0 _1] =μ1+μ12−μ2 = _1+ _1^2- _2 For the covariance terms, we have [Y0Y1|a0=C] [Y_0Y_1|a_0=C] =[Y0[Y1|Y0,a0=C]] =E[Y_0E[Y_1|Y_0,a_0=C]] =[Y0μ1] =E[Y_0 _1] =μ12 = _1^2 [Y0Y1|a0=D] [Y_0Y_1|a_0=D] =[Y0[Y1|Y0,a0=D]] =E[Y_0E[Y_1|Y_0,a_0=D]] =[Y0(Y0(1−Y0)+Y0μ1)] =E[Y_0(Y_0(1-Y_0)+Y_0 _1)] =μ2−μ3+μ1μ2 = _2- _3+ _1 _2 The variance terms are then: Var(r1|a0=C)=b2μ1(1−μ1)+c2x(1−x),Var(r1|a0=D)=b2(μ1+μ12−μ2)(1−μ1−μ12+μ2)+c2x(1−x) (r_1|a_0=C)=b^2 _1(1- _1)+c^2x(1-x),\;\;Var(r_1|a_0=D)=b^2( _1+ _1^2- _2)(1- _1- _1^2+ _2)+c^2x(1-x) The covariance terms are Cov(r0,r1|a0=C) (r_0,r_1|a_0=C) =b2Cov(Y0,Y1|a0=C) =b^2Cov(Y_0,Y_1|a_0=C) =b2([Y0Y1|a0=C]−[Y0|a0=C][Y1|a0=C]) =b^2(E[Y_0Y_1|a_0=C]-E[Y_0|a_0=C]E[Y_1|a_0=C]) =b2(μ12−μ12)=0 =b^2( _1^2- _1^2)=0 Cov(r0,r1|a0=D) (r_0,r_1|a_0=D) =b2Cov(Y0,Y1|a0=D) =b^2Cov(Y_0,Y_1|a_0=D) =b2([Y0Y1|a0=D]−[Y0|a0=D][Y1|D]) =b^2(E[Y_0Y_1|a_0=D]-E[Y_0|a_0=D]E[Y_1|D]) =b2(μ2−μ3+μ1μ2−μ1(μ1+μ12−μ2)) =b^2( _2- _3+ _1 _2- _1( _1+ _1^2- _2)) =b2(μ2−μ12+2μ1μ2−μ13−μ3) =b^2( _2- _1^2+2 _1 _2- _1^3- _3) It remains to condition on the second step, h=1h=1. Explicitly, the terms are given by [r1|a1=C] [r_1|a_1=C] =bm1(x) =bm^1(x) =b∫01y[(1−x)(1−y)ρ(y)+ρ(y)(1−(1−x)(1−μ1))]y =b _0^1y[(1-x)(1-y)ρ(y)+ρ(y)(1-(1-x)(1- _1))]dy =b(μ1−(1−x)(μ2−μ12)) =b( _1-(1-x)( _2- _1^2)) [r1|a1=D] [r_1|a_1=D] =bm1(x)+c =bm^1(x)+c =b(μ1−(1−x)(μ2−μ12))+c =b( _1-(1-x)( _2- _1^2))+c Var(r1|a1=a) (r_1|a_1=a) =b2m1(x)(1−m1(x)) =b^2m^1(x)(1-m^1(x)) =b2(μ1−(1−x)(μ2−μ12))(1−μ1+(1−x)(μ2−μ12))) =b^2( _1-(1-x)( _2- _1^2))(1- _1+(1-x)( _2- _1^2))) Collating these terms, we can summarise the elements of the second moment in Table 7. Table 7: Second moments (SahS^h_a) of episodic reward when H=2H=2 under ROFT. h a Var(Rh|ah=a)Var(R^h|a_h=a) (Gah−β)2(G^h_a-β)^2 0 C 2b2μ1(1−μ1)+c2x(1−x)2b^2 _1(1- _1)+c^2x(1-x) [2bμ1+c(1−x)−β]2[2b _1+c(1-x)-β]^2 D b2[(μ1+μ12−μ2)(1−μ1−μ12+μ2)+μ1(1−μ1)+2(μ2−μ12+2μ1μ2−μ13−μ3)]+c2x(1−x)b^2[( _1+ _1^2- _2)(1- _1- _1^2+ _2)+ _1(1- _1) +2( _2- _1^2+2 _1 _2- _1^3- _3)]+c^2x(1-x) [b(2μ1+μ12−μ2)+c(2−x)−β]2[b(2 _1+ _1^2- _2)+c(2-x)-β]^2 11 C b2(μ1−(1−x)(μ2−μ12))(1−μ1+(1−x)(μ2−μ12))b^2( _1-(1-x)( _2- _1^2))(1- _1+(1-x)( _2- _1^2)) [b(μ1−(1−x)(μ2−μ12))−β]2[b( _1-(1-x)( _2- _1^2))-β]^2 D b2(μ1−(1−x)(μ2−μ12))(1−μ1+(1−x)(μ2−μ12))b^2( _1-(1-x)( _2- _1^2))(1- _1+(1-x)( _2- _1^2)) [b(μ1−(1−x)(μ2−μ12))+c−β]2[b( _1-(1-x)( _2- _1^2))+c-β]^2 Always Stay Since the transition is independent of the action, for both actions a∈C,Da∈\C,D\ the expected cooperation rate of the opponent in the next step given the current opponent and action is [Y1|Y0,a0=a]=Y0. [Y_1|Y_0,a_0=a]=Y_0. Then computing the conditional mean, [Y1|a0=a] [Y_1|a_0=a] =[[Y1|Y0,a0=a]]=[Y0]=μ1 =E[E[Y_1|Y_0,a_0=a]]=E[Y_0]= _1 For the covariance terms, we have [Y0Y1|a0=a] [Y_0Y_1|a_0=a] =[Y0[Y1|Y0,a0=a]] =E[Y_0E[Y_1|Y_0,a_0=a]] =[Y0Y0] =E[Y_0Y_0] =μ2 = _2 The variance terms are then: Var(r1|a0=a)=b2μ1(1−μ1)+c2x(1−x) (r_1|a_0=a)=b^2 _1(1- _1)+c^2x(1-x) The covariance terms are Cov(r0,r1|a0=a) (r_0,r_1|a_0=a) =b2Cov(Y0,Y1|a0=a) =b^2Cov(Y_0,Y_1|a_0=a) =b2([Y0Y1|a0=a]−[Y0|a][Y1|a0=a]) =b^2(E[Y_0Y_1|a_0=a]-E[Y_0|a]E[Y_1|a_0=a]) =b2(μ2−μ12) =b^2( _2- _1^2) =b2Var(ρ) =b^2Var(ρ) It remains to condition on the second step, h=1h=1. Explicitly, the terms are given by [r1|a1=C] [r_1|a_1=C] =bm1(x) =bm^1(x) =b∫01yρ(y)y =b _0^1yρ(y)dy =bμ1 =b _1 [r1|a1=D] [r_1|a_1=D] =bm1(x)+c =bm^1(x)+c =bμ1+c =b _1+c Var(r1|a1=a) (r_1|a_1=a) =b2m1(x)(1−m1(x)) =b^2m^1(x)(1-m^1(x)) =b2μ1(1−μ1) =b^2 _1(1- _1) Collating these terms, we can summarise the elements of the second moment in Table 8. Table 8: Second moments (SahS^h_a) of episodic reward when H=2H=2 under Always Stay. h a Var(Rh|ah=a)Var(R^h|a_h=a) (Gah−β)2(G^h_a-β)^2 0 C 2b2(μ1+μ2−2μ12)+c2x(1−x)2b^2( _1+ _2-2 _1^2)+c^2x(1-x) [2bμ1+c(1−x)−β]2[2b _1+c(1-x)-β]^2 D 2b2(μ1+μ2−2μ12)+c2x(1−x)2b^2( _1+ _2-2 _1^2)+c^2x(1-x) [2bμ1+c(2−x)−β]2[2b _1+c(2-x)-β]^2 11 C b2μ1(1−μ1)b^2 _1(1- _1) [bμ1−β]2[b _1-β]^2 D b2μ1(1−μ1)b^2 _1(1- _1) [bμ1+c−β]2[b _1+c-β]^2 Always Switch Since the transition is independent of the action, for both actions a∈C,Da∈\C,D\ the expected cooperation rate of the opponent in the next step given the current opponent and action is [Y1|Y0,a0=a]=μ1. [Y_1|Y_0,a_0=a]= _1. Then computing the conditional mean, [Y1|a] [Y_1|a] =[[Y1|Y0,a0=a]]=[μ1]=μ1 =E[E[Y_1|Y_0,a_0=a]]=E[ _1]= _1 For the covariance terms, we have [Y0Y1|a0=a] [Y_0Y_1|a_0=a] =[Y0[Y1|Y0,a0=a]] =E[Y_0E[Y_1|Y_0,a_0=a]] =[Y0μ1] =E[Y_0 _1] =μ12 = _1^2 The variance terms are then: Var(r1|a0=a)=b2μ1(1−μ1)+c2x(1−x) (r_1|a_0=a)=b^2 _1(1- _1)+c^2x(1-x) The covariance terms are Cov(r0,r1|a0=a) (r_0,r_1|a_0=a) =b2Cov(Y0,Y1|a0=a) =b^2Cov(Y_0,Y_1|a_0=a) =b2([Y0Y1|a0=a]−[Y0|a0=a][Y1|a0=a]) =b^2(E[Y_0Y_1|a_0=a]-E[Y_0|a_0=a]E[Y_1|a_0=a]) =b2(μ12−μ12)=0 =b^2( _1^2- _1^2)=0 Note that for conditioning at step h=1h=1, the same derivation as Always Stay holds. Collating these terms, we can summarise the elements of the second moment in Table 9. Table 9: Second moments (SahS^h_a) of episodic reward when H=2H=2 under Always Switch. h a Var(Rh|ah=a)Var(R^h|a_h=a) (Gah−β)2(G^h_a-β)^2 0 C 2b2μ1(1−μ1)+c2x(1−x)2b^2 _1(1- _1)+c^2x(1-x) [2bμ1+c(1−x)−β]2[2b _1+c(1-x)-β]^2 D 2b2μ1(1−μ1)+c2x(1−x)2b^2 _1(1- _1)+c^2x(1-x) [2bμ1+c(2−x)−β]2[2b _1+c(2-x)-β]^2 11 C b2μ1(1−μ1)b^2 _1(1- _1) [bμ1−β]2[b _1-β]^2 D b2μ1(1−μ1)b^2 _1(1- _1) [bμ1+c−β]2[b _1+c-β]^2 Appendix B Missing Proofs B.1 Conditional Rewards We begin by analysing the base case of k=0k=0. Defining Δρ0,h+1(y)=ρh+1(y|x,C0)−ρh+1(y|x,D0) ρ^0,h+1(y)=ρ^h+1(y|x,C^0)-ρ^h+1(y|x,D^0), the system reduces to analysing the following recursion: Δρ0,h+1(y) ρ^0,h+1(y) =xyΔρ0,h(y)−xρ(y)Δm0,h =xy ρ^0,h(y)-xρ(y) m^0,h Δm0,h m^0,h =∫01yΔρ0,h(y)y = _0^1y ρ^0,h(y)dy Define the operator T as the recursive update on a function g. Formally, g:=Yg−[Yg] :=Yg-E[Yg] then gh+1=ghg^h+1=Tg^h. For h>1h>1 gh(y) g^h(y) =ygh−1(y)−[Ygh−1(Y)],g1(y)=y−μ1. =yg^h-1(y)-E[Yg^h-1(Y)], g^1(y)=y- _1. (12) Then we can isolate the partner dependence y, and from the policy x. Lemma B.1. Let ρ(y)ρ(y) be the initial distribution, and ghg^h be defined by Equation (12). Then for all h≥1h≥ 1, the difference in distribution and mean of the opponent is given by Δρ0,h(y) ρ^0,h(y) =xh−1gh(y)ρ(y),Δm0,h=xh−1[Ygh(Y)]. =x^h-1g^h(y)ρ(y), m^0,h=x^h-1E[Yg^h(Y)]. Proof. We proceed by induction. Let h=1h=1, then Δρ0,1(y) ρ^0,1(y) =(y−μ1)ρ(y)=x0g1(y)ρ(y) =(y- _1)ρ(y)=x^0g^1(y)ρ(y) Δm0,1 m^0,1 =Var(ρ)=[Y2−μY]=x0[Yg1(Y)] =Var(ρ)=E[Y^2-μ Y]=x^0E[Yg^1(Y)] Assume the form holds at time step h, then for h+1h+1, Δρ0,h+1(y) ρ^0,h+1(y) =xyΔρ0,h(y)−xρ(y)Δm0,h =xy ρ^0,h(y)-xρ(y) m^0,h =xy(xh−1ghρ(y))−xρ(y)xh−1[Ygh(Y)] =xy (x^h-1g^hρ(y) )-xρ(y)x^h-1E[Yg^h(Y)] =xh(ygh−[Ygh(Y)])ρ(y) =x^h (yg^h-E[Yg^h(Y)] )ρ(y) =xhgh+1(y)ρ(y) =x^hg^h+1(y)ρ(y) Δm0,h+1 m^0,h+1 =xh∫01ygh+1(y)ρ(y)y =x^h _0^1yg^h+1(y)ρ(y)dy =xh[Ygh+1(Y)] =x^hE[Yg^h+1(Y)] ∎ This form enables us to provide a simple condition for the future rewards to be non-negative. Specifically, for Δm0,h≥0 m^0,h≥ 0 we need [Ygh(Y)]≥0E[Yg^h(Y)]≥ 0. This can be analysed through the properties of the operator T, by splitting the steps of the update. Define the inner product of two functions as ⟨f,g⟩:=[Yf(Y)g(Y)]. f,g :=E[Yf(Y)g(Y)]. (13) Lemma B.2. The operator T is self-adjoint with respect to the inner product in Equation (13). Proof. ⟨f,g⟩ ,g =⟨Yf−[Yf],g⟩ = Yf-E[Yf],g =[Y(Yf−[Yf])g] =E[Y(Yf-E[Yf])g] =[Y2fg]−[Yf][Yg] =E[Y^2fg]-E[Yf]E[Yg] =[Yf(Yg−[Yg])] =E[Yf(Yg-E[Yg])] =⟨f,Yg−[Yg]⟩ = f,Yg-E[Yg] =⟨f,g⟩. = f,Tg . ∎ Lemma B.3. Let ghg^h be defined by Equation (12). For any m,n∈ℕm,n , [Ygmgn]=[Ygm+n]. [Yg^mg^n]=E[Yg^m+n]. Proof. By Lemma B.2, we have ⟨gm,gn⟩ g^m,g^n =⟨m−1g1,n−1g1⟩ = ^m-1g^1,T^n-1g^1 =⟨g1,m+n−2g1⟩ = g^1,T^m+n-2g^1 =⟨g1,gm+n−1⟩ = g^1,g^m+n-1 and using the definition of the inner product, for any k∈ℕk ⟨g1,gk⟩ g^1,g^k =[Yg1gk] =E[Yg^1g^k] =[Y(Y−μ1)gk] =E[Y(Y- _1)g^k] =[Y2gk]−μ1[Ygk] =E[Y^2g^k]- _1E[Yg^k] From the recursion, we can take expectation over Ygk+1Yg^k+1, giving the identity [Ygk+1] [Yg^k+1] =[Y2gk]−μ1[Ygk] =E[Y^2g^k]- _1E[Yg^k] Therefore, ⟨g1,gk⟩ g^1,g^k =[Yg1gk]=[Ygk+1] =E[Yg^1g^k]=E[Yg^k+1] Finally, [Ygmgn] [Yg^mg^n] =⟨gm,gn⟩=⟨g1,gm+n−1⟩=[Ygm+n]. = g^m,g^n = g^1,g^m+n-1 =E[Yg^m+n]. ∎ Lemma B.4. Under the OFT rule, for all h≥1h≥ 1 the reward difference Δr0,h=bΔm0,h≥0 r^0,h=b m^0,h≥ 0. Proof. We will separate into odd and even cases. First, let m=n=lm=n=l in Lemma B.3, then Δm0,2l m^0,2l =x2l−1[Yg2l]=x2l−1[Yglgl]=x2l−1[Y(gl)2]≥0 =x^2l-1E[Yg^2l]=x^2l-1E[Yg^lg^l]=x^2l-1E[Y(g^l)^2]≥ 0 Likewise, let m=lm=l and n=l+1n=l+1. By Lemma B.3 Δm0,2l+1 m^0,2l+1 =x2l[Yg2l+1] =x^2lE[Yg^2l+1] =x2l[Yglgl+1] =x^2lE[Yg^lg^l+1] =x2l[Ygl(Ygl−[Ygl])] =x^2lE[Yg^l(Yg^l-E[Yg^l])] =x2l([(Ygl)2]−[Ygl]2)≥0 =x^2l (E[(Yg^l)^2]-E[Yg^l]^2 )≥ 0 Hence, since b>0b>0, Δr0,h≥0 r^0,h≥ 0 for all h≥1h≥ 1. ∎ Here, we show how the same result applies to ROFT. Lemma B.5. Under the ROFT rule, for all h≥1h≥ 1 the reward difference Δr0,h=bΔm0,h≥0 r^0,h=b m^0,h≥ 0. Proof. Under the ROFT partner selection rule, Δρ0,h+1(y)=(1−x)(1−y)Δρ0,h(y)+(1−x)ρ(y)Δm0,h(x) ρ^0,h+1(y)=(1-x)(1-y) ρ^0,h(y)+(1-x)ρ(y) m^0,h(x) and using the same inductive proof in Lemma B.1, we get Δρ0,h(y) ρ^0,h(y) =(1−x)h−1fh(y)ρ(y) =(1-x)^h-1f^h(y)ρ(y) Δm0,h m^0,h =(1−x)h−1[Yfh(Y)] =(1-x)^h-1E[Yf^h(Y)] where fh(y)=(1−y)fh−1(y)+[Yfh−1(Y)]f^h(y)=(1-y)f^h-1(y)+E[Yf^h-1(Y)], with f1(y)=y−μ1f^1(y)=y- _1 and [fh(Y)]=0E[f^h(Y)]=0. This means we can write fh(y) f^h(y) =(1−y)fh−1(y)−[(1−Y)fh−1(Y)]. =(1-y)f^h-1(y)-E[(1-Y)f^h-1(Y)]. Now consider the transformation Z=1−YZ=1-Y, and the function ϕh(z)=−fh(1−z)φ^h(z)=-f^h(1-z). Then, ϕh(z)=zϕh−1(z)−[Zϕh−1(Z)],ϕ1=z−[Z] φ^h(z)=zφ^h-1(z)-E[Zφ^h-1(Z)], φ^1=z-E[Z] this is the exact same recursion defined for the OFT partner selection rule, and therefore the result of Lemma B.4 follows. Concretely, Δm0,h m^0,h =(1−x)h−1[Yfh(Y)] =(1-x)^h-1E[Yf^h(Y)] =(1−x)h−1[Zϕh(Z)]≥0. =(1-x)^h-1E[Zφ^h(Z)]≥ 0. Since b>0b>0, the result follows. ∎ We now consider the general case of conditioning on the action at step k. Note that by dividing through by ρ(y)ρ(y), the one-step update for OFT becomes qh+1(y|x) q^h+1(y|x) =xyqh(y|x)+(1−xmh(x)),mh(x)=[Yqh(Y|x)] =xyq^h(y|x)+(1-xm^h(x)), m^h(x)=E[Yq^h(Y|x)] This distribution can then be expressed in terms of the recursive function, gh(y)g^h(y). Lemma B.6. For all h, the distribution qh(y|x)q^h(y|x) is given by qh(y|x) q^h(y|x) =1+∑j=1hxjgj(y) =1+Σ^h_j=1x^jg^j(y) Proof. We proceed by induction. Let h=0h=0, then q0(y|x) q^0(y|x) =ρ0(y|x)ρ(y)=1. = ρ^0(y|x)ρ(y)=1. Assume that the statement holds for h=kh=k. Then for h=k+1h=k+1, qk+1(y|x) q^k+1(y|x) =xyqk(y|x)+1−xmk(x) =xyq^k(y|x)+1-xm^k(x) =xy(1+∑j=1kxjgj(y))+1−x[Y(1+∑j=1kxjgj(Y))] =xy (1+ _j=1^kx^jg^j(y) )+1-xE [Y (1+ _j=1^kx^jg^j(Y) ) ] =1+x(y−μ1)+∑j=1kxj+1(ygj(y)−[Ygj(Y)]) =1+x(y- _1)+ _j=1^kx^j+1 (yg^j(y)-E[Yg^j(Y)] ) =1+xg1(y)+∑j=1kxj+1gj+1(y) =1+xg^1(y)+ _j=1^kx^j+1g^j+1(y) =1+∑j=1k+1xjgj(y) =1+ _j=1^k+1x^jg^j(y) ∎ Proof of Proposition 3.1 Proof. Consider the kkth step of the REINFORCE algorithm. Then the subsequent update under the OFT rule is, Δρk,k+1 ρ^k,k+1 =yρk(y|x)−mk(x)ρ(y) =yρ^k(y|x)-m^k(x)ρ(y) gk,k+1 g^k,k+1 =yqk(y|x)−mk(x) =yq^k(y|x)-m^k(x) =y(1+∑j=1kxjgj(y))−(μ1+∑j=1kxj[Ygj(Y)]) =y (1+Σ^k_j=1x^jg^j(y) )- ( _1+Σ^k_j=1x^jE[Yg^j(Y)] ) =y−μ1+∑j=1kxj(ygj(y)−[Ygj(Y)]) =y- _1+Σ^k_j=1x^j (yg^j(y)-E[Yg^j(Y)] ) =y−μ1+∑j=1kxjgj+1(y) =y- _1+Σ^k_j=1x^jg^j+1(y) =∑j=0kxjgj+1(y) = _j=0^kx^jg^j+1(y) Then for h>k+1h>k+1, we have Δρk,h+1(y) ρ^k,h+1(y) =xyΔρk,h(y)−xρ(y)Δmk,h =xy ρ^k,h(y)-xρ(y) m^k,h and therefore the same operator T applies. Consequently, gk,h(y) g^k,h(y) =h−k−1gk,k+1(y) =T^h-k-1g^k,k+1(y) =∑j=0kxjh−k−1gj+1(y) = _j=0^kx^jT^h-k-1g^j+1(y) =∑j=0kxjgh−k+j(y) = _j=0^kx^jg^h-k+j(y) The difference in the conditional means at this step is then Δmk,h m^k,h =xh−k−1[Ygk,h(Y)] =x^h-k-1E[Yg^k,h(Y)] =xh−k−1[Y∑j=0kxjgh−k+j(Y)] =x^h-k-1E[Y _j=0^kx^jg^h-k+j(Y)] =∑j=0kxh−k+j−1[Ygh−k+j(Y)] = _j=0^kx^h-k+j-1E[Yg^h-k+j(Y)] =∑j=0kΔm0,h−k+j≥0 = _j=0^k m^0,h-k+j≥ 0 where the final line follows by Lemma B.4. Note that the proof has an identical structure for the ROFT rule, invoking Lemma B.5 in the final step. ∎ Proof of Corollary 3.1 Proof. For each step k of the REINFORCE algorithm, we obtain a reward difference of Δrk,k=−c r^k,k=-c. The future rewards must compensate for this. Consider the actions conditioned at step k. By the proof of Proposition 3.1, at step k+1k+1 Δrk,k+1(x) r^k,k+1(x) =bΔmk,k+1(x) =b m^k,k+1(x) =b∑j=0kΔm0,j+1 =b _j=0^k m^0,j+1 ≥bΔm0,1. ≥ b m^0,1. Since we have Δm0,1=Var(ρ) m^0,1=Var(ρ), and for all h>k+1h>k+1 the reward difference Δrk,h(x)≥0 r^k,h(x)≥ 0, then ∑h=0H−1ΔGh[ρ] Σ^H-1_h=0 G^h[ρ] =∑k=0H−1∑h=kH−1Δrk,h =Σ^H-1_k=0 _h=k^H-1 r^k,h ≥∑k=0H−1(Δrk,k+Δrk,k+1) ≥Σ^H-1_k=0 ( r^k,k+ r^k,k+1 ) ≥∑k=0H−1(−c+bVar(ρ)) ≥Σ^H-1_k=0 (-c+bVar(ρ) ) ≥(H−1)(bVar(ρ)−c)−c. ≥(H-1)(bVar(ρ)-c)-c. which is increasing in H when ΔG[ρ]>0 G[ρ]>0 for H=2H=2. ∎ B.2 Mean Dynamics Proof of Theorem 3.1 Proof. The characteristic flow is generated by dtXt(x)=−4αcXt2(x)(1−Xt(x))2,X0(x)=x0, ddtX_t(x)=-4α cX^2_t(x)(1-X_t(x))^2, X_0(x)=x_0, This is a separable ODE, which integrating gives ∫1X2(1−X)2X 1X^2(1-X)^2dX =∫−4αcdt = -4α c\;dt 11−Xt−1Xt+2ln|Xt1−Xt| 11-X_t- 1X_t+2 | X_t1-X_t| =−4αct+C =-4α ct+C Let F(x) F(x) =11−x−1x+2ln|x1−x| = 11-x- 1x+2 | x1-x| then the solution at time t is given by F(Xt)=F(x0)−4αct. F(X_t)=F(x_0)-4α ct. For c>0c>0 and α>0α>0, then as t→∞t→∞ we have F(Xt)→−∞F(X_t)→-∞ for every x0∈(0,1)x_0∈(0,1). Since F(x)→−∞F(x)→-∞ as x→0x→ 0, we have that Xt→0X_t→ 0. Since the velocity field u(x)=−4αcx2(1−x)2∈C1([0,1])u(x)=-4α cx^2(1-x)^2∈ C^1([0,1]), the solution to the continuity equation is given by the pushforward of the initial point under the characteristic flow: ρ(t,⋅)=(Xt)#ρ0. ρ(t,·)=(X_t)_\# _0. Then for any bounded test function φ∈C([0,1]) ∈ C([0,1]), ∫01φ(y)ρ(t,y)y=∫01φ(X(t,x))ρ0(x)x _0^1 (y)ρ(t,y)\;dy= _0^1 (X(t,x)) _0(x)dx Since X(t,x)→0X(t,x)→ 0 pointwise for all x∈(0,1)x∈(0,1), then φ(X(t,x))→φ(0) (X(t,x))→ (0). Moreover, |φ(X(t,x))ρ0(x)| | (X(t,x)) _0(x)| ≤‖φ‖∞|ρ0(x)| ≤\| \|_∞| _0(x)| By the dominated convergence theorem, ∫01φ(X(t,x))ρ0(x)x→φ(0)∫01ρ0(x)x=φ(0) _0^1 (X(t,x)) _0(x)dx→ (0) _0^1 _0(x)dx= (0) which means ∫01φ(y)ρ(t,y)y→φ(0) _0^1 (y)ρ(t,y)\;dy→ (0) and therefore ρ(t)⇀δ0ρ(t) _0. ∎ Lemma B.7. The non-local advection is spatially Lipschitz and sub-linear. That is, there exists an L1,M∈ℝL_1,M such that |u[ρ](x)−u[ρ](y)|≤L1|x−y|,|u[ρ](x)|≤M(1+|x|). |u[ρ](x)-u[ρ](y)|≤ L_1|x-y|, |u[ρ](x)|≤ M(1+|x|). Proof. The spatial derivative is ∂xu[ρ](x)=4αΔG[ρ]⋅x(1−x)(1−2x) _xu[ρ](x)=4α G[ρ]· x(1-x)(1-2x) which is uniformly bounded by ΔG[ρ] G[ρ] on [0,1][0,1]. Therefore by the mean value theorem, |u[ρ](x)−u[ρ](y)| |u[ρ](x)-u[ρ](y)| ≤supz|∂xu[ρ](z)||x−y| _z| _xu[ρ](z)|\;|x-y| ≤2α|ΔG[ρ]||x−y| ≤ 2α| G[ρ]||x-y| Since the episode length and per-game reward are bounded, |ΔG[ρ]|| G[ρ]| is uniformly bounded. Therefore exists an L1L_1 such that u[ρ](y)u[ρ](y) is spatially Lipschitz with constant L1L_1. The condition for sub-linearity follows from |u[ρ](x)| |u[ρ](x)| =|2αΔG[ρ]x2(1−x)2| =|2α G[ρ]x^2(1-x)^2| ≤2α|ΔG[ρ]||x2(1−x)2| ≤ 2α| G[ρ]||x^2(1-x)^2| ≤α8|ΔG[ρ]| ≤ α8| G[ρ]| which is uniformly bounded. ∎ Lemma B.8. For every ρ1,ρ2∈ _1, _2 , there exists a positive constant L2∈ℝL_2 such that ‖u[ρ1](⋅)−u[ρ2](⋅)‖([0,1])≤L2W1(ρ1,ρ2). \|u[ _1](·)-u[ _2](·)\|_C([0,1])≤ L_2W_1( _1, _2). Proof. |u[ρ1](x)−u[ρ2](x)| |u[ _1](x)-u[ _2](x)| =|2αΔG[ρ1]x2(1−x)2−2αΔG[ρ2]x2(1−x)2| =|2α G[ _1]x^2(1-x)^2-2α G[ _2]x^2(1-x)^2| =2α|x2(1−x)2(ΔG[ρ1]−ΔG[ρ2])| =2α|x^2(1-x)^2( G[ _1]- G[ _2])| ≤α8|ΔG[ρ1]−ΔG[ρ2]| ≤ α8| G[ _1]- G[ _2]| Therefore, ‖u[ρ1](⋅)−u[ρ2](⋅)‖C0([0,1])≤α8|ΔG[ρ1]−ΔG[ρ2]| \|u[ _1](·)-u[ _2](·)\|_C^0([0,1])≤ α8| G[ _1]- G[ _2]| It suffices to show that |ΔG[ρ1]−ΔG[ρ2]| | G[ _1]- G[ _2]| ≤CW1(ρ1,ρ2) ≤ CW_1( _1, _2) for some constant C∈ℝC . Using the Kantorovich-Rubinstein duality [villani], the Wasserstein distance with p=1p=1 can be expressed as W1(ρ1,ρ2)=sup∫[0,1]f(x)d(ρ1−ρ2)|f∈1,f:[0,1]→ℝ,Lip(f)≤1. W_1( _1, _2)= \ _[0,1]f(x)d( _1- _2)\; |\;f ^1,\;f:[0,1] ,\;Lip(f)≤ 1 \. For any integer i, the iith moment is mi m_i =∫01xiρ(x) = _0^1x^idρ(x) Let f(x)=xiif(x)= x^ii, then f is continuous and has Lipschitz constant 1. By the Wasserstein metric, 1i|mi(ρ1)−mi(ρ2)|≤W1(ρ1,ρ2). 1i|m_i( _1)-m_i( _2)|≤ W_1( _1, _2). Therefore, |ΔG[ρ1]−ΔG[ρ2]| | G[ _1]- G[ _2]| =b|Var(ρ1)−Var(ρ2)| =b|Var( _1)-Var( _2)| ≤b|m2(ρ1)−m2(ρ2)|+b|m1(ρ1)2−m1(ρ2)2| ≤ b|m_2( _1)-m_2( _2)|+b|m_1( _1)^2-m_1( _2)^2| ≤2bW1(ρ1,ρ2)+b|m1(ρ1)+m1(ρ2)||m1(ρ1)−m1(ρ2)| ≤ 2bW_1( _1, _2)+b|m_1( _1)+m_1( _2)||m_1( _1)-m_1( _2)| ≤2bW1(ρ1,ρ2)+b⋅2⋅W1(ρ1,ρ2) ≤ 2bW_1( _1, _2)+b· 2· W_1( _1, _2) =4bW1(ρ1,ρ2) =4bW_1( _1, _2) Therefore the result holds with L2=b4L_2= b4. ∎ B.2.1 Proof of Proposition 3.2 Proof. By Theorem 2 in [bonnet_pontryagin_2019], it suffices to show the conditions of Lemma B.7 and Lemma B.8 are satisfied. These require that ΔG[ρ] G[ρ] is bounded, which is always satisfied for ρ0∈([0,1)] _0 ([0,1)]. ∎ The characteristic ODE admits a solution XtX_t dtXt ddtX_t =2αΔG[ρ]Xt2(1−Xt)2 =2α G[ρ]X_t^2(1-X_t)^2 ∫1X2(1−X)2X 1X^2(1-X)^2dX =2α∫ΔG[ρ]t =2α G[ρ]dt F(Xt) F(X_t) =F(x0)+2αK(t) =F(x_0)+2α K(t) where F(x)F(x) is defined in the proof of Theorem 3.1. Clearly, the time evolution is uniquely determined by K. Since F′(x)>0F (x)>0 for all x∈(0,1)x∈(0,1), F is strictly increasing and therefore invertible on the domain. Therefore the solution can be expressed by XK(x)=F−1(F(x0)+2αK). X_K(x)=F^-1(F(x_0)+2α K). The next result proves that under this flow, the cooperation levels increase when the initial population has sufficient variance. Proof of Theorem 3.2 Proof. Define the autonomous equations K′ K =h(K), =h(K), h(K) h(K) =bVar((XK)#ρ0)−2c =bVar ((X_K)_\# _0 )-2c The initial condition means h(0) h(0) =bVar(ρ0)−2c>0 =bVar( _0)-2c>0 Note that for every x>0x>0 limK→∞XK(x)→1,∂KXK=2αXK2(1−XK)2≤α/8, _K→∞X_K(x)→ 1, _KX_K=2α X_K^2(1-X_K)^2≤α/8, therefore XKX_K is uniformly Lipschitz in K and bounded on [0,1][0,1]. In particular, |XK1(x)−XK2(x)|≤α8|K1−K2|,∀x∈[0,1]. |X_K_1(x)-X_K_2(x)|≤ α8|K_1-K_2|, ∀ x∈[0,1]. For X∼ρ0X _0, the first moment is |m1(K1)−m1(K2)| |m_1(K_1)-m_1(K_2)| =|[XK1(X)]−[XK2(X)]| =|E[X_K_1(X)]-E[X_K_2(X)]| ≤[|XK1(X)−XK2(X)]| [|X_K_1(X)-X_K_2(X)]| ≤α8|K1−K2|. ≤ α8|K_1-K_2|. It follows that the second moment is Lipschitz: |m2(K1)−m2(K2)| |m_2(K_1)-m_2(K_2)| =|[XK12(X)]−[XK22(X)]| =|E[X_K_1^2(X)]-E[X^2_K_2(X)]| =|[XK12(X)−XK22(X)]| =|E[X_K_1^2(X)-X^2_K_2(X)]| ≤[|XK12(X)−XK22(X)|] [|X_K_1^2(X)-X^2_K_2(X)|] ≤[2|XK1(X)−XK2(X)|] [2|X_K_1(X)-X_K_2(X)|] ≤α4|K1−K2| ≤ α4|K_1-K_2| Therefore the mean and variance are Lipschitz continuous in K, which means there is a unique solution to the IVP K′=h(K),K(0)=0K =h(K),K(0)=0. Under the limit as K→∞K→∞, [XK2(X)]→1E[X^2_K(X)]→ 1, so Var(K)=[XK2(X)]−[XK(X)]2→0Var(K)=E[X^2_K(X)]-E[X_K(X)]^2→ 0. The left and right limits of h are h(0)>0,limK→∞h(K)=−2c<0, h(0)>0, _K→∞h(K)=-2c<0, so by the intermediate value theorem (h is Lipschitz continuous), there exists a K∗=infK>0|h(K)=0K^*= \K>0\;|h(K)=0\. For all K<K∗,h(K)>0K<K^*,\;h(K)>0 therefore K is increasing. Therefore, K→K∗K→ K^*. ∎ Proof of Proposition 4.1 Proof. Let f∈C2([0,1])f∈ C^2([0,1]) be an arbitrary test function, then dt[f(Xt)]=[ℒf(Xt)] ddtE[f(X_t)]=E[Lf(X_t)] where ℒL is the infinitesimal generator given by ℒf(x)=Aρ(t,x)f′(x)+12Bρ2(t,x)f′(x). (x)=A_ρ(t,x)f (x)+ 12B_ρ^2(t,x)f (x). In particular, letting f(x)=xf(x)=x gives the evolution of the mean as dt[Xt]=[Aρ(t,x)] ddtE[X_t]=E[A_ρ(t,x)] which gives the ODE dtm1(t)=∫01Aρ(t,x)ρ(t,x)x ddtm_1(t)= _0^1A_ρ(t,x)ρ(t,x)dx =2α∫01x2(1−x)2ΔG[ρ]ρ(t,x)x+O(α2). =2α _0^1x^2(1-x)^2 G[ρ]ρ(t,x)dx+O(α^2). To show that this is increasing, we begin by showing the O(α2)O(α^2) term is uniformly bounded by a constant C. Each variance term is bounded by Var(Uh) (U^h) ≤[(Uh)2]≤(H(b+c)+|β|)2 [(U^h)^2]≤(H(b+c)+|β|)^2 and the covariance is bounded by |Cov(Uh,hk)| |Cov(U^h,h^k)| ≤Var(Uh)Var(Uk)≤(H(b+c)+|β|)2. ≤ Var(U^h)Var(U^k)≤(H(b+c)+|β|)^2. Therefore, we have |ΣCC| | _C| =|α2∑h=0H−1Var(Uh)+2α2∑h<kH−1Cov(Uh,Uk)| = |α^2 _h=0^H-1Var(U^h)+2α^2 _h<k^H-1Cov(U^h,U^k) | ≤α2(H+H(H−1))(H(b+c)+|β|)2 ≤α^2(H+H(H-1))(H(b+c)+|β|)^2 =α2H2(H(b+c)+|β|)2 =α^2H^2(H(b+c)+|β|)^2 This give a simple bound on the O(α2)O(α^2) term as |2x(1−x)(1−2x)ΣCC| |2x(1-x)(1-2x) _C| ≤14α2H2(H(b+c)+|β|)2. ≤ 14α^2H^2(H(b+c)+|β|)^2. Define I(t) I(t) =∫01x2(1−x)2ΔG[ρ]ρ(t,x)x = _0^1x^2(1-x)^2 G[ρ]ρ(t,x)dx and let α∗=4I(0)H2(H(b+c)+|β|)2α^*= 4I(0)H^2(H(b+c)+|β|)^2. Since I(t)I(t) is continuous and I(0)>0I(0)>0, then there exists a T such that I(t)>I(0)/2>0∀t∈[0,T]. I(t)>I(0)/2>0 ∀ t∈[0,T]. Hence on the interval [0,T][0,T], dtm1(t) ddtm_1(t) =2αI(t)+O(α2) =2α I(t)+O(α^2) ≥αI(0)−14α2H2(H(b+c)+|β|)2 ≥α I(0)- 14α^2H^2(H(b+c)+|β|)^2 >0 >0 for all α<α∗α<α^*. ∎ To show the regularity conditions required for a well-posed solution to the steady-state equation, we provide the following Lipschitz continuity result on moments |mi(η)−mi(ν)| |m_i(η)-m_i(ν)| =|∫01yid(η−ν)(y)|≤iW1(η,ν). =| _0^1y^id(η-ν)(y)|≤ iW_1(η,ν). Next, we use this moment regularity to show the conditional reward terms are also Lipschitz. Lemma B.9. The drift AηεA _η and diffusion Bη2,εB^2, _η are Lipschitz continuous in both x and η. Specifically, ‖Aηε−Aνε‖∞≤LAW1(η,ν),‖Bη2,ε−Bν2,ε‖∞≤LBW1(η,ν), \|A _η-A _ν\|_∞≤ L_AW_1(η,ν), \|B^2, _η-B^2, _ν\|_∞≤ L_BW_1(η,ν), and |Aηε(x)−Aηε(y)|≤LA|x−y|,|Bη2,ε(x)−Bη2,ε(y)|≤LB|x−y|, |A _η(x)-A _η(y)|≤ L_A|x-y|, |B^2, _η(x)-B^2, _η(y)|≤ L_B|x-y|, where LA,LB∈ℝ+L_A,L_B ^+ and are independent of ε . Proof. First note that ΔG[ρ] G[ρ] is a polynomial in m1m_1 and m2m_2, and Σ is a polynomial in x,m1,m2,m3x,m_1,m_2,m_3. Since the product of polynomials is a polynomial, we have Aηε=f(x,m1,m2,m3)A_η =f(x,m_1,m_2,m_3) is a polynomial such that x,m1,m2,m3∈[0,1]x,m_1,m_2,m_3∈[0,1]. Therefore, f is Lipschitz in each of its arguments such that Lf=sup[0,1]4‖∇f(x,m1,m2,m3)‖∞<∞. L_f= _[0,1]^4\|∇ f(x,m_1,m_2,m_3)\|_∞<∞. Hence, |f(x,m1(η),m2(η),m3(η))−f(x,m1(ν),m2(ν),m3(ν))| |f(x,m_1(η),m_2(η),m_3(η))-f(x,m_1(ν),m_2(ν),m_3(ν))| ≤Lf(|m1(η)−m1(ν)|+|m2(η)−m2(ν)| ≤ L_f(|m_1(η)-m_1(ν)|+|m_2(η)-m_2(ν)| +|m3(η)−m3(ν)|) +|m_3(η)-m_3(ν)|) ≤Lf(W1(η,ν)+2W1(η,ν)+3W1(η,ν)) ≤ L_f(W_1(η,ν)+2W_1(η,ν)+3W_1(η,ν)) =6LfW1(η,ν) =6L_fW_1(η,ν) For Lipschitz continuity in x, note that since f is a polynomial, its derivative is bounded. In particular, sup[0,1]4|∂xf|<∞ _[0,1]^4| _xf|<∞ and therefore for any η, |Aη(x)−Aη(y)|≤sup[0,1]4|∂xf||x−y| |A_η(x)-A_η(y)|≤ _[0,1]^4| _xf||x-y| so AηA_η is Lipschitz in x. Taking LAL_A as the maximum of each Lipschitz bound gives the result. The proof for Bη2,εB^2, _η is identical, since this is also a polynomial. ∎ Note that since the drift and diffusion are Lipschitz and the domain is bounded, then both the drift and diffusion are bounded above such that Aηε(x)≤MAA _η(x)≤ M_A and Bη2,ε(x)≤MBB^2, _η(x)≤ M_B. The next result shows the iterative operator is well-defined. Lemma B.10. The operator ℱε[η]:Sε→SεF [η]:S → S is well-defined. Proof. Note that |2Aηε(x)Bη2,ε(x)| | 2A _η(x)B^2, _η(x) | ≤|2Aηε(x)ε|≤2MAε ≤ | 2A _η(x) |≤ 2M_A so for every x∈[0,1]x∈[0,1] |∫0x2Aηε(y)Bη2,ε(y)y| | _0^x 2A _η(y)B^2, _η(y)dy | ≤2MAε ≤ 2M_A and e−2MAε≤exp(∫0x2Aηε(y)Bη2,ε(y)y)≤e2MAε e^- 2M_A ≤ ( _0^x 2A _η(y)B^2, _η(y)dy )≤ e 2M_A This means 0<wηε(x)≤e2MAε 0<w _η(x)≤ e 2M_A We have that wηε(x)∈C[0,1]w _η(x)∈ C[0,1], wηε(x)>0w _η(x)>0 and the integral across the domain is finite and strictly positive. Moreover, ∫01ℱε[η](x)x=1. _0^1F [η](x)dx=1. Therefore ℱε:Sε→SεF :S → S is a well-defined mapping. Since Bη2,ε(x)≤MB+εB^2, _η(x)≤ M_B+ , we also have the bound 1MB+εe−2MAε 1M_B+ e^- 2M_A ≤wηε(x)≤e2MAε ≤ w _η(x)≤ e 2M_A 1MB+εe−2MAε 1M_B+ e^- 2M_A ≤∫01wηε(x)x≤e2MAε ≤ _0^1w _η(x)dx≤ e 2M_A hence 0≤ℱε[η](x)≤Mε0 [η](x)≤ M for some constant MεM which is independent of η. ∎ Lemma B.11. The set ℱε(S):=ℱε[η]:η∈Sε⊂C([0,1])F (S):=\F [η]:η∈ S \⊂ C([0,1]) is relatively compact. Proof. From the Arzela-Ascoli Theorem [rudin_principles], it suffices to show that the set ℱε(Sε)F (S ) is equibounded and equicontinuous. Firstly, by Lemma B.10, the set is uniformly bounded by MεM . Now let ψηε(x) ψ _η(x) =∫0x2Aηε(y)Bη2,ε(y)y = _0^x 2A _η(y)B^2, _η(y)dy Then, |∂xψηε(x)| | _xψ _η(x)| =|2Aηε(x)Bη2,ε(x)| = | 2A _η(x)B^2, _η(x) | ≤2MAε ≤ 2M_A To pass the bound to the operator, we use that Bη2,ε(x)B^2, _η(x) is uniformly Lipschitz in x and has a lower bound ε to get |∂xlog(Bη2,ε)| | _x (B^2, _η)| =|∂xBη2,ε(x)Bη2,ε(x)|≤LBε. =| _xB^2, _η(x)B^2, _η(x)|≤ L_B . then since wηε(x)=eψηε(x)/Bη2,ε(x)w _η(x)=e^ψ _η(x)/B^2, _η(x), logwηε(x) w _η(x) =ψηε(x)−log(Bη2,ε(x)) =ψ _η(x)- (B^2, _η(x)) has a uniform Lipschitz bound. By Lemma B.10, the normalisation is bounded away from zero. Therefore, ℱε[η](x)F [η](x) has a uniform Lipschitz bound. Hence, ℱε[η](x)F [η](x) is uniformly bounded and equicontinuous and therefore relatively compact in C([0,1])C([0,1]). ∎ Proposition B.1. The mapping ℱε:Sε→SεF :S → S is Lipschitz continuous. Proof. Consider the integrand |2AηεBη2+ε−2AνεBν2+ε| | 2A _ηB^2_η+ - 2A _νB^2_ν+ | ≤|2(Aηε−Aνε)Bη2+ε|+2|Aνε||Bη2−Bν2||(Bη2+ε)(Bν2+ε)| ≤ | 2(A _η-A _ν)B^2_η+ |+ 2|A _ν||B^2_η-B^2_ν||(B^2_η+ )(B^2_ν+ )| ≤2|Aηε−Aνε|ε+2MA|Bη2−Bν2|ε2 ≤ 2|A _η-A _ν| + 2M_A|B^2_η-B^2_ν| ^2 ≤(2LAε+2MALBε2)W1(η,ν) ≤( 2L_A + 2M_AL_B ^2)W_1(η,ν) where in the final line we take the supremum over x. This implies ‖ψηε−ψνε‖∞ \|ψ _η-ψ _ν\|_∞ ≤(2LAε+2MALBε2)W1(η,ν) ≤( 2L_A + 2M_AL_B ^2)W_1(η,ν) Now take |wηε(x)−wνε(x)| |w _η(x)-w _ν(x) | =|exp(ψηε)Bη2,ε−exp(ψνε)Bν2,ε| = | (ψ _η)B^2, _η- (ψ _ν)B^2, _ν | ≤|exp(ψηε)−exp(ψνε)Bη2,ε|+|exp(ψνε)||Bη2,ε−Bν2,ε||(Bη2,ε)(Bν2,ε)| ≤ | (ψ _η)- (ψ _ν)B^2, _η |+ | (ψ _ν)||B^2, _η-B^2, _ν||(B^2, _η)(B^2, _ν)| ≤exp(maxψηε,ψνε)|ψηε−ψνε|ε+exp(ψηε)LBW1(η,ν)ε2 ≤ ( \ψ _η,ψ _ν\)|ψ _η-ψ _ν| + (ψ _η)L_BW_1(η,ν) ^2 ≤e2MAε(2LAε2+2MALBε3+LBε2)W1(η,ν) ≤ e 2M_A ( 2L_A ^2+ 2M_AL_B ^3+ L_B ^2 )W_1(η,ν) where in the penultimate line we use mean-value theorem, and in the final line the uniform upper bound on ψηεψ _η and Lipschitz bound above. Denote this value by CwC_w, then normalisation has the same Lipschitz bound. Hence |ℱε[η]−ℱε[ν]| |F [η]-F [ν] | =|wηε(x)∫01wηε(x)x−wνε(x)∫01wνε(x)x| = | w _η(x) _0^1w _η(x)dx- w _ν(x) _0^1w _ν(x)dx | ≤|wηε(x)−wνε(x)||∫01wηε(x)x|+|wνε(x)||∫01wηε(x)−wνε(x)dx||(∫01wηε(x)x)(∫01wνε(x)x)| ≤ |w _η(x)-w _ν(x) | | _0^1w _η(x)dx |+ |w _ν(x)|| _0^1w _η(x)-w _ν(x)dx|| ( _0^1w _η(x)dx )( _0^1w _ν(x)dx)| ‖ℱε[η]−ℱε[ν]‖∞ \|F [η]-F [ν] \|_∞ ≤(MB+ε)e2MAεCwW1(η,ν)+(MB+ε)2εe6MAεCwW1(η,ν) ≤(M_B+ )e 2M_A C_wW_1(η,ν)+ (M_B+ )^2 e 6M_A C_wW_1(η,ν) :=LℱW1(η,ν) :=L_FW_1(η,ν) ∎ Proof of Proposition 4.2 Proof. First note SεS is a non-empty, closed, bounded and convex subset of a Banach space. A solution to the steady state equation is uniquely expressed by the fixed point of the mapping ℱε:Sε→SεF :S → S . By Lemma B.10, this operator is well-defined and an invariant mapping. By Proposition B.1 it is continuous and by Lemma B.11 its relatively compact in C([0,1])C([0,1]). Applying Schauder’s Fixed Point Theorem [shapiro_fixed_point], there exists at least one fixed point of the mapping ℱε:Sε→SεF :S → S , and therefore a solution to the regularised steady state equation. ∎ Proof of Theorem 4.1 Proof. Let με(dx)=ρε(x)dxμ (dx)=ρ (x)dx. Since [0,1][0,1] is compact, the space of probability measures ([0,1])P([0,1]) is sequentially compact under weak convergence. Along any sequence εn→0 _n→ 0, there exists a subsequence and measure μ∈([0,1])μ ([0,1]) such that [billing_probability_measures] μεnk⇀μ. μ _n_k μ. The weak form of the regularised steady-state equation at the stationary distribution μεμ is ∫01[Aμε(x)φ′(x)+12(Bμε2(x)+ε)φ′(x)]με(x) _0^1 [A_μ (x) (x)+ 12(B^2_μ (x)+ ) (x) ]dμ (x) =0 =0 for all φ∈C2([0,1]) ∈ C^2([0,1]) such that φ′(0)=φ′(1)=0 (0)= (1)=0. Along the subsequence, we have ∫01[Aμεnk(x)φ′(x)+12(Bμεnk2(x)+εnk)φ′(x)]μεnk(x) _0^1 [A_μ _n_k(x) (x)+ 12(B^2_μ _n_k(x)+ _n_k) (x) ]dμ _n_k(x) =0. =0. Since the drift, Aμ(x)A_μ(x), and diffusion terms Bμ2(x)B^2_μ(x) are continuous and only depend on finitely many bounded moments, weak convergence implies Aμεnk(x)φ′(x)+12(Bμεnk2(x)+εnk)φ′(x)→Aμ(x)φ′(x)+12Bμ2(x)φ′(x) A_μ _n_k(x) (x)+ 12(B^2_μ _n_k(x)+ _n_k) (x)→ A_μ(x) (x)+ 12B^2_μ(x) (x) uniformly on [0,1][0,1]. We can take the limit as k→∞k→∞ and obtain ∫01[Aμ(x)φ′(x)+12Bμ2(x)φ′(x)]μ(x) _0^1 [A_μ(x) (x)+ 12B^2_μ(x) (x) ]dμ(x) =0 =0 which is the weak form of the stationary unregularised problem. ∎ Appendix C Empirical Study All simulations and numerical solutions have been run on a laptop with 12th Gen Intel(R) Core(TM) i7-1260P CPU. The total compute time for the experiments was approximately 40 hours. For the learning rate time scaling, the simulation and numerical solution for α=0.001α=0.001 was run for ten times longer than for the α=0.01α=0.01 case. Figure 3: Evolution of the strategy distribution where the population is initialised with X∼Uni(0,1)X (0,1) for the 4 rules. The theoretical solution (solid line) matches the simulations (histogram). Figure 4: Evolution of the strategy distribution where the population is initialised with X∼X Beta(3,3) for the 4 rules. The theoretical solution (solid line) matches the simulations (histogram).