Preface↩︎

Reinforcement learning (RL) has become a central paradigm for learning decisions directly from interaction with an uncertain environment. Its successes in games, robotics, operations, economics, and finance illustrate both the flexibility of the approach and the breadth of its possible applications. Although most of its early accomplishments concerned single-agent problems, many important environments are not single-agent systems. Traffic networks, energy markets, communication platforms, financial markets, and large online communities all involve many decision makers whose choices affect one another. A direct multi-agent description can quickly become too large to analyze or compute with, while a single-agent description misses the feedback loop between individual decisions and the collective behavior they generate.

Mean field reinforcement learning (MFRL) addresses this difficulty by combining RL, stochastic control, and mean field models. Its guiding idea is to replace a large population of weakly interacting agents by a representative agent coupled to the distribution of the population. This is more than a dimension reduction device. It changes the state of the decision problem: probability measures become part of the dynamics, conditional laws may have to be tracked, and common sources of randomness can keep the population correlated even in the infinite-population limit. These features are what make the subject both powerful and mathematically delicate.

There is now a substantial literature on multi-agent reinforcement learning (MARL), mean field games (MFGs), mean field control (MFC), and mean field Markov decision processes (MFMDPs). Yet these areas are often presented with different conventions, different levels of mathematical detail, and different algorithmic priorities. This monograph is meant to help close that gap. It develops a self-contained discrete-time route from rigorous MFC models to concrete RL algorithms, with particular care for the distinction between finite populations and their mean field limits, the role of common noise and common randomization, and the relation between open-loop, closed-loop, and randomized controls.

The exposition is mathematical, but it is guided by computation. We do not seek the greatest mathematical generality. Instead, the emphasis is on formulations that are precise enough to derive dynamic programming principles, propagation of chaos arguments, and convergence statements, while remaining close to implementable learning procedures.

The monograph is organized as follows. It first places MFRL in the broader RL and MARL context, and recalls the probabilistic and control-theoretic tools needed for their analysis. We then study two complementary models that cover different aspects of the framework, and analyze them at a rigorous mathematical level. A first abstract model with common noise is used to clarify the different roles of open and closed loop policies and the relation between finite-population systems and representative-agent formulations. It also connects optimization by a central planner with the lifted MFMDP viewpoint. Then, a linear quadratic model provides an explicit benchmark in which stability, optimality, and policy gradient methods can be analyzed in detail. The final part proposes numerical implementations of these formulations: tabular and discretized algorithms for the abstract model, and policy gradient methods with several simulator viewpoints for the linear quadratic model.

The intended reader can be a graduate student or an early researcher in applied mathematics coming from probability, stochastic control, reinforcement learning, or multi-agent systems. The guiding message is that MFRL is not simply RL with a larger state space, nor MFC with numerical experiments appended to it. It is a framework in which large population limits, information structures, and learning procedures have to be designed together.

This monograph is an expanded version of the lecture notes prepared by the authors for a set of tutorial lectures at the University of Chicago NSF Institute for Mathematics, Statistics and Innovation (IMSI) in March 2026, in preparation for the week-long workshop Theoretical Foundations of Multi-agent and Mean Field Reinforcement Learning held at IMSI May 18-22, 2026.³

Acknowledgements:

We would like to thank IMSI for a superb job at hosting the workshop. Our gratitude also goes to Zongjun Tan with whom we had numerous interactions during his PhD studies at Princeton, and with whom we co-authored several papers on RL and MFRL; some of the results of these papers are used in this monograph. R.C. was partially supported by AFOSR grant FA-9550-23-1-0324.

1 Stochastic Control and Reinforcement Learning Landscape↩︎

1.1 Multi Agent Reinforcement Learning versus Mean Field Reinforcement Learning↩︎

1.1.1 Introduction↩︎

The relevance and economic importance of applications of sequential decision making by multiple interacting agents has emerged as a major impetus for the development of the theory of Multi-Agent Reinforcement Learning (MARL). The survey-style papers [1] and [2] give an overview of the challenges and some of the existing algorithms involved. Noticeably, [3] offers an extensive bibliography of \(320\) references. The applications come from domains as diverse as autonomous vehicle coordination, robotic swarms, financial markets, smart grids, and distributed resource allocation.

While the classical theory of single-agent Reinforcement Learning (RL) focuses on the approximation of agent optimal policies that depend only on the state of the agent, their own actions and their subsequent potential rewards, and possibly some random shocks, the search for optimal policies in multi-agent systems is much more involved. Indeed, each agent’s decision influences the time evolution (classically prescribed by transition probabilities) of the states and the rewards of the other agents, creating intricate dependencies that compound the complexity of MARL models. In fact, the theoretical complexity of the model is not the only roadblock to a straightforward generalization of the theory and algorithms developed for single agent RL. The formulation of the optimization problem becomes a delicate matter as the computational complexity tends to scale with the number of agents. As modeling practical applications increasingly requires the control of tens, hundreds, or even thousands of individual agents, the need for sound models of the interactions between the agents, and for scalable and theoretically grounded algorithms becomes of critical importance.

The goal of this chapter is to review the pros and cons of the theoretical models and their implementations as they appear in the published literature for finitely many agents RL. Highlighting some of their limitations will serve as a motivation for the introduction of Mean Field Reinforcement Learning (MFRL) models, to which this monograph is devoted. Our contention is that these models provide tractable solutions to some of the main difficulties encountered in the analysis of large populations of agents.

We first recall separately the main definitions and properties of MARL and MFRL. We then highlight the main differences in a comparison intended to firmly establish the basis of our mathematical treatment of MFRL in the rest of the monograph.

Note to the reader: while most economists and engineers set up models so that the optimization part is the maximization of some form of utility or reward, most of the models treated by mathematicians involve the minimization of costs. Since these two forms of models are mathematically equivalent, we choose, because of our own backgrounds and biases, the latter approach. Throughout the monograph, we shall minimize costs instead of maximizing rewards, unless this choice is detrimental to comprehension of the application.

1.1.2 Multi-Agent Reinforcement Learning↩︎

This section is intended to set the stage for a rigorous mathematical discussion of the foundations and paradigms of MARL as found in the published literature on the subject.

1.1.2.1 Problem Formulation

The premise of MARL is to extend the technology developed in classical single agent reinforcement learning to \(N\ge 2\) agents interacting within a shared environment. The natural generalization is to replace the standard mathematical setting of a Markov Decision Process (MDP) by some sort of joint MDP, and because of the multivariate nature of the optimizations, the most convenient way to formulate the model is to treat it as a stochastic game with \(N\) players. Although attractive at first, this option requires specifying the kind of information available to each agent when making their decisions. The mathematical set-up we use to describe the models includes the following components:

• A joint state space \(S\) that captures the global system configuration. It will be assumed to be a nice measurable space, typically a Polish space.

• Individual action spaces \(A^1,\cdots, A^N\), one for each agent, which we shall assume to be closed convex subsets of a finite dimensional Euclidean space, say \(\RR^k\). We shall also use the space of action profiles, that is the product space \(A = A^1\times\cdots\times A^N\).

• A state transition probability kernel \(p(\,\cdot\,|s,\ba)=\PP[\cdot\,|\,s, a^1, \cdots, a^N]\) which for each state \(s\in S\) and joint action profile \(\ba=(a^1,\cdots,a^N)\) is a probability measure on \(S\). We denote by \(\cP(S)\) the space of probability measures on \(S\). Notice that when the state space \(S\) is finite, say with \(n\) elements, \(\cP(S)\) can be identified with the \(n-1\) dimensional simplex \[\Delta_{n-1}=\{(p^1,\cdots,p^n)\,|\, \sum_{i=1}^np^i=1, \; p^i\ge 0,\,i=1,\cdots,n\}.\]

  To conform with the vast majority of existing articles on the subject, we define state dynamics with a transition probability. Still, as we explain in the following subsection, the dynamics of the state can alternatively be given by a system function. We explain below that these two formulations are essentially equivalent and that we shall choose one over the other when convenient.

• Individual reward functions \(S\times A\ni (s,a^1,\cdots,a^N)\mapsto r^i(s, a^1,\cdots, a^N)\) that are typically assumed to be bounded measurable functions of the state and the entire action profile (i.e. the actions of all the agents);

  Most of the MARL literature relates to models in which rewards are maximized, hence our choice to introduce reward functions in this section. But as explained at the end of the introduction, despite the fact that the early RL literature maximizes rewards, we shall switch to the mathematical tradition to minimize cost functions, and starting from the next chapter, we shall use cost functions \(f^i\) instead of reward functions \(r^i\).

• Individual policies \(\pi^i: S \to \cP(A^i)\) that each agent uses to choose which action to take in order to optimize their objective. A reader familiar with the stochastic control or MDP literature could expect that the policies are strict, i.e. measurable functions \(\pi^i: S \to A^i\), however, as it is standard in the theoretical models and the implementation algorithms of RL, one privileges mixed policies, also called randomized policies for which, given the state \(s\in S\) of the system, agent \(i\) will choose an action \(a^i\in A^i\) drawn at random from the probability distribution \(\pi^i(s)\in\cP(A^i)\).

Given this set-up, each agent \(i\) aims to maximize (resp. minimize) their expected cumulative discounted reward (resp. cost) \[\label{fo:ith95objective} J^i(\pi)=\EE\Bigl[\sum_{n\in\cN} \gamma^n r^i(X_n,a^1_n,\cdots,a^N_n)\Bigr]\tag{1}\] while other agents simultaneously optimize their own objectives. This formulation encompasses cooperative settings (where agents share a common reward), competitive scenarios captured by stochastic game models as well as mixed-motive environments. We shall give more details and references for the differences between these models. Here

• \(\cN\) is the time horizon over which the optimization is performed. \(\cN=\{0,1,\cdots,T\}\) for a finite integer \(T\) for finite horizon models, and \(\cN=\{0,1,\ldots\}\) for infinite horizon models.

• \(\gamma\in(0,1]\) is a discount factor that is used to actualize the values of the reward or the cost over time. It is often chosen to be \(1\) for finite horizon models, and it is most often chosen to be strictly smaller than \(1\) for infinite horizon models, helping with the convergence of the infinite summation.

• for each \(i\in\{1,\cdots,N\}\), \(a^i_n\sim\pi^i(X_n)\) is the action taken by agent \(i\) at time \(n\) by sampling from the probability distribution \(\pi^i(X_n)\), and \(X_n\) represents the random realization of the state of the system at time \(n\).

According to our introduction of the transition probability \(p(\,\cdot\,|s,\ba)\), the time evolution of the state of the system is given by a stochastic process \(\bX=(X_n)_{n\in\cN}\) whose dynamics are given by \[X_{n+1}\sim \int_A p(\,\cdot\,|\,X_n,\ba)\,\pi^1(da^1|X_n)\cdots\pi^N(da^N|X_n)\] where throughout the monograph we use the notation \(Y\sim\mu\) to state that the random variable \(Y\) is a sample realization from the probability distribution \(\mu\), implying that \(\mu\) is the distribution of \(Y\).

1.1.2.2 Transition Probability or System Function?

It is often convenient to define the time evolution of the state of the system using a system function \(F:S\times A\times E\to S\) and a dynamical equation \[\label{fo:system} X_{n+1}=F(X_n,\ba_n,\epsilon_{n+1})\tag{2}\] for an independent identically distributed ( i.i.d. for short) sequence \((\epsilon_n)_{n\ge 1}\), of innovation random variables taking values in a measurable space \(E\). This is the approach used in [4], [5] and [6] to analyze the optimization of mean field MDPs with general open loop and closed loop Markovian control processes, even in the presence of common noise. Because of their ease of implementation, and because of their ability to preserve the desirable Markov property, closed loop policies have been studied and implemented to a larger extent than the more general open loop policies whose implementation in practical applications is problematic. We shall limit ourselves to closed loop policies in this introductory tutorial. We shall provide a comparative study of the optimizations with respect to open and closed loop policies in Chapter 2.2.

If the dynamics of the state are prescribed via a system function like in 2 , it is straightforward to retrieve the transition probability. Indeed, in this case \[\label{fo:transition} p(\,\cdot\,|\,x,\ba)=\PP[F(x,\ba,\epsilon_{n+1})\in\,\cdot\,],\tag{3}\] in other words, the probability \(p(\,\cdot\,|\,x,\ba)\) is merely the image (i.e. push-forward) of the common law of the random shocks \(\epsilon_{n+1}\) under the mapping \(E\ni\epsilon\mapsto F(x,\ba,\epsilon)\).

On the other hand, if the dynamics of the state are given by a transition probability and we wish to study the model as if it were given by a system function, the most convenient way to introduce the appropriate system function is to use the Blackwell-Dubins theorem [7]. This important result from measure theory is a generalization of the famous Skorohod representation theorem, which is particularly well suited to our needs in developing a rigorous mathematical theory of MDPs and MFRL. We will give a precise statement in Lemma 1.1 for later references.

1.1.2.3 Key Differences between MARL and Single Agent RL

MARL differs fundamentally from single-agent RL in several critical aspects which we briefly highlight below. The main differences arise from the multitude of possible information structures for individual agents. Specifically, the main determining factor in the construction of a model and the development of optimization algorithms is the answer to the following question:what kind of information is available to an individual agent making a decision concerning the best action they should take? Usually, individual agents are cognizant of their own individual states, but do they know anything about the states of the other agents and of the common sources of random shocks? In order to use the information available to them, do agents collaborate? form coalitions? compete among themselves? Also, when it comes to deciding what action to take, should the individual agent perform their own optimization? or should they implement policies determined by a central planner? These are some of the typical issues of multivariate optimization theory and multi-objective decision theory. Not surprisingly, they complicate the development of a uniform theory for MARL based on a straightforward generalization of the successes of the single agent theory. Here we mention a small number of obvious consequences of the plethora of differences with the single agent theory.

• Growth of the Size of the Joint State-Action Space: The effective state-action space grows exponentially with the number of agents. Even if we assume that all the agents can choose their actions from a common finite space, say that all the action spaces \(A^i\) are equal and have \(k\) elements, the size of the set of action profiles, i.e. the joint action space, is \(k^N\), leading to exponential growth of the complexity of the value functions and the policies.

• Non-Stationarity: In most MARL models, the optimization problem faced by each individual agent is more complex than in single agent RL models. This extra complexity comes in no small part from the non-stationary learning environment created by the other agents concurrently updating their own policies during learning. From each individual agent’s perspective, the time evolution of the environment and reward structure they are experiencing shifts as other agents learn, violating the stationarity assumption underlying classical convergence guarantees in single agent RL.

• Decentralization and Partial Observability: Many MARL models impose decentralized execution, requiring individual agents to choose their actions on the sole basis of the information they have access to, typically local observations without access to global state information or other agents’ private observations. This serious restriction on the nature of the admissible individual policies is akin to the problem of partial observation in stochastic control, and it introduces additional complexity.

Beyond these structural differences, MARL also faces challenges in credit assignment (i.e., determining each agent’s contribution to the aggregate reward) which becomes prohibitive in large cooperative games, as exemplified by the computation of Shapley indices. Models requiring communication between the agents add yet another layer of complexity; see, e.g., [8] and [9].

1.1.2.4 Major MARL Paradigms

As suggested earlier, MARL frameworks are similar to those of multiple player games. Indeed, a crucial component of the definition of a MARL model is the description of the interactions between the agents. Deciding which types of interaction are admissible is of crucial importance because the nature of these interactions plays a fundamental role in the theoretical analysis of these models, as well as the development of algorithmic solutions to address the challenges raised by the multitude of agents trying to optimize simultaneously their respective objectives. The following bullet points describe in broad strokes some of the most frequently found types of interaction in the published literature.

• Independent Learning: [10] is an early example of decentralized RL in cooperative multi agent systems where a team of learning robots tries to coordinate their individual behavior to reach a coherent joint behavior. The authors of the paper propose an extension of Q-learning which they call hysteretic Q-learning, algorithm that does not require any additional communication between robots who can learn independently of each other and maximize an aggregate reward. More generally, independent learning covers multi agent models in which each agent treats the actions and rewards of the other agents as part of the environment, and learns independently using single-agent RL algorithms. While simple to implement, this strategy suffers from non-stationarity and falls short of providing expected results when the interactions between the agents are strong.

• Centralized Training with Decentralized Execution (CTDE): This paradigm leverages centralized information (such as global state or other agents’ actions) while deploying decentralized execution policies. [11] and [12] are typical examples of this strategy. However, this approach does not always resolve the scalability issue when the number \(N\) is large. Indeed, algorithms still depend on joint actions and global state information.

• Communication-Based Methods: These approaches enable agents to exchange messages or learn communication protocols to facilitate coordination under partial observability. Although effective for small teams, communication overhead and the challenge of learning meaningful protocols limit scalability to large populations. See, for example, [8] for algorithms to learn communication protocols and still fit into the framework of centralized learning and decentralized execution. In a separate line of attack, [13] develops an algorithm to help agents learn when to communicate based on scenarios, and profitability in cooperative and competitive MARL settings. See also [14] for a policy gradient algorithm skipping some steps to avoid communication overheads.

The underpinnings of the first two bullet points will resurface in our discussion of Mean Field RL strategies, though in a slightly different form, allowing us to use an important extra symmetry assumption.

1.1.2.5 The Scalability Bottleneck

MARL algorithms face fundamental scalability barriers when the number of agents becomes too large. Even if individual action spaces have the same finite cardinality, say \(k=|A^i|\) for all \(i\in\{1,\cdots,N\}\), and even if the size of the state space does not grow with \(N\), the centralized policies and the value functions are mappings over spaces of size at least \(k^N\). Even with function approximation, using, for example, neural networks, representing and learning such functions becomes prohibitively expensive. Indeed, the number of samples required to learn accurate value function estimates or policies grows exponentially with the dimension of the joint state-action space. Also, storing value tables or neural network parameters scales poorly with the number of agents, and memory becomes a serious issue. This is a major impediment for applications such as traffic control, robotic swarms, or crowd simulation, where hundreds or thousands of agents interact, rendering conventional MARL approaches impractical.

1.1.3 Mean Field RL as a Scalable Alternative↩︎

The idea of mean field approximation of a large system of agents acting similarly (at least in a statistical sense) comes from statistical physics, and especially the study of large systems of identical particles with symmetric interactions. In MFRL, the way each agent is affected by, and interacts with, the other agents is captured by aggregate quantities defining the behavior of the other agents, these aggregate quantities being quantified using the statistical distributions of the state of the system and the other agent actions. So instead of a complicated network of interactions, each individual agent interacts with a probability distribution, and the symmetry assumption is used to guarantee that this distribution does not depend upon the particular agent that we are tagging for the analysis.

The way we envision MFRL is captured by the statement centralized learning leading to decentralized execution. The policies of the \(N\) agents are coordinated by a central controller. See [15], [16]

1.1.3.1 Core Concept and Theoretical Foundation

As argued in [3], [17]–[19], MFRL addresses the scalability challenge of MARL by replacing explicit modeling of \(N\) individual agents with a representative agent that interacts with a probability distribution—the mean field. Such a clean formulation of the mathematical set-up comes from the theory of mean field games [20]–[22].

At an intuitive level, one should argue that in large homogeneous populations in which interactions between agents are either weak or asymptotically symmetric, the effect of all other agents on a given agent can be well-approximated by sufficient statistics of the population distribution. Instead of keeping track of the state of the system and the behavior of the other \(N-1\) agents, a good approximation of the system is to assume that each individual agent interacts with a mean field \(\mu_n\in\cP(S \times A)\), which represents the empirical distribution of the states and the actions of the agents at time \(n\). This reformulation is justified by the propagation of chaos which has been known for a long time in statistical physics. See for example [22] for a detailed proof for system dynamics driven by stochastic differential equations, or Proposition [pr:chaos] in Chapter 2.1 for a version tailored to the abstract RL frameworks of this monograph.

At a minimum, the formulation of a MFRL should include:

• A Representative Agent: A single agent whose state dynamics and rewards depend on their own state-action pair \((s, a)\) and the mean field \(\mu\);

• The Mean Field Dynamics: The time evolution of the population distribution which we can formalize as \(\mu_{n+1} = \cT(\mu_n,\pi)\), where \(\cT\) is a form of transition operator and \(\pi\) is a policy:

• The Optimization of a Mean Field MDP: The representative agent optimizes a single-agent MDP with state space \(S\) and action space \(A\), parameterized by the mean field \(\mu\).

When the solution of a MFRL reformulated in this way is used as an approximation of an \(N\)-agent MARL, the \(O(N^2)\) complexity of the pairwise interactions is dramatically reduced into a single agent–mean-field interaction, drastically reducing the computational burden [15], [18].

1.1.3.2 Key Assumptions

In order to reinforce the highlights of the previous subsection, we emphasize once more that the limiting argument at the core of the MFRL’s tractability relies on the following assumptions:

• Agent Homogeneity and Symmetry: in a statistical sense, agents are assumed to be exchangeable or statistically permutation-invariant, implying that the identity of individual agents is irrelevant in the sense that only the probability distribution of the state-action couples matters. This is the basis on which we can characterize the agent population by a single distribution \(\mu\) rather than by the state-action pairs of \(N\) individual agents.

• Weak Coupling: Implied in part by the above assumption is the fact that the influence of any single agent on the others is sufficiently small that aggregate effects dominate. This assumption justifies rejection of individual agent identities and focus on population-level statistics.

• Large Population Asymptotics: A limiting theorem of the type known as propagation of chaos must hold. While the de Finetti extension to exchangeable systems of the classical law of large numbers ensures that the empirical distribution of the state of the system and the actions of the agents converges to the mean field distribution \(\mu\), and individual agent fluctuations become negligible, the propagation of chaos theoretical result describes how the dynamical equation of the state of the finite agent system needs to be modified to describe the dynamics of the state interacting with the limiting mean field.

1.1.3.3 Mean Field Approximation

In the mathematical mean field theory of stochastic systems, whether we work on a Mean Field Game (MFG) model or a Mean Field Control (MFC) problem, the basic approximation result is a form of propagation of chaos formalization and quantification. See, for example, [21] or [22] for precise statements. Typically, the mean field limit justifies the replacement of the original \(N\)-agent system with a limiting system where each agent faces an evolving mean field which can be deterministic or stochastic in the presence of a common noise, i.e. random shocks affecting all the agents. The rate of convergence in the various forms of the propagation of chaos theorems provides a quantification of the approximation error based on the difference between characteristics of the \(N\)-agent system and their corresponding mean field limits. These guaranteed convergence rates are typically of the order of \(O(1/\sqrt{N})\) under suitable regularity conditions. See, for example, [23] or [24] for specific examples.

As expected, attempts to prove that this rate of convergence still holds in the added presence of learning appeared in the MFRL literature. For example [16] proposes a mathematical set up to approximate cooperative MARL models using a MFC approach, and shows that the approximation error is still of the order of \(O(1/\sqrt{N})\). This type of analysis requires detailed control of the various approximations used in the algorithmic implementations of the theoretical limiting MFRL model to connect it to the original \(N\)-agent cooperative MARL. Typically, successful analyses do estimate the mean field:

• by approximating the empirical distributions by histograms or kernel density estimates of the agent states and actions;

• using subsampling to approximate the mean field by appropriately selected subsets of agents, with explicit approximation error bounds.

See, for example, [6], [25], [26] for instances of implementations of this strategy in specific models and the second part of this monograph devoted to numerical implementations.

1.1.4 MARL vs. MFRL: Key Differences and Trade-offs↩︎

Faced with the task of solving a practical problem, the choice between MARL and MFRL depends on the characteristics of the problem and modeling priorities:

MARL offers higher modeling freedom for problems with a very small number of agents, especially when the population is heterogeneous and interactions are agent specific. However, it faces insurmountable computational barriers at scale. On the other hand, MFRL may sacrifice some faithfulness to some of the population characteristics by assuming exchangeability (and weak coupling), but achieves tractable learning and stronger theoretical guarantees in large-population settings. For limitations of mean-field approximations under weak assumptions, see e.g. [27]; for partial observability of the mean field, see e.g. [28].

1.1.4.1 Application Domains

MFRL emerged to make learning and control feasible in settings with large populations of agents by trading exact replication of multi-agent interactions for a tractable population-level description. This trade-off enables both algorithmic scalability and cleaner theoretical analysis, addressing the fundamental limitations of conventional MARL in large-population regimes.

MFRL is particularly well-suited to models for which assuming that agents are exchangeable (in a statistical sense) is reasonable. When agents are in large numbers (and as a result MARL algorithms fail to deliver reasonable approximations), modeling the interactions with aggregate characteristics provides realistic results. Such models have been used for the following applications:

• Robotic Swarms: Coordinate hundreds of simple robots for coverage, formation control, or collective transport. See, for example, [29] for an instance of MFRL application to flocking, possibly with multiple populations, or the more recent [30] providing numerical comparisons with applications of state-of-the-art MARL algorithms to Kuramoto and Vicsek swarming models; see also [31]–[35].

• Traffic and Transportation: Managing vehicle routing, traffic signal control, and ride-sharing with thousands of participants. See, for example, [36], [37] for an application to the vehicle positioning problem; related mean field models for routing and traffic include [38]–[42]. See also [43] for the use of off-line algorithms for crowd navigation purposes.

• Smart Grids and Resource Allocation: Distributed energy management and load balancing across many devices. See, for example, [44] for a hybrid model with an MFC and MFG framework for integrating Distributed Energy Resources (DER) aggregators into wholesale electricity markets, or [45] for a note containing an application to demand response in smart grids. Other mean field game/control applications to energy systems include [46]–[51].

• Epidemiology and Social Dynamics: Modeling disease spread, opinion dynamics, or crowd behavior. See, for example, [52]–[61].

• Financial Markets: High-frequency trading and market making with many competing agents. See, for example, [62] for a model of dealer market, or [63] for an example of energy trading; see also [64]–[69].

• Security and communication networks: botnet defense, wireless networks, distributed caching, and interference control. See, for example, [70]–[75].

1.1.4.2 Theoretical Guarantees

On the theoretical front, MFRL takes full advantage of two decades of remarkable developments in stochastic analysis, optimal transport and analysis on Wasserstein spaces of probability measures, as well as the probabilistic approach to MFGs and MFC. In the following, we highlight some of these developments and their relevance to progress in MFRL.

• Dynamic Programming and Optimality Characterization: MFC formulations lead to dynamic programming principles and characterizations of value functions and optimal policies via Hamilton-Jacobi-Bellman equations and Pontryagin maximum principles in Wasserstein spaces of probability distributions. See for example [76] and [22] for a textbook form exposé.

• Approximation Bounds: Finite-population approximation results quantify the error between N-agent systems and mean field limits, typically showing \(O(1/\sqrt{N})\) convergence. See, for example, [23] or [24], and [22] for a complete proof of the sharp rate of convergence in the propagation of chaos.

• Convergence of Learning Algorithms: Rigorous bounds on learning algorithm complexity and convergence guarantees are highly technical and notoriously difficult to prove in classical MARL settings, and these difficulties are even greater in the MFRL setting. [6], [15], [16], [25], [77]–[86] are recent works establishing sample complexity bounds and convergence guarantees for mean field Q-learning, policy gradient methods, and actor-critic algorithms in the mean field setting.

• Extensions Beyond Homogeneity: The rare models of MFRL for which the fundamental assumption of homogeneity has been relaxed draw inspiration from the developments in the theory of MFGs. They include several sub-population (typically \(2\)) models, and major-minor agent formulations. See, for example, [87]–[91]. For graphon-based cooperative models with heterogeneous interactions, see [92]; for approximate-symmetry approaches, see [93].

• Inclusion of Constraints: Despite the presence of constraints, the error in the MFC approximation of a MARL can still be controlled in some specific models. See, for example, [94] for an example.

2 Mathematical Models↩︎

2.1 Mathematical Preliminaries and Model Foundations↩︎

2.1.1 Probability Spaces and Measure Theory↩︎

We first introduce the notations and the definitions used to define and characterize the stochastic dynamical systems we plan to optimize.

2.1.1.1 Basic Probability Structure

Throughout the monograph, unless stated otherwise, all the random quantities are assumed to be defined on a probability space \((\Omega, \cF, \PP)\), where:

• \(\Omega\) is the sample space,

• \(\cF\) is a \(\sigma\)-algebra of events,

• \(\PP\) is a probability measure on \((\Omega, \cF)\).

We denote by \(\EE[X]\) the expectation with respect to \(\PP\) of a random element \(X\), by \(\PP_X\) or \(\cL(X)\) its distribution (also called its law), and by \(\EE[X|\cG]\) and \(\PP_{X|\cG}\) or \(\cL(X|\cG)\) its conditional expectation and its conditional distribution given a sub-\(\sigma\)-algebra \(\cG \subseteq \cF\).

2.1.1.2 Borel and Polish Spaces

Our goal is to provide a rigorous mathematical presentation of concepts and theories that are too often discussed only at an intuitive level. For that reason, most of our arguments will be grounded in abstract measure theory. Throughout the monograph, a standard Borel space means a measurable space isomorphic to the Borel sigma-field of a Borel subset of a Polish space. For simplicity, whenever possible, we shall work directly with Polish spaces, namely complete separable metrizable spaces. They offer tractable analysis tools and are widely adopted as a convenient setting for stochastic control problems. For a Polish space \(S\), we denote by:

• \(\cB_S\) its Borel \(\sigma\)-field,

• \(\cP(S)\) the space of probability measures on \((S, \cB_S)\).

We implicitly assume that the space \(\cP(S)\) is equipped with the topology of the weak convergence of probability measures, according to which a sequence \((\mu_n)_{n\ge 0}\) in \(\cP(S)\) converges toward \(\mu\in\cP(S)\) if and only if \[\lim_{n\to\infty}\int_S \varphi \, d\mu_n = \int_S \varphi \, d\mu\] for every bounded continuous function \(\varphi: S \to \RR\). With this uniform structure, the space \(\cP(S)\) is itself a Polish space. We shall denote by \(d_w\) a distance for this Polish space structure, in other words, a distance for the weak convergence of probability measures.

2.1.1.3 Probability Kernels and Disintegration

In most reinforcement learning studies, the dynamics of the states are prescribed by probability kernels.

Given a probability measure \(\mu\) on \(X\) and a kernel \(\kappa\) from \(X\) to \(Y\), their composition \(\mu \circ \kappa\) is the probability measure on \(Y\) defined by \[(\mu \circ \kappa)(B) = \int_X \kappa(x, B) \, \mu(dx), \quad B \in \cY,\] and their product \(\mu \measprod \kappa\) is the probability measure on \(X\times Y\) defined by \[\int_{X\times Y}\varphi(x,y)(\mu \measprod \kappa)(dx,dy) = \int_X \mu(dx)\int_Y \varphi(x,y)\kappa(x, dy)\] for all bounded measurable functions \(\varphi\) on \(X\times Y\).

Conditional expectations and conditional probabilities are best understood and manipulated when they can be handled with probability kernels. If \(\cA\subset\cF\) is a sub \(\sigma\)-field of \(\cF\) and \(X\) is a random element in a measurable space \((S,\cS)\), by a regular conditional distribution of \(X\) given \(\cA\) we mean a version of the function \(\PP[X\in\cdot|\cA]\) on \(\Omega\times S\) which is a probability kernel from \((\Omega,\cA)\) to \((S,\cS)\), hence a \(\cA\)-measurable random probability measure on \((S,\cS)\).

We use the notation \(\PP[\cdot|Y]\) for the conditional probability given the \(\sigma\)-field \(\cF_Y\) generated by \(Y\), namely the subset of events depending only on \(Y\). Similarly \(\EE[\cdot|Y]\) stands for \(\EE[\cdot|\cF_Y]\). The following result should be interpreted as a generalization of the classical Fubini theorem of integration theory.

The following theorem is the most general form of disintegration of probability measures which we shall use.

We use these two forms of disintegration in the sequel.

2.1.1.4 Skorohod and Blackwell-Dubins Tools

The weak convergence of probability measures whose definition we recalled earlier is a fundamental tool in the analysis of the regularity and the asymptotic properties of stochastic structures. In the late fifties, Skorohod made the following discovery: if one is willing or able to change the probability space on which random elements are defined, weak convergence of random elements is not any different from almost sure convergence. This remarkable result turned out to be a powerful tool allowing all sorts of results from real analysis to be applied to random elements whose distributions converge weakly. To be specific:

This result is instrumental in the proof of convergence of value functions in stochastic control and reinforcement learning, and of continuity properties in the weak topology.

Generalizing some of the ideas at the heart of the proof of Skorohod’s representation theorem, Blackwell and Dubins developed a fundamental tool for constructing random variables with prescribed distributions. We shall use it repeatedly in the sequel.

We shall use such a function \(\rho_S\) and call it "the" Blackwell-Dubins function of the space \(S\) even if it is not unique. The atomless randomization space \(([0,1],\mathcal{B}([0,1]),\mathrm{Leb})\) can be replaced by another atomless probability space when needed.

This last point implies that if the state dynamics are given via a transition probability \(p(\,\cdot\,|\,x,\ba)\), the function \(S\times A\times E\ni (x,\ba,\epsilon)\mapsto F(x,\ba,\epsilon)=\rho_S(p(\,\cdot\,|\,x,\ba),\epsilon)\) if used as a system function would provide the same dynamical system.

This strategy was repeatedly used in [5] in mean field reinforcement learning with common noise.

2.1.2 Functional Analysis and Function Spaces↩︎

2.1.2.1 Spaces of Bounded Functions

For a Polish space \(S\), we denote by:

• \(\cB(S)\) the vector space of bounded Borel measurable functions \(\varphi: S \to \RR\). We shall implicitly equip this space with the supremum norm \(\left\|\varphi\right\|_\infty = \sup_{x\in S} \left|\varphi(x)\right|\), even though this space is not complete in general.

• \(\cC_b(S)\) the subspace of \(\cB(S)\) formed by bounded continuous functions. This is a Banach space for the sup norm.

• \(\cL_{sc}(S)\) the space of bounded lower semi-continuous functions.

Lower semi-continuity is preserved under pointwise suprema, and is crucial for establishing existence of optimal policies.

2.1.2.2 Contraction Mappings and Fixed Points

When the Dynamic Programming Principle (DPP) holds, optimal value functions are often characterized as solutions of an equation called the Bellman equation. Their existence and uniqueness are typically proven by identifying those solutions as fixed points of a contraction map. So if one can show that a certain Bellman operator is in fact a contraction on a Banach space of functions, existence and uniqueness of an optimal value function will follow. We shall use this strategy in Chapter 2.2 and we describe its steps in broad brushstrokes in Section 2.1.5 below.

2.1.2.3 Matrix Spaces and Norms

Linear Quadratic (LQ) models play an essential role in stochastic control, the theory of MDPs, and as a result in Reinforcement Learning (RL). The main reason is their simplicity and the fact that their optimization is solved via explicit (or at least semi-explicit) formulas. Their popularity spread from theoretical circles to real-world applications and practitioners. For this reason, one of the two main classes of models we consider in this monograph will be the class of LQ models. In such a framework, spaces of finite-dimensional matrices play an important role, so in preparation for their analyses, we introduce notations and definitions for some of the matrix spaces we use in the sequel. For convenience, we shall use several notions of norm. For a matrix \(U \in \RR^{d \times d}\), we shall use the:

• Spectral norm: \(\left\|U\right\| = \sigma_{\max}(U)\), the largest singular value,

• Frobenius norm: \(\left\|U\right\|_F = \left(\sum_{i,j} U_{ij}^2\right)^{1/2} = \left(\mathop{\mathrm{Tr}}(U^T U)\right)^{1/2}\).

We shall also use the Trace inner product: \(\left\langle U, V \right\rangle_{\text{tr}} = \mathop{\mathrm{Tr}}(UV^T)\).

2.1.3 Stochastic Dynamical Systems and Information Structures↩︎

We aim at optimizing systems evolving in time and for the purpose of implementation we shall assume that time is discrete. In other words, observations, measurements, actions, \(\ldots\) , are performed at regularly spaced times labelled by natural integers \(n=0,1,\cdots\). A decision at time \(n\) can only be made on the basis of information available at that time. Typically, the information structure is given by a filtration \(\FF=(\cF_n)_{n\ge 0}\), namely a non-decreasing sequence of sub \(\sigma\)-fields of \(\cF\): \(\cF_n\) contains the events known to have happened (or not) at time \(n\) or before. Since we investigate dynamical systems whose states are random, we rely on models based on discrete time stochastic processes.

We shall say that the process is a noise process if the random variables \(X_n\) are independent of each other. A central feature of both models subsequently studied in this monograph is the presence of two types of noise processes.

In models including both idiosyncratic and common noise processes, we shall assume that the sequences \((\varepsilon^i_n)_{n\geq 0, i=1,\cdots,N}\) and \((\varepsilon^0_n)_{n\geq 0}\) are mutually independent.

The goal of this dual noise structure is to capture both individual heterogeneity and aggregate uncertainty, which are essential features of many economic and engineering systems.

Filtrations and Information Structures

The information available to the agents evolves over time as captured by the choice of filtrations for information structure. These filtrations vary according to the model and the type of decision maker.

• Open loop policies: When the model is set up with open-loop controls, the \(\sigma\)-field \(\cF_n\) will typically include the realizations of all past random noise terms (both idiosyncratic and common, when appropriate) and the sources of randomness used to randomize past actions.

• Agent-level filtration: In the case of finitely many agents, the \(\sigma\)-field \(\cF^i_n\) representing the information available to agent \(i\) at time \(n\), includes their own state history (or the current value of their state in the case of Markovian policies), and possibly observations of the population distribution. In the mean field limit, the \(\sigma\)-field \(\cF^c_n\), representing the information available to a generic agent at time \(n\) includes the history up to time \(n\) of a generic state and observations of the population distribution.

For each of the specific models studied in the sequel, we shall specify explicitly the various sources of randomness, including noise sequences and randomization devices, used in the definitions of the filtration we use.

2.1.3.1 Dynamics Prescription and Objective Optimization

If the state of the system is given by a generic element \(x\) of a Polish space \(S\) and the actions agents are allowed to take are elements of a Polish space \(A\), the stochastic dynamics of the state of the system are given by a stochastic process \((X_n)_{n\ge 0}\) in \(S\). Historically, and especially in the existing RL literature, the rules governing the time evolution of the state are given by an initial distribution \(\mu_0\in\cP(S)\) and a sequence \((p_n)_{n\ge 0}\) of probability kernels from \(S\times A\) into \(S\). Intuitively speaking, \(\mu_0=\PP_{X_0}\) is the initial distribution of the state, and for each integer \(n\ge 0\), state-action couple \((x,a)\) and measurable set \(B\in\cB_S\), the quantity \(p_n\bigl((x,a),B\bigr)\) gives the probability that the state \(X_{n+1}\) is in the set \(B\) whenever (i.e. conditioned by the fact that) the state at time \(n\) is \(x\), i.e. \(X_n=x\) and the action taken at that time was \(a\in A\).

Mostly for mathematical convenience, some form of stationarity in time is assumed, so that the probability kernel \(p_n\) is independent of \(n\).

In order to define the dynamics of the state using a system function instead of a transition probability kernel, we use the Blackwell-Dubins function \(\rho_S\) of the Polish space \(S\) (recall Lemma 1.1) and define, independently of the time \(n\ge 0\), the function \(F:S\times A\times [0,1]\mapsto S\) by \[F(x,a,u)=\rho_S\bigl(p_n\bigl((x,a),\,\cdot\,\bigr),u\bigr),\] so that,

• if \((a_n)_{n\ge 0}\) is a prescribed sequence of actions in \(A\),

• if \((U_n)_{n\ge 0}\) is an i.i.d. sequence of uniform random variables in \([0,1]\),

• the stochastic process \((X_n)_{n\ge 0}\) defined by \[X_0=\rho_S(\mu_0,U_0)\quad\text{and}\quad X_{n+1}=F(X_n,a_n,U_{n+1}),\qquad n=0,1,\cdots\]

satisfies all the requirements prescribed above in the description of the dynamics of the state of the system. Since the original description of the state dynamics did not specify the information structure, the above use of the Blackwell-Dubins map allows us to specify the state dynamics by means of a sequence of independent random shocks (which we shall sometimes call innovations) and a (deterministic) function mapping a triple \((x,a,u)\) of a state, an action and a random shock into another state. As explained in the previous chapter, we call such a function \(F\) a system function.

When the state of the system is affected by only one controller, the goal of the latter is to choose a sequence \((a_n)_{n\ge 0}\) of actions in order to minimize an expected cost aggregated over time. In practice, the sequence \((a_n)_{n\ge 0}\) is random since at each time \(n\), the individual action \(a_n\) chosen by the controller, is picked as a function of the information available to the controller at that time. Specifying the exact nature of the available information will force the measurability of the \(A\)-valued random variables \(\alpha_n\) from which the actions \(a_n\) are sampled. So the objective function can be viewed as being of the form \[J = \EE\left[\sum_{n=0}^\infty \gamma^n f(X_n, \alpha_n, \epsilon_n)\right],\] where

• \(f:S\times A\times E \mapsto \RR\) is a bounded measurable function, \(f(x,a,\epsilon)\) giving the one stage cost incurred when the state is \(x\in S\), the action taken by the controller is \(a\in A\) and the realization of the noise (random shock) is \(\epsilon\in E\);

• \(\gamma\in (0,1)\) is a discount factor guaranteeing that the infinite sum converges,

the expectation being over the random values of the states \(X_n\), all the sources of randomness: initial conditions, the realizations \(\epsilon_n\) of the random shocks, and the random actions \(\alpha_n\) from which the actions \(a_n\) are sampled, \(\ldots\).

Multi-Agent Model Definitions

When the system involves \(N\) agents (\(N\) is a finite integer \(N\ge 1\)), each agent \(i\) chooses an action \(a^i_n\in A^i\) at time \(n\), the action space in which agent \(i\) chooses their action being possibly different from agent to agent. Accordingly, each agent may have a different one stage cost function and agent \(i\) aims at minimizing the cost \[J^i(\balpha) = \EE\left[\sum_{n=0}^\infty \gamma^n f^i(X_n, \balpha_n, \epsilon^i_n)\right],\] where

• \((\epsilon^i_n)_{n\ge 0, i=1,\cdots,N}\) is an idiosyncratic noise process,

• \(\bA=A^1\times\cdots\times A^N\) is the set of action profiles,

• \(\balpha=(\balpha_n)_{n\ge 0}\) is a stochastic process in the space \(\bA\) of action profiles,

• \(f^i:S\times \bA\times E \to \RR\) is a bounded measurable function, \(f^i(x,\ba,\epsilon)\) giving the one stage cost incurred by agent \(i\) when the state is \(x\in S\), the actions taken by all the agents are \(\ba=(a^1,\cdots,a^N)\in \bA\) and the realization of the noise (random shock) is \(\epsilon^i\in E\);

• and as before, \(\gamma\in (0,1)\) is a discount factor guaranteeing that the infinite sum converges.

Notice that the interactions between the agents enter the costs \(J^i\) through the values \(X_n\) of the state of the system and the actions \(a^j_n\) of all the agents. We mentioned earlier that we were assuming some form of stationarity in time allowing us to use one stage cost functions independent of the stage index \(n\). Now, and again purely for mathematical convenience, we shall assume a form of symmetry among the agents and use the same one-stage cost function \(f\) for all the agents. As emphasized earlier, this exchangeability between the agents is crucial in the models we want to study, and it will be especially useful when we consider the mean field limit \(N\to\infty\).

In many \(N\)-agent models, the state \(X_n\) of the system at time \(n\) is best described by an \((N+1)\)-tuple \(X_n=(X^0_n,X^1_n,\cdots,X^N_n)\) where the stochastic process \((X^0_n)_{n\ge 0}\) with values in a Polish space \(S^0\) captures stochastic factors common to all the agents, and for each \(i=1,\cdots,N\), the stochastic process \((X^i_n)_{n\ge 0}\) with values in a Polish space \(S^i\) captures stochastic factors specific to agent \(i\). When this is the case, \(X^i_n\) can be interpreted as the value at time \(n\) of the private state of agent \(i\), and the state space of the system can be factorized in the form \(S=S^0\times S^1\times\cdots\times S^N\). Accordingly, the system function \(F\) splits in a form \(F=(F^0, F^1,\cdots, F^N)\), the components being interpreted as individual system functions \(F^j\) for \(j=0,1,\cdots,N\) describing the time evolutions of the common factor and the individual private states \[X^j_{n+1}=F^j(X_n,a_n,\epsilon^j_n),\qquad n=0,1,\cdots,\quad j=0,1,\cdots,N.\] For the desirable assumption of symmetry among the agents to hold, we shall often assume that the individual system functions \(F^i\) are the same for \(i=1,\cdots,N\).

Competition or Cooperation?

Even though we are not yet ready to give the precise details of how the agents choose their actions, in other words the precise definitions of the filtrations \(\FF^i=(\cF^i_n)_{n\ge 0}\) with respect to which the control processes \(\balpha^i=(\alpha^i_n)_{n\ge 0}\) need to be adapted for the sequences of sampled actions \(\ba^i=(a^i_n)_{n\ge 0}\) to be admissible, we can still describe the major two forms of optimization.

Whether or not the individual system functions are identical, each agent \(i\) has their own objective cost \(J^i\) and this form of multivariate optimization requires formalizing the behavior of the agents.

If the agents compete among themselves and optimize their objectives independently, the most frequently searched outcome is given by the notion of Nash equilibrium.

Here (and in the following) we use the notation \((\balpha^{-i},\bbeta)\) for the control process profile where the \(i\)-th component of \(\balpha\) is replaced by \(\bbeta\). In other words, when the \(N\) agents use the control processes \(\balpha=(\balpha^i)_{i=1,\cdots,N}\), no single agent \(i\) can be better off by deviating alone from using their control \(\balpha^i=(\alpha^i_n)_{n\ge 0}\).

Another alternative is to assume that all the agents cooperate and aim at minimizing an aggregate or average cost \[\tilde{J}(\balpha)=\frac{1}{N}\sum_{i=1}^NJ^i(\balpha),\] which is often referred to as social cost. The factor \(1/N\) does not play any role when the number of agents is fixed. However, it will facilitate the control of the limit \(N\to\infty\) when we study large populations of agents.

For all the reasons already articulated in Chapter 1.1, the thrust of this monograph is the study of models of large populations of cooperative agents. As a result, we shall concentrate on cooperative equilibrium, and only mention Nash equilibria when appropriate, essentially to provide insightful comparisons.

Including a Common Noise

Many realistic situations cannot be modeled with idiosyncratic noise as the only source of random shocks. In order to add random factors and random shocks affecting simultaneously all the agents, we can use the flexibility of the system function to add an extra variable from another Polish space, say \(E^0\), and use an i.i.d. sequence of random variables \(\bepsilon^0=(\epsilon^0_n)_{n\ge 0}\) to include these common random effects in the dynamics of the state of the system. In this way, the dynamics of the private state of the \(i\)-th agent is now given by: \[\label{fo:dynamics95N} X^i_{n+1}=F^i(X_n,\alpha_n,\epsilon^i_n,\epsilon^0_n),\qquad n=0,1,\cdots.\tag{4}\] Note that the important Proposition [pr:chaos] below, referred to as the propagation of chaos is stated in this very setting including a common noise. Accordingly, the one-stage cost function of agent \(i\) can also be made dependent upon the common noise and read: \[\label{fo:cost95N} f^i(X_n,\alpha_n,\epsilon^i_n,\epsilon^0_n).\tag{5}\]

Randomization and Mixed Strategies

Whether the motivation is to help the theoreticians have more tools to prove existence of optimal strategies, or to help the individual agents explore better the state space in hope of finding more favorable regions with lower costs, randomization is a time honored trick used very frequently in stochastic optimization. In the context of continuous time stochastic control, it relies on the introduction of relaxed controls , altering the dynamics while keeping the same optimal values. In the discrete time setting considered in this monograph, the reinforcement learning practice is to choose actions by independently sampling from probability distributions on the action space. So instead of using adapted stochastic processes with values in the action spaces \(A^i\), the controls processes \((\alpha^i_n)_{n\ge 0}\) are replaced by adapted stochastic processes \((\pi^i_n)_{n\ge 0}\) with values in the spaces \(\cP(A^i)\) of probability measures on the action spaces \(A^i\), which are called mixed controls. Consequently, the action values \(\alpha^i_n\) actually used in the state update formula 4 and the cost evaluation formula 5 are samples of the probability distributions \(\pi^i_n\). Including sampling from probability distributions in the theoretical models naturally calls for the Blackwell-Dubins lemma, and indeed, the implementation of this idea of randomization in the theoretical analysis of the model will rely on this fundamental result.

2.1.4 Mean Field Models and their Limiting Behavior↩︎

In this section, we are still considering a multi-agent model with \(N\) agents, and we assume stationarity in time and symmetry among the agents as we mentioned earlier. First, we review the main properties which we require for the class of Mean Field models.

2.1.4.1 Symmetry Requirements

A first step in enforcing the symmetry among the agents is to assume that all the individual state spaces are the same, say \(S^i=S\), that all the action spaces are the same, say \(A^i=A\), that all the system functions are the same, i.e. \(F^i=F\), and all the one-stage cost functions are the same, i.e. \(f^i=f\) for all \(i\in[N]\). In fact, we require a stronger form of symmetry by assuming that the interactions between the individuals are symmetric in the following sense: for each agent \(i\in[N]\), we assume that the dynamics of the state \(X^i_n\) is given by equation 4 with a common system function \(F\), the dependence of the right hand side upon \(X_n=(X^1_n,\cdots,X^N_n)\) being given by a function of \(X^i_n\) and a symmetric function of the \(X^j_n\) for \(j\ne i\). Following the intuition that symmetric functions of a large number of variables are essentially functions of their empirical distribution (see for example [22] for a mathematical justification), we can imagine that because of the symmetry in the interactions between the states of the agents, the dependence of the system function \(F\) on \(X_n=(X^1_n,\cdots,X^N_n)\) could be a function of \(X^i_n\) and the empirical measure \[\mu^X_n=\frac{1}{N}\sum_{j=1}^N\delta_{X^j_n}.\] Similarly, one can imagine that the dependence of \(F^i\) on \(\alpha_n=(\alpha^1_n,\cdots,\alpha^N_n)\) could be through a function of \(\alpha^i_n\) and the empirical measure \[\mu^\alpha_n=\frac{1}{N}\sum_{j=1}^N\delta_{\alpha^j_n}.\] In fact, when we say that the agents are in mean field interaction, we shall mean that the dynamics of the private states of the agents are given an equation of the form \[X^i_{n+1}=F^i\Bigl(X^i_n, \alpha^i_n, \frac{1}{N}\sum_{j=1}^N\delta_{(X^j_n,\alpha^j_n)}, \varepsilon^i_{n+1}, \varepsilon^0_{n+1}\Bigr)\] for a system function \(F\) defined on the product space \(S\times A\times \cP(S\times A)\times E\times E^0\). De facto, the interactions between the agents are modeled via the empirical distribution of the state-action couples \((X^j_n,\alpha^j_n)\). Similarly, we can assume that the common individual one-stage cost functions \(f\) have the same structure.

2.1.4.2 Limit \(N\to\infty\) and Propagation of Chaos

Even though we added significant sources of coupling between the agents experiences, exchangeability among the agents persists, and in the limit \(N\to\infty\), a conditional form of the propagation of chaos should lead to conditional independence of individual dynamics. This asymptotic conditional independence is quite intuitive if one is familiar with the standard theory of the propagation of chaos for Itô stochastic differential equations. However, stating the discrete time conditional form of this result precisely with full mathematical rigor is rather difficult, and it requires new notation.

Special Notation: For any sequence \((\zeta_n)_{n\ge 0}\) we shall use the notation \(\underline{\zeta}_n=(\zeta_0,\zeta_1,\cdots,\zeta_n)\).

Even if we do not call them by their names, in what follows, we introduce the elements of the definition of open-loop policies which is given later on in Definition 1.13.

We fix a sequence \(\bpi=(\pi_n)_{n\ge 0}\) of deterministic measurable functions \(\pi_n:\Xi_n\to \cP(A)\) where \(\Xi_n = \cY\times (\Theta\times\Theta^0\times E\times E^0)^n\times \Theta^0\) where \(\cY\), \(\Theta\) and \(\Theta^0\) are Polish spaces from which samples are drawn to randomize the initial condition, the individual action choices and the common mixed action. Next, for each integer \(i\ge 1\), which we may want to think of as the index of an agent in an infinite population, we consider a random variable \(Y^i\) in \(\cY\) to randomize the initial condition, and i.i.d. sequences \((\theta^i_n)_{n\ge 0}\) to randomize the action choices of the individual agents, \((\theta^0_n)_{n\ge 0}\) to randomize the action choices of the central unit, \((\varepsilon^i_n)_{n\ge 0}\) as realizations of the idiosyncratic random shocks, and \((\varepsilon^0_n)_{n\ge 0}\) as realizations of the common random shocks. The families indexed by \(i\) are independent across agents and independent of the shared common sequences.

Then, again for each \(i\ge 1\) fixed, we consider the stochastic process \((X^{\bpi,i}_n)_{n\ge 0}\) describing the time evolution of the \(i\)-th state given by:

\[\label{fo:MFC95individual95dynamics} X^{\bpi,i}_{n+1}=F(X^{\bpi,i}_n,\alpha^{\bpi,i}_n,\PP^0_{(X^{\bpi,i}_n,\alpha^{\bpi,i}_n)},\varepsilon^i_{n+1},\varepsilon^0_{n+1}),\quad n\ge 0,\tag{6}\] where the random action \(\alpha^{\bpi,i}_n\) is chosen as \[\label{fo:alpha95n} \alpha^{\bpi,i}_n=\pi_n\bigl(Y^i,\underline\theta^i_n,\underline\theta^0_n,\underline\varepsilon^i_n,\underline\varepsilon^0_{n+1}\bigr),\qquad n\ge 0.\tag{7}\] Note that all the processes \((X^{\bpi,i}_n)_{n\ge 0}\) are conditionally independent given the shared common sequences and they follow the same dynamical equation \[\label{fo:MFC95dynamics} X_{n+1}=F(X_n,\alpha_n,\PP^0_{(X_n,\alpha_n)},\varepsilon_{n+1},\varepsilon^0_{n+1}),\quad n\ge 0,\tag{8}\] where the notation \(\PP^0_{(X_n,\alpha_n)}\) stands for the conditional law of the state-action couple \((X_n,\alpha_n)\) given the common noise \((\varepsilon^0_n)_{n\ge 1}\), and where the action process \((\alpha_n)_{n\ge 0}\) is obtained from \(\bpi=(\pi_n)_{n\ge 0}\) and a sequence of randomization random variables and idiosyncratic and common noise sequences as in 7 .

Accordingly, each agent \(i\) of this infinite population of conditionally independent agents, should solve the following optimization problem: \[\label{fo:MFC95optimization} \inf_{\balpha=(\alpha_n)_{n\ge 0}} \EE\Bigl[\sum_{n=0}^\infty \gamma^n f\bigl(X_n,\alpha_n,\PP^0_{(X_n,\alpha_n)}\bigr)\Bigr].\tag{9}\] The intuition that such set-up is the proper asymptotic regime for a symmetric population of \(N\) agents when \(N\to\infty\) is substantiated by the following mathematical result. By constructing the \(N\)-agent model on the same probability structure (i.e. with the same sequences of randomization random variables and idiosyncratic and common random shocks), it is possible to prove that this convergence is actually in the almost sure sense.

Recall that the distance \(d_w\) is the metric of the weak convergence of probability measures introduced in Subsection 2.1.1.2.

More Measurability Considerations

The form of Mean Field discrete time Control problem with common noise given by 8 and 9 is now very close to the model we investigate in the next chapter. The last question we would like to address before turning to the theoretical analysis of the model is:

Could the conditional distributions \(\PP^0_{(X_n,\alpha_n)}\) depend upon some extra sources of randomness beyond the randomness of the common noise?

Indeed, if individual agents and the central unit performing the social cost optimization are allowed to use mixed strategies, they need independent sources of randomness to randomize their actions and decisions at each time. So since the cost optimization and the choice of action strategy are performed by such a central unit,

• the central unit needs a source of randomness to randomize the choice of a mixed policy to be dispatched to the individual agents, and

• the individual agents need to sample their actions from the policy sampled for them by the central unit.

In the limit \(N\to\infty\), we expect that the idiosyncratic randomizations of the individual agents average out. On the other hand, we do not see why the central unit randomization should also average out, so we should expect that it remains. The motivation for the common randomization is two-fold. First, in full analogy with single-agent RL, we expect the common randomization to help exploring the state space and the action space during learning. This point will be of crucial importance when we show that the mean field Q-learning algorithm converges under the assumption that the state-action sequence covers repeatedly the whole state-action space. Moreover, this form of randomization will also prove useful for technical reasons.

As a result, the conditioning in \(\PP^0_{(X_n,\alpha_n)}\) should be with respect to the past values of the common noise \(\varepsilon^0_1,\cdots,\varepsilon^0_n\) as well as all the sources of randomness used by the central unit to randomize the policy handed out to the individual agents for implementation. We shall provide rigorous mathematical definitions for this conditioning in the next chapter.

2.1.5 A General Abstract MFC Model↩︎

The next two chapters are devoted to the complete mathematical analyses of two fundamental models whose definitions we make explicit in this section for the sake of future reference. Both models will be set in infinite horizon for the sake of definiteness, and we shall use the notation \(\gamma\in(0,1)\) for the discount factor in both cases. We start with the model that we study in Chapter 2.2.

The model aims at describing the time evolution of a homogeneous population of agents, the dynamics of the state of a generic agent being given by: \[X_{n+1} = F(X_n, a_n, \mu_n, \varepsilon_{n+1}, \varepsilon^0_{n+1}),\] where at each time \(n\):

• \(X_n \in S\) is the state of the generic agent,

• \(a_n \in A\) is the action taken by the agent,

• \(\mu_n \in \cP(S \times A)\) is the joint distribution of states and actions across the population,

• \((\epsilon_n)_{n\ge 1}\) is an independent identically distributed (i.i.d. for short) sequence of random elements in \(E\) whose common distribution we denote by \(\nu\), \(\varepsilon_{n+1} \in E\) representing the next value of the idiosyncratic random shock,

• \((\epsilon^0_n)_{n\ge 1}\) is an i.i.d. sequence of random elements in \(E^0\) whose common distribution we denote by \(\nu^0\), \(\varepsilon^0_{n+1} \in E^0\) is the next value of the common noise.

The one-stage cost incurred at time \(n\) by the generic agent is \(f(X_n, a_n, \mu_n)\). Notice that it may depend on the population distribution \(\mu_n\), capturing mean field interactions.

The goal of the analysis is the minimization of the expected discounted social cost \[\label{fo:J} J = \EE\left[\sum_{n=0}^\infty \gamma^n f(X_n, a_n, \mu_n)\right],\tag{11}\] where the expectation is over all the sources of randomness: initial conditions, noise sequences, and sources of randomness used for policy randomization. Notice that we did not specify yet the space of actions - controls - policies over which the optimization is performed. We shall discuss the admissible information structures later on.

Assumptions for MFC with Common Noise

The results stated and proven in Chapter 2.2 hold under various sets of assumptions. Still, for future reference, we list here the assumptions that will be in force throughout the analysis of the model. In some instances, more restrictive assumptions will be added in order to derive specific results.

Since \(S\) and \(A\) are assumed to be Polish spaces, \(S\times A\) and \(\cP(S\times A)\) are also Polish spaces, and the notion of continuity alluded to above is understood with respect to the corresponding weak topology on probability measures.

2.1.5.1 Policies and Policy Spaces

Policies are functions with values in the action space \(A\), or possibly in the space \(\cP(A)\) of probability distributions on \(A\) when using randomized actions i.e. mixed strategies. Their measurability properties offer a convenient way to specify the information available to the agent making a decision. Indeed, this information is coded by the arguments of the policy function. Because we shall use several classes of policies, we state their definitions for future reference. Recall the special notation introduced earlier: for any given sequence \((\zeta_n)_{n\ge 0}\) we use \(\underline{\zeta}_n=(\zeta_0,\zeta_1,\cdots,\zeta_n)\). We want to think of the minimization of the objective 11 as the job of a central unit. As suggested by economics applications, we shall sometimes use the terminology central planner. Still, the decisions made by the central unit will have to be implemented by the individual agents, typically the generic agent in the mean field limit. Both the central unit and the generic agent should be able to randomize their choices, so for this purpose, we introduce spaces \(\Theta^0\) and \(\Theta\) and i.i.d. sequences \((\theta^0_n)_{n\ge 0}\) and \((\theta_n)_{n\ge 0}\) of random variables in \(\Theta^0\) and \(\Theta\) respectively which they will use to randomize their action choices.

We denote by \(\bPi^{OL}\) the set of open loop policies. Typically, an action taken from an open loop Markov policy is a sample from a distribution on the action space \(A\), which depends on the the initial condition, the past sources of randomization of the central unit and the generic agent, the past values of the idiosyncratic noise and the common noise observed up to time \(n\) and the current common randomization.

We denote by \(\bPi^{CL}\) the set of closed loop Markov policies. In other words, an action taken from a closed-loop Markov policy is a sample from a distribution on the action space \(A\), which depends on the current state \(X_n\), the current population distribution \(\mu_n\), and the current common randomization \(\theta^0_n\).

As expected, we shall say that a policy \(\pi = (\pi_n)_{n\geq 0}\) is stationary if \(\pi_n = \hat{\pi}\) for all \(n\), for some fixed function \(\hat{\pi}\) independent of \(n\).

2.1.5.2 The Lifted Mean Field MDP

The analysis of McKean-Vlasov dynamics, is typically based on stochastic differential equations whose coefficients depend upon the distribution of the solution, sometimes called nonlinear stochastic differential equations. These equations are also called mean field equations as they appear in the limit of the description of the time evolution of particle systems interacting in a mean field manner. Most theoretical analyzes of these equations rely on a lifting procedure replacing the original problem by a (deterministic) dynamical system that describes the time evolution of the probability distribution of the original state of the system. As can be seen in the monograph [95], this idea is already present in early works on the control of Markov Decision Processes (MDPs). However, due to the presence of two independent sources of random shocks, a simple lifting of our model will not lead to the optimization of a deterministic dynamical system. It will lead to a simpler MDP, which is still stochastic, though driven only by the common noise. So in order to end up with a deterministic dynamical system, we will have to go through two successive lifting operations.

The lifted MFMDP is a standard MDP in the space \(\hat{S}\) of probability measure on the agent state space \(S\), the new state of which evolves according to the recursion \[\mu_{n+1} \sim \hat{P}(\cdot|\mu_n, \bar{a}_n),\quad\text{or equivalently}\quad \mu_{n+1} = \hat{F}(\mu_n, \bar{a}_n,\varepsilon^0),\quad \bar{a}_n \in \hat{U}(\mu_n),\] with one-stage cost \(\hat{f}(\mu_n, \bar{a}_n)\). We shall use tools from the standard theory of MDPs, relying heavily on [95], to analyze this lifted MFMDP in the next chapter. Under suitable assumptions, the value function of this lifted MDP should satisfy the Bellman equation \[J(\mu) = \inf_{\bar{a} \in \hat{U}(\mu)} \left\{\hat{f}(\mu, \bar{a}) + \gamma \, \EE[J(\hat{F}(\mu, \bar{a}, \varepsilon^0))]\right\},\] where the expectation is over the common noise \(\varepsilon^0 \sim \nu^0\). Recall that \(\hat{F}\) is the system function of the lifted MDP, namely \(\hat{F}(\mu, \bar{a}, \varepsilon^0)\) is the distribution of \(F(X, a, \bar{a}, \varepsilon, \varepsilon^0)\) with \((X, a) \sim \bar{a}\), \(\varepsilon \sim \nu\). In order to guarantee solvability of the Bellman equation, we shall prove:

Consequently, by the Banach fixed point theorem (Theorem 1.3), \(T\) has a unique fixed point \(J^* = TJ^*\), which is the optimal value function. Moreover, \(J^*\) can be obtained as the limit of value iterations: \[J^* = \lim_{n \to \infty} T^n J_0,\] for any initial bounded function \(J_0\).

State-Action Value Functions and Existence of Optimal Policies

Because of our interest in RL, we shall use the state-action value function also known as the Q-function defined by \[Q^*(\mu, \bar{a}) = \hat{f}(\mu, \bar{a}) + \gamma \, \EE[J^*(\hat{F}(\mu, \bar{a}, \varepsilon^0))],\] With its help, one will characterize optimal states for the lifted MDP:

Finally, using the fact that \(Q^*\) is the unique solution of the Bellman equation for the Q-functions, i.e. \[Q^*(\mu, \bar{a}) = \hat{f}(\mu, \bar{a}) + \gamma \, \EE_{\mu' \sim \hat{P}(\cdot|\mu,\bar{a})} \left[\min_{\bar{a}' \in \hat{U}(\mu')} Q^*(\mu', \bar{a}')\right].\] one will derive existence of pure stationary optimal policies for the lifted MDP model:

2.1.6 The Linear-Quadratic Mean Field Control Model↩︎

We now switch to the second MFC model studied in this monograph. In some sense, the LQ MFC model which we analyze in detail in Chapter 2.3 can be viewed as a particular case of the model of Chapter 2.2 for which the system function is linear and the one-stage cost function is quadratic. To be more specific, we assume that the state and action spaces are Euclidean spaces, say \(S= \RR^d\) and \(A= \RR^\ell\) \[\label{eq:lq95sys95function95F} F(x,a,\mu,\epsilon,\epsilon^0)=Ax +Ba+\bar A\bar\mu^x+\bar B \bar\mu^a+\epsilon+\epsilon^0,\tag{12}\] where the matrices \(A \in \RR^{d \times d}\) and \(B \in \RR^{d \times \ell}\) drive the dynamics and the matrices \(\bar{A} \in \RR^{d \times d}\) and \(\bar{B} \in \RR^{d \times \ell}\) capture the mean field coupling as we denote by \(\bar\mu^x\) the mean of the first marginal of \(\mu\), namely \(\bar\mu^x=\int_{S\times A}x\mu(dx,da)\) and by \(\bar\mu^a\) the mean of the second marginal of \(\mu\), namely \(\bar\mu^a=\int_{S\times A}a\mu(dx,da)\). Obviously, the marginal distributions \(\nu\) and \(\nu^0\) of the idiosyncratic and common noises are now probability measures on the Euclidean space \(\RR^d\). They are most commonly assumed to be mean-zero Gaussian distributions, but we shall show that many results still apply to more general distributions.

Similarly, we assume that the one-stage cost function is of the form \[\label{fo:LQ95cost} f(x,a,\mu) = (x-\bar{\mu}^x)^T Q (x-\bar{\mu}^x) + (a-\bar{\mu}^a)^T R (a-\bar{\mu}^a) + (\bar{\mu}^x)^T (Q+\bar{Q}) \bar{\mu}^x + (\bar{\mu}^a)^T(R+\bar{R}) (\bar{\mu}^a),\tag{13}\] where we use the exponent \({}^T\) to denote the transpose of a matrix or a vector and where:

• \(Q, \bar{Q} \in \text{SM}_d\) are state cost matrices with \(Q \geq 0\), \(Q + \bar{Q} > 0\),

• \(R, \bar{R} \in \text{SM}_\ell\) are control cost matrices with \(R > 0\), \(R + \bar{R} > 0\).

While it would be tempting to apply the results of Chapter 2.2 to the present model, most of the proofs rely on the compactness of some spaces or the boundedness of the system and cost functions, hypotheses which are obviously not satisfied here. Still the linear-quadratic structure of the model makes the problems amenable to explicit or semi-explicit solutions and we shall take full advantage of this possibility.

2.1.6.1 Finite-Agent LQ Model

In the \(N\) agent version of the model, each agent state \(X^i_n \in \RR^d\) at time \(n\) evolves according to the recursive equation \[\label{fo:N95LQ95agent95dynamics} X^i_{n+1} = A X^i_n + B a^i_n + \bar{A} \bar{X}_n + \bar{B} \bar{a}_n + \varepsilon^i_{n+1} + \varepsilon^0_{n+1}, \quad i = 1, \ldots, N,\tag{14}\] where:

• \(\bar{X}_n = \frac{1}{N} \sum_{i=1}^N X^i_n\) and \(\bar{a}_n = \frac{1}{N} \sum_{i=1}^N a^i_n\) are the empirical mean state and control,

• \(\varepsilon^i_{n+1} \sim \nu\) are i.i.d.idiosyncratic noise (mean zero, covariance \(\Sigma_\varepsilon\)),

• \(\varepsilon^0_{n+1} \sim \nu^0\) is common noise (mean zero, covariance \(\Sigma_{\varepsilon^0}\)), independent of all \((\varepsilon^i_{n+1})\).

At time \(n\), the one-stage cost for agent \(i\) is \[\label{eq:finite-agent-cost} f^i(X_n,a_n)= (X^i_n-\bar X_n)^T Q (X^i_n-\bar X_n) + (a^i_n-\bar a_n)^T R (a^i_n -\bar a_n)+ \bar{X}_n^T (Q+\bar{Q}) \bar{X}_n + \bar{a}_n^T (R+\bar{R}) \bar{a}_n,\tag{15}\] where the cost matrices \(Q, \bar{Q}, R, \bar{R}\) satisfy the same conditions as in 13 .

Since we are restricting ourselves to cooperative agents, the objective is to minimize the expected discounted social cost \[J^N = \EE\left[\sum_{n=0}^\infty \gamma^n \frac{1}{N} \sum_i f^i(X_n,a_n)\right].\]

2.1.6.2 LQ MFC Model (Mean Field Limit)

Given what we have seen so far, it is not surprising that in the asymptotic regime of large homogeneous populations of agents, i.e. when \(N \to \infty\), the finite-agent model converges to the LQ MFC model, where the state of a generic agent evolves as \[\label{eq:mfc-dynamics} X_{n+1} = A X_n + B \alpha_n + \bar{A} \bar{X}_n + \bar{B} \bar{\alpha}_n + \varepsilon_{n+1} + \varepsilon^0_{n+1},\tag{16}\] with:

• \(\bar{X}_n = \EE[X_n\mid\sigma\{\underline\varepsilon^0_n\}]\) the mean field state (population average),

• \(\bar{\alpha}_n = \EE[\alpha_n\mid\sigma\{\underline\varepsilon^0_n\}]\) the mean field control.

The one-stage cost is \[\label{eq:mfc-cost} f(X_n,\alpha_n,\bar X_n,\bar\alpha_n)= (X_n-\bar X_n)^T Q (X_n-\bar X_n) + (\alpha_n-\bar\alpha_n)^T R (\alpha_n-\bar\alpha_n)+ \bar{X}_n^T (Q+\bar Q) \bar{X}_n + \bar{\alpha}_n^T (R+\bar R) \bar{\alpha}_n,\tag{17}\] and the objective is to minimize \[J = \EE\left[\sum_{n=0}^\infty \gamma^n f(X_n, \alpha_n, \bar{X}_n, \bar{\alpha}_n)\right].\]

2.1.6.3 Assumptions for LQ MFC

We state the main assumptions we shall use in the proofs in Chapter 2.3 for future reference.

We shall also need an extra assumption, used in a technical proof, ensuring that a specific set of coupled matrix Riccati equations admit unique positive definite solutions.

2.1.6.4 Admissible Control Parameters and Optimality Characterization

The mainstay of our analysis of the optimization of the LQ MFC problem is to show that one can reduce the search for optimal control to the special class of controls in linear closed-loop Markov (feedback) form: \[\label{eq:linear-feedback} \alpha_n = -K (X_n-\bar X_n) - L \bar{X}_n,\tag{18}\] for a couple \((K,L)\) of (deterministic) matrices:

• \(K \in \RR^{\ell \times d}\) representing the individual feedback gain,

• \(L \in \RR^{\ell \times d}\) representing the mean field feedback gain.

The rationale for these admissibility conditions is that they ensure that both the individual dynamics (\(X_{n+1} \approx (A - BK)X_n + \text{noise}\)) and the mean field dynamics (\(\bar{X}_{n+1} \approx (A + \bar{A} - (B + \bar{B})L)\bar{X}_n + \text{noise}\)) are contractive under discounting.

Subsequently, the main result allowing us to characterize optimality in this class of admissible linear controls will be stated at this stage in an informal form as

Note that in the case of the LQMFC model, the Bellman equation reads:

2.1.7 Reinforcement Learning Framework↩︎

2.1.7.1 Model-Free Learning

In many applications, the system function \(F\) and/or the cost function \(f\) are not explicitly known. Reinforcement learning (RL) provides a framework for learning approximate optimal policies through interaction with the environment, treating it as a black box.

2.1.7.2 Q-Learning

In the particular case of the lifted MDP, Q-learning proposes a model-free algorithm that estimates \(Q^*\) by iteratively updating estimates based on observed transitions: \[\label{eq:q-learning} Q_{k+1}(\mu, \bar{a}) = (1 - \alpha_k) Q_k(\mu, \bar{a}) + \alpha_k \left[r + \gamma \min_{\bar{a}'} Q_k(\mu', \bar{a}')\right],\tag{19}\] where:

• \(\alpha_k \in (0, 1)\) is the learning rate,

• \(r = \hat{f}(\mu, \bar{a})\) is the observed cost,

• \(\mu' \sim \hat{P}(\cdot|\mu, \bar{a})\) is the next state.

The popularity of Q-learning is due to the fact that under appropriate conditions (diminishing step sizes, sufficient exploration), \(Q_k \to Q^*\) almost surely. Such an algorithmic search for an approximation of an optimal solution works surprisingly well when the state space and the action space are finite with a very small number of elements. As argued in Chapter 1.1, unfortunately, Q-learning performance deteriorates very quickly when the cardinality of any of these spaces increases. So when the state or action spaces are large or continuous, \(Q\) is approximated by a parametric function \(Q_\theta\) (e.g., a neural network), and parameters \(\theta\) are learned via updates provided by stochastic gradient descent.

2.1.7.3 Mean Field Reinforcement Learning (MFRL)

Common approaches include:

• Finite-dimensional approximations: we discretize \(\cP(S)\), or represent distributions via finite-dimensional sufficient statistics (e.g., moments, histograms).

• Particle-based methods: we approximate \(\mu\) by empirical distributions of finite particle systems.

• Neural network function approximation: we use deep networks to represent value functions or policies over distribution spaces.

We shall provide examples of applications of all these methods.

Simplex Discretization

For state or action spaces that are probability simplexes (e.g., discrete distributions), we approximate continuous distributions by finite-dimensional simplexes.

Histogram Approximations

For continuous state spaces, histogram approximations partition the state space into bins and represent distributions by histograms (piecewise constant densities). This reduces the infinite-dimensional problem to a finite-dimensional one, enabling standard RL algorithms.

Monte Carlo Simulation and Particle Systems

To implement environments for RL algorithms, we simulate the mean field dynamics using particle systems: a large but finite collection of \(N\) agents whose empirical distribution approximates the population distribution \(\mu\).

1. Initialize \(N\) agents with states \(\{x^1_0, \ldots, x^N_0\}\) sampled from the initial distribution.

2. At each time step, compute the empirical distribution \(\mu_N = \frac{1}{N} \sum_i \delta_{x^i}\).

3. Each agent takes action \(a^i\) according to the policy, and transitions to \(x^i_+\) via the dynamics \(F(x^i, a^i, \mu_N, \varepsilon^i, \varepsilon^0)\).

4. Repeat for \(T\) time steps, collecting cost observations.

Neural Network Function Approximation

For complex, high-dimensional problems, deep neural networks are used to approximate:

• the Value functions \(J_\theta(\mu)\) or \(Q_\theta(\mu, \bar{a})\),

• the Policies: \(\pi_\theta(a \mid x, \mu)\).

Common architectures include:

• Feedforward fully connected networks: Standard multilayer perceptrons.

• Convolutional networks: For spatial or image-based state representations.

• Recurrent networks: For temporal dependencies (though less common in our discrete-time Markov settings).

2.1.7.4 Policy Gradient Methods for LQ MFC

In the particular case of the LQ framework considered here, the control parameters \(\theta = (K, L)\) can be optimized via policy gradient descent: \[\label{eq:policy-gradient} \theta_{k+1} = \theta_k - \beta_k \nabla_\theta J(\theta_k),\tag{20}\] where \(\beta_k > 0\) is the step size and \(\nabla_\theta J(\theta)\) is the gradient of the expected social cost with respect to \(\theta\).

In our subsequent analysis in Chapter 2.3, we derive gradient estimates from finite differences by perturbing \(\theta\) and estimating gradients from cost differences. Our control of the convergence depends upon conditions such as \(\gamma \left\|A - BK\right\|^2 < 1\) which ensure that the closed-loop system is exponentially stable, so that perturbations decay over time and value functions remain bounded.

Sub-Gaussian Concentration Inequalities

LQ models are tailor-made for the theory of Gaussian processes. However, while convenient mathematically, assuming that random shocks are Gaussian is often unrealistic. So in order to improve the realm of application of theoretical results, a valuable endeavor is to extend their reach to sub-Gaussian and sub-exponential random variables. Accordingly, we employ sub-Gaussian and sub-exponential norms to control tail behavior.

We also use the standard definition and notation of \(L^p\)-norms for a random vector \(\xi \in \RR^d\): \[\left\|\xi\right\|_{L^p(\RR^d)} = \left(\EE[\left\|\xi\right\|^p]\right)^{1/p},\] where \(\left\|\cdot\right\|\) denotes the Euclidean norm on \(\RR^d\). For a random vector \(\xi \in \RR^d\), we use the vector sub-Gaussian norm \[\left\|\xi\right\|_{\psi_2} = \sup_{\substack{u\in\RR^d\\ \left\|u\right\|=1}} \left\|\langle u,\xi\rangle\right\|_{\psi_2},\] where the norm on the right-hand side is the scalar norm of Definition 1.20.

In the LQ framework, sub-Gaussian concentration inequalities turn out to be a powerful tool providing tail bounds for sums of independent random vectors. A standard form, see e.g. [96], is the following.

This inequality controls deviations of sample averages from their expectations, with the expected dependence on the dimension.

2.1.8 Summary and Roadmap↩︎

This chapter has laid the mathematical groundwork for the MFC problems and learning methods studied in the subsequent chapters:

1. Chapter 2.2: Mean Field Control with Common Noise

   • Focuses on the abstract MFC model with both idiosyncratic and common noise.

   • Analyzes the lifted MFMDP and establishes equivalence between open-loop and closed-loop formulations.

2. Chapter 2.3: Linear-Quadratic Mean Field Control

   • Specializes to the LQ setting with quadratic costs and linear dynamics.

   • Derives explicit solutions and studies finite-agent approximations.

3. Chapter 3.1: Numerical Methods for Mean Field Reinforcement Learning

   • Develops Q-learning and deep RL algorithms for model-free learning.

   • Presents approximation and discretization methods for computational implementation.

4. Chapter 3.2: Policy Gradient Methods for LQ Mean Field Control

   • Develops policy-gradient methods for learning feedback gains.

   • Analyzes sample complexity, approximation errors, and numerical experiments.

2.1.9 Bibliographical Notes and Complements↩︎

The notations and the definitions from the measure theoretical part at the beginning of the chapter are generic in stochastic analysis. For consistency, we tried to follow [97] and [98]. In particular, a proof of the disintegration result stated in Proposition [pr:disintegration] can be found in [97], and the universal disintegration result stated in Theorem 1.1 is from [98].

A proof of the original Blackwell-Dubins lemma 1.1 can be found in the original paper [7] by the authors. See also [22] for a self-contained explicit construction. A detailed discussion of the randomization practice in reinforcement learning can be found in [99]. There, the similarities and the differences with the use of relaxed controls in continuous time control problems are identified and put in perspective.

The proof of propagation of chaos for continuous-time dynamical systems driven by stochastic differential equations can be found in many textbooks on stochastic analysis. For convenience, we refer the interested reader to [22] and the references given therein, including classical propagation of chaos results in statistical physics, see [100]. Discrete-time versions are few and far between. Even though it is only proven for compact state and action spaces, the result of Motte and Pham [4], quoted in Proposition [pr:chaos], is particularly well suited to our setting because it treats mean field Markov decision processes with common noise and open-loop controls. For a recent discrete-time propagation-of-chaos result for nonlinear Markov chains, see also [101].

Randomization via the use of mixed controls is ubiquitous in the practice of reinforcement learning. A recent wave of interest in the connection between the theory of continuous-time stochastic control and the applications of discrete-time reinforcement learning has been the source of numerous publications and some confusion. While speculation on these connections can be enlightening in some cases, it has some limitations, especially when it comes to handling the similarities and the differences between the randomization procedures used in both fields as the use of mixed controls in discrete time reinforcement learning, and relaxed controls in continuous time stochastic controls have very different impacts on the system being controlled. The reader interested in this issue may want to consult the review article [99] for a discussion of the parallels and the differences between the two approaches in the algorithmic implementations of these two classes of models.

2.2 A General Mathematical MFRL Model with Common Noise↩︎

2.2.1 Introduction↩︎

This chapter provides a rigorous mathematical analysis of the first model introduced in Chapter 2.1. After a quick review of the model, we study the lifted mean-field MDP (MFMDP) and its connections with the original mean field control (MFC) problem (Section 2.2.4), and we introduce the state-action value function (Q-function) of the MFMDP (Section 2.2.5).

The thrust of this chapter is twofold. First, we study a setting which allows for common randomization, i.e., mixed policies at the mean-field level. To the best of our knowledge, this is a departure from the standard literature on mean field games (MFGs) and MFC problems. We believe that common randomization is of fundamental importance in the RL setup to ensure exploration in the space of the central planner’s actions. Second, we prove equality of the closed-loop and open-loop value functions for any given policy, and not just for the optimal ones: in Section 2.2.4.2, we prove that for any open-loop policy, there is a closed-loop policy achieving the same value, and vice versa.

In the second part of the monograph, we describe the implementation of several RL methods based on the theoretical results of this chapter, and we discuss the results: a tabular Q-learning algorithm relying on a discretization of the mean-field state simplex for which we prove convergence (Theorem 1.1), and a deep RL method to deal with continuous state and action spaces, which allows us to avoid simplex discretization and to deal with randomized actions.

2.2.2 Model Description and Notations↩︎

We consider the model introduced in Section 2.1.5 of Chapter 2.1. We shall repeat the notations and the definitions already stated there only when necessary for the clarity of our statements.

Notice that because we work with a homogeneous population of agents, we use a single action space \(A\) common to all the agents, and a common cost function \(f\) for all the agents. This is in contrast with our original discussion of MARL in Chapter 1.1 where the action spaces and the cost functions were specific to the agent.

As we argued in Chapter 1.1, the dynamics of a MARL or MFRL model can be given equivalently by a set of transition probabilities or a system function. Since we plan to devote a significant part of this chapter to dependencies with respect to the whole past via open loop policies, we find it more convenient to work with the stochastic dynamical equations of the systems as defined by a system function and innovation processes, both idiosyncratic and common. Our system function \(F\) is used to describe the evolution of the state process based on a state, an action, a mean-field interaction term⁴ and two noise terms.

Even if we are willing to postpone the regularity conditions on \(F\) and \(f\), the definition of the model cannot be complete without explaining how the actions are taken and how these decisions affect the evolution of the system over time. The following features of our model are important:

1. The mean-field interactions are conditioned on all shared information.

2. The actions are randomized by additional sources of randomness.

were clearly motivated in Chapter 2.1 by our discussion of the propagation of chaos and how taking the limit when the number of agents goes to infinity could help identify the right mathematical features of the mean field limit.

As explained earlier in Chapter 2.1, the formulation of the model is motivated by the limit of an \(N\)-agent model when the size of the population \(N\) tends to infinity. Still, because we plan to use randomized controls and policies, we need to introduce the various sources of randomization at play. Beyond the i.i.d. sequences \((\varepsilon_{n+1})_{n\ge 0}\) with distribution \(\nu \in \cP(E)\) modeling the idiosyncratic random shocks and \((\varepsilon^0_{n+1})_{n\ge 0}\) with distribution \(\nu^0 \in \cP(E^0)\) modeling the common noise, we shall assume the existence of the following sources of randomness on the basic probability space \((\Omega,\cF,\PP)\):

1. a random variable \(\mathscr{U}\) with distribution \(\PP_{\mathscr{U}}\) in a Polish space \((\Upsilon, \cB_{\Upsilon})\) providing the randomization for the initial state;

2. an i.i.d. sequence \((\vartheta_n)_{n\ge 0}\) of random variables in a Polish space \((\Theta, \cB_{\Theta})\) with distribution \(\PP_{\vartheta}\) providing the generic agent with a source of randomization for their action choices;

3. an i.i.d. sequence \((\vartheta^0_n)_{n\ge 0}\) of random variables in a Polish space \((\Theta^0, \cB_{\Theta^0})\) with distribution \(\PP_{\vartheta^0}\) providing the central unit with a source of randomization for the choices of policies.

We assume that all these random sequences are independent of each other. We also assume that \(\PP_{\vartheta}\) and \(\PP_{\vartheta^0}\) are both atomless. This guarantees the existence of Borel measurable functions \(h:\Theta\to[0,1]\) and \(h^0:\Theta^0\to[0,1]\) which are uniformly distributed when viewed as random variables on the probability spaces \((\Theta,\cB_\Theta,\PP_\vartheta)\) and \((\Theta^0,\cB_{\Theta^0},\PP_{\vartheta^0})\) respectively.

The uniform random variables constructed with the functions \(h\) and \(h^0\) will be used repeatedly with the Blackwell-Dubins Lemma.

2.2.2.1 Probabilistic Framework and Classes of Policies

We now introduce the framework that will be used to rigorously study the MFC problem before we implement the theory so developed through RL algorithms. As explained intuitively in the introduction, we will distinguish several types of randomness. We will also distinguish between actions (elements of \(A\)), and controls (probability measures on \(A\)). They are the building blocks needed to define later the notion of policy (see § 2.2.2.1.2).

We recall the following useful notation which was already introduced in Chapter 2.1. For a given sequence \((\zeta_n)_{n\ge 0}\) we use the notation \(\underline{\zeta}_n=(\zeta_0,\zeta_1,\cdots,\zeta_n)\).

Next, we introduce the terminology which is going to help us characterize the information available to a generic agent and the central planner to make their choices of actions and policies. Incidentally, we start referring to everything related to agents (e.g. actions, controls, policies, costs, etc) as level-0. This is in contrast with the lifted stochastic optimization model which we introduce in Section 2.2.3 whose elements will be referred to as level-1. For reasons which will become clear later, the central unit controlling the agents will be called the level-1 controller. This new terminology frees our presentation of the theoretical results from the details of the application we used to motivate the set-up.

For later use, set \[\cF_n^0=\sigma\{\underline{\varepsilon}_n^0,\underline{\vartheta}_{n-1}^0\}, \qquad \cG_n^0=\sigma\{\underline{\varepsilon}_n^0,\underline{\vartheta}_n^0\},\] with the convention that \(\underline{\vartheta}_{-1}^0\) is empty. The sigma-field \(\cF_n^0\) represents the common information before sampling the current common randomization \(\vartheta_n^0\), while \(\cG_n^0\) includes it.

Intuitively, an action process is the realization of a control process, where the sampling is done using \((\vartheta_n)_{n \ge 0}\). Generally speaking, we use the term mixed for probability measures and randomized for the corresponding random variable. Typically, a level-0 action is denoted by \(a\) and a level-0 control is denoted by \(\fa\). Unless specified otherwise, all the controls we consider at the level-0 are implicitly assumed to be mixed.

For notational convenience, we use both notations \(\PP_\xi\) and \(\cL(\xi)\) interchangeably for the distribution of a random element \(\xi\). We also use natural extensions to denote conditional distributions. The following simple measure theory result will prove to be useful in the sequel.

2.2.2.1.1 Conditional distribution and state process

We are now in a position to describe precisely the mean-field interactions in the system function, and provide a clear definition of the state process driven by a mixed control process in the mean-field model with common noise.

Notice that for \(n\ge 1\), the state \(X^{\balpha, \mu_0}_n\) is measurable with respect to \(\sigma\{ \sU, (\vartheta_k, \vartheta_k^0, \varepsilon_{k+1}, \varepsilon_{k+1}^0 )_{k = 0, \ldots, n-1} \}\) which does not include \(\vartheta_n\) and \(\vartheta_n^0\) as \(X_n^{\balpha, \mu_0}\) does not depend on \(\vartheta_n\) and \(\vartheta_n^0\). These random variables are used to define the control at time \(n\), which will in turn influence the state at the next time step, namely, \(X_{n+1}^{\balpha, \mu_0}\).

For each level-\(0\) action process \(\balpha\) and each \(n \geq 0\), we denote by \(\PP^0_{X_n^{\balpha, \mu_0}}\) a regular version of the conditional distribution \(\cL(X_n^{\balpha, \mu_0}\,|\,\cF_n^0)\). It holds: \[\label{fo:P0X95n} \PP^0_{X_n^{\balpha, \mu_0}} =\cL(X_n^{\balpha, \mu_0}\,|\,\cG_n^0), \qquad \mathbb{P}-a.s..\tag{21}\] Indeed, conditioned on \(\cF_n^0\), \(X_n^{\balpha,\mu_0}\) is independent of \(\vartheta_n^0\).

2.2.2.1.2 Open-loop and closed-loop policies

The definitions of this subsection were already introduced in Subsection 2.1.5.1 of Chapter 2.1. We repeat them here because of their technicalities and their importance for the subsequent developments. We first consider open-loop policies. For each \(n\ge 0\), let \(\Xi_n = \Upsilon \times (\Theta \times \Theta^0 \times E \times E^0)^n \times \Theta^0\), with the convention that \(\Xi_0 = \Upsilon \times \Theta^0\), and let us define \(\bxi=(\xi_n)_{n\ge 0}\) by: \[\xi_0 = (\sU, \vartheta_0^0), \qquad\text{and}\qquad \xi_n = (\sU, (\vartheta_k, \vartheta_k^0, \varepsilon_{k+1}, \varepsilon_{k+1}^0)_{k = 0, \ldots, n-1}, \vartheta_n^0), \qquad n \ge 1.\]

Because of measurability restrictions, the correspondence between level-0 control processes \(\bfa\) and open-loop policies \(\bpi\) is understood up to \(\PP\)-a.s. equality of the induced control processes, after choosing measurable versions in ?? . We shall use one object or the other depending on whether we want to emphasize the stochastic process or the sequence of deterministic functions defining the process. We shall often short-circuit the control process \(\bfa\) and say that an action process \(\balpha = (\alpha_n)_{n \geq 0}\) is (a realization of the control process) \(\bpi\) if \(\cL( \alpha_n\, | \, \sigma\{\sU,\underline{\vartheta}_{n-1},\underline{\vartheta}^0_n,\underline{\epsilon}_n,\underline{\epsilon}^0_n\} ) = \pi_n(\xi_n),\) \(\PP-a.s.\), for all \(n \geq 0\). We now consider closed-loop policies.

The choice of this form of closed-loop Markov strategy function is suggested by the mean-field nature of the dynamics ?? and the form of the costs ?? . Taking values in \(\cP(A)\) instead of \(A\) indicates that the strategy functions are mixed. Typically, they suggest that at each time \(n \geq 0\), the action \(\alpha_n \in A\) taken according to such a policy should be sampled from a probability measure depending directly on the values of \(X_n^{\balpha, \mu_0}\) and \(\PP^0_{X_n^{\balpha, \mu_0}}\), and the random variable \(\vartheta_n^0\) used by the level-\(1\) controller to randomize their choice. We formalize this procedure in Definition 2.5 below.

In Definition 2.5, we can view the state and action processes as constructed simultaneously by alternatively invoking the system dynamics ?? and the sampling procedure consistent with ?? . A convenient way to construct an action process \(\balpha\) satisfying ?? is to use the Blackwell-Dubins Lemma  1.1. Indeed, if \(\rho_A\) is the Blackwell-Dubins function of \(A\) and the uniformly distributed random variables \(U_n\) is given by \(U_n =h(\vartheta_n)\), we can choose \(\alpha_n= \rho_A \big( \pi_n \big( X_n, \, \PP^0_{X_n} ,\, \vartheta_n^0 \big), U_n \big)\), \(\PP\)-a.s., \(n \geq 0\).

2.2.2.2 Optimization and value functions

We now complete the definition of the agent level (level-0) MFC problem by stating the optimization problems with open-loop and closed-loop policies. First, we introduce the value function associated to an action process.

The value of \(J^\balpha(\mu)\) is well-defined for every \(\mu \in \cP(S)\) because \(f\) is measurable and bounded. Furthermore, this value depends only on the sequence of time-marginal laws of the \(\cP(S\times A)\)-valued process \(\Bigl(\PP^0_{(X_n^{\balpha, \mu_0},\alpha_n)} \Bigr)_{n\ge 0}\).

For any open-loop or closed-loop policy \(\bpi\), we can show that \(\mathbb{P}^0_{(X_n,\alpha_n)}\) depends on the action process only through the policy \(\bpi\) provided \((\bX,\balpha)\) is generated by \(\bpi\). See Lemma 2.2 in the appendix for a proof for open-loop policy, and a similar proof can be derived for closed-loop policy. As a consequence, we define, for any level-\(0\) action process \(\balpha\) generated by \(\bpi\), \[J^{\bpi}(\mu) := J^{\balpha}(\mu), \qquad \mu \in \cP(S).\] Accordingly, we define the optimal open-loop value function and the optimal closed-loop value function as: \[J^{\tinyol,*}(\mu) = \inf_{ \bpi \in \bPi^\tinyol } J^{\bpi}( \mu ), \qquad J^{\tinycl,*}(\mu) = \inf_{ \bpi \in \bPi^\tinycl } J^{\bpi}( \mu ), \qquad \mu \in \cP(S),\] which are finite because we assume that the one-stage cost function \(f\) is bounded and \(\gamma \in (0,1)\).

2.2.3 Mean-Field MDP↩︎

In this section, we introduce the Markov Decision Process (MDP) which we use to identify the optimal closed-loop Markov policies for our original MFC model. This MDP is obtained by lifting the agent level (level-0) dynamical system to a level-1 dynamical system on a space of probability measures.

2.2.3.1 Mean-field MDP framework

The key observation is that, for a closed-loop Markov policy \(\bpi \in \bPi^\tinycl\), the associated value function \(J^{\bpi}\) can be viewed as the value function of an MDP with a state process \((\PP^0_{X_n})_{n \geq 0}\) on state space \(\cP(S)\), and an action process \(\bigl(\PP^0_{(X_n,\alpha_n)}\bigr)_{n \geq 0}\) with values in a new action space \(\cP(S \times A)\). Actions need to be consistent with the state in the sense that the first marginal of any action which can be taken while in a given state, has to be equal to the state itself. We provide a rigorous definition of this MDP, and we show its connection to the original MFC model.

The lifted consists of a six-tuple \((\bar S, \bar A, \bar \Gamma, P, \bar f, \gamma)\) as described below:

• The state space is the Polish space \(\bar S := \cP(S)\); a generic element in \(\bar S\) is denoted by \(\mu\).

• The action space is the Polish space \(\bar A := \cP(S \times A)\); a generic element in \(\bar A\) is denoted by \(\bar a\).

• The control constraint is a set-valued function \({\bar U}\) from \({\bar S}\) into the set of non-empty subsets of \({\bar A}\) defined by: \[\label{eq:Ubar} {\bar U}(\mu) :=\{{\bar a}\in{\bar A};\; \text{pr}_1({\bar a})=\mu\}, \qquad \forall \, \mu\in{\bar S};\tag{26}\] where \(\pr_1: \bar A \to \bar S\) is the projection function that maps \(\bar a \in \bar A\) onto its first marginal distribution on \(S\). We shall also use the notation \[\label{eq:bar95Gamma} \bar\Gamma :=\{(\mu,{\bar a})\in{\bar S}\times{\bar A};\; {\bar a}\in{\bar U}(\mu)\}.\tag{27}\] for the graph of the constraint.

• A transition probability kernel \(P: \bar \Gamma \to \cP( \bar S)\), which is a measurable function.

• A one-stage cost function \(\bar f: \bar \Gamma \to \RR\), which is a bounded measurable function.

• The discount coefficient \(\gamma \in (0,1)\).

Here and in the following we denote by the push-forward of a measure \(\nu\) by a measurable function \(g\). We can check that \(\bar F\) is Borel measurable; see e.g. [95].We chose to characterize the dynamics by the transition kernel \(P\) (and not only by the system function \(\bar F\)) for two reasons: 1) to conform with the standard literature on MDPs which typically uses transition kernels to system functions; 2) we shall restrict ourselves to Markovian policies and control processes given by feedback functions of the state for the analysis of the lifted MDP.

To highlight the tight connections between the lifted MFMDP and the original MFC, we define the notion of mixed strategy and mixed Markov policy for MFMDP.

We shall identify every \(\varphi \in\overline{\Pi}^{p}\) with the corresponding level-1 mixed strategy \(\bar\pi\) defined by \(\bar \pi(\mu) = \delta_{ \varphi(\mu) }\) for \(\mu \in \bar S\). For the sake of brevity, we sometimes omit the term “mixed" or”randomized” when \(\bar \pi \in \overline{\Pi}\), but we always keep the term “pure" or”non-randomized” when \(\bar \pi \in \overline{\Pi}^{p}\). We now define policies as sequences of strategy functions.

We restrict ourselves to these policies and refrain from using history dependent policies because we are mostly interested in optimizing value functions, and we know that for MDPs like our lifted MFMDP, for each history dependent policy, there exists a Markov policy with the same value function (as defined below in Definition 2.11). See for example [95].

2.2.3.2 Assumptions and optimization problem for MFMDP

Some of our proofs rely on the following assumptions which we state here for the sake of later reference.

To show existence of an optimal policy, we will make use of the following extra assumption:

As a consequence of Lemma 2.4, the value function \(\bar J^{\bar \bpi}\) given in ?? is well defined because the expectation in ?? does not depend upon the particular choice of the pair of state action processes \((\bmu,\boldsymbol{\bar a})\) generated by \((\bar \bpi, \mu)\).

We will show below in Theorem 2.1 that the value function \(\bar J^{\bar \bpi}\) is lower semi-continuous when the policy is an optimal policy \(\bar\bpi = \bar\bpi^*\).

With the value function for MFMDP at hand, we define the as: \[\bar J^{*}(\mu) = \inf_{ \bar \bpi \in \bar\bPi}{\bar J}^{\bar \bpi}(\mu), \qquad \mu \in \bar S.\]

2.2.3.3 Dynamic Programming principle for MFMDP

We state and prove the DPP for the optimal value function with Borel measurable mixed Markov policies. We stress the fact that we are working in the lifted MFMDP setting by using overlines in the notations of all the quantities relative to the MFMDP. Still, it is worth noticing that, due to the presence of the common noise, our lifted model is still a stochastic model while in general, one lifting is enough to get to a deterministic dynamical model. This is for example the case in the monograph [95] to which we refer frequently for basic results.

The DPP ?? is known to hold for universally measurable policies and lower semi-analytic one-stage cost function. See for example [95]. The gist of the above theorem is to show that it also holds for Borel measurable policies.

2.2.4 Relations between the models↩︎

In the previous section, we defined the level-1 lifted MFMDP and we proved some of its main properties. We now reconnect this model with the original agent level MFC problem and derive properties of the latter from what we just found out about the lifted model

2.2.4.1 Relations between MFC closed-loop policies and MFMDP policies

We start by highlighting that, intuitively, a closed-loop policy for the MFC can be viewed as sampled from a policy for the MFMDP by picking the common randomness. The following definition formalizes this idea.

Note that if \(\bpi\) and \(\bar\bpi\) correspond to each other and if \(\bpi\) is stationary, then so is \(\bar\bpi\). Conversely, if \(\bpi\) and \(\bar\bpi\) correspond to each other and if \(\bar\bpi\) is stationary, then \(\bpi\) can be chosen to be stationary. We will provide a specific construction in equation 30 below. The main result of this section is the following.

In preparation for the proof of this result, we start with two useful technical lemmas describing the properties of conditional distributions of the level-0 state and action processes.

Intuitively, \((\kappa_n)_{n \ge 0}\) plays the role of the conditional law of the action process. We will come back to this interpretation in the proof of Lemma 2.8. We are now ready to prove Theorem 2.2.

At this stage, we need to emphasize the crucial role played by the common randomization provided by the sequence \((\vartheta^0_n)_{n\ge 0}\). Its presence is what allowed us to prove that the value of a policy for the lifted MDP can always be achieved by a closed loop policy of the original MFC.

Notice that the above result proves equality of the value functions, policy by policy, which is much stronger that equality of the value functions after taking the optimum values.

Saying that the policy does not depend upon the common randomization amounts to say that (by an abuse of notation) for every \(n\ge 0\), the strategy function takes the form \(\pi_n: S\times \cP(S)\to \cP(A)\).

2.2.4.2 Relations between MFC closed-loop and open-loop policies

We first prove existence of optimal closed-loop Markov policies and then, equality of the open loop and closed loop value functions.

We recall that by “pure” we mean that the policy does not depend upon the common randomization. However, the representative agent actions can still be randomized.

We now show equality of the open and closed loop optimal value functions of the original agent level MFC optimization problem.

Note that when both Assumption 2.1 and Assumption 2.2 hold, then combining Theorem 2.3 with Proposition [pr:existence95opt95policy95mean95field] yields the existence of a pure stationary closed loop Markov policy for the original MFC that achieves the value \(J^{\tinyol,*}\).

The first part of Theorem 2.3 is given by Lemma 2.7 below. For the converse direction, the result is obtained by combining Lemma 2.8 below and Theorem 2.2. Under the additional compactness assumption, the existence of an optimal closed-loop policy stems from Theorem 2.2, which entails the existence of an optimal open-loop policy by Lemma 2.7.

Again, while this type of equality between the optimal value functions could be expected to hold under different assumptions and without the central randomization, we prove it here by leveraging the equalities proven in Lemma 2.7 and Lemma 2.8 policy by policy, before computing optima over sets of policies.

First, it is expected that every closed-loop policy for the MFC can be viewed as an open-loop policy for the MFC, which leads to the following result in terms of value functions.

Next, every open-loop policy for the MFC corresponds to a policy for the MFMDP, as we show in the following result.

2.2.5 Mean-Field Q-Function↩︎

In preparation for our discussion of learning in the third part of the monograph, and in line with our analysis of the first theoretical model of MFC, we introduce the state-action value function, for the lifted MFMDP. Using standard arguments, we derive its first theoretical properties. This will pave the way to the algorithmic search for approximations of the solution of the MFC problem in a model-free setting, i.e. assuming the model is unknown.

2.2.5.1 State-action value function

In order to take full advantage of the strongest results proven so far, from now on we assume that both assumptions 2.1 and 2.2 hold. Under these assumptions, recall that the DPP given by Theorem 2.1 holds. In this section, we restrict ourselves to pure stationary policies. Recall that when \(\bar \bpi = (\bar \pi, \bar \pi, \ldots) \in \bar \bPi\) is stationary, we use the notation \(\bar J^{\bar \pi} = \bar J^{\bar \bpi}\). Also, recall that we use the notation \(\bar \Pi^p\) for the set of pure stationary strategy functions.

For each \(\bar \pi \in \bar \Pi^p\), the mapping \(\bar S \ni \mu \mapsto \delta_{\bar \pi(\mu)} \in \cP(\bar A)\) which assigns to each \(\mu \in \bar S\) the Dirac point mass at the point \(\bar \pi(\mu) \in \bar A\) is measurable by definition of the Borel \(\sigma\)-field of \(\cP(\bar A)\).

Now, for each \(\bar\pi \in \bar \Pi^p\), we introduce the state-action value function \(\bar Q^{\bar \pi} : \bar\Gamma \to \RR\) defined by: \[\label{eq:Q95function95stationary95control} \bar Q^{\bar \pi} (\mu, \bar a) := \bar f(\mu, \bar a) + \sum_{n \geq 1} \gamma^{n} \EE [\bar f(\mu_n, \bar\pi(\mu_n))], \qquad (\mu, \bar a) \in \bar\Gamma,\tag{34}\] where the process \((\mu_n)_{n\geq 0}\) starting at \(\mu_0 = \mu\) satisfies \(\mu_1 = \bar F(\mu, \bar a, \varepsilon^0)\) with \(\cL(\varepsilon^0) = \nu^0\), and for every \(n \geq 1\): \(\mu_{n+1} = \bar F \bigl( \mu_n, \bar \pi(\mu_{n}), \varepsilon_{n+1}^0 \bigr).\) Next we define the optimal state-action value function by: \[\label{fo:optimal95Q} \bar Q^* (\mu, \bar a) := \inf_{\bar\pi\in\bar \Pi^p}\bar Q^{\bar \pi} (\mu, \bar a), \qquad (\mu,\bar a)\in \bar\Gamma.\tag{35}\]

The goal of this section is to prove the following dynamic programming principle for \(\bar Q^*\).

We will prove this result by showing that \(\bar Q^*\) is the unique fixed point of the state-action Bellman operator \(T\) defined on the set \(\Lsc(\bar\Gamma)\) of bounded lower semi-continuous functions on \(\bar\Gamma\), by: \[\label{fo:state-action95Belman} [T \bar Q](\mu,\bar a) := \bar f(\mu,\bar a) +\gamma\EE \Bigg[ \inf_{\bar a'\in \bar U \big( \bar F(\mu, \bar a, \varepsilon^0) \big) } \bar Q\bigl(\bar F(\mu,\bar a,\varepsilon^0),\bar a' \bigr) \Bigg],\qquad (\mu,\bar a)\in \bar\Gamma.\tag{36}\] We first justify in Lemma 2.9 below the fact that the operator \(T\) is well-defined. Since \(\bar\Gamma\) is a closed subset of the Polish space \(\bar S \times \bar A\), it is a Borel space on its own, and the space \(\cBD_u(\bar\Gamma)\) of bounded universally measurable real-valued functions on \(\bar\Gamma\) endowed with the sup norm \(\| f \|_{\infty} = \sup_{(\mu, \bar a) \in \bar\Gamma} | f(\mu, \bar a) |\) is a Banach space. While the set \(\Lsc(\bar\Gamma)\) is not a vector space, it is a closed subset of \(\cBD_u(\bar\Gamma)\), hence a complete metric space for the metric \(d_\infty(f, f') = \| f - f' \|_{\infty}\).

Using the Markov property, we can rewrite the state-action value function \(\bar Q^{\bar \pi}\) in terms of the state value function \(\bar J^{\bar \pi}\): \[\label{eq:connection95J95Q} \bar J^{\bar \pi}(\mu) = \bar f(\mu, \bar \pi(\mu) ) + \gamma \EE [ \bar J^{\bar \pi}( \bar F(\mu, \bar \pi(\mu), \varepsilon^0) ) ] = \bar Q^{\bar \pi}\bigl( \mu, \bar \pi(\mu)\bigr), \qquad \bar \pi \in \bar \Pi^p, \mu \in \bar S.\tag{37}\] Now, for the optimal value functions, we have the following.

We can now complete the proof of Theorem 2.4.

2.2.6 Notes and Complements↩︎

MFMDPs have been studied at least since the work of Gast et al. [102], [103] who studied the convergence of a finite-agent MDP towards a MFMDP. More recently, Gu et al. studied in [15], [16] dynamic programming for the state-action value function arising in MFC problems without common noise, and they proposed RL methods for this setting. Following the posting on arxiv of the original version of our paper [5], Motte and Pham tackled in [4] the open loop case in the presence of common noise, in a model very similar to the model considered in this chapter. They proved a propagation of chaos result in the spirit of Proposition [pr:chaos] of Chapter 2.1, a DPP for the MFMDP, as well as the equality of optimal value functions for open-loop and closed-loop policies in their model. In this chapter, we proposed a finer model by paying close attention to all the sources of randomization, and by relying on fine measure theoretical arguments. Comparing to [4], while they prove equality of the optimal value functions without the central randomization, the latter allows us to prove the identity of the value functions, policy by policy, even before taking the optimum values, allowing for a deeper understanding of the agent level MFC, as did for example in Corollary 2.1.

Note that when the state space \(S\) is discrete, assuming continuity of \(F(\cdot, \cdot, \cdot, e, e^0)\) implies that the dependence on the mean-field term is basically constant when the other inputs are fixed. This is very much in line with the fact that many MFG and MFC models appearing in the literature do not have mean-field interactions in the dynamics of the states, the interactions appearing explicitly only in the cost functions. See for example the original models in [20] as well as the applications to Bertrand and Cournot competitions [104], [105], and the models of crowd motion [31] or flocking [106], [107]. To cover models with interactions in the dynamics of the states in a finite state space setting, one should use, instead of the assumptions mentioned in Remark [rem:assumptions-mfc], other assumptions that are sufficient to guarantee Assumption 2.1. In particular, it is expected that the idiosyncratic noise has a regularizing effect that helps ensuring continuity of \(\bar F\) even when \(F\) is discontinuous. For instance [4] does not require continuity of \(F\) but only some kind of almost everywhere continuity.

2.3 A Linear Quadratic Mean Field Model with Common Noise↩︎

2.3.1 Introduction↩︎

The thrust of the chapter is to properly introduce the notations and the definitions of the linear-quadratic (LQ) mean field control (MFC), or LQMFC, stochastic optimization problem which was introduced in broad strokes in Chapter 2.1. After arguing that the model appears naturally as the limit of a large multi agent model with a common noise, we prove that the optimal admissible controls exist, are unique, and linear in the state and its conditional mean processes. We formulate and prove the necessary and the sufficient conditions for a discrete-time form of the Pontryagin maximum principle for our LQMFC model. This result implies existence and uniqueness of the optimal control process. Incidentally, this provides a crucial link between the MFC problem and learning by a finite population of agents, as it gives an approximately optimal control for the problem with a finite number of learners.

2.3.2 Models↩︎

We work with the model introduced at the beginning of Section 2.1.6 of Chapter 2.1 and we repeat the meaning of the notations only when necessary. Because the original probability space \((\Omega, \cF,\PP)\) is assumed to be rich enough, to facilitate the analysis for LQ models, we can introduce two initial random variables, denoted by \(\varepsilon_0^0\) and \(\varepsilon_0\), called initial perturbations, that are independent of each other and independent of the noise sequences \(\bvarepsilon=(\varepsilon_n)_{n\ge 1}\) and \(\bvarepsilon^0=(\varepsilon^0_n)_{n\ge 1}\). The random variable \(\varepsilon_0^0\) represents the common shock at initial time from the environment. We distinguish the two initial perturbations \(\varepsilon_0\) and \(\varepsilon_0^0\) from the aforementioned noise terms. They are not assumed to follow the mean-zero distributions of the noise sequences. Also, they can be constant, i.e. random variables with a Dirac measure as distribution. This is the case when we want the LQ model to start with a deterministic initial condition.

In this chapter, we assume that the initial perturbations \(\varepsilon_0\) and \(\varepsilon_0^0\), and all noise terms \((\varepsilon_{n+1})_{n \geq 0}\) and \((\varepsilon_{n+1}^0)_{n \geq 0}\) are sub-Gaussian random vectors. Recall Definition 1.20 at the end of Chapter 2.1.

We shall extensively use the filtrations \((\sigma\{\underline\epsilon_n,\underline\epsilon^0_n\})_{n \geq 0}\) and \((\sigma\{\underline\epsilon^0_n\})_{ n \geq 0}\) that did not play much of a role in the analysis of Chapter 2.2. Recall the meaning of the underline, namely \(\underline\varepsilon_n=(\varepsilon_0,\varepsilon_1,\cdots,\varepsilon_n)\) and similarly for \(\underline\varepsilon_0^0\).

2.3.2.1 Model with a finite number of agents

As stated above, we justify the model by first considering the corresponding finite agent model and taking the limit as the number of agents goes to \(\infty\).

We assume that the state of the system is given by the \(N\)-tuple of individual agent states, and that the impact of the common environment is captured by a common noise included in the system functions by a transformation à la Blackwell-Dubins. Our desire for linear dynamics and mean field interactions dramatically restricts the form of the system function. The individual state dynamics in the stochastic optimization problem for \(N\) agents with linear dynamics and symmetric interactions are given by the system of equations: \[\begin{align} \label{fo:N-multi95state} X_{n+1}^i = {\mathrm A} X_n^i + \bar{\mathrm A} \bar{X}_n^{N} + {\mathrm B} a^i_n + \bar{{\mathrm B} } \bar{a}_n^N + \comnoise{n+1} + \idynoise{n+1}^{i},\qquad i=1,\cdots,N,\quad n\ge 0 \end{align}\tag{38}\] with initial conditions \(X_0^{i} = \varepsilon_0^0 + \varepsilon_0^{i}\), for \(i = 1, \ldots, N\), where \[\bar{X}^N_n = \frac{1}{N} \sum_{i=1}^N X_n^{i} \qquad\text{and}\qquad \bar{a}^N_n = \frac{1}{N}\sum_{i=1}^N a_n^{i}\] denote the sample averages of the individual states and actions, and \(( \varepsilon_0^{i}, (\idynoise{n+1}^{i})_{n \geq 0})_{i=1,\ldots, N}\) are \(N\) identical and independent copies of the initial perturbation and the idiosyncratic noise sequence \(( \varepsilon_0, (\varepsilon_{n+1})_{n \geq 0} )\).

Averaging the \(N\) dynamical equations above gives the time evolution of the sample average of the state as: \[\label{fo:N-mean95state95dynamics} \bar{X}^N_{n+1} = ({\mathrm A} + \bar{\mathrm A} ) \bar{X}^N_n + ({\mathrm B} + \bar {\mathrm B} ) \bar{a}^N_n + \varepsilon^0_{n+1} + \frac{1}{N}\sum_{i=1}^N \varepsilon^{i}_{n+1},\tag{39}\] with initial value \(\bar{X}^N_0 = \comnoise{0} + \frac{1}{N} \sum_{i=1}^N \idynoise{0}^{i}\). Note that for \(n\) fixed, the classical law of large numbers says that the last term in 39 should converge when \(N\to\infty\) toward the common expectation of the \(\varepsilon^{i}_{n+1}\) which is \(0\). On the other hand, the random shock \(\varepsilon^0_{n+1}\) affects all the agents and its impact should not disappear in the limit \(N\to\infty\).

Letting \(\underline{X}_n= Vect\big(X_n^{1}, \cdots, X_n^{N} \big)\) and \(\underline{a}_n= Vect \big( a_n^{1}, \cdots, a_n^{N} \big)\) for \(n\ge 0\), we can rewrite the system 38 of \(N\) dynamical equations in a single dynamical equation: \[\label{fo:N-vector95state} \underline{X}_{n+1} = {\mathrm A} ^N \underline{X}_n + {\mathrm B} ^N \underline{a}_n + \underline{\varepsilon}^0_{n+1} + \underline{\varepsilon}_{n+1}\tag{40}\] with \(\underline{\varepsilon}_{n+1} = Vect( \idynoise{n+1}^{1}, \ldots ,\idynoise{n+1}^{N} )\), \(\underline{\varepsilon}^0_{n+1} = Vect( \comnoise{n+1},\ldots,\comnoise{n+1} )\), \({\mathrm A} ^{N} = I_N \otimes {\mathrm A} + \frac{1}{N} \bone_N \otimes \bar{{\mathrm A} }\) and \({\mathrm B} ^{N} = I_N \otimes {\mathrm B} + \frac{1}{N} \bone_N \otimes \bar{{\mathrm B} }\). In other words, the matrix \({\mathrm A} ^{N}\) is the sum of the block-diagonal matrix with \(N\) blocks identically equal to \({\mathrm A}\), and of the matrix with \(N \times N\) blocks \(\frac{1}{N}\bar{{\mathrm A} }\), and similarly for the matrix \({\mathrm B} ^{N}\).

Accordingly, since we assume that the agents cooperate with each other (instead of competing against each other), the goal is to minimize the social cost of the population, defined as \[\label{fo:social95cost95of95population} J^N( \underline{\ba} ) = \EE \Big[ \sum_{n \geq 0} \gamma^{n} \bar f^N(\underline{X}_n, \underline{a}_n) \Big]\tag{41}\] over the set of admissible control processes \(\underline{\ba} = (\underline{a}_n)_{n \geq 0}\) consisting of \(\RR^{\ell N}\)-valued processes adapted to the filtration \((\sigma\{\underline\epsilon^0_n,\underline\epsilon_n\})_{n \geq 0}\) for which \(J^N(\underline{\ba}) < \infty\). Here, the one-stage social cost function \(\bar f^N : \RR^{dN} \times \RR^{\ell N} \to \RR\) is \[\label{fo:N-cost} \bar f^{N}( \underline{X}_n, \underline{a}_n) := \frac{1}{N} \sum_{i=1}^N c \big( X_n^i, \bar{X}_n^N, a_n^i, \bar{a}_n^N \big) = \underline{X}_n^\top Q^N \underline{X}_n + \underline{a}_n^\top R^N \underline{a}_n,\tag{42}\] where \(Q^N = \frac{1}{N} I_N \otimes Q + \frac{1}{N^2} \boldsymbol{1}_{N} \otimes \bar{Q} \in \RR^{dN \times dN}\) and \(R^N = \frac{1}{N} I_N \otimes R + \frac{1}{N^2} \boldsymbol{1}_{N} \otimes \bar{R} \in \RR^{\ell N \times \ell N}\).

This is now a classical infinite-horizon discounted-cost linear-quadratic control problem with \((d N)-\)dimensional state process \(\underline{\bX} = (\underline{X}_n)_{n \geq 0}\) following dynamics 40 . Under Assumptions 3.1 and 3.2 stated below, an optimal control \(\underline{\ba}^{*,N}\) exists and is necessarily of the form \(\underline{a}^{*,N}_n = \Phi^{*,N} \underline{X}_n\) for a deterministic constant matrix \(\Phi^{*,N} \in \RR^{(\ell N) \times (d N)}\).

2.3.2.2 The mean-field control model

Motivated by the analysis of the large \(N\) limit of the above model, we envision a generic agent the state of which, say \(X_n\in\RR^d\), evolves over time according to the equation: \[\begin{align} X_{n+1} & = {\mathrm A} X_n + \bar{{\mathrm A} }\bar{X}_n + {\mathrm B} a_n + \bar{{\mathrm B} }\bar{a}_n + \varepsilon^0_{n+1} + \varepsilon_{n+1}, \label{fo:MKV-state} \end{align}\tag{43}\] for \(n \geq 0\), with an initial state \(X_0 =\varepsilon^0_{0}+\varepsilon_{0}\). Let \(\tilde{\mu}^0_0 = \cL(\varepsilon_0^0)\) and \(\tilde{\mu}_0 = \cL(\varepsilon_0)\) be their respective distributions. The terms \(\bar{X}_n = \EE[X_n|\sigma\{\underline\epsilon^0_n\}]\) and \(\bar{a}_n = \EE[a_n|\sigma\{\underline\epsilon^0_n\}]\) appearing in the dynamics 43 of the state are the conditional means of \(X_n\) and \(a_n\) given \(\sigma\{\underline\epsilon^0_n\}\), the past of the common noise up to and including time \(n\). The Mean Field cost function is \[\label{fo:MKV-discounted95cost} J(\ba)=\EE \Big[ \sum_{n \geq 0} \gamma^n f(X_n, \bar{X}_n, a_n, \bar{a}_n) \Big].\tag{44}\] We choose a one-stage cost function \(f\) which is quadratic in the state and the action. However, in order to make the derivations in the proofs of this chapter slightly easier to follow, we depart slightly from the notation of the previous section on \(N\) agents and we choose: \[\label{fo:quadratic95cost} f(x, \bar{x}, a, \bar{a}) = (x-\bar x)^\top Q (x-\bar x) + \bar{x}^\top (Q+\bar Q) \bar{x} + (a-\bar a)^\top R (a-\bar a) + \bar{a}^\top (R+\bar R) \bar{a}\tag{45}\] where we use the superscript \({}^\top\) to denote the transpose of a vector or a matrix. Notice that the expected one stage cost at time \(n\) reads \[\label{fo:expression95expected95instateneous95cost} \begin{align} \EE[ f(X_n, \bar{X}_n, a_n, \bar{a}_n)] &= \EE[ (X_n-\bar X_n)^\top Q (X_n-\bar X_n) + \bar{X}_n^\top (Q+\bar Q) \bar{X}_n\\ &\hskip 45pt + (a_n-\bar a_n)^\top R (a_n-\bar a_n) + \bar{a}_n^\top (R+\bar R) \bar{a}_n ]. \end{align}\tag{46}\] The set of admissible control processes is denoted by \(\cA_{ad}\) and defined as the collection of sequences of random variables on \((\Omega,\cF,\PP)\) adapted to the filtration \((\sigma\{\underline\epsilon_n,\underline\epsilon^0_n\})_{n \geq 0}\) and square integrable (a.k.a. \(L^2-\)discounted integrable):

\[\label{def:admissible95control} \cA_{ad} := \left\{ (a_n)_{n \geq 0} \left|\, a_n \in \RR^{\ell} \text{ is } \sigma\{\underline\epsilon_n,\underline\epsilon_n^0\}-\text{measurable}, \quad \EE \Big[ \sum_{n \geq 0} \gamma^n \Vert a_n \Vert^2 \Big] <\infty \right. \right\}.\tag{47}\] We also consider the set \(\cX\) of \(L^2-\)discounted-integrable processes in \(\RR^{d}\) defined by: \[\label{def:L295discoutned95integratble95process} \cX := \left\{ (X_n)_{n \geq 0} \left|\, X_n \in \RR^{d} \text{ is } \sigma\{\underline\epsilon_n,\underline\epsilon_n^0\}-\text{measurable}, \quad \EE \Big[ \sum_{n \geq 0} \gamma^n \Vert X_n \Vert^2 \Big] <\infty \right. \right\}.\tag{48}\] Note that both \(\cA_{ad}\) and \(\cX\) are linear spaces. The MFC problem is then: \[\label{pb:MFC95L295admissible95control} \inf_{\bu \in \cA_{ad}} J(\bu).\tag{49}\]

The above stochastic control problem has some unique characteristics: 1) it is set in infinite horizon; 2) it includes a common noise \(\bvarepsilon^0\), 3) The optimization is set over open-loop controls, and 4) the interaction is not only through the conditional mean of the state, but also through the conditional mean of the control.

2.3.2.3 Assumptions and the well-definedness of the problem

Throughout the whole chapter, we impose the following assumptions which we already stated in the introduction of the model in Chapter 2.1.

As we shall see in the following sections, Assumption 3.1 is required to show that the state process is well-defined and the discounted mean field cost is finite. Taking advantage of the strong convexity in the control variables provided by \(R\succ0\) and \(R+\bar R\succ0\), we show in Section 2.3.3 the existence and uniqueness of the optimal control process for the MFC problem.

We will sometimes use the notation \(\bX^\ba\) if necessary to emphasize the fact that the state process following dynamics 43 is controlled by the admissible control process \(\ba \in \cA_{ad}\). Lemma 3.1 and Proposition [pr:admissible95of95X] below justify the consistency of considering the set \(\cA_{ad}\) of admissible control processes and the set \(\cX\) for controlled state processes, and that the discounted mean field cost \(J(\ba)\) for admissible control process \(\ba \in \cA_{ad}\) is finite.

The following result provides a simple procedure to identify \(L^2-\)discounted-integrable state processes.

In Section 2.3.3, we will show existence and uniqueness of an optimal admissible control process for the MFC problem 49 using a discrete-time version of Pontryagin’s maximal principle. Moreover, the optimal control process \(\ba\) will be shown to be in a feedback form and to be linear in \((\bX^\ba, \bar \bX^\ba)\), with coefficients independent of time. To rigorously establish these results with the aforementioned admissible controls, we define the following set of admissible control parameters.

The set \(\Theta\) is a convex subset of \(\RR^{\ell \times d} \times \RR^{\ell \times d}\) because the spectral norm \(\| \cdot \|\) is convex. It is worth noting that in discrete-time control problems, contrary to their continuous-time analogs, there is no issue concerning the existence of the state process controlled by a control process given in feedback form. The following lemma states that when the feedback form is defined through a policy function \(\pi_\theta : (x, \bar x) \mapsto \pi_\theta(x, \bar x)\) defined below, and the control parameter \(\theta\) involved within \(\pi_\theta\) is admissible, the control process \(\bu\) is also admissible.

We shall from time to time use the notation \(\ba^\theta\) and \(\bX^\theta\) to stress the dependence on \(\theta\). The collection of (admissible) parameterized control processes is denoted by \(\cA^\Theta_{ad} = \{ \ba^\theta \, | \, \theta \in \Theta \} \subset \cA_{ad}\).

2.3.3 Analysis of the MFC Optimization Problem↩︎

In this section, we provide an analysis of the MFC problem. To simplify the exposition we use the following notations: \[\tilde{\mathrm A} = {\mathrm A} + \bar {\mathrm A},\;\; \tilde{\mathrm B} = {\mathrm B} + \bar {\mathrm B}, \;\; \tilde{Q} = Q + \bar Q, \;\;\text{and}\;\; \tilde{R} = R + \bar R.\]

We will use the following assumption.

The main result of this section is the following.

2.3.3.1 Admissibility of the optimal control

Before proving optimality, we record the algebraic identities relating the Riccati construction to the feedback matrices in Theorem 3.1. The membership \(\theta^*\in\Theta\) is imposed in Assumption 3.3. To start, we construct \((P, \bar P)\) based on the solutions \((P^{*, y}, P^{*, z})\) of another pair of discrete-time algebraic Riccati equations (DAREs): \[\label{fo:DARE95closeloop95feedback95excat95PG} \left\{ \begin{array}{rcl} P^{*, y} &=& Q + \gamma {\mathrm A} ^\top P^{*, y} {\mathrm A} - \gamma^2 {\mathrm A} ^\top P^{*, y} {\mathrm B} ( R + \gamma {\mathrm B} ^\top P^{*, y} {\mathrm B} )^{-1} {\mathrm B} ^\top P^{*, y} {\mathrm A} , \\ P^{*, z} &=& \tilde{Q} + \gamma \tilde{{\mathrm A} }^\top P^{*, z} \tilde{{\mathrm A} } - \gamma^2 \tilde{{\mathrm A} }^\top P^{*, z} \tilde{{\mathrm B} } ( \tilde{R} + \gamma \tilde{{\mathrm B} }^\top P^{*, z} \tilde{{\mathrm B} } )^{-1} \tilde{{\mathrm B} }^\top P^{*, z} \tilde{{\mathrm A} }. \end{array} \right.\tag{53}\]

2.3.3.2 Optimality of the control process

In this subsection, we show existence and uniqueness of an optimal control using a discrete-time version of Pontryagin’s maximal principle. For convenience, we use the notation \(\check {\mathrm A} = {\mathrm A} - I_d\) and \(\zeta = (x, \bar x, a, \bar a) \in \RR^{d} \times \RR^{d} \times \RR^{\ell} \times \RR^{\ell}\). We define the drift function \(b\) by \[\label{fo:drift95function95b} b(\zeta) = b(x, \bar x, a, \bar a) = \check{{\mathrm A} } x + \bar{{\mathrm A} } \bar x + {\mathrm B} a + \bar{{\mathrm B} } \bar{a}\tag{54}\] and the Hamiltonian function \(h\) by: \[\label{fo:Hamilton95function95h} h(\zeta, p) = h(x, \bar x, a, \bar a, p) = b(\zeta) \cdot p + f(\zeta) - \delta\; x \cdot p\tag{55}\] where \(p \in \RR^{d}\), \(\delta = (1 - \gamma) / \gamma\), and \(f(\zeta)=f(x, \bar x, a, \bar a)\) is the one-stage cost function defined by 45 . Then, for an admissible control process \((a_n)_{n \geq 0}\), the state dynamics 43 of the controlled state process \((X_n)_{n \geq 0}\) can be rewritten as \[\label{fo:delta95X95t95with95b} X_{n+1} - X_n = b(X_n, \bar X_n, a_n, \bar a_n) + \varepsilon_{n+1}^0 + \varepsilon_{n+1} = b(\zeta_n) + \varepsilon_{n+1}^0 + \varepsilon_{n+1}.\tag{56}\] To compute the partial derivatives of \(h\) with respect to its arguments \((x, \bar x, a, \bar a, p)\), we treat \(\zeta\) as a \((2d + 2\ell)\) dimension vector, \(Vect(\zeta) = (x^\top, \bar x^\top, a^\top, \bar a^\top)^\top \in \RR^{(2d + 2 \ell)}\), and by a slight abuse of notation we use the same notation \(\zeta\). Now, for every fixed \(p \in \RR^{d}\), we have \[\label{fo:gradient95h} \nabla_{\zeta} h (\zeta, p) = \begin{pmatrix} \partial_x h(\zeta,p)\\ \partial_{\bar{x}} h(\zeta,p)\\ \partial_{a} h(\zeta,p)\\ \partial_{\bar{a}} h(\zeta,p) \end{pmatrix} = \begin{pmatrix} (\check{{\mathrm A} }-\delta I_d)^\top p + 2Q (x-\bar{x}) \\ \bar{{\mathrm A} }^\top p - 2 Q (x - \bar x) + 2 \tilde{Q} \bar x \\ {\mathrm B} ^\top p + 2 R (a - \bar{a})\\ \bar{{\mathrm B} }^\top p - 2 R (a - \bar a) + 2\tilde{R} \bar a \end{pmatrix},\tag{57}\] and the Hessian of \(h\) with respect to \((a, \bar a)\) is given by \[\label{fo:hession95h95u} \nabla^2_{(a,\bar{a}),(a,\bar{a})} h(\zeta,p)= \begin{pmatrix} 2R & -2R \\ -2R & 2(R+\tilde{R}) \end{pmatrix},\tag{58}\] which is positive definite because \[\begin{pmatrix} a & \bar{a} \end{pmatrix} \begin{pmatrix} R & -R \\ -R & R+\tilde{R} \end{pmatrix} \begin{pmatrix} a \\ \bar{a} \end{pmatrix} =(a - \bar a) R(a - \bar a) + \bar a \tilde{R} \bar a \ge 0\] and equal to \(0\) only when \(a=\bar a=0\). Indeed, under Assumption 3.2, \(R\) and \(\tilde{R}\) are assumed to be (strictly) positive definite.

We now introduce the notion of adjoint process associated to a state process generated by an admissible control process.

The transversality condition should be viewed as a terminal condition controlling the possible growth of \(p_n\) when \(n\to\infty\). Notice that equation ?? can equivalently be written as: \[\label{fo:simpler95adjoint} p_n = \gamma \EE [ {\mathrm A} ^\top p_{n+1} + 2 Q X_{n+1} + \bar {\mathrm A} ^\top \bar p_{n+1} + 2 \bar Q \bar X_{n+1} \;| \; \sigma\{\underline\epsilon_n,\underline\epsilon_n^0\} ].\tag{59}\]

Using the adjoint process, we compute the Gateaux derivative of \(J\), and derive a necessary condition for optimality in the spirit of the classical Pontryagin’s maximum principle.

2.3.3.3 Identification of the optimal control

Let us now assume that \(\ba\) is an admissible optimal control and show that \(a_n\) can be written as a linear combination of the process \((p_n, \bar p_n)_{n\ge 0}\) where \(p_n\) is the adjoint process of \(\bX^{\ba}\), and as usual we use a bar on the top of a random variable to denote its conditional expectation with respect to \(\sigma\{\underline\epsilon^0_n\}\). By taking conditional expectation \(\EE[\,\cdot\,|\,\sigma\{\underline\epsilon^0_n\}]\) in equation 66 , we get \[({\mathrm B} + \bar {\mathrm B} )^\top \bar p_n + 2 (R + \bar R) \bar a_n = 0\] from which we derive: \[\bar a_n=-\frac{1}{2} (R +\bar R)^{-1}({\mathrm B} +\bar {\mathrm B} )^\top\bar p_n,\] and going back to 66 we finally get: \[a_n = -\frac{1}{2} R^{-1} {\mathrm B} ^\top p_n - \frac{1}{2} R^{-1} \big( \bar {\mathrm B}^\top - \bar R (R + \bar R)^{-1} ({\mathrm B} + \bar {\mathrm B} )^\top \big) \bar p_n \label{fo:control95ut95expression951}\tag{67}\] Using the notations ?? we get formula 52 for the process \(\ba = (a_n)_{n \geq 0}\). Given the well-known fact that at the optimum of standard LQ models, the adjoint processes are affine functions of the state, one may wonder if optimal controls can also be shown to be linear in the controlled state process \(\bX^\ba\) and the conditional mean process \(\bar{\bX}^\ba\). In order to prove that it is indeed the case, we introduce the following system of forward-backward stochastic equations for two \(L^2-\)discounted integrable processes \((\bX, \bp) \in \cX \times \cX\): \[\label{fo:FBDSE95state95adjoint} \left\{ \begin{array}{rcl} X_{n+1} - X_n &=& \check {\mathrm A} X_n + \bar {\mathrm A} \bar{X}_n + {\mathrm B} \Gamma (p_n - \bar{p}_n) + ({\mathrm B} + \bar {\mathrm B} ) \Lambda \bar{p}_n + \varepsilon_{n+1}^0 + \varepsilon_{n+1} \\ p_{n+1} - p_{n} &=& - \gamma \big[ (\check {\mathrm A} - \delta I_d)^\top p_{n+1} + 2Q X_{n+1} + \bar {\mathrm A} ^\top \bar{p}_{n+1} + 2 \bar Q \bar{X}_{n+1} \big] \\ && \hskip 180pt + Z_{n+1}^0 \varepsilon_{n+1}^0 + Z_{n+1} \varepsilon_{n+1}, \end{array} \right.\tag{68}\] where the adapted processes \((Z_n)_{n \geq 0}\) and \((Z_n^0)_{n \geq 0}\) are part of the solution, and where \(X_n\) satisfies the initial condition at time \(n=0\), and \(p_n\) satisfies the transversality condition which we view as a terminal condition at \(\infty\). This justifies our terminology of forward-backward system, abbreviated as FBDTSE (Forward-Backward Discrete-Time Stochastic Equations). Notice that the first equation in 68 is nothing but the state equation for a linear control of the form we want to restrict ourselves to. Also, by taking conditional expectations \(\EE[\cdot | \sigma\{\underline\epsilon_n,\underline\epsilon_n^0\}]\) of both sides of the second equation in 68 , we recover the fact that \(\bp\) is necessarily the adjoint process of \(\bX\).

The following proposition provides a solution to the forward-backward system 68 .

The following theorem shows that we can construct an optimal control process \(\ba\) for the MFC problem 49 based on a solution of the FBDTSE system 68 .

2.3.3.4 Proof of Theorem 3.1

We now provide the proof of the main result of this section, namely Theorem 3.1.

2.3.4 Notes and Complements↩︎

Discrete-time LQ mean field control problems without common noise are treated in [111], [112], where the analysis leads to decoupled Riccati equations. Continuous-time LQ mean field control problems and related McKean–Vlasov control formulations are studied, e.g., in [113]–[115]; see also [22] for a systematic treatment of the continuous-time theory.

Classical LQR problems without mean-field interactions have also been studied from the reinforcement learning perspective, in particular for learning linear feedback gains; see [25], [108], [116]. The policy-gradient analysis for the LQ mean field model considered here is developed in Chapter 3.2.

3 Numerical Implementations↩︎

3.1 Implementation of the First Abstract Model↩︎

3.1.1 Introduction↩︎

This chapter is the first of two dedicated to the numerical resolution of Mean Field Control (MFC) through reinforcement learning. It bridges the gap between the rigorous mathematical models developed in Chapter 2.2 and their practical implementation. As discussed in Chapter 2.2, finding optimal policies for mean field Markov decision processes (MFMDPs) requires overcoming significant computational hurdles, primarily due to the potential infinite-dimensional nature of the lifted state space consisting of probability measures. We revisit the concept of the “lifted” (level-1) problem, where the central planner controls the flow of distributions, and contrast it with the level-0 dynamics of individual agents. To address this challenge in the general setting, we present two distinct approaches.

First, we discuss a method based on space discretization, leading to a tabular version of the Q-learning algorithm. By projecting the mean-field distributions onto a finite grid (a discretized simplex), we can rigorously establish the convergence of the algorithm (Theorem 1.1). However, as we shall see, this approach suffers from the curse of dimensionality, making it impractical for problems where the level-0 state and action spaces are large.

Second, in order to circumvent the limitations of discretizations, we introduce function approximation techniques, specifically adapting the Deep Deterministic Policy Gradient (DDPG) algorithm. This approach operates directly on continuous state and action spaces by relying on artificial neural networks to approximate the Q-function (the critic) and the policy (the actor). While theoretical convergence guarantees for deep RL in the mean field setting remain an active area of research, the numerical experiments reported in Section 3.1.4 demonstrate the practical efficacy of this approach for models ranging from cybersecurity to distribution planning.

By exploring both a theoretically grounded discrete method and a scalable continuous method, this chapter equips the reader with the tools necessary to deploy Mean Field Reinforcement Learning algorithms in diverse and complex environments. These methods prepare the way for Chapter 3.2, where we specifically focus on Linear-Quadratic MFC problems for which we develop a more refined policy gradient analysis and provide stronger convergence guarantees.

3.1.2 Background on Reinforcement Learning Methods↩︎

Before delving into the specific algorithms tailored to MFC, we briefly review the core concepts of Reinforcement Learning (RL) that underpin the methodologies presented in this and the subsequent chapter.

3.1.2.1 The Reinforcement Learning paradigm

The core objective of RL is to solve an optimal control problem formulated as a Markov Decision Process (MDP). An MDP is defined by a state space \(S\), an action space \(A\), a transition probability \(P(x' | x, a)\), a one stage cost function \(f(x, a)\), and a discount factor \(\gamma \in (0,1)\). Given a policy \(\pi\) that determines the choice of actions \(a_n\in A\) based on the current value of the state \(x_n\in S\), the goal is to minimize the expected discounted cumulative cost: \[\label{fo:single95agent95objective} J(\pi) = \mathbb{E}^\pi \left[ \sum_{n=0}^\infty \gamma^nf(X_n, A_n) \Big| X_0=x_0 \right],\tag{76}\] where we use upper case characters to denote random variables and lower cases for specific values or realizations of these random variables. Thus typically \(x_{n}\) will denote a realization or a sample of the random variable \(X_n\) describing the state as a random variable in \(S\). So the description of the MDP stipulates that \(X_{n+1}\) is a random variable in \(S\) with distribution \(P(\cdot | x_n, a_n)\) when \(X_n=x_n\) and \(A_n=a_n\), and its realizations \(x_{n+1}\) are sampled from \(P(\cdot | x_n, a_n)\). The objective \(J(\pi)\) is analogous to the multi-agent reward maximization described in Chapter 1.1 (see for instance equation 1 ), but adapted to our cost-minimization framework.

A distinguishing feature of RL is its model-free approach. In many practical scenarios, the agent does not have access to the internal dynamics of the environment (the transition kernel \(P\)). Instead, they learn by interacting with the environment, observing states, taking actions, and incurring costs without a prior model of the environment’s transition or cost structure. This paradigm shift requires the agent to balance the exploration of possibly new parts of the state and action spaces, with the exploitation of known low-cost regions, a challenge that lies at the heart of the algorithms discussed in this chapter.

3.1.2.2 Value-based methods and Q-learning

Value-based methods rely on the evaluation of the state-action value function, or Q-function, the value \(Q(x,a)\) representing the expected cumulative cost incurred by starting in a state \(x\), taking an initial action \(a\), and then following a policy \(\pi\) thereafter: \[Q^\pi(x, a) = \mathbb{E}^\pi \left[ \sum_{n=0}^\infty \gamma^n f(X_n, A_n) \Big| X_0 = x, A_0 = a \right].\] The optimal Q-function, denoted \(Q^*\), satisfies the Bellman optimality equation: \[Q^*(x,a) = f(x, a) + \gamma \mathbb{E}_{X' \sim P(\cdot|x,a)} \left[ \inf_{a' \in A} Q^*(X', a') \right].\] A fundamental value-based algorithm is Q-learning [117], which iteratively updates an estimate of \(Q^*\) using temporal difference learning. During training, upon observing a transition \((x_n, a_n, c_n, x_{n+1})\), the Q-value is updated via: \[Q(x_n, a_n) \leftarrow (1-\eta_n) Q(x_n, a_n) + \eta_n \left( c_n + \gamma \inf_{a'} Q(x_{n+1}, a') \right),\] where \(\eta_n\) is the learning rate. For environments with finite state and action spaces, these Q-values can be stored in a table, justifying the terminology tabular Q-learning. A critical aspect of such algorithms is the exploration-exploitation trade-off, often managed via an \(\epsilon\)-greedy strategy, where the agent selects the action prescribed by the current policy with probability \(1-\epsilon\), and with probability \(\epsilon\), explores the action space by choosing an action uniformly at random. This tabular approach will be the foundation for the discretized simplex Mean Field Q-learning algorithm in Section [se:tabularQ-discrete].

3.1.2.3 Policy gradient and zeroth-order optimization

When dealing with continuous action spaces or highly parameterized policies (such as neural networks), value-based methods become computationally expensive due to the computation of the infimum \(\inf_{a'} Q(x', a')\) for each state \(x'\). Policy gradient methods address this by directly parameterizing the policy \(\pi_\theta(x)\) with weights \(\theta\) and performing stochastic gradient descent on the expected total cost \(J(\theta)\).

However, in many RL applications, calculating exact gradients is not possible because the environment acts as a black box. Zeroth-order optimization (also known as derivative-free optimization) overcomes this obstacle by estimating gradients using only cost evaluations. By applying small random perturbations \(u\) to the parameters and measuring the resulting cost \(J(\theta + u)\), one can construct an unbiased estimate of the gradient of a smoothed objective function. This approach will be exploited in the case of LQ MFC models in Chapter 3.2.

3.1.2.4 Actor-Critic and Deep RL

To combine the benefits of value-based and policy-based methods, Actor-Critic architectures maintain two separate structures: the actor, representing the policy parameterized by \(\theta^\pi\), and the critic, estimating the value function parameterized by \(\theta^Q\). The critic evaluates the current policy’s performance, providing a low-variance critique that the actor uses to guide its parameter updates.

When state and action spaces are large or continuous, tabular methods are infeasible due to the curse of dimensionality. Deep RL solves this quandary by using artificial neural networks as function approximators for the actor and the critic. A prominent algorithm in this class is DDPG [118], an actor-critic method specifically designed for continuous action spaces. DDPG utilizes experience replay buffers and target networks to stabilize the neural network training. This methodology underpins the Mean Field DDPG algorithm presented in Section [sec:DDPG-algo], enabling scalable numerical solutions applied to complex problems in Section 3.1.4.

3.1.2.5 Mean Field MDP: notations and environment

In the context of MFRL, we adapt the standard MDP paradigm to the lifted problem described in Chapter 2.2. We recall that the population distribution \(\mu_n \in \mathcal{P}(S)\) (where \(S\) is the individual state space) serves as the state of the MFMDP, often referred to as a level-1 state. The action, or level-1 action, is a joint distribution \(\bar a_n \in \mathcal{P}(S \times A)\) consistent with \(\mu_n\), meaning its first marginal is \(\mu_n = \text{pr}_1(\bar a_n)\). The environment in this setting is a mean field simulator that implements the lifted dynamics \(\bar F\): given the current distribution \(\mu_n\) and the lifted action \(\bar a_n\), it returns the next distribution \(\mu_{n+1}\) and the aggregate social cost \(\bar f(\mu_n, \bar a_n)\).

Ideally, this environment provides a perfect simulation of the infinite-population dynamics. However, in practice, simulating a continuous distribution \(\mu_n\) is often intractable, especially with high-dimensional state spaces. Consequently, numerical implementations frequently adopt a particle system simulator. In this approach, the environment simulates a system of \(N\) interacting agents \(\underline{X}_n = (X^1_n, \dots, X^N_n)\) and provides the empirical distribution \(\hat{\mu}^N_n = \frac{1}{N} \sum_{i=1}^N \delta_{X^i_n}\) as an observation. Under suitable propagation-of-chaos assumptions, and conditionally on the common information when common noise is present, such empirical distributions approximate the corresponding mean-field distribution as \(N \to \infty\).

In the numerical examples of Section 3.1.4 in this chapter, we primarily focus on simulating the mean field dynamics directly. In contrast, Chapter 3.2 investigates more granular configurations. Specifically, Section 3.2.1.4 discusses an idealized McKean-Vlasov simulator, while Section 3.2.1.5 introduces a finite-population simulator. The latter allows for a theoretical study and numerical illustration of the additional error induced by the size of the finite population.

3.1.3 Mean Field Q-Learning↩︎

3.1.3.1 Controls for finite state and action spaces

 

In the rest of this section, we assume that \(S\) and \(A\) are finite, we denote their cardinalities by \(|S|\) and \(|A|\) respectively, and we denote by \(x^{(1)},\dots,x^{(|S|)}\) and \(\alpha^{(1)},\dots,\alpha^{(|A|)}\) their elements. We first revisit the description of the action space and then propose two RL methods for this setting. Subsequently, we explain in Section [sec:DDPG-algo] how to adapt these RL techniques to continuous spaces.

Before introducing the mean-field Q-learning algorithm, we first provide a representation of the set \(\bar \Gamma \subseteq \bar S \times \bar A = \cP(S) \times \cP(S \times A)\) on which the \(\bar Q^*\) function is defined. Recall that the lifted state-action space \(\bar \Gamma\), defined in Equation 27 of Chapter 2.2, captures the requirement that the marginal of a lifted action must be consistent with the current mean field state.

Since we assume that \(S\) is finite, its lifted space \(\cP(S)\) can be identified with a simplex \(\mathfrak{S}\) in \(\RR^{| S |}\). In other words, we treat a distribution \(\mu \in \cP(S)\) as an \(|S|\)-dimensional vector \((\mu^{(i)})_{i=1,\dots,|S|}\) whose non-negative coordinates sum up to one. Similarly, since \(A\) is finite, we identify \(\cP(A)\) with a simplex \(\mathfrak{A}\) in \(\RR^{|A|}\). However, representing admissible actions \(\bar a \in \bar U(\mu) \subseteq \cP(S \times A)\) of the lifted MDP requires careful handling due to the underlying constraints. A first approach is to identify \(\cP(S \times A)\) with a simplex in \(\RR^{|S| \times |A|}\) and to view a lifted action \(\bar a\) as a \(|S|\times |A|\) matrix \(\big( \bar a( x^{(i)}, \alpha^{(j)}) \big)_{1 \leq i \leq |S|, 1 \leq j \leq |A|}\) of non-negative numbers summing up to \(1\). Then a pair \((\mu, \bar a) \in \bar S \times \bar A\) is in \(\bar \Gamma\) if and only if the following linear constraint is satisfied: \[\sum_{j=1}^{|A|} \bar a( x^{(i)}, \alpha^{(j)} ) = \mu^{(i)}, \qquad i = 1, \cdots, |S|.\] While this transformation is straightforward, it is insufficient for our purposes. It only provides a representation of the actions and controls from the central planner’s perspective. It fails to capture the strategy functions of non-randomized stationary mixed Markovian closed-loop policies for an individual agent in our original optimization problem.

For any pair \((\mu, \bar a) \in \bar \Gamma\), we can define a kernel \(k_{\mu}: S \to \cP(A)\) by setting, for \(i=1,\dots,|S|\) and \(j=1,\dots,|A|\), \[k_\mu(x^{(i)})(\alpha^{(j)}) = \frac{\bar a( x^{(i)}, \alpha^{(j)} )}{\mu(x^{(i)})}, \qquadif\mu(x^{(i)})>0,\] and by choosing an arbitrary probability measure on \(A\) when \(\mu(x^{(i)})=0\).

Note that common randomization is absent here. As proved above (see Theorem 2.1), there exists a non-randomized stationary policy for the lifted MDP. So the central planner can look for strategy functions within the set: \[\label{eq:def-calA-S-PA} \tilde{\cA} := \{ \tilde{a}: S \to \cP(A) \; | \; \tilde{a} \text{ Borel measurable} \}.\tag{77}\] So, because of the constraint, we can identify a probability measure \(\bar a\) on the product space \(S\times A\) whose first marginal is given and equal to \(\mu\), with the probability kernel \(\tilde{a}\) appearing in its disintegration with respect to its first projection. In this way, we can view the Q-function \(\bar Q^*\) as a function of \(\mu\) and \(\tilde{a}\) instead of \(\mu\) and \(\bar a\), avoiding the redundancy imposed by the constraint. So if we introduce the function \(\tilde{Q}^*: \cP(S) \times \tilde{\cA} \to \RR\) defined by: \[\tilde{Q}^*(\mu, \tilde{a}) := \bar Q^*(\mu, \mu \measprod \tilde{a}),\] the Bellman equation ?? becomes: \[\label{eq:optimal95Q95kernel95version} \tilde{Q}^*(\mu, \tilde{a}) = \int_{S \times A} f(x, \alpha, \mu \measprod \tilde{a} ) \tilde{a}(x, d\alpha) \mu( dx)+ \gamma \EE \left[ \inf_{\tilde{a}'\in \tilde{\cA}} \tilde{Q}^*( \mu_1, \tilde{a}')\right],\tag{78}\] for \((\mu, \tilde{a}) \in \cP(S) \times \tilde{\cA}\), where \(\mu_1 = \bar F( \mu, \mu \measprod \tilde{a}, \varepsilon^0)\). Here \(\bar F\) denotes the lifted system function introduced in equation ?? of Chapter 2.2, representing how the population distribution evolves given the current state and the chosen level-1 action. An important observation from Chapter 2.2 is that the discrete-time evolution remains stochastic in the presence of common noise \(\varepsilon^0\), as shown in Equation ?? . Although \(S\) and \(A\) are finite, equation 78 is still characterized as a fixed point in the space of bounded lower semi-continuous functions on a closed subset of a finite-dimensional Euclidean space, consistent with the measurability considerations addressed in deriving equation ?? . We also introduce the function \(\tilde{f}: \cP(S) \times \tilde{\cA} \to \RR\) defined by: \[\tilde{f}(\mu, \tilde{a}) := \bar f(\mu, \mu \measprod \tilde{a}) = \int_{S \times A} f(x, \alpha, \mu \measprod \tilde{a}) \tilde{a}(x, d\alpha) \mu( dx), \qquad (\mu, \tilde{a}) \in \cP(S) \times \tilde{\cA}.\] The function \(\bar f\) relates back to the lifted cost function defined in Equation ?? of Chapter 2.2, which aggregates individual agent costs into a single scalar value for the central planner. In the remainder of this section, we propose two model-free algorithms based on the optimal state-action value function \(\bar Q^* : \bar \Gamma \to \RR\) or, equivalently, \(\tilde{Q}^*: \cP(S) \times \tilde{\cA} \to \RR\).

3.1.3.2 Simplex discretization and tabular MFQ-learning

 

We consider two settings, depending on whether the controls at level-0 are mixed or pure. In both cases, we discretize the simplexes, which leads to a finite-state finite-action MDP for which the tabular Q-learning algorithm can be used. We establish its convergence. When using pure controls, we can prove the convergence of the value function as well as that of the associated optimal control.

3.1.3.2.1 Q-learning with level-0 mixed controls

We first consider the situation where the controls are mixed at the level-0, meaning that the generic agent can use randomized actions. Since the simplexes \(\mathfrak{S}\) and \(\mathfrak{A}\) are infinite, a tabular version of the Q-learning algorithm cannot be directly applied to approximate \(\tilde{Q}^*\). To address this, we first approximate these simplexes using finite subsets \(\check{\mathfrak{S}} \subset \mathfrak{S}\) and \(\check{\mathfrak{A}} \subset \mathfrak{A}\). Let \(\check{\cA} = \{ \check{a} : S \to \check{\mathfrak{A}} \}\). In particular, \(|\check{\cA}| = | \check{\mathfrak{A}} |^{|S|}\) because we identify functions in \(\check{\cA}\) with \(|S|\)-dimensional vectors whose entries take values in the finite set \(\check{\mathfrak{A}}\). To ensure that the mean-field term takes values in the finite set \(\check{\mathfrak{S}}\), we use a projection: At time \(n\), given \(\mu_n \in \check{\mathfrak{S}}\), we compute \(\mu_{n+1} = \bar F( \mu_n, \mu_n \measprod \check{a}, \varepsilon_{n+1}^0)\), and then we project \(\mu_{n+1}\) back on \(\check{\mathfrak{S}}\) using a projection operator \(\proj_{\check{\mathfrak{S}}}: \cP(S) \to \check{\mathfrak{S}}\). Precise definitions of the discretization and the projection are provided below, after introducing a discrete version of the original MFC problem.

More precisely, we consider the : \[\inf_{ \check \pi \in \check \Pi} \check J^{\check \pi}(\mu_0),\qquad \mu_0 \in \check{\mathfrak{S}},\] where \(\check \Pi = \{ \check \pi: S \times \check{\mathfrak{S}} \to \check{\mathfrak{A}} \}\) is a (finite) set of policies, and for every strategy function \(\check \pi : S \times \check{\mathfrak{S}} \to \check{\mathfrak{A}}\), \(\check J^{\check \pi}: \check{\mathfrak{S}} \to \RR\) is defined by: \[\label{eq:generic-MFC-fctmeasure-cost-proj} \check J^{\check \pi}(\mu_0) = \EE \left[ \sum_{n \geq 0} \gamma^n \tilde{f} \Big( \mu^{\mu_0, \check \pi}_n, \check \pi(\cdot, \mu^{\mu_0, \check \pi}_n) \Big) \right]\tag{79}\] where \[\label{eq:generic-MFC-fctmeasure-dyn-proj} \mu^{\mu_0, \check \pi}_{n+1} = \proj_{\check{\mathfrak{S}}} \circ \bar F \Big( \mu_n^{\mu_0, \check \pi}, \mu_n^{\mu_0, \check \pi} \measprod \check \pi( \cdot, \mu_n^{\mu_0, \check \pi}), \varepsilon_{n+1}^0 \Big) =: \check \Phi^{\check \pi, \varepsilon_{n+1}^0}(\mu_n^{\mu_0, \check \pi}).\tag{80}\] We will denote by \(\check J^*\) and \(\check Q^*\) respectively the optimal state and state-action value functions of this projected MFC problem. Here \(\check Q: \check{\mathfrak{S}} \times \check{\cA} \to \RR\) can be represented by a matrix (also called a table) in \(\RR^{ |\check{\mathfrak{S}}| \times | \check{\cA}| }\) and is viewed as an approximation of \(\tilde{Q}^* : \cP(S) \times \tilde{\cA} \to \RR\) of the original MFC problem.

This problem can be viewed as an MDP with finite state and action spaces. In this case, one can adapt in a straightforward way the classical tabular Q-learning algorithm. The pseudo-code is provided in Algorithm 1. Note that, even in the absence of common noise, this algorithm is possibly stochastic since at each episode, the order in which the state-action pairs are picked is potentially random. In practice, the order could be fixed in advance or stem from a sampled trajectory, which we use in the numerical examples provided in Section 3.1.4.

Algorithm 1 returns the table \(\check Q_{N_{\mathrm{epi}}}\) after \(N_{\mathrm{epi}}\) episodes. We prove below that this table converges to the optimal Q-function \(\tilde{Q}^*\) in a suitable sense.

To keep the chapter at a reasonable length, we will make the following simplifying assumptions.

We endow the simplexes \(\mathfrak{S}\) and \(\mathfrak{A}\) respectively with the Euclidean distances \(d_{\mathfrak{S}}\) and \(d_{\mathfrak{A}}\) of the spaces \(\RR^{|S|}\) and \(\RR^{|A|}\). Because \(S\) is finite, we can identify \(\tilde{\cA}\) defined in 77 with \(\cP(A)^{|S|}\) with the distance \(d_{\tilde{\cA}}(\tilde{a}, \tilde{a}') = \sup_{x \in S} d_{\mathfrak{A}}( \tilde{a}(x), \tilde{a}'(x))\) for \(\tilde{a}, \tilde{a}' \in \tilde{\cA}\). Furthermore, we consider the following discretizations of the simplexes. Let \(\varepsilon_{\mathfrak{S}}>0\) satisfying: for all \(\mu \in \mathfrak{S},\) there exists \(\check \mu \in \check{\mathfrak{S}}\) s.t. \(d_{\mathfrak{S}}(\mu, \check \mu) \le \varepsilon_{\mathfrak{S}}\). Similarly, let \(\varepsilon_{\mathfrak{A}} > 0\) satisfying: for all \(\nu \in \mathfrak{A}\), there exists \(\check \nu \in \check{\mathfrak{A}}\) such that \(d_{\mathfrak{A}}( \nu, \check \nu ) \leq \varepsilon_{\mathfrak{A}}\). Because \(S\) is finite and the definition of the distance \(d_{\tilde{\cA}}\), we have for every \(\tilde{a} \in \tilde{\cA}\), there exists \(\check{a} \in \check{\cA}\), such that \(d_{\tilde{\cA}}( \tilde{a}, \check{a}) \leq \varepsilon_{\mathfrak{A}}\).

Note that the boundedness of the one-stage cost \(\tilde{f}\) from Assumption 1.1 together with the fact that \(\gamma \in (0,1)\) ensures the existence of a finite bound \(\check J_{bound}\) for the state value function of the projected MFC problem. The regularity of \(\bar J^*\) in 1.2 can typically be ensured through suitable conditions on the data of the problem, as e.g. in [119]–[121]. Assumption 1.3 is similar to the covering time assumption in [122]. One way to satisfy this assumption is to fix an ordering of the state-action pairs and loop through it repeatedly. In this case, \(T_{cov} = |\check{\mathfrak{S}} \times \check{\cA}|\). In practice, it is common to take \(\check{a}_{n}\) using a randomized policy, such as an \(\epsilon\)-greedy policy, which amounts to taking action \(\argmin \check{Q}_n(\check{\mu}_n, \cdot)\) with probability \(1-\epsilon\), and otherwise to pick an action uniformly at random in \(\check{\cA}\). In standard Q-learning, this is interpreted as a randomized policy. In our mean-field context, this is a form of common randomization. The main advantage of this approach is that, as the Q-function is learned, the sequence of actions is mostly concentrated on the optimal choices instead of spending many iterations on irrelevant actions. However, when the Q-function is not yet well estimated, using an \(\epsilon\)-greedy policy with a small \(\epsilon\) can lead to poor decisions. This is an illustration of the exploration and exploitation trade-off. The proof of convergence of the Q-learning algorithm can be extended to stochastic state-action sequences, provided the corresponding covering time is bounded by a constant with high probability. For instance, it is enough to assume that the covering time is smaller than \(L\) with probability \(1/2\). In practice, exploration of the state-action space is often guaranteed by a combination of action randomization and exploring starts at different episodes (provided the learner can query an oracle which simulates transitions from any \((\mu, \tilde{a})\)).

We consider projection operators \(\proj_{\check{\mathfrak{S}}}: \mathfrak{S} \to \check{\mathfrak{S}}\) and \(\proj_{\check{\cA}}: \tilde{\cA} \to \check{\cA}\) such that \[(\proj_{\check{\mathfrak{S}}}(\mu), \proj_{\check{\cA}}(\tilde{a}) ) := (\check \mu, \check{a})\] for every \((\mu, \tilde{a}) \in \mathfrak{S} \times \tilde{\cA}\) where \((\check \mu, \check{a})\) is the closest point (or one of the closest points, in case of ties) in \(\check{\mathfrak{S}} \times \check{\cA}\) with respect to \(d_{\mathfrak{S}}\) and \(d_{\tilde{\cA}}\). Based on simplex discretizations, this point satisfies \(\| \mu - \check \mu\|_{d_{\mathfrak{S}}} \leq \varepsilon_{\mathfrak{S}}\) and \(d_{\tilde{\cA}}( \tilde{a} , \check{a} ) \leq \varepsilon_{\mathfrak{A}}\). We set \(\beta = (1-\gamma)/2\), and for \(\delta \in (0,1)\), we let \(T_{cov}(\delta) = \lceil T_{cov} \log_2(1/(2\delta))\rceil\).

As an example, if we cover the simplexes with balls of radius \(\varepsilon_{\mathfrak{S}}\) and \(\varepsilon_{\mathfrak{A}}\) respectively, the number of points in the covering would be of order \(\Theta(\varepsilon_{\mathfrak{S}}^{1-|S|})\) and \(\Theta(\varepsilon_{\mathfrak{A}}^{1-|A|})\) respectively. Then in the above result the number of episodes is of order: \[\Omega\left( \left(\frac{(1-|S|)\ln \left(\varepsilon_{\mathfrak{S}} \right) + (1-|A|)|S|\ln \left(\varepsilon_{\mathfrak{A}}\right) + \ln \left( 1 / \varepsilon \right)}{ \varepsilon^2}\right)^{\frac{1}{\kappa}} \hskip -6pt + \left(\ln \left( \frac{1}{\varepsilon} \right)\right)^{\frac{1}{1-\kappa}} \right),\] while the error bound \(\varepsilon'\) is given by: \(\varepsilon' = \varepsilon + C (\varepsilon_{\mathfrak{S}} + \varepsilon_{\mathfrak{A}}).\)

Note that \(\varepsilon\) can be chosen as small as desired provided \(N_{\mathrm{epi}}\) is large enough. The second and third terms in the error \(\varepsilon'\) are proportional to \(\varepsilon_{\mathfrak{S}}\) and \(\varepsilon_{\mathfrak{A}}\), a result that is generally unavoidable due to the projection onto the finite sets \(\check{\mathfrak{S}}\) and \(\check{\mathfrak{A}}\). However, this error vanishes as \(\varepsilon_{\mathfrak{S}} \to 0\) and \(\varepsilon_{\mathfrak{A}} \to 0\), i.e., as \(\check{\mathfrak{S}}\) and \(\check{\mathfrak{A}}\) become increasingly accurate approximations of \(\cP(S)\) and \(\cP(A)\), respectively.

We prove this result below. The proof can be summarized in the following three steps: (1) For \(N_{\mathrm{epi}}\) large enough, we have \(\check Q_{N_{\mathrm{epi}}} \approx \check Q^*\) on \(\check{\mathfrak{S}} \times \check{\cA}\); (2) \(\check Q^* \approx \tilde{Q}^*\) on \(\check{\mathfrak{S}} \times \check{\cA}\); (3) For every \((\mu, \tilde{a}) \in \mathfrak{S} \times \tilde{\cA}\), \(\tilde{Q}^*( \proj_{\check{\mathfrak{S}}}(\mu), \proj_{\check{\cA} }(\tilde{a} ) ) \approx \tilde{Q}^*(\mu, \tilde{a})\). The first step leverages standard Q-learning convergence results [122], while the other two steps arise from the regularity assumptions and the approximation of \((\mathfrak{S}, \tilde{\cA})\) by \((\check{\mathfrak{S}}, \check{\cA})\).

3.1.3.2.2 Q-learning with pure level-0 controls

The approach described above is designed for cases where we seek optimal actions that may be randomized at the individual level. Searching in the space \(\cP(A)\) comes with a computational cost that is reflected in the bounds through the cardinality of the discrete simplex \(\check{\mathfrak{A}}\). However, in some situations it may be interesting to directly search for actions that are pure at the individual level. In this paragraph, we restrict our attention to pure level-0 policies, i.e., without idiosyncratic randomization at the agent level. For such policies, given the recommendation of the central planner (which might depend on a common randomization), each agent chooses their action in a deterministic way, based purely on their own state. This amounts to removing the random variables \(\vartheta_n\) from the probabilistic framework presented in Section 2.2.2.1. For simplicity we will use the same notations for the value functions \(\bar{J}^*\), \(\bar{Q}^*\) and \(\check{Q}\). The goal is now to approximate \(\bar{Q}^*\), which is the Q-function optimized over the set of policies that are pure at the level \(0\). In this case, instead of 77 , the set of strategy functions is (for simplicity we keep the notation \(\tilde{\cA}\)): \[\tilde{\cA} := \{ \tilde{a}: S \to A \} = A^S.\] In Algorithm 1, we replace \(\check{a} \in \check{\cA}\) by \(\tilde{a} \in \tilde{\cA}\).

A similar result can be obtained if \(\bar J^*\) is only Hölder continuous, consistent with Remark [ref:holderJ].

The result of the above theorem provides convergence guarantees for the Q-function. Let us now derive a consequence in terms of the optimizer. To this end, we employ the following additional assumption regarding the gap between the values of the best and second-best actions, a standard requirement in approximation algorithms based on tabular Q-functions [123], [124].

The above assumption is based on the simplex discretizations discussed previously. In particular, the constant \(K_A\) may depend on \(\check{\mathfrak{S}}\) and \(\check{\mathfrak{A}}\). This assumption is satisfied, for example, if \(\tilde{Q}^*( \check \mu, \cdot)\) is strictly convex for every \(\check \mu \in \check{\mathfrak{S}}\).

To recover minimizers or approximate minimizers, it will be convenient to work with the following operators. In general, they are defined on the vector space \(\RR^{\texttt{m}}\). For \(\tau>0\) and \(x = (x_1,\dots,x_{\texttt{m}}) \in \RR^{\texttt{m}}\), we define \(\softmin_\tau : \RR^{\texttt{m}} \to \RR^{\texttt{m}}\) by \[\softmin_\tau(x) = (e^{- \tau x_1}, \dots, e^{- \tau x_{\texttt{m}}}) / \sum_{j} e^{- \tau x_j}.\] For \(x \in \RR^{\texttt{m}}\), we define \(\argmine: \RR^{\texttt{m}} \to [0,1]^{\texttt{m}}\) by \[\argmine(x) = \left( \mathbf{1}_{i \in \argmin(x)} \right)_{i=1}^{\texttt{m}} / |\argmin(x)|.\] where \(\argmin(x) = \{ j \in \{1,\ldots, {\texttt{m}} \} \, : \, x_j = \min\{x_1,\ldots, x_{\texttt{m}} \} \}\). In the sequel, we use these operators with the dimension \({\texttt{m}} = | \tilde{\cA} | = | A |^{|S|}\). For any function \(q: \tilde{\cA} \to \RR\), we identify \(q\) with the vector \((q(\tilde{a}))_{\tilde{a} \in \tilde{\cA}}\).

The proof is provided below. The \(\argmine\) in the second term is here in case there are several optimal controls. The \(\softmin\) operator regularizes the best action predicted by the estimation \(\check Q_{N_{\mathrm{epi}}}\) of the function \(\tilde{Q}^*\).

Notice that, contrary to Theorem 1.1, here \(\cP(A)\) is not involved hence the error bound in Theorem 1.2 does not depend on any discretization parameter \(\varepsilon_{\mathfrak{A}}\). Consequently, the error in Corollary 1.1 can be made arbitrarily small by choosing a sufficiently large \(\tau\) and a sufficiently small \(\varepsilon\).

3.1.3.3 Deep reinforcement learning for MFMDP

 

The tabular Q-learning method of the previous section is sufficiently straightforward to allow for a detailed analysis. However, it cannot be used in practice for large state or action spaces because of the computational cost of discretizing the simplexes. An alternative is to work directly with continuous spaces, in which case the policies and value functions cannot be represented in a tabular way. Instead, we can rely on function approximation. To this end, we now propose methods from deep RL which are more suitable for continuous spaces. The motivations are twofold.

First, if \(S\) and \(A\) are finite but we want to learn an optimal policy that is potentially randomized at level-0, the discretization approach proposed in § 3.1.3.2.1 has a complexity that increases with the number of points in the discretization of \(\cP(A)\), which itself increases exponentially quickly with the cardinality of \(A\). For this reason, it can be interesting to tackle directly \(\cP(A)\) as a continuous action space and to use deep RL methods for continuous action space MDPs.

Second, some MFC problems are naturally posed with a continuous state space \(S\). In this case, under mild conditions, the optimal policy is in fact non-randomized not only at the level-1 but even at the level-0 (see, e.g., [22], [127] and the references therein). However, the state of the MFMDP is an element of the infinite dimensional space \(\cP(S)\). As noted in Chapter 2.2, this measure-valued state allows the central planner to treat the collective behavior of an infinite population as a single Markovian process on the simplex \(\bar S = \cP(S)\). From a numerical viewpoint, we need two ingredients: (1) a finite-dimensional approximation of the mean field and (2) a parameterized approximation of the value function or the policy taking this finite-dimensional representation of the mean field state as an input. For the second point, we will again use deep neural networks. For the first point, for the sake of definiteness, we choose to simply replace \(\cP(S)\) by \(\cP(\check S)\) where \(\check S\) is a discretization of \(S\) with a finite number of points. We assume that, given \(\check\mu \in \cP(\check S)\) and a level-1 action \(\check a: S \to A\), one can get from the environment a sample of the next state and the associated cost \(\tilde{f}(\check\mu, \check a)\). The problem thus reduces to an MDP with finite dimensional (but potentially continuous) state and action spaces. Such MDPs can be solved with a variety of deep RL algorithms. In the following, we provide numerical illustrations based on the proposed in [118]. It relies on two neural networks, one for the Q-function (the critic) and one for the policy (the actor). The core of the algorithm involves alternating updates: the critic is updated by minimizing an empirical square error, while the actor is updated through gradient descent. To improve exploration, a Gaussian noise \(\epsilon^a_{n+1}\) is added to the action prescribed by the actor. Furthermore, for more stability, target networks are also added. We refer to [118] for more details. We write the algorithm below in cost-minimization form: the critic approximates cumulative cost and the actor is updated by descending the estimated cost. The algorithm is summarized in Algorithm 2.

3.1.4 Numerical Examples↩︎

In this section, we illustrate the above RL algorithms on two different examples. In each example, the representative agent’s state space and action space are finite. The first example is a cybersecurity model with 4 states, on which we illustrate both the tabular Q-learning algorithm and the DDPG algorithm. The second example features a larger number of states, leading us to apply only the DDPG algorithm, which offers the advantage of not requiring simplex discretization.

The code is available at: https://github.com/mlauriere/mfrl-tutorial-code.

3.1.4.1 Example 1: Cybersecurity model

 

For a first testbed, we start with a finite state problem. We revisit the cybersecurity example introduced in [74], now adopting the perspective of a central planner (such as a large company or a state) aiming to protect its system against hacker attacks. The situation can be phrased as an MFC problem.

In this model, the population consists of a large group of computers which can be either defended (D) or undefended (U), and either infected (I) or susceptible (S) of infection. Hence the set \(S\) has four elements corresponding to the four possible combinations: DI, DS, UI, US. The action set is \(A = \{0,1\}\), where \(0\) is interpreted as the fact that the central planner is satisfied with the current protection level (D or U) of the computer under consideration, whereas \(1\) indicates a desire to modify this level. In the latter case, the update occurs at a (fixed) rate \(\lambda >0\). If the controls are pure at level-0, the central planner assigns a single action to each of the four states, which is then applied to all computers at that state. If the controls are mixed at level-0, the planner chooses a distribution over actions for each state, and each computer in that state then independently selects an action according to the chosen distribution. When infected, each computer may recover at rate \(q_{rec}^D\) or \(q_{rec}^U\) depending on whether it is defended or not. On the other hand, a computer may be infected either directly by a hacker, at rate \(v_H q_{inf}^D\) (resp. \(v_H q_{inf}^U\)) if it is defended (resp. undefended), or by undefended infected computers, at rate \(\beta_{UU}\mu(\{UI\})\) (resp. \(\beta_{UD}\mu(\{UI\})\)) if it is undefended (resp. defended), or by defended infected computers, at rate \(\beta_{DU}\mu(\{DI\})\) (resp. \(\beta_{DD}\mu(\{DI\})\)) if it is undefended (resp. defended). Here \(v_H\) can be interpreted as the attack intensity parameter.

In short, the infinitesimal generator is given by: \[\label{eq:cybersecurity-MFC-generator} Q^{\mu,a} = \begin{pmatrix} \dots & Q^{\mu,a}_{DS \rightarrow DI} & \lambda a & 0 \\ q_{rec}^D & \dots & 0 & \lambda a \\ \lambda a & 0 & \dots & Q^{\mu,a}_{US \rightarrow UI} \\ 0 & \lambda a & q_{rec}^U & \dots \end{pmatrix}\tag{81}\] where \[\begin{align} & Q^{\mu,a}_{DS \rightarrow DI} = v_H q_{inf}^D + \beta_{DD} \mu(\{DI\}) + \beta_{UD} \mu(\{UI\}) , \\ & Q^{\mu,a}_{US \rightarrow UI} = v_H q_{inf}^U + \beta_{UU} \mu(\{UI\}) + \beta_{DU} \mu(\{DI\}), \end{align}\] and all instances of \(\dots\) are chosen so that each row sums to zero, i.e., each diagonal entry is the negative of the sum of the off-diagonal rates in its row. At each time step, the central planner pays a protection cost \(k_D>0\) for each defended computer, and a penalty \(k_I>0\) for each infected computer. The instantaneous cost in the MFMDP is thus defined as: \[\bar f(\mu, \bar a) = k_D \mu(\{DI, DS\}) + k_I \mu(\{DI, UI\}), \qquad (\mu, \bar a) \in \bar S \times \bar A.\] The optimal control and optimal flow of distributions can be characterized by a forward-backward ODE system which can be obtained in a way similar to what is done in the MFG setting e.g. in  [22]. We will use this solution as a benchmark.

Tabular Q-learning. For the sake of illustration, we present results obtained by tabular Q-learning with simplex discretization as described in § [se:tabularQ-discrete]. The state space for the population distribution is \(\bar S\), which is identified with the simplex \(\mathfrak{S} = \{(\mu^{(i)})_{i=1,\dots,4} \in [0,1]^4\,:\, \sum_{i} \mu^{(i)} = 1\}\). In this example, we consider that the actions are pure not only at level 1 (the central planner picks a single control from \(\bar A\)) but also at level 0 (individual agents do not randomize). This means that for each state \(i \in \{1,\dots,4\}\), the action \(a^{(i)}\) belongs to \(\{0, 1\}\). Thus, the set of actions \(\bar A\) for the population is the finite set \(\{0,1\}^4\) of cardinality 16. We replace \(\mathfrak{S}\) by the finite set: \[\label{eq:grid-mu-simplex} \check{\mathfrak{S}} = \Big\{ (\mu^{(i)})_{i=1,\dots,4} \in \{0,1/N_m,\dots,1-1/N_m, 1\}^4\,:\, \sum_{i} \mu^{(i)} = 1 \Big\},\tag{82}\] where \(N_m+1\) is the number of points in the discretization of each dimension. We consider three different training schemes for the Q-function:

1. SyncFull: A synchronous approach where the Q-function is updated for all state-action pairs \((\check{\mu}, \check{a}) \in \check{\mathfrak{S}} \times \check{\cA}\) at each iteration using the exact expectation of the next state’s value. This serves as a benchmark for the discrete problem.

2. SyncGreedy: A synchronous approach where at each iteration, for every state \(\check{\mu}\), only the Q-value for the current \(\epsilon\)-greedy action is updated based on a single sample transition.

3. AsyncTraj: An asynchronous approach where the Q-values are updated along trajectories sampled from the environment, using an \(\epsilon\)-greedy policy to explore the state space. This corresponds to the standard Q-learning setting.

For these simulations, we used the parameters summarized in Table 2. The time step is set to \(\Delta t = 0.2\) and the horizon to \(T=10.0\) (50 steps per episode). The continuous discount rate is set by the convention \(\rho=-\log(0.5)/0.1\), so the per-step discount used in these tabular simulations is \(\gamma=\exp(-\rho\Delta t)=0.25\). The discretization mesh size is \(N_m=50\) for the case without common noise and \(N_m=30\) for the case with common noise.

We first present results without common noise. In this case, the MFMDP is deterministic. Figure 3 shows the evaluation trajectories for the three solvers, which all successfully learn to drive the system toward the same stationary configuration although they do not match perfectly the reference ones because of the discretization of the state space. Figure 4 displays the training progress in terms of costs.

Next, we consider the case with common noise, where the attack intensity \(v_H\) oscillates. The results are shown in Figures 5, 6, and 7. Despite the stochasticity, the solvers manage to learn robust policies that stabilize the system. Figure 7 specifically shows slices of the learned value function on the simplex, comparing it with the benchmark solution.

DDPG. The solution can also be learned without discretizing the simplex \(\mathfrak{S}\), by directly considering the value function on the continuous space. We use the DDPG algorithm [118] to learn a deterministic policy \(\pi: \mathfrak{S} \to \bar A\). For this continuous solver, we consider the action space to be the continuous hypercube \(\bar A = [0,1]^4\), where \(a^{(i)} \in [0,1]\) represents the defense intensity in state \(i\). Notice that this is a relaxed action space compared to the tabular setting, where \(\bar A = \{0,1\}^4\).

The actor and critic networks both consist of two hidden layers of 32 neurons each, with sigmoid activations. We trained the networks for 500 episodes with a learning rate of \(10^{-3}\) for the actor and \(2 \times 10^{-3}\) for the critic. The soft update parameter for the target networks was set to \(\rho_{\mathrm{tar}}=0.005\). Here the DDPG implementation uses \(\Delta t=0.1\), so the same continuous discount convention \(\rho=-\log(0.5)/0.1\approx 6.93\) gives the DDPG per-step discount \(\gamma=\exp(-\rho\,0.1)=0.5\).

The results for both deterministic and stochastic environments are displayed in Figure 8 for the evaluation trajectories, Figure 9 for the value function simplex slices, and Figure 10 for the training progress in terms of costs. DDPG leads to trajectories that are closer to the reference ones than with tabular Q-learning because the simplex has not been discretized. We also observe that in the presence of common noise, the trajectories remain close to the reference ones, although there are small fluctuations.

3.1.4.2 Example 2: Discrete distribution planning

 

We now apply the MFC framework to a problem where the objective is to match a target distribution. We consider a model with parameters summarized in Table 3, where the agents’ actions correspond to moving left, staying, or moving right. At the boundaries, attempted moves outside the grid are folded back to the boundary. The dynamics follow \(F(x,a,\mu,e,e^0)=\Pi_X(\Pi_X(x+a)+e^0)\). Here \(\Pi_X\) denotes projection onto the finite grid \(X=\{1,\dots,11\}\); the outer projection applies the same boundary convention after the common-noise shift. A representative agent incurs a one-step cost of \(1\) for any movement (left or right), plus the distribution penalty \(2\|\mu-\mu_{\texttt{target}}\|_2\). Thus, \(f(x,a,\mu)=\mathbf{1}_{\{a\neq0\}}+2\|\mu-\mu_{\texttt{target}}\|_2\). For this example, the precise target distribution \(\mu_{\texttt{target}}\) is given in Table 3. This model does not include idiosyncratic noise.

Importantly, except for specific pairs of initial and target distributions, the population generally cannot match the target distribution unless agents are permitted to randomize their actions at the individual level. We therefore use \(\cP(A)^X\) as the level-1 action space. For a randomized action kernel \(\nu\in\cP(A)^X\), the movement-cost term is \(\sum_{x\in X}\mu(x)\sum_{a\in A}\mathbf{1}_{\{a\neq0\}}\nu_x(a)\). The resulting continuous action space justifies the application of the DDPG method.

We present results for the DDPG algorithm both with and without common noise. In the stochastic case, the common noise follows the distribution specified in Table 3. The networks are trained for 3000 episodes using the hyperparameters detailed in the same table.

Figures 11 and 12 display the evaluation results for four different initial distributions, respectively without and with common noise. Each row corresponds to a single test case and contains five subplots: (1) Initial state distribution vs target, (2) Final/average state distribution vs target, (3) Initial action distribution across states, (4) Terminal action distribution, and (5) The realized trajectory of the common noise (which is zero in the non-stochastic case). Without common noise, the trajectories are deterministic and the final state distribution is very close to the target distribution. With common noise, the trajectories are stochastic and the final state distribution is slightly different from the target distribution. But we observe that in both scenarios, the DDPG agent successfully learns to move the population towards the target distribution. At the initial time, the control predominantly favors movement towards the center of the target distribution, while at the terminal time, the "stay" action becomes dominant as the target has been reached. Figure 13 shows the training progress in terms of the social cost for both cases.

3.1.5 Notes and Complements↩︎

RL has been applied to solve MFC and MFMDPs using Q-learning as presented in this chapter in [5], as well as in [15], [16]; see also [128]. The deep RL approach presented in this chapter was first introduced in [5]. Extensions include applications to mean-field type games, namely, Nash equilibria between a finite number of central planners solving MFC problems: [129] studied convergence of Nash Q-learning and application of DDPG to such problems, and an adaptation of the PPO algorithm to MFTGs has been proposed in [130].

Other algorithms include a two-timescale method introduced in [131] and whose convergence has been established in [132]. This approach has then been extended to study actor-critic methods for MFC in [132]. Alternative actor-critic methods have been proposed in [83], [133]. [134] proposed a learning algorithm based on linear function approximation. [135] proposed a method based on an adaptation of the classical policy gradient method. In the line of continuous-time RL developed by [136], [137] studied q-learning for MFC in continuous time. For model-free continuous-time mean-field control from discrete-time transition data, see also [138]. RL methods for partially observable MFC models has been proposed in [30], [139]. Related deep-learning methods for MFC with common noise, based on population-dependent controls and neural approximations of the mean field, have been studied in [140].

The model for cybersecurity considered in this chapter was introduced, in the context of mean field games, by [74]; related MFG and MFRL references include [78]–[80], [82], [126], [141]. Other examples of applications of RL for MFC include economics and finance [142], robotic swarm control [143], and UAVs [144], [145].

We conclude by mentioning that while this chapter focuses on a model-free perspective, an alternative approach is to rely on model-based methods such as in [37], [146].

3.2 Implementation of the LQ Mean Field Model↩︎

Following the general numerical methods discussed in Chapter 3.1, this chapter focuses on the Linear-Quadratic Mean Field Control (LQMFC) setting. As introduced in Chapter 2.3, the family of Linear-Quadratic (LQ) models is central to optimal control due to its mathematical tractability and wide-ranging applicability. In reinforcement learning (RL), these models serve as an important benchmark for analyzing the convergence of iterative algorithms [108], [116].

Mean-field type control, or McKean-Vlasov control, extends the LQ framework by incorporating state and action distributions. While Chapter 3.1 explored tabular and deep RL for general MFMDPs, the specific structure of LQ models allows us to establish global convergence guarantees for Policy Gradient (PG) methods, even in model-free settings. Our analysis generalizes the results of Fazel et al. [108] to the mean-field case, proving that the optimal parameters of the feedback control can be effectively learned.

A key feature of our model is the presence of common noise, which captures random shocks affecting all agents simultaneously (see Chapter 2.3). While it complicates the analysis, it also makes the limiting population distribution stochastic. In our model-free investigations, we study two types of simulators: an idealized McKean–Vlasov (MKV) simulator in § 3.2.1.4 and a more realistic finite-population simulator in § 3.2.1.5. Proving convergence of the latter relies on the propagation of chaos. Numerical experiments illustrate these methods in § 3.2.5.

3.2.1 Policy gradient methods↩︎

Besides the notation introduced for the LQMFC problem in Section 2.1.6 and Chapter 2.3, it will sometimes be convenient to use the following notation. We will use \(u\) for the action. The one-step cost function will be denoted by \(c: S \times S \times A \times A \to \RR\), with model parameters \((Q, \bar Q, R, \bar R) \in \RR^{d \times d} \times \RR^{d \times d} \times \RR^{\ell \times \ell} \times \RR^{\ell \times \ell}\) given by \[\label{fo:lq95one95step95cost} c(x, \bar{x}, \ctrl, \bar{\ctrl}) = (x-\bar{x})^\top Q (x-\bar{x}) + \bar{x}^\top (Q + \bar{Q}) \bar{x} + (\ctrl-\bar{\ctrl})^\top R (\ctrl-\bar{\ctrl}) + \bar{\ctrl}^\top (R + \bar{R}) \bar{\ctrl}.\tag{83}\] Thus, using the notation \(f\) introduced in 13 , we have \(f(x, u, \mu) = c(x, \bar\mu^x, u, \bar\mu^a)\), where \(\bar\mu^x\) and \(\bar\mu^a\) denote the means of \(\mu\) for states and actions, respectively. We also recall that the system function \(F\) is defined in 12 .

In this section, we search for the optimal control process using PG algorithms. Theorem 3.1 states that the optimal control process has feedback coefficients \(\theta^* = (K^*, L^*)\). Inspired by formula ?? , we limit our search for the optimal control process to the set of controls \(\cU_{ad}^\Theta\) that can be parameterized by an element in \(\Theta\). We present model-based and model-free PG algorithms to find the best admissible control parameter \(\theta^* \in \Theta\).

3.2.1.1 Reformulation of the MFC problem

We reformulate the MF problem 49 with the set of admissible control parameters \(\Theta\). For a parameterized admissible control process \(\bctrl^\theta \in \cU_{ad}^\Theta\), we denote the discounted MF cost 44 by \[\label{eq:def-cost-C-theta} C(\theta) = J(\bctrl^\theta).\tag{84}\] Lemma 3.2 states that for a control parameter \(\theta \in \Theta\), the corresponding parameterized control process \(\bctrl^\theta\) belongs to \(\cU_{ad}\) (and thus \(\bctrl^\theta \in \cU_{ad}^\Theta\)), and that the controlled state process \(\bX^{\theta} \in \cX\). For this reformulation, we use Assumption 3.3, which ensures that the optimal feedback coefficients belong to the admissible parameter class. Proposition [pr:reformulation95of95MF95problem] below justifies reformulating the search for the optimal control using parameters in \(\Theta\).

Henceforth, we focus on searching for a minimizer of \(C\) in the set \(\Theta\) to find the optimal control process for the MFC problem. We will sometimes adopt a slight abuse of terminology and refer to \(\theta \in \Theta\) as a “control" or say that the state process is”controlled” by \(\theta\).

In the following, we cast the problem of minimizing the LQMFC cost \(C(\theta)\) into two standard LQ problems depending only on the parameters \(K\) and \(L\), respectively. For a state process \((X_{n}^{\theta})_{n\geq 0}\) controlled by \(\theta = (K, L) \in \Theta\) following dynamics 43 and starting from the initial position \(X_0^\theta = \varepsilon_0^0 + \varepsilon_0\), we introduce two auxiliary processes, \(\bY = (Y_{n})_{n\geq 0}\) and \(\bZ = (Z_{n})_{n\geq 0}\), defined by: \[\label{eq:def-Y-Z-fct-X} Y_{n} := X^{\theta}_{n} - \bar{X}^{\theta}_{n}, \qquad\text{and}\qquad Z_{n} := \bar{X}^{\theta}_{n},\tag{86}\] for \(n \geq 0\), where the dynamics of processes \(\bY\) and \(\bZ\) satisfy \[\begin{align} Y_{n+1} & = ({\mathrm A} - {\mathrm B} K) Y_{n} + \varepsilon_{n+1}, & & Y_{0} = \varepsilon_0 - \EE[\varepsilon_0]; \tag{87} \\ Z_{n+1} & = ( \tilde{\mathrm A} - \tilde{\mathrm B} L) Z_{n} + \varepsilon^0_{n+1}, & & Z_{0} = \varepsilon^0_0 + \EE[\varepsilon_0]. \tag{88} \end{align}\] The processes \(\bY\) and \(\bZ\) depend on the parameters \(K\) and \(L\), respectively. The noise processes and their initial states involved in 87 and 88 are independent. As a result, the processes \(\bY\) and \(\bZ\) are independent. We introduce the corresponding cost functions \(C_y, C_z: \RR^{\ell \times d} \to \RR\) as: \[\label{eq:def95Cy95Cz} C_y(K) = \mathbb{E} \big[ \sum_{n \geq 0} \gamma^{n} (Y_{n})^\top (Q + K^\top R K) Y_{n} \big] , \quad C_z(L) = \mathbb{E} \big[ \sum_{n \geq 0} \gamma^{n} (Z_{n})^\top (\tilde{Q} + L^\top \tilde{R} L) Z_{n} \big].\tag{89}\] We may sometimes use \(\bY^K = (Y_{n}^K)_{n \geq 0}\) and \(\bZ^L = (Z_{n}^L)_{n \geq 0}\) to emphasize the dependence on \(K\) and \(L\) for the processes \(\bY\) and \(\bZ\).

Let \(\Theta_K = \{ K \, | \, \gamma \| {\mathrm A} - {\mathrm B} K \|^2 < 1 \}\) and \(\Theta_L = \{ L \, | \, \gamma \| \tilde{\mathrm A} - \tilde{\mathrm B} L \|^2 < 1 \}\). The following lemma shows the benefit of the re-parametrization 86 inspired by Theorem 3.1.

3.2.1.2 Additional notations

For an admissible control parameter \(\theta = (K, L) \in \Theta\), let \(\bY\) and \(\bZ\) be the two auxiliary processes following dynamics 87 –88 . The symbols \(\Sigma_K^y\) and \(\Sigma_L^z\) stand for the infinite discounted sums of variance/covariance matrices given by \[\label{eq:variance95matrices95y95and95z} \Sigma_K^y := \EE \Big[ \sum \nolimits_{n \geq 0} \gamma^{n} Y_{n} (Y_{n})^\top \Big], \qquad \Sigma_L^z := \EE \Big[ \sum \nolimits_{n \geq 0} \gamma^{n} Z_{n} (Z_{n})^\top \Big].\tag{90}\] By a slight abuse of terminology, we call the matrices \(\Sigma_K^y\) and \(\Sigma_L^z\) the variance matrices of the auxiliary processes \(\bY\) and \(\bZ\). To alleviate notation, we may omit the superscript in \(\Sigma_K^y\) and \(\Sigma_L^z\) and use \(\Sigma_K, \Sigma_L\) when the context is clear. With the expression of \(C_y(K)\) and \(C_z(L)\) in 89 , we have \[\label{eq:expression95cost95sigma} C_y(K) = \langle Q + K^\top R K, \, \Sigma_K \rangle_{tr}, \qquad C_z(L) = \langle \tilde{Q} + L^\top \tilde{R} L ,\, \Sigma_L \rangle_{tr}.\tag{91}\] We denote the exact gradients of the cost \(C(\theta)\) with respect to \(K\) and \(L\) by \(\nabla_K C(\theta)\) and \(\nabla_L C(\theta)\), respectively. Lemma 2.1 implies immediately that \[\nabla_K C(\theta) = \nabla_K C_y(K), \quad \nabla_L C(\theta) = \nabla_L C_z(L).\] For the quantitative PG estimates below, we use the following stronger state-cost condition than Assumption 3.2.

We also make the following non-degeneracy assumption.

Under the non-degenerate conditions ?? in Assumption 2.2, we have \[\label{eq:min95eigenvalue95of95variance95matrices} \lambda_y^1 := \lambda_{min} \big( \Sigma_{Y_{0}} + \frac{\gamma}{1-\gamma} \Sigma^1 \big) > 0 , \quad \lambda_z^0 := \lambda_{min}\big( \Sigma_{Z_{0}} + \frac{\gamma}{1-\gamma} \Sigma^0 \big) > 0.\tag{92}\] The strict positivity of the above two quantities \(\lambda_y^1\) and \(\lambda_z^0\) is crucial for establishing a uniform upper bound for the norms of the gradients \(\nabla_K C_y(K)\) and \(\nabla_L C_z(L)\) within the PG algorithms, and thus they are essential for the global convergence of the PG algorithms.

Since we assume sub-Gaussian distributions for the initial perturbations in the LQ problem, the random variables \(Y_{0}\) and \(Z_{0}\) are also sub-Gaussian. We let \(C_{init, noise}\) denote an upper bound on the sub-Gaussian norms of \(Y_{0}\), \(Z_{0}\), and the i.i.d. noise terms \((\varepsilon_{n+1})_{n \geq 0}\) and \((\varepsilon_{n+1}^0)_{n \geq 0}\). That is, for every \(n \geq 0\): \[\max\big\{ \| Y_{0} \|_{\psi_2}, \| Z_{0} \|_{\psi_2}, \| \varepsilon_{n+1} \|_{\psi_2}, \| \varepsilon_{n+1}^0 \|_{\psi_2} \big\} \leq C_{init, noise}.\]

3.2.1.3 Exact PG for MFC

Assuming full knowledge of the model, we compute the optimal parameters using a standard gradient descent algorithm. With a fixed learning rate \(\eta>0\) and initial parameter \(\theta^{(0)} = (K_0,L_0)\), we update the parameter from \(\theta^{(k)}\) to \(\theta^{(k+1)} = (K^{(k+1)}, L^{(k+1)})\) via : \[\label{eq:one-step-update-exact-pg} \left\{ \begin{array}{rcl} K^{(k+1)} & = & K^{(k)} - \eta \nabla_K C_y(K^{(k)}), \\ L^{(k+1)} & = & L^{(k)} - \eta \nabla_L C_z(L^{(k)}). \end{array} \right.\tag{93}\] To show the convergence of \(C(\theta^{(k)})\) towards the optimal cost \(C(\theta^*)\) with the scheme 93 , we adapt the gradient descent methods with the Polyak-Lojasiewicz condition [147]–[149] for the nonconvex functions \(C_y\) and \(C_z\). Below, we first state the Polyak-Lojasiewicz condition for the auxiliary costs \(C_y\) and \(C_z\) in Proposition [pr:PL95cond95Cy95Cz], followed by a local smoothness condition in the direction of the gradients in Proposition [pr:local95smoothness95Cy95Cz]. The proofs of these results can be found in Sections 3.2.2.7 and 3.2.2.8, respectively.

We recall that a differentiable function \(f: D \to \RR\) on an open set \(D \subseteq \RR^m\) is said to satisfy the \(\nu\)-Polyak-Lojasiewicz (\(\nu\)-PL for short) condition on \(D\) for a constant \(\nu > 0\) if \[\label{eq:PL95condition95definition} \nu \big( f(x) - f(x^*) \big) \leq \| \nabla_x f(x) \|^2, \qquad \forall \, x \in D\tag{94}\] where \(x^* \in D\) is a global minimum in \(D\), i.e., \(f(x) \geq f(x^*)\) for all \(x \in D\).

Now, we consider four technical constants \(h_{var, y}(C_0)\), \(h_{hess, y}(C_0)\), \(h_{var, z}(C_0)\), and \(h_{hess, z}(C_0)\) defined in Lemma 2.13 in Section 3.2.2, where we conduct a detailed perturbation analysis of the control parameters. These constants are polynomials of \(( C_0, 1/\lambda_y^1, 1/\lambda_z^0)\) and other model parameters. The following Lemma 2.2 states that when the learning rate \(\eta\) is small enough, the perturbed control parameter \(\theta' = (K', L')\) in Proposition [pr:local95smoothness95Cy95Cz] is admissible and the coefficients related to local smoothness in equation 96 are bounded. Its proof is found in Lemma 2.13.

We are now ready to state and prove a global convergence result based on the PL condition and a local smoothness condition on the cost functions. Consider the sequence \((\theta^{(k)})_{k \geq 0}\) generated by the PG updates 93 .

3.2.1.4 Model-free PG for MFC with MKV simulator

We now turn to a modification of the PG updates 93 , in which the gradient will be approximated in a model-free way. The new update scheme takes the form of : \[\label{eq:one-step-update-model-free-MKV} \left\{ \begin{array}{rcl} K^{(k+1)} & = & K^{(k)} - \eta \tilde{\nabla}^{T,M,\tau}_K (\theta^{(k)}), \\ L^{(k+1)} & = & L^{(k)} - \eta \tilde{\nabla}^{T,M,\tau}_L (\theta^{(k)}), \end{array} \right.\tag{97}\] where \(\tilde{\nabla}^{T, M, \tau}(\theta) = \big( \tilde{\nabla}^{T, M,\tau}_K(\theta), \tilde{\nabla}^{T, M,\tau}_L(\theta) \big)\) is defined below, based on samples generated by a simulator.

Let us assume that we have access to the following (stochastic) simulator, called MKV simulator \(\cS^T_{MKV}\): given an admissible control parameter \(\theta \in \Theta\), \(\cS^T_{MKV}(\theta)\) returns a sample of the discounted MF cost (see Section 2.3.2.2 and equation 44 ) for the MKV dynamics 43 controlled by the process \((\ctrl_{n}^\theta)_{n \geq 0}\) and truncated at a finite horizon \(T\). In other words, it returns a realization of the truncated MF cost \[\label{eq:MKV95simulator95123n125runcated95MF95cost} \tilde{C}^T(\theta) := \sum_{n=0}^{T-1} \gamma^{n} c_{n}(\theta),\tag{98}\] where \(c_{n}(\theta) = c\big(X_{n}^{\theta}, \bar{X}_{n}^{\theta}, \ctrl_{n}^\theta, \bar{\ctrl}_{n}^\theta \big)\) is the instantaneous cost at time \(n\) (see equation 83 ).

While the simulator uses the model to generate samples, we aim to develop a model-free technique to compute the optimal control parameter. Without full knowledge of the model, we use derivative-free methods (see Section 3.2.3) to estimate the cost gradient. The main idea is to evaluate the function to minimize on a set of perturbed parameters so that one can move towards a weighted average direction that improves the value of the target function [150]. It has been shown that a careful implementation of ES can provide a performance improvement compared with classical RL algorithms based on back-propagation. Moreover, the method can be easily scaled up and adapted to distributed computer systems [151], [152]. We define some notation for the gradient estimation. Let \(\mathbb{B}_\tau\subset \mathbb{R}^{\ell \times d}\) be the ball of radius \(\tau\) under the Frobenius norm centered at the origin, and \(\mathbb{S}_\tau = \partial \mathbb{B}_\tau\) be its boundary. The uniform distributions on \(\mathbb{B}_\tau\) and \(\mathbb{S}_\tau\) are denoted by \(\mu_{\mathbb{B}_\tau}\) and \(\mu_{\mathbb{S}_\tau}\) respectively. We should think of \(v = (v^{(idy)}, v^{(com)}) \in \SS_\tau \times \SS_\tau\) as a direction for a perturbation of a parameter \(\theta = (K, L) \in \Theta\). For \(M\) perturbation directions \(\underline{v} = (v_i)_{i=1}^M = \big( (v_i^{(idy)}, v_i^{(com)}) \big)_{i=1}^M\) in which \(v_i^{(idy)}\) and \(v_i^{(com)}\) are independently sampled from \(\mu_{\SS_\tau}\), we denote by \(\theta_i = (K + v^{(idy)}_i, L + v^{(com)}_i)\) the control parameter \(\theta\) perturbed with \(v_i\). Then, the approximated gradient based on MKV simulator \(\cS^T_{MKV}\) and \(M\) perturbation directions \(\underline{v}\) on \(\SS_{\tau} \times \SS_{\tau}\) is defined by \[\label{eq:definition95approx95gradient95MKV95simulator} \tilde{\nabla}^{T, M,\tau}_K (\theta) = \frac{\ell d}{\tau^2}\frac{1}{M} \sum_{i=1}^M \tilde{C}^T(\theta_i) v^{(idy)}_i, \quad \tilde{\nabla}^{T, M,\tau}_L (\theta) = \frac{\ell d}{\tau^2}\frac{1}{M} \sum_{i=1}^M \tilde{C}^T(\theta_i) v^{(com)}_i.\tag{99}\] We omit the perturbation directions \(\underline v\) in equation 99 for clarity here.

Algorithm 14 provides a biased estimator of the true PG \(\nabla C(\theta)\). This method follows the spirit of [108], except that here we have two components playing different roles for the state and the conditional mean. The approximate gradient \(\tilde{\nabla}^{T, M, \tau}(\theta)\) thus depends on three additional hyperparameters, \((T, M, \tau)\): the cost truncation horizon, the number of perturbation directions, and the perturbation radius. In Proposition [pr:gradient95approx95with95MKV95simulator] below, we provide conditions on \((T, M, \tau)\) so that Algorithm 14 yields a sufficiently good gradient estimator to ensure the convergence of the approximate PG.

The proof of Proposition [pr:gradient95approx95with95MKV95simulator] is deferred to Section 3.2.3.5 after we demonstrate the approximation of the exact PG using derivative-free techniques based on parameter perturbations (Lemma 2.16), the approximation of the infinite-horizon cost with a truncated finite-horizon cost (Lemma 2.18), and the approximation of the expected truncated cost using sampled costs from a simulator (Lemma 2.20). It is worth noting that the bounds on \((\tau, T, M)\) in 100 , 101 , and 102 are independent of the control parameter \(\theta \in \Theta\) as long as \(C(\theta) \leq C_0\).

Now, we are ready to state the convergence results in the model-free setting with the MKV simulator \(\cS_{MKV}^T\). We consider a bound for the learning rate: \[\begin{align} &\phi_{lrate, MKV}(C_0) \\ &=\min \Big\{ \min \{ h_{var, y}^{-1}(C_0), h_{var, z}^{-1}(C_0), h_{hess, y}^{-1}(C_0), h_{hess, z}^{-1}(C_0), \\ & \qquad\qquad\qquad m_y(C_0)/(\|{\mathrm B}\|h_{grad,y}(C_0)), m_z(C_0)/(\|\tilde{\mathrm B}\|h_{grad,z}(C_0)) \} / 4, \\ & \qquad\qquad \min\big\{ h_{small-pert, y}^{-1}(C_0), h_{small-pert, z}^{-1}(C_0) \big\}\, (2 h_{cost}(C_0)) / \nu_{pl}, m_y(C_0)h_{cost}(C_0)/(\nu_{pl}\|{\mathrm B}\|), m_z(C_0)h_{cost}(C_0)/(\nu_{pl}\|\tilde{\mathrm B}\|) \Big\}, \end{align}\] where the constants \(h_{var, y}(C_0)\), \(h_{var, z}(C_0)\), \(h_{hess, y}(C_0)\), \(h_{hess, z}(C_0)\), \(h_{small-pert,y}(C_0)\), \(h_{small-pert, z}(C_0)\) depending on \(C_0\) are defined in Lemma 2.12 for small perturbation of control parameters, and \(h_{cost}(C_0)\) is defined in Lemma 2.14. The margins \(m_y(C_0)\) and \(m_z(C_0)\) are defined as in Theorem 2.1. We consider a sequence of control parameters \((\theta^{(k)})_{k \geq 0}\) generated with model-free PG update scheme 97 with the approximated gradients \(\big(\tilde{\nabla}^{T, M, \tau} (\theta^{(k)}) \big)_{k \geq 0}\) from Algorithm 14 using an MKV simulator \(\cS_{MKV}^T\).

3.2.1.5 Model-free PG for MFC with population simulator

We now turn our attention to a more realistic setting where one does not have access to an oracle for simulating the MKV dynamics, but only to an oracle capable of simulating the evolution of \(N\) agents. We then use the state sample average instead of the theoretical conditional mean, and for the social cost, the empirical average instead of the mean-field cost provided by the MKV simulator.

We rely on the following population simulator \(\cS^{T,N}_{pop}\): given a control parameter \(\theta\), \(\cS^{T,N}_{pop}(\theta)\) returns a sample of the social cost obtained by generating realizations of the \(N\)-agent state trajectories controlled by \(\theta\) and computing the associated cost for the population, see 38 and 42 . In other words, it returns a realization of \[\tilde{C}^{T,N}(\theta) = \sum_{n=0}^{T-1} \gamma^{n} \frac{1}{N} \sum_{j=1}^N c \big( X^{(j), \theta}_{n}, \bar{X}^{N, \theta}_{n}, \ctrl^{(j), \theta}_{n}, \bar \ctrl^{N, \theta}_{n} \big)\] where \(c\) is the one-step cost function 83 , and \((X_{n}^{(j), \theta}, \ctrl_{n}^{(j), \theta})_{j=1}^N\) are the state-action pairs at time \(n\) for the \(N\) agents who adopt the same policy function so that \(\ctrl_{n}^{(j), \theta} = -K \big( X_{n}^{(j), \theta} - \bar{X}_{n}^{N, \theta} \big) - L \bar{X}^{N, \theta}_{n}\) for every \(n \geq 0\) and \(j=1, \ldots, N\), and where \(\bar{X}_{n}^{N, \theta}\) and \(\bar{\ctrl}_{n}^{N, \theta}\) are the average state and average control at time \(n\). Let \(C^N(\theta) = J^N(\underline{\bU}^\theta)\) in equation 41 when all agents adopt the same control parameter \(\theta \in \Theta\).

Notice that this population simulator \(\cS^{T,N}_{pop}\) is arguably more realistic, though less powerful, than the previous MKV simulator \(\cS^T_{MKV}\). The former uses only a noisy approximation of the true mean processes, while the MKV simulator generates the exact means. We emphasize that, in this population simulator, all agents adopt the same control and therefore use the same perturbed version of the control parameter in Algorithm 15. This design aligns with the idea that the problem corresponds to an optimization problem for a central planner or a group of cooperative agents using the same control rule to minimize the social cost.

For the one-step gradient update scheme with population simulator \(\cS_{pop}^{T, N}\), we replace the gradient \(\tilde{\nabla}^{T, M, \tau}(\theta)\) obtained from Algorithm 14 with an MKV simulator by another zero-th order approximation \(\tilde{\nabla}^{T, N, M,\tau}(\theta)\) of the gradient based on \(\cS_{pop}^{T, N}\) generated by Algorithm 15. The term \(\tilde{\nabla}^{T, N, M, \tau}(\theta) = \big( \tilde{\nabla}^{T, N, M, \tau}_K( \theta), \tilde{\nabla}^{T, N, M, \tau}_L(\theta) \big)\) is called the sampled population policy gradient at \(\theta \in \Theta\). It is defined with \(M\) perturbation directions \((v_i)_{i=1}^M = (v_i^{(idy)}, v_i^{(com)})_{i=1}^M\) on \(\SS_\tau \times \SS_\tau\): \[\label{eq:sampled95population95gradient} \tilde{\nabla}_K^{T, N, M, \tau}(\theta) = \frac{\ell d}{\tau^2} \frac{1}{M} \sum_{i=1}^M \tilde{C}^{T, N}(\theta_i) v_i^{(idy)}, \quad \tilde{\nabla}_L^{T, N, M, \tau}(\theta) = \frac{\ell d}{\tau^2} \frac{1}{M} \sum_{i=1}^M \tilde{C}^{T, N}(\theta_i) v_i^{(com)},\tag{103}\] where \(\theta_i = \theta + v_i\) is the perturbed control parameter with \(v_i\). The at iteration step \(k\) for a parameter \(\theta^{(k), pop} = (K^{(k), pop}, L^{(k), pop} )\) is defined by \[\label{eq:one-step-update-model-free-PoP} \left\{ \begin{array}{rcl} K^{(k+1), pop} & = & K^{(k), pop} - \eta \tilde{\nabla}^{T, N, M,\tau}_K (\theta^{(k), pop}), \\ L^{(k+1), pop} & = & L^{(k), pop} - \eta \tilde{\nabla}^{T, N, M,\tau}_L (\theta^{(k), pop}). \end{array} \right.\tag{104}\]

The main result of this section is Theorem 2.3. It establishes the convergence of the above learning scheme toward the optimal control parameter \(\theta^* = (K^*, L^*)\) of the MF problem, using a population simulator \(\cS_{pop}^{T, N}\) of the social cost of \(N\) homogeneous agents. Before proving this result, we state Proposition [pr:approx95with95social95cost95no95MF95cost] and Proposition [pr:modelfree95population95gradient95approx], which provide crucial approximation estimates for the social cost by its \(N\)-agent equivalent, as well as the approximation of the sampled population policy gradient. Their proofs are given in Section 3.2.4.2.

Now, we show that Algorithm 15 successfully learns the optimal control parameter of the MF problem with a global linear convergence rate in the social cost. To state the theorem, we consider a sequence of control parameters \((\theta^{(k), pop})_{k \geq 0}\) generated by the model-free PG update scheme with the sampled population policy gradient 104 .

3.2.2 Technical results in Section 3.2.1.3↩︎

This section collects the technical lemmas used for proving the statements in Section 3.2.1.3. These lemmas are adapted from [108] with necessary modifications to cope with the idiosyncratic noise and the common noise processes present in the state dynamics. To make the paper self-contained, we provide their proofs here. Because there is an analogy between the processes \(\bY\) and \(\bZ\) in their definitions and usages, we will state and prove results related to the process \(\bY\) and the parameter \(K\) only.

3.2.2.1 Cost expressions

For an admissible control \(\theta = (K, L) \in \Theta\), consider two \(L^2\)-discounted integrable processes in \(\bY^{K, \xi_y}, \bZ^{L, \xi_z} \in S\) following dynamics 87 , 88 with initial values \(Y_{0} = Y_{0}^{K, \xi_y} = \xi_y\) and \(Z_{0} = Z_{0}^{L, \xi_z}= \xi_z\). We defined the corresponding value functions \(V_y\) and \(V_z\) on \(\RR^{\ell \times d} \times \RR^d\) by \[V_y(K, \xi_y) = \EE_{\boldsymbol{\varepsilon}} \Big[ \sum_{n \geq 0} \gamma^{n} f(Y_{n}^{K, \xi_y}, K, Q, R) \Big], \quad V_z(L, \xi_z) = \EE_{\boldsymbol{\varepsilon^0}} \Big[ \sum_{n \geq 0} \gamma^{n} f(Z_{n}^{L, \xi_z}, L, \tilde{Q}, \tilde{R}) \Big]\] where \(f(\xi, K, Q, R) := \xi^\top( Q + K^\top R K ) \xi\) and the expectations \(\EE_{\boldsymbol{\varepsilon}}, \EE_{\boldsymbol{\varepsilon^0}}\) are integrations with respect to the idiosyncratic and the common noise processes.

These value functions can be computed simply with the help of solutions \(P_y^K\) and \(P_z^L\) to the two discrete-time Lyapunov equations (DLEs): \[\label{eq:lyapunov95eq95theta} \left\{ \begin{array}{rcl} P^y_K & = & Q + K^\top R K + \gamma ({\mathrm A} - {\mathrm B} K)^\top P^y_K ({\mathrm A} - {\mathrm B} K) \\ P^z_L & = & \tilde{Q} + L^\top \tilde{R} L +\gamma ( \tilde{{\mathrm A}}- \tilde{\mathrm B} L)^\top P^z_L ( \tilde{{\mathrm A}} - \tilde{\mathrm B} L). \end{array} \right.\tag{105}\] The existence of solutions to 105 is deferred to the next Section 3.2.2.4. Let us also define \[\label{eq:value95function95noise95constant} \alpha^y_K = \EE \big[ \sum \nolimits_{n \geq 1} \gamma^{n} (\varepsilon_{n})^\top P^y_K \varepsilon_{n} \big], \qquad \alpha^z_L = \EE \big[ \sum \nolimits_{n \geq 1} \gamma^{n} (\varepsilon^0_{n})^\top P^z_L \varepsilon^0_{n} \big].\tag{106}\] Then, using the dynamics of \(\bY^{K, \xi_y}\) 87 and the discrete-time Lyapunov equation 105 for \(P_K^y\), we have \(V_y(K, \xi_y) = \xi_y^\top P_K^y \xi_y + \alpha_K^y.\) When the initial states are random, e.g., \(Y_{0} = X_{0} - \bar{X}_0\), we have \[\label{eq:cost95expression95with95PK95y} C_y(K) = \EE[ V_y(K, Y_{0}) ] = \langle P_K^y, \, \Sigma_{Y_{0}} + \frac{\gamma}{1-\gamma} \Sigma^1 \rangle_{tr} = \langle Q + K^\top R K, \, \Sigma_K \rangle_{tr}.\tag{107}\] Because \(\lambda_y^1 = \lambda_{min}( \Sigma_{Y_{0}} + \frac{\gamma}{1-\gamma} \Sigma^1 )\), we derive the following bounds: \[\label{eq:upper95bound95Riccati95and95Variance95matrix} \|P^y_K\| \leq C_y(K) / \lambda_y^1, \quad \|\Sigma^y_K\| \leq C_y(K) / \lambda_{min}(Q).\tag{108}\] Similar expressions for the cost \(C_z(L)\), value \(V_z(L, \xi_z)\), and the bounds for \(P_L^z\) and \(\Sigma_L^z\) follow with \((P_L^z, \Sigma_{Z_{0}}, \Sigma^0, \lambda_{z}^0, \tilde{Q}, \tilde{R})\).

For agent \(j\), her control at time \(n\) takes the form \(\ctrl_{n}^{j} = - K (X_{n}^{j} - \bar{X}_{n}^{N} ) - L \bar{X}_{n}^N\). And the corresponding auxiliary processes \(Y_{n}^{K, j, N} := X_{n}^{j} - \bar{X}^N_{n}\) and \(Z_{n}^{L, N} := \bar{X}^N_{n}\) follow the dynamics similar to the ones defined in equations 87 –88 , but by replacing the noise terms with \(\varepsilon_{n+1}^{1, j, N} := \varepsilon_{n+1}^{1, (j)} - \frac{1}{N} \sum_{i=1}^N \varepsilon_{n+1}^{1, (i)}\) and \(\varepsilon_{n+1}^{0, N} := \varepsilon_{n+1}^0 + \frac{1}{N} \sum_{i=1}^N \varepsilon_{n+1}^{1, (i)}\). Then, the cost for agent \(n\) can be expressed into By consequence, the social cost of the population given by 41 is bounded by

When these agents adopt the same control parameter \(\theta\), with the same arguments used in the homogeneous case, the difference between the auxiliary cost \(C^{(j)}_y(K)\) for agent \(j\) and the auxiliary cost in the mean-field setting \(C_y(K)\) satisfies The gap between the social cost and the MF cost tends towards 0 when we impose a bounded domain for the control parameters and when \(h_{\delta Q} \to 0\) and \(N \to \infty\).

3.2.2.2 Gradient expression

We express the gradient of cost \(C_y(K)\) in terms of the variance matrix \(\Sigma_K\). This gradient expression is crucial in the whole paper. Consequently, the approximation of gradient \(\nabla_K C_y(K)\) could hinge on an estimation of variance matrix \(\Sigma_K\). For completeness, we provide the proof here, which is analogous to [108].

3.2.2.3 Cost variation expression

The following lemma shows the variation in cost \(C_y\).

3.2.2.4 Solution to discrete-time Lyapunov equations

For any admissible parameters \(\theta \in \Theta\), we can construct solutions to the DLEs 105 .

The following lemma completes the connections between the solution matrices \((P, \bar P)\) of the Riccati equations ?? and the solution matrices of (DLEs) 105 with the specific parameters \(K^*\) and \(L^*\) defined in Theorem 3.1.

3.2.2.5 Perturbation analysis

3.2.2.5.1 Perturbation on variance

In this subsection, we show that when \(\Delta K = K' - K\) is small enough, the norm between variance matrices is controlled by \(\|\Delta K \|\). To do so, we introduce a function \(f_{var}\) on matrix \(K\) defined by \[\label{eq:f95var95K} f_{var}(K) = 4 ( 1 + \| {\mathrm A} - {\mathrm B} K \| ) \| {\mathrm B} \| \frac{ C_y(K)}{ \lambda_y^1 \lambda_{min}(Q) }.\tag{110}\] We also introduce the following operators: for any \(\Sigma \in \cM\cS\), \[\label{eq:def-ope-F95K95y} \cF_K(\Sigma) = \gamma ({\mathrm A}- {\mathrm B} K) \Sigma ({\mathrm A}-{\mathrm B} K)^\top, \, \cT_K(\Sigma) = \sum_{n \geq 0} \gamma^{n} ({\mathrm A} - {\mathrm B} K)^{n} \Sigma \big( ({\mathrm A} - {\mathrm B} K)^{n} \big)^\top.\tag{111}\] Consider the operator norms on \(\cM\cS\) for \(\cF_K\) and \(\cT_K\), and denote them by \(\vertiii{\cF_K}\) and \(\vertiii{\cT_K}\). Because \(\vertiii{\cF_K} = \sup_{\| \Sigma \| = 1} \cF_K(\Sigma) \leq \gamma \| {\mathrm A} - {\mathrm B} K \|^2 < 1\) for a control parameter \(K\) in \(\theta = (K, L ) \in \Theta\), we then have \[\label{eq:cT} \cT_K(\Sigma) = \sum \nolimits_{n \geq 0} (\cF_K)^{n}(\Sigma) = ( \bI - \cF_K)^{-1} (\Sigma).\tag{112}\]

The following lemma is about the bounds with operators \(\cF_K\) and \(\cT_K\).

Now, we are ready to show a perturbation result for the variance matrix.

3.2.2.5.2 Perturbation on the DLE solution matrix

Consider the solutions \(P_K^y\) and \(P_{K'}^y\) of the first DLE equation in 105 for parameters \(K\) and \(K'\). We introduce a function \(f_{riccati}\) on the matrix \(K\) defined by \[\label{eq:f95riccati95K} f_{riccati}(K, \Delta K) = \frac{ C_y(K) }{\lambda_y^1} \bigg( \gamma f_{var}(K) + \frac{4\|R\|}{\lambda_{min}(Q) } \| K \| + \frac{ 2 \| R \|}{ \lambda_{min}(Q)} \| \Delta K \| \bigg)\tag{113}\]

Perturbation on gradient

We show in this subsection that the difference between gradients \(\nabla_K C_y\) at \(K\) and \(K'\) is controlled by \(\| K' - K \|\). To do so, we consider a function \[\label{eq:f95grad95K} f_{grad}(K) = \Big( 6 \| R \| + \frac{ 4 \| {\mathrm B}\|^2 C_y(K)}{\lambda_y^1} + 6 f_{var}(K) \| R \| \| K \| + \frac{3}{2} ( f_{var}(K))^2 \lambda_{min}(Q) \Big)\frac{C_y(K)}{\lambda_{min}(Q)}.\tag{114}\]

3.2.2.6 Admissibility of perturbed control parameter

This subsection is to show that when the perturbations \(\| \Delta K \| = \| K' - K \|\) and \(\| \Delta L \| = \| L' - L \|\) are small enough, the perturbed parameter \(\theta' = (K', L')\) is still admissible, i.e. \(\theta' \in \Theta\).

3.2.2.7 Proof of Proposition [pr:PL95cond95Cy95Cz] for PL condition

We can also derive a lower bound for \(C_y(K) - C_y(K^*)\). Consider \(K' = K - (R + {\mathrm B}^\top P_K^y {\mathrm B})^{-1} E_K\) so that the advantage quantities in Lemma 2.4 attain the lower bounds. Consequently, with \(R \succ 0\) and \(P_K^y \succ 0\) (Assumption 3.2 and Proposition [prop:existence95of95lyapunov95equation]), we have \[\begin{align} C_y(K) - C_y(K^*) & \geq C_y(K) - C_y(K') \\ & = \langle (R + \gamma {\mathrm B}^\top P_K^y {\mathrm B})^{-1} E_K, \, E_K \Sigma_{K'} \rangle_{tr} \nonumber \\ & \geq \lambda_y^1 \, \| E_K \|_F^2 / \| R + \gamma {\mathrm B}^\top P_K^y {\mathrm B} \|. \label{eq:lower95bound95for95Cy95minus95Copt} \end{align}\tag{115}\]

3.2.2.8 Proof of Proposition [pr:local95smoothness95Cy95Cz] for local smoothness

3.2.2.9 Bounds in perturbation analysis

To demonstrate the convergence results for the exact PG descent algorithm in Theorem 2.1, we first establish several bounds on the coefficients \(\lambda_{var, y}\) and \(\lambda_{hess,y}\) introduced in Proposition [pr:local95smoothness95Cy95Cz].

3.2.2.9.1 Perturbation on cost function

3.2.3 Technical lemma for Section 3.2.1.4↩︎

3.2.3.1 Intuition for zeroth-order optimization

Derivative-free optimization (see e.g. [150], [154]) tries to find an optimizer of a function \(f: x \mapsto f(x)\) using only pointwise values of \(f\), without having access to its gradient. The basic idea is to express the gradient of \(f\) in a point \(x\) by using values of \(f(\cdot)\) at a Gaussian perturbation region around point \(x\). This can be achieved by introducing an approximation of \(f\) using a Gaussian smoothing. Formally, if we define the perturbed function for some \(\sigma >0\) as \(f_{\sigma^2}(x) = \mathbb{E}[f(x+ \varepsilon)],\) where \(\varepsilon \sim \mathcal{N}(0,\sigma^2 I)\), then its gradient can be written as \(\frac{d}{dx}f_{\sigma^2}(x) = \frac{1}{\sigma^2} \mathbb{E}[f(x+\varepsilon) \varepsilon].\) We can replace the Gaussian smoothing term \(\varepsilon\) by a uniform random variable with bounded norms.

Consider \(U = (U^{(idy)}, U^{(com)})\) with \(U^{(idy)},U^{(com)}\) being two independent random matrices in \(\RR^{\ell \times d}\) uniformly distributed on the ball with radius \(\tau\): \(\mathbb{B}_{\tau} = \{ W \in \RR^{\ell \times d} : \| W \|_F \leq \tau \}\). The unit ball under Frobenius norm for matrices in \(\RR^{\ell \times d}\) can be identified as a unit ball for vectors in \(\RR^{\ell d}\). Let \(\mathbb{S}_\tau = \{ W \in \RR^{\ell \times d}: \| W \|_F = \tau \}\) be the sphere of \(\mathbb{B}_\tau\) of dimension \(\ell d -1\). The ratio between the sphere area of a ball with radius \(\tau\) in \(\RR^{\ell d}\) to its volume equals \((\ell d) / \tau\) [155]. Furthermore, the spectral norm of \(W \in \mathbb{B}_\tau\) satisfies that \(\| W \| \leq \| W \|_F = \tau\).

For any \(\theta = (K,L)\) and \(\tau > 0\), we introduce the following smoothed version \(C_\tau\) of \(C\), defined by: \[\label{eq:formula95perturbation95definition95C} C_\tau(\theta) = \mathbb{E}_U[C(\theta + U)] = \mathbb{E}_{(U^{(idy)},U^{(com)})}[C_y(K+U^{(idy)}) + C_z(L+U^{(com)})],\tag{116}\] where \(\mathbb{E}_U\) denotes the expectation over \(U\), and the costs \(C_y(K+U^{(idy)})\) and \(C_z(L + U^{(com)})\) are evaluated using the auxiliary processes \(\bY^{K + U^{(idy)}}\) and \(\bZ^{L + U^{(com)}}\) as defined in equations 89 . Let \(V^{(idy)} \sim \tau U^{(idy)} / \| U^{(idy)} \|_F\) and \(V^{(com)} \sim \tau U^{(com)} / \| U^{(com)} \|_F\) be two random matrices independently and uniformly distributed on the sphere \(\SS_\tau\); then \(V = (V^{(idy)}, V^{(com)}) \in \SS_\tau \times \SS_\tau\). From [155] or [108], we get the following expressions for the gradients of \(C_\tau(\theta)\): \[\begin{align} \tag{117} \nabla_K C_\tau(\theta) & = \frac{\ell d}{\tau^2} \mathbb{E}_{V^{(idy)}} \big[ C_y(K + V^{(idy)}) V^{(idy)} \big] = \frac{\ell d}{\tau^2}\mathbb{E}_V \big[C(\theta + V) V^{(idy)} \big], \\ \tag{118} \nabla_L C_\tau(\theta) & = \frac{\ell d}{\tau^2} \mathbb{E}_{V^{(com)}} \big[ C_z(L + V^{(com)}) V^{(com)} \big] = \frac{\ell d}{\tau^2} \mathbb{E}_V \big[C(\theta+V) V^{(com)} \big]. \end{align}\] The last equalities in 117 and 118 are justified by Lemma 2.1 and the fact that \(V^{(idy)}\) and \(V^{(com)}\) are independent zero-mean random variables.

To show that the gradient estimator from Algorithm 14 provides a good approximation of the exact PG, we introduce two more gradient estimators before we dive into the perturbation analysis in this section. Let \(M\) be the number of perturbation directions, \(\tau > 0\) be the perturbation radius, and \(T\) be the truncation horizon for the infinite discounted MF cost. The following four notations of gradient will be used in the next subsections:

1. [Policy gradient] \(\nabla C(\theta) = \big(\nabla_K C(\theta), \nabla_L C(\theta) \big)\) represents the exact gradient of cost \(C(\theta)\) with respect to \(\theta = (K, L)\).

2. [Perturbed policy gradient] \(\hat{\nabla}^{M,\tau}(\theta, \underline v)\) defined in 119 below is an approximation of exact PG based on the zeroth-order approximation technique with \(M\) perturbation directions \(\underline v = (v_i)_{i=1}^M\) where \(v_i \sim \mu_{\mathbb{S}_\tau} \otimes \mu_{\mathbb{S}_\tau}\).

3. [Truncated policy gradient] \(\hat{\nabla}^{T, M, \tau}(\theta, \underline v)\) defined in 126 below is an approximation of perturbed policy gradient using a cost truncated to a finite horizon \(T\).

4. [Sampled policy gradient] \(\tilde{\nabla}^{T, M,\tau} \big(\theta, \underline v \big)\) defined in equation 99 is the output of Algorithm 14.

Before detailing the exact bounds, we sketch the asymptotic requirements to achieve an \(\varepsilon\)-approximation of the exact PG. We will show that the perturbation radius must scale as \(\tau = \mathcal{O}(\varepsilon)\), the truncation horizon as \(T = \mathcal{O}(\log(1/\varepsilon))\), and the sample size as \(M = \tilde{\mathcal{O}}(\varepsilon^{-4})\). The following sections make these dependencies explicit through threshold functions \(\phi_{pert, radius}, \phi_{trunc, T},\) and \(\phi_{sample, size}\).

Our proof strategy bounds the total approximation error via the triangle inequality: \[\begin{align} \| \tilde{\nabla}^{T, M, \tau} - \nabla C \| & \leq \underbrace{\| \tilde{\nabla}^{T, M, \tau} - \hat{\nabla}^{T, M, \tau} \|}_{\text{Sampling Error}} + \underbrace{\| \hat{\nabla}^{T, M, \tau} - \hat{\nabla}^{M, \tau} \|}_{\text{Truncation Error}} + \underbrace{\| \hat{\nabla}^{M, \tau} - \nabla C \|}_{\text{Zeroth-order Error}}. \end{align}\] The following subsections bound each of these distinct sources of error.

3.2.3.2 Approximation with smoothed cost

In this subsection, we show an approximation result of the exact PG \(\nabla C(\theta)\) using the perturbed policy gradient \(\hat{\nabla}^{M, \tau}(\theta, \underline v)\) defined in 119 below. To start, we consider the gradient of a smooth cost \(C_\tau(\theta)\) defined in 116 with respect to \(\theta = (K, L)\) as \(\nabla C_\tau(\theta) = \big( \nabla_K C_\tau(\theta) , \nabla_L C_\tau(\theta) \big)\). Equations 117 and 118 provide the expression for \(\nabla C_\tau(\theta)\) without taking the gradients. Now, we approximate the expectations \(\EE_V\) therein using \(M\) independent samples \((v_i^{(idy)})_{i=1}^M\) of \(V^{(idy)}\) and \(M\) independent samples \((v_i^{(com)})_{i=1}^M\) of \(V^{(com)}\) respectively. Let \(\underline v = (v_i)_{i=1}^M\) with \(v_i = (v_{i}^{(idy)}, v_{i}^{(com)}) \in \RR^{\ell \times d} \times \RR^{\ell \times d}\) and \(v_i \sim \mu_{\mathbb{S}_\tau} \otimes \mu_{\mathbb{S}_\tau}\). We define the “perturbed policy gradient" in \(K\) and in \(L\) by \[\label{eq:def95perturbed95PG95MKV} \hat{\nabla}_K^{M, \tau}(\theta, \underline v) = \frac{\ell d}{\tau^2} \frac{1}{M} \sum_{i=1}^M C(\theta + v_i) v_i^{(idy)}, \quad \hat{\nabla}_L^{M, \tau}(\theta, \underline v) = \frac{\ell d}{\tau^2} \frac{1}{M} \sum_{i=1}^M C(\theta + v_i) v_i^{(com)}.\tag{119}\] Let \(\hat{\nabla}^{M, \tau}(\theta, \underline v) = \big( \hat{\nabla}_K^{M, \tau}(\theta, \underline v), \hat{\nabla}_L^{M, \tau}(\theta, \underline v) \big)\). In 119 , we consider \(\tau\) small enough so that by the stability Lemma 2.11, we have \(\theta + v_i \in \Theta\) for all \(i=1, \ldots, M\), and thus \(\hat{\nabla}^{M, \tau}(\theta, \underline v)\) is well-defined. It is clear that the randomness of the perturbed policy gradient comes from the random directions \(\underline v \in (\mathbb{S}_{\tau} \times \mathbb{S}_\tau )^{M}\).

Based on the conditions in Lemma 2.15, we define the following two functions for the perturbation radius and the number of perturbation directions: \[\begin{align} & \phi_{pert, radius}(\varepsilon, C_0) = \max \big\{ h_{cost}(C_0) / C_0, h_{small-pert, y}(C_0), h_{small-pert, z}(C_0), \nonumber \\ &4 h_{f\_grad\_K}(C_0) / \varepsilon, 4 h_{f\_grad\_L}(C_0) / \varepsilon \big\} \tag{120} \\ & \phi_{pert, size}(\varepsilon, \tau, C_0, \delta_{pert} ) = \log \Big( \frac{ 2(\ell + d) }{\delta_{pert}} \Big) \Big( \frac{16 \ell d C_0 }{ \tau \varepsilon } + 1 \Big)^2 \tag{121} \end{align}\] where \(h_{cost}(C_0)\), \(h_{small-pert, y}(C_0)\), \(h_{small-pert, z}(C_0)\), \(h_{f\_grad\_K}(C_0)\), \(h_{f\_grad\_L}(C_0)\) are defined in Lemma 2.12 and Lemma 2.14.

3.2.3.3 Approximation with truncated cost

In this subsection, we first provide an approximation of the MF cost \(C(\theta)\) with a truncated cost, then we show a result for the approximation of truncated policy gradient \(\hat{\nabla}^{T, M, \tau}\) defined later in the subsection. For \(T \geq 0\), let us define the truncated variances and the truncated costs by: \[\begin{align} \tag{122} \Sigma_{K}^{y,T} & = \mathbb{E}\big[\sum_{n=0}^{T-1} \gamma^{n} (Y_{n}^{K}) (Y_{n}^K)^\top \big], \quad \Sigma_{L}^{z,T} = \mathbb{E} \big[ \sum_{n=0}^{T-1} \gamma^{n} (Z_{n}^L) (Z_{n}^L)^\top \big], \\ \tag{123} C^T_y(K) & = \mathbb{E} \big[ \sum_{n=0}^{T-1} \gamma^{n} f(Y_{n}^{K}, K, Q, R) \big] , \quad C^T_z(L) = \mathbb{E} \big[ \sum_{n=0}^{T-1} \gamma^{n} f(Z_{n}^{L}, L, \tilde{Q}, \tilde{R} ) \big], \end{align}\] where \(f(\xi, \phi, q, r) = \xi^{\top} (q + \phi^\top r \phi)\xi\). We define \(C^T(\theta) = \EE[ \sum_{n=0}^{T-1} \gamma^{n} c_{n}(\theta) ]\) where \(c_{n}(\theta) = c(X_{n}^\theta, \bar{X}_{n}^\theta, \ctrl_{n}^\theta, \bar \ctrl_{n}^\theta)\) is the instantaneous cost at time \(n\). With the help of auxiliary process \((\bY^K, \bZ^L)\) from 87 and 88 , we have \(C^T(\theta) = C_y^T(K) + C_z^T(L)\).

For \(\theta = (K, L) \in \Theta\) with \(C(\theta) \leq C_0\), assume that the sublevel set has the positive Euclidean margins \(m_y(C_0)\) and \(m_z(C_0)\) defined in Theorem 2.1. Then \[\gamma \| {\mathrm A} - {\mathrm B} K \|^2 \leq \gamma (1/\sqrt{\gamma}-m_y(C_0))^2, \quad \gamma \| \tilde{\mathrm A} - \tilde{\mathrm B} L \|^2 \leq \gamma (1/\sqrt{\gamma}-m_z(C_0))^2.\] Let us define a constant \(\gamma_{pert}\) and a function \(\phi_{truncation}(\varepsilon)\) defining the lower bound on the truncation time horizon: \[\begin{align} & \gamma_{pert} = \max \big\{ \gamma, \gamma (1/\sqrt{\gamma}-m_y(C_0))^2,\, \gamma (1/\sqrt{\gamma}-m_z(C_0))^2 \big\}, \tag{124} \\ & \phi_{truncation}(\varepsilon) = \Big( \frac{\log( 2d C_{init, noise}^2) - \log(\varepsilon (1 - \gamma_{pert})^2 ) + 1 }{\log(1/\gamma_{pert})} \Big)^2. \tag{125} \end{align}\] We have the following result between the truncated cost \(C^T(\theta)\) and the MF cost \(C(\theta)\).

We define the truncated policy gradients in \(K\) and in \(L\) at time \(T\) with perturbation size, perturbation radius, and perturbation directions \((M, \tau, \underline v)\) by \[\label{eq:def95123n125runcated95gradient95MKV} \hat{\nabla}_K^{T, M, \tau}(\theta, \underline v) = \frac{\ell d}{M \tau^2} \sum_{i=1}^M C^T(\theta + v_i) v_i^{(idy)}, \quad \hat{\nabla}_L^{T, M, \tau}(\theta, \underline v) = \frac{\ell d}{M \tau^2} \sum_{i=1}^M C^T(\theta + v_i) v_i^{(com)}.\tag{126}\] Let \(\hat{\nabla}^{T, M, \tau}(\theta, \underline v) = \big( \hat{\nabla}_K^{T, M, \tau}(\theta, \underline v), \, \hat{\nabla}_L^{T, M, \tau}(\theta, \underline v) \big)\). Here, we assume that the perturbation radius \(\tau\) is small enough so that \(\theta + v_i \in \Theta\) for all \(i=1, \ldots, M\), and thus the truncated policy gradient \(\hat{\nabla}^{T, M, \tau}(\theta, \underline v)\) is well-defined.

To show the approximation of perturbed policy gradient \(\hat{\nabla}^{M, \tau}\) with a truncated policy gradient \(\hat{\nabla}^{T, M, \tau}\), based on Lemma 2.17, we define \[\begin{align} h_T(C_0) & = d.(\| Q \| + \| \tilde{Q}\| + \|R \| h_K(C_0)^2 + \| \tilde{R} \| h_L(C_0)^2 ), \tag{127} \\ \phi_{trunc, T}(\varepsilon, \tau, C_0) & = \phi_{truncation}\big( \frac{\tau \varepsilon }{2 \ell d. h_T(C_0)} \big) \nonumber \\ & = \Big( \frac{ \log( 4 \ell d^2 C_{init, noise}^2 h_T(C_0) ) - \log( \varepsilon \tau (1 - \gamma_{pert})^2 ) + 1 }{\log(1/\gamma_{pert})} \Big)^2. \tag{128} \end{align}\] where \(h_{K}(C_0)\) and \(h_L(C_0)\) are defined in Lemma 2.12, and \(\gamma_{pert}\), \(\phi_{truncation}\) are defined in equations 124 and 125 .

3.2.3.4 Sampled costs with MKV simulator

In this section, we will first show a matrix concentration inequality (Lemma 2.19) for a quadratic form \((Y_{n}^K)^\top ( Q + K^\top R K) Y_{n}^K\) where \(Y_{n}^K \in \RR^{d}\) is a term of the auxiliary process \(\bY^{K}\) at time \(n\), specifically a sub-Gaussian random vector with dependent coordinate entries. These quadratic forms from the auxiliary process \(\bY^K\) as well as those from \(\bZ^L\) are closely related to the sample costs simulated from an MKV simulator \(\cS_{MKV}^T\). Based on this result, we will quantify the approximation error between the sampled policy gradients as outputs of Algorithm 14 and the truncated policy gradient defined in equation 126 . To start, we define a constant \[\label{eq:h95psi9524432y} h_{\psi_2, y}(\tau, C_0) = \frac{2 \tau d^2 C_{init, noise}^2 \Big( \| Q \| + \|R\| ( \tau + h_K(C_0) )^2 \Big) }{(1 - \sqrt{\gamma})^2}.\tag{129}\] A similar constant \(h_{\psi_2, z}(\tau, C_0)\) is defined with \((\tilde{Q}, \tilde{R}, h_L(C_0))\) in the above equation. To abbreviate the notation, we also consider a real-value function \(f\) on \(\RR^{d} \times \RR^{\ell \times d} \times \RR^{d \times d} \times \RR^{\ell \times \ell}\) given by \(f(\xi, \varphi, q, r) \mapsto \xi^\top ( q + \varphi^\top r \varphi ) \xi\).

Now, we show the approximation between the sampled policy gradients from Algorithm 14 and truncated policy gradients defined in 126 . We recall that the sampled policy gradients at \(\theta = (K, L)\) with \(M\) perturbation directions \(\underline v = ( (v_i^{(idy)}, v_i^{(com)} ) )_{i=1}^M\) on \(\SS_{\tau} \times \SS_{\tau}\) with radius \(\tau > 0\) and from an MKV simulator \(\cS_{MKV}^{T}\) is given by \[\label{eq:def95sampeld95policy95gradient95MKV95restate} \tilde{\nabla}^{T, M,\tau}_K (\theta, \underline v) = \frac{\ell d}{\tau^2}\frac{1}{M} \sum_{i=1}^M \tilde{C}^T(\theta_i) v^{(idy)}_i, \quad \tilde{\nabla}^{T, M,\tau}_L (\theta, \underline v) = \frac{\ell d}{\tau^2}\frac{1}{M} \sum_{i=1}^M \tilde{C}^T(\theta_i) v^{(com)}_i,\tag{130}\] where \(\theta_i = (K + v_i^{(idy)}, L + v_i^{(com)})\) for \(i=1, \ldots, M\). We also define the following constants to bound the perturbation size \(M\) in the sampled costs: \[\begin{align} \phi_{sample,size}(\varepsilon, \tau, T, C_0, \delta_{sample}) = \log \Big( \frac{4T (\ell + d)}{\delta_{sample}} \Big) \Big( \frac{8 \ell d T}{\varepsilon \tau^2 } h_{\psi_2}(\tau, C_0) + \frac{1}{4} \Big)^2 \label{eq:phi95sample95size} \end{align}\tag{131}\] where \(h_{\psi_2}(\tau, C_0) = \max \big\{ h_{\psi_2, y}(\tau, C_0),\, h_{\psi_2, z}(\tau, C_0) \big\}\) with \(h_{\psi_2, y}(\tau, C_0)\) and \(h_{\psi_2, z}(\tau, C_0)\) defined in 129 .

3.2.3.5 Proof of Proposition [pr:gradient95approx95with95MKV95simulator]

In this subsection, we show that when the truncation horizon \(T\) is large enough for an MKV simulator \(\cS_{MKV}^T\), the perturbation radius \(\tau\) is small enough, and the perturbation size \(M\) is large enough in Algorithm 14, the output of the approximated policy gradient \(\tilde{\nabla}^{T, M, \tau}\) from the algorithm is close to the exact PG \(\nabla C\) with high probability. To do so, we define \[\begin{align} \phi_{pert, radius, MKV}( \varepsilon, C_0) & = \phi_{pert, radius}( \varepsilon / 3, C_0) \tag{132} \\ \phi_{trunc, T, MKV}(\varepsilon, \tau, C_0) & = \phi_{trunc, T}( \varepsilon / 3, \tau, C_0) + 2 \tag{133} \\ \phi_{sample, size, MKV}(\varepsilon, \tau, T, C_0, \delta_{approx}) & = \phi_{pert, size} \big(\varepsilon / 3, \tau, C_0, \delta_{approx} / 2 \big) \nonumber \\ &+ \phi_{sample, size} \big(\varepsilon / 3, \tau, T, C_0, \delta_{approx} / 2\big) \tag{134} \end{align}\] where \(\phi_{pert, radius}\) 120 , \(\phi_{pert, size}\) 121 , \(\phi_{trunc, T}\) 128 , \(\phi_{sample, size}\) 131 are polynomial in \((d, \ell, C_0, (\lambda_y^1)^{-1}, (\lambda_z^0)^{-1}, C_{init, noise}, \varepsilon^{-1}, \tau^{-1})\) and model parameters defined in sections 3.2.3.2, 3.2.3.3, and 3.2.3.4.
Let \(\varepsilon > 0\) be a target precision. We assume that

3.2.3.6 Additional lemmas for sub-Gaussian random vectors

The following lemma is similar to [157].

The following lemma can be found in [96] to connect the \(L^p\) norms with the sub-exponential norm and the sub-Gaussian norm.

The following lemma is about the dilatation technique used for matrix concentration inequalities with non-symmetric matrices; see [156].

3.2.4 Technical proofs for Section 3.2.1.5↩︎

3.2.4.1 Auxiliary costs for \(N\)-agent problem

To show the model-free population-based gradient estimator \(\tilde{\nabla}^{T,N,M, \tau}\) from Algorithm 15 provides a good approximation for the exact PG \(\nabla C\), similar to the approach in Section 3.2.3 for the MKV simulator, we define the perturbed population policy gradient \(\hat{\nabla}^{N, M, \tau}\) and the truncated population policy gradient \(\hat{\nabla}^{T, N, M, \tau}\) by replacing the perturbed MF cost \(C(\theta + v_i)\) in 119 and the truncated perturbed cost \(C^T(\theta + v_i)\) in 126 with the corresponding \(N\)-agent social costs adapted from equation 41 .

Let \([N] = \{1, \ldots, N\}\). For \(\theta = (K, L) \in \Theta\), we introduce auxiliary processes \((\bY^{K, j, N})_{j \in [N]}\) and \(\bZ^{L, N}\) for the \(N\)-agents by \(Y_{n}^{K, j, N} = X_{n}^{(j), \theta} - \bar X_{n}^{N, \theta}\) and \(Z_{n}^{L, N} = \bar X_{n}^{N, \theta} = \frac{1}{N} \sum_{j=1}^N X_{n}^{(j), \theta}\) for all \(n \geq 0\) and \(j \in [N]\), where \((X_{n}^{(j), \theta})_{n\geq 0}\) is the state process for agent \(j\) following dynamics 38 with control at time \(n\) given by \(\ctrl_{n}^{(j), \theta} = - K (X_{n}^{(j), \theta} - \bar{X}_{n}^{N, \theta}) - L \bar{X}_{n}^{N, \theta}\). The dynamics of \((\bY^{K, j, N})_{j \in [N]}\) and \(\bZ^{L, N}\) then follow \[\begin{align} & Y_{n+1}^{K, j, N} = ( {\mathrm A} - {\mathrm B} K ) Y_{n}^{K, j, N} + \tilde{\varepsilon}_{n+1}^{(j)}, \quad & & Y_{0}^{j, N} = \varepsilon_0^{(j)} - \tfrac{1}{N}\sum \nolimits_{i \in [N]} \varepsilon_0^{(i)}, \tag{135} \\ & Z_{n+1}^{L, N} = ( \tilde{\mathrm A} - \tilde{\mathrm B} L ) Z_{n}^{L, N} + \varepsilon_{n+1}^0 + \overline{{\varepsilon}_{n+1}^{N}}, \quad & & Z_{0}^{N} = \varepsilon_0^0 + \tfrac{1}{N}\sum \nolimits_{i \in [N]} \varepsilon_0^{(i)}, \tag{136} \end{align}\] where \(\overline{{\varepsilon}_{n+1}^{N}} = \frac{1}{N} \sum_{j=1}^N \varepsilon_{n+1}^{(j)}\) and \(\tilde{\varepsilon}_{n+1}^{(j)} = \varepsilon_{n+1}^{(j)} - \overline{{\varepsilon}_{n+1}^{N}}\) for \(n \geq 0\). We define the variance matrices \(\Sigma_K^{y, j, N}\) and \(\Sigma_L^{z, N}\) by replacing \(Y_{n}\) and \(Z_{n}\) in 90 with \(Y_{n}^{K, j, N}\) and \(Z_{n}^{L, N}\), and define \(\Sigma_K^{y, N} = \frac{1}{N} \sum_{j=1}^N \Sigma_K^{y, j, N}\). Similar modifications are applied to the initial perturbation variances \(\Sigma_{Y_{0}}^N = \frac{1}{N} \sum_{j=1}^N \EE[ Y_{0}^{j, N} (Y_{0}^{j,N})^\top ]\) and \(\Sigma_{Z_{0}}^N = \EE[ Z_{0}^N (Z_{0}^N)^\top ]\), and for the noise variances \(\Sigma^{1, N} = \frac{1}{N} \sum_{j=1}^N \EE\big[ \tilde{\varepsilon}_{n}^{(j)} \big( \tilde{\varepsilon}_{n}^{(j)} \big)^\top \big]\) and \(\Sigma^{0, N} = \EE\big[ \big( \varepsilon_{n}^0 + \overline{{\varepsilon}_{n}^{N}} \big) \big( \varepsilon_{n}^0 + \overline{{\varepsilon}_{n}^{N}} \big)^\top \big]\). For the auxiliary costs in the \(N-\)agent problem, we adjust equation 89 with \((\bY^{K, j, N})_{j\in [N]}\) and \(\bZ^{L, N}\) so that \[\begin{align} C_y^{N}(K) &= \frac{1}{N} \sum \nolimits_{j \in [N]} \EE \big[ \sum \nolimits_{n\geq 0} \gamma^{n} (Y_{n}^{K, j, N})^\top (Q + K^\top R K ) Y_{n}^{K, j, N} \big] = \langle Q + K^\top R K, \Sigma_{K}^{y, N} \rangle_{tr} \\ C_z^{N}(L) &= \EE \big[ \sum \nolimits_{n\geq 0} \gamma^{n} (Z_{n}^{L, N})^\top ( \tilde{Q} + L^\top \tilde{R} L ) Z_{n}^{L, N} \big] = \langle \tilde{Q} + L^\top \tilde{R} L, \Sigma_{L}^{z, N} \rangle_{tr}. \end{align}\] For \(\theta = (K, L) \in \Theta\), \(C^N(\theta) = C_y^N(K) + C_z^N(L)\). Because all agents adopt the same control parameter \(\theta\), the \(N\)-agent auxiliary costs \(C_y^N(K)\) and \(C_z^{N}(L)\) can then be expressed with solution matrices \(P_K^y\) and \(P_L^z\) to the DLEs 105 : \[\label{eq:N95agent95Cy95Cz95expression95PK95PL} C_y^{N}(K) = \langle P_K^y, \, \Sigma_{Y_{0}}^N + \frac{\gamma}{1 - \gamma} \Sigma^{1, N} \rangle_{tr}, \quad C_z^N(L) = \langle P_L^z, \, \Sigma_{Z_{0}}^N + \frac{\gamma}{1 - \gamma} \Sigma^{0, N} \rangle_{tr}.\tag{137}\] The truncated versions at horizon \(T\) for the variances and the costs \(\Sigma_K^{y, N, T}\), \(\Sigma_{L}^{z, N, T}\), \(C^{N, T}_y(K)\), \(C^{N, T}_z(L)\), \(C^{N, T}(\theta)\), and the sampled truncated costs with population simulator \(\tilde{C}_y^{N, T}(K),\) \(\tilde{C}_z^{N, T}(L),\) \(\tilde{C}^{N, T}(\theta)\) are defined accordingly.

By definition of the matrices \(\Sigma^{1,N}\), \(\Sigma^{0, N}\), \(\Sigma_{Y_{0}}^{N}\) and \(\Sigma_{Z_{0}}^N\), we have \[\begin{align} \Sigma^{1, N} & = (1 - \tfrac{1}{N} ) \Sigma^1, & \Sigma^{0, N} & = \Sigma^0 + \tfrac{1}{N}\Sigma^1, \tag{138} \\ \Sigma_{Y_{0}}^{N} & = ( 1 - \tfrac{1}{N})\Sigma_{Y_{0}}, & \Sigma_{Z_{0}}^{N} & = \Sigma_{Z_{0}} + \tfrac{1}{N}\Sigma_{Y_{0}}. \tag{139} \end{align}\] Moreover, because the sub-Gaussian norms of \(Y_{0}\), \(Z_{0}\) and the noise processes are bounded by \(C_{init,noise}\), we have for all \(n \geq 0\) and \(j \in [N]\), \[\label{eq:bound95phi295y095z095N95agent} \max \Big\{ \| Y_{0}^{j, N} \|_{\psi_2}, \| Z_{0}^{N} \|_{\psi_2}, \| \tilde{\varepsilon}_{n+1}^{(j)} \|_{\psi_2}, \big\| \varepsilon_{n+1}^0 + \overline{{\varepsilon}_{n+1}^{N}} \big\|_{\psi_2} \Big\} \leq 2 C_{init, noise}\tag{140}\] We define the minimum eigenvalues associated with the initial and noise variance matrices: \[\lambda_y^{1,N} = \lambda_{min}\big( \Sigma_{Y_{0}}^N + \gamma \Sigma^{1, N } / (1 - \gamma) \big), \quad \lambda_z^{0,N} = \lambda_{min}\big( \Sigma_{Z_{0}}^N + \gamma \Sigma^{0, N} / ( 1- \gamma) \big).\]

3.2.4.2 Proof of Proposition [pr:approx95with95social95cost95no95MF95cost] and [pr:modelfree95population95gradient95approx]

Proof of Proposition [pr:approx95with95social95cost95no95MF95cost]

We show in the following that the output of Algorithm 15, \(\tilde{\nabla}^{T, N, M, \tau}\), is close enough to the exact PG \(\nabla C\). We use the same arguments as in Section 3.2.3. In the following, we specify only the differences in the corresponding proof arguments when we use a social cost instead of a mean field cost. To begin with, we define a few bounds for the parameters \((T, N, M)\). \[\begin{align} & \phi_{agent, size, pop}(\varepsilon, \tau, C_0) = \big( \ell d / (\varepsilon \tau) \big) \phi_{social-cost, factor}(C_0) \nonumber \\ & = \frac{2 \ell d^2 C_{init, noise}^2 C_0}{\tau \varepsilon ( 1- \gamma)} \Big( \frac{1}{\lambda_y^1} + \frac{1}{\lambda_z^0} \Big) \tag{141} \\ & \phi_{sample,size, pop}(\varepsilon, \tau, T, C_0, \delta_{sample}) = \phi_{sample, size}(\varepsilon / 4, \tau, T, C_0, \delta_{sample}) \nonumber \\ & = \log \Big( \frac{4T (\ell + d)}{\delta_{sample}} \Big) \Big( \frac{32 \ell d T}{\varepsilon \tau^2} ( h_{\psi_2}(\tau, C_0) ) + \frac{1}{4} \Big)^2 \tag{142} \\ & \phi_{trunc, T, pop}(\varepsilon, \tau, C_0, N) = \phi_{trunc, T}(\varepsilon / (1 + 1/N), \tau, C_0) \nonumber \\ & = \Big( \frac{1}{\log(1 / \gamma_{pert})} \Big)^2 \Big[ \log \Big( \frac{4 \ell d C_{init, noise}^2 h_T(C_0)}{\varepsilon \tau ( 1- \gamma_{pert})^2} (1 + \tfrac{1}{N}) \Big) + 1 \Big]^2 \tag{143} \end{align}\] where \(\phi_{social-cost, factor}(C_0)\) is defined in Proposition [pr:approx95with95social95cost95no95MF95cost] , and \(\phi_{sample, size}\), \(\phi_{trunc, T}\) are defined in equations 131 and 128 .

Proof of Proposition [pr:modelfree95population95gradient95approx]

3.2.5 Numerical Results↩︎

In this section, we provide numerical illustrations of the RL algorithms developed for the LQMFC problem. We compare the performance of the exact PG (model-based) with the model-free approach using both the MKV simulator and the finite-population simulator.

The code is available at: https://github.com/mlauriere/mfrl-tutorial-code.

3.2.5.1 Experimental Setup

We consider a scalar LQMFC model where the dynamics and costs are symmetric across the agent population. The representative agent’s objective is to minimize the socially optimal cost by learning the parameters \(\theta = (K, L)\), which determine the decentralized control \(\ctrl_{n} = -K(X_{n} - \bar{X}_{n}) - L\bar{X}_{n}\).

The model parameters and hyperparameters used for the simulation are summarized in Table 4. The dynamics are characterized by coefficients \(A, \bar{A}, B, \bar{B}\), and the cost function by \(Q, \bar{Q}, R, \bar{R}\). The noise in the dynamics is modeled as zero-mean Gaussian with variances \(\sigma_\epsilon^2\) and \(\sigma_{\epsilon_0}^2\) for the idiosyncratic and common noise, respectively.

We note that for the numerical illustrations, we use the Adam optimizer to adapt the learning rate during training. While the theoretical analysis presented in the previous sections relies on a fixed learning rate for simplicity and to derive clear convergence rates, Adam is employed here for its practical efficiency and superior convergence properties in stochastic environments. This choice serves to illustrate the robust performance of the developed algorithms in a more standard machine learning implementation.

We study numerically the three primary configurations analyzed above:

1. Exact PG: Uses model-based exact gradients as a baseline.

2. Adam MF-MKV: Model-free PG using the idealized MKV simulator.

3. Adam MF-Pop (\(N\)): Model-free policy gradient using the finite-population simulator with \(N\) agents. We compare \(N=1\) and \(N=10\) to study the impact of population size on the estimation error.

3.2.5.2 Discussion of Results

The results are aggregated over several independent seeds (see Table 4) and smoothed using a running average. In the figures below, the solid lines represent the mean value across seeds, while the shaded areas represent the standard deviation.

Figure 16 illustrates the convergence of the social cost. The left plot shows the cumulative cost for each method approaching the theoretical optimal value (the black dashed line). The right plot displays the logarithm of the distance to the optimal cost, highlighting the efficiency of the gradient estimation. We observe that while the MKV simulator closely tracks the exact gradient performance, the population simulator exhibits higher variance and a larger gap, which decreases significantly as \(N\) increases from 1 to 10.

Figure 17 provides the convergence of the policy parameters. The left plot shows the Euclidean distance between the learned parameters \((K, L)\) and the optimal values \((K^*, L^*)\). The middle and right plots track the individual evolution of \(K\) and \(L\), respectively. Consistent with the cost results, the exact-gradient, MKV, and large-population simulators move toward the optimal parameters. When \(N=1\), the parameter \(K\) is not identifiable. Indeed, for \(N=1\), \(X_n^1-\bar X_n^N=0\), so the centered feedback term involving \(K\) is identically zero.

3.2.6 Notes and Complements↩︎

The convergence proof in this chapter builds upon the approach proposed by [108] for policy-gradient methods in classical LQR, and extends it by including stochastic dynamics and mean-field interactions.

PG methods have also been applied to linear-quadratic mean field type games, i.e., non-cooperative games between central planners solving MFC problems, see [158]–[160]. Further related work on LQ mean field learning includes [81], [161], [162].

Other RL algorithms for LQ MFC problems include the two-timescale method of [131] and the policy iteration method studied in [163]. [164], [165] proposed an approach based on integral RL for the case of multiplicative noise and robust mean field social control.

References↩︎

---

1. Department of Operations Research and Financial Engineering & Program in Applied and Computational Mathematics, Princeton NJ 08544, USA,
   email: rcarmona@princeton.edu↩︎

2. Shanghai Frontiers Science Center of Artificial Intelligence and Deep Learning; NYU-ECNU Institute of Mathematical Sciences at NYU Shanghai; NYU Shanghai, 567 West Yangsi Road, Shanghai, 200126, People’s Republic of China
   email: mathieu.lauriere@nyu.edu↩︎

3. See https://www.imsi.institute/activities/theoretical-advances-in-reinforcement-learning-and-control/foundations-of-multi-agent-and-mean-field-reinforcement-learning/↩︎

4. The interactions are through the joint state-action distribution, which is sometimes referred to as “extended” MFC or MFC of controls.↩︎