ACTOR-CRITIC PROVABLY FINDS NASH EQUILIBRIA OF LINEAR-QUADRATIC MEAN-FIELD GAMES
Year: 2020

Abstract:
We study discrete-time mean-field Markov games with infinite numbers of agents where each agent aims to minimize its ergodic cost. We consider the setting where the agents have identical linear state transitions and quadratic cost functions, while the aggregated effect of the agents is captured by the population mean of their states, namely, the mean-field state. For such a game, based on the Nash certainty equivalence principle, we provide sufficient conditions for the existence and uniqueness of its Nash equilibrium. Moreover, to find the Nash equilibrium, we propose a mean-field actor-critic algorithm with linear function approximation, which does not require knowing the model of dynamics. Specifically, at each iteration of our algorithm, we use the single-agent actor-critic algorithm to approximately obtain the optimal policy of the each agent given the current meanfield state, and then update the mean-field state. In particular, we prove that our algorithm converges to the Nash equilibrium at a linear rate. To the best of our knowledge, this is the first success of applying model-free reinforcement learning with function approximation to discrete-time mean-field Markov games with provable non-asymptotic global convergence guarantees.

heading: 1 INTRODUCTION
text: In reinforcement learning (RL) (Sutton and Barto, 2018), an agent learns to make decisions that minimize its expected total cost through sequential interactions with the environment. Multi-agent reinforcement learning (MARL) (Shoham et al., 2003; 2007; Busoniu et al., 2008) aims to extend RL to sequential decision-making problems involving multiple agents. In a non-cooperative game, we are interested in the Nash equilibrium (Nash, 1951), which is a joint policy of all the agents such that each agent cannot decrease its expected total cost by unilaterally deviating from its Nash policy. The Nash equilibrium plays a critical role in understanding the social dynamics of self-interested agents (Ash, 2000; Axtell, 2002) and constructing the optimal policy of a particular agent via fictitious selfplay (Bowling and Veloso, 2000; Ganzfried and Sandholm, 2009). With the recent development in deep learning (LeCun et al., 2015), MARL with function approximation achieves tremendous empirical successes in applications, including Go (Silver et al., 2016; 2017), Poker (Heinrich and Silver, 2016; Moravčı́k et al., 2017), Star Craft (Vinyals et al., 2019), Dota (OpenAI, 2018), autonomous driving (Shalev-Shwartz et al., 2016), multi-robotic systems (Yang and Gu, 2004), and solving social dilemmas (de Cote et al., 2006; Leibo et al., 2017; Hughes et al., 2018). However, since the capacity of the joint state and action spaces grows exponentially in the number of agents, such MARL approaches become computationally intractable when the number of agents is large, which is common in real-world applications (Sandholm, 2010; Calderone, 2017; Wang et al., 2017a).
Mean-field game is proposed by Huang et al. (2003; 2006); Lasry and Lions (2006a;b; 2007) with the idea of utilizing mean-field approximation to model the strategic interactions within a large population. In a mean-field game, each agent has the same cost function and state transition, which depend on the other agents only through their aggregated effect. As a result, the optimal policy of
each agent depends solely on its own state and the aggregated effect of the population, and such an optimal policy is symmetric across all the agents. Moreover, if the aggregated effect of the population corresponds to the Nash equilibrium, then the optimal policy of each agent jointly constitutes a Nash equilibrium. Although such a Nash equilibrium corresponds to an infinite number of agents, it well approximates the Nash equilibrium for a sufficiently large number of agents (Bensoussan et al., 2016). Also, as the aggregated effect of the population abstracts away the strategic interactions between individual agents, it circumvents the computational intractability of the MARL approaches that do not exploit symmetry.
However, most existing work on mean-field games focuses on characterizing the existence and uniqueness of the Nash equilibrium rather than designing provably efficient algorithms. In particular, most existing work considers the continuous-time setting, which requires solving a pair of Hamilton-Jacobi-Bellman (HJB) and Fokker-Planck (FP) equations, whereas the discrete-time setting is more common in practice, e.g., in the aforementioned applications. Moreover, most existing approaches, including the ones based on solving the HJB and FP equations, require knowing the model of dynamics (Bardi and Priuli, 2014), or having the access to a simulator, which generates the next state given any state-action pair and aggregated effect of the population (Guo et al., 2019), which is often unavailable in practice.
To address these challenges, we develop an efficient model-free RL approach to mean-field game, which provably attains the Nash equilibrium. In particular, we focus on discrete-time mean-field games with linear state transitions and quadratic cost functions, where the aggregated effect of the population is quantified by the mean-field state. Such games capture the fundamental difficulties of general mean-field games and well approximates a variety of real-world systems such as power grids (Minciardi and Sacile, 2011), swarm robots (Fang, 2014; Araki et al., 2017; Doerr et al., 2018), and financial systems (Zhou and Li, 2000; Huang and Li, 2018). In detail, based on the Nash certainty equivalence (NCE) principle (Huang et al., 2006; 2007), we propose a mean-field actor-critic algorithm which, at each iteration, given the mean-field state µ, approximately attains the optimal policy π∗µ of each agent, and then updates the mean-field state µ assuming that all the agents follow π∗µ. We parametrize the actor and critic by linear and quadratic functions, respectively, and prove that such a parameterization encompasses the optimal policy of each agent. Specifically, we update the actor parameter using policy gradient (Sutton et al., 2000) and natural policy gradient (Kakade, 2002; Peters and Schaal, 2008; Bhatnagar et al., 2009) and update the critic parameter using primal-dual gradient temporal difference (Sutton et al., 2009a;b). In particular, we prove that given the mean-field state µ, the sequence of policies generated by the actor converges linearly to the optimal policy π∗µ. Moreover, when alternatingly update the policy and mean-field state, we prove that the sequence of policies and its corresponding sequence of mean-field states converge to the unique Nash equilibrium at a linear rate. Our approach can be interpreted from both “passive” and “active” perspectives: (i) Assuming that each self-interested agent employs the single-agent actorcritic algorithm, the policy of each agent converges to the unique Nash policy, which characterizes the social dynamics of a large population of model-free RL agents. (ii) For a particular agent, our approach serves as a fictitious self-play method for it to find its Nash policy, assuming the other agents give their best responses. To the best of our knowledge, our work establishes the first efficient model-free RL approach with function approximation that provably attains the Nash equilibrium of a discrete-time mean-field game. As a byproduct, we also show that the sequence of policies generated by the single-agent actor-critic algorithm converges at a linear rate to the optimal policy of a linear-quadratic regulator (LQR) problem in the presence of drift, which may be of independent interest.
Related Work. Mean-field game is first introduced in Huang et al. (2003; 2006); Lasry and Lions (2006a;b; 2007). In the last decade, there is growing interest in understanding continuous-time mean-field games. See, e.g., Guéant et al. (2011); Bensoussan et al. (2013); Gomes et al. (2014); Carmona and Delarue (2013; 2018) and the references therein. Due to their simple structures, continuous-time linear-quadratic mean-field games are extensively studied under various model assumptions. See Li and Zhang (2008); Bardi (2011); Wang and Zhang (2012); Bardi and Priuli (2014); Huang et al. (2016a;b); Bensoussan et al. (2016; 2017); Caines and Kizilkale (2017); Huang and Huang (2017); Moon and Başar (2018); Huang and Zhou (2019) for examples of this line of work. Meanwhile, the literature on discrete-time linear-quadratic mean-field games remains relatively scarce. Most of this line of work focuses on characterizing the existence of a Nash equilibrium and the behavior of such a Nash equilibrium when the number of agents goes to infinity
(Gomes et al., 2010; Tembine and Huang, 2011; Moon and Başar, 2014; Biswas, 2015; Saldi et al., 2018a;b; 2019). See also Yang et al. (2018a), which applies maximum entropy inverse RL (Ziebart et al., 2008) to infer the cost function and social dynamics of discrete-time mean-field games with finite state and action spaces. Our work is most related to Guo et al. (2019), where they propose a mean-field Q-learning algorithm (Watkins and Dayan, 1992) for discrete-time mean-field games with finite state and action spaces. Such an algorithm requires the access to a simulator, which, given any state-action pair and mean-field state, outputs the next state. In contrast, both our state and action spaces are infinite, and we do not require such a simulator but only observations of trajectories under given mean-field state. Correspondingly, we study the mean-field actor-critic algorithm with linear function approximation, whereas their algorithm is tailored to the tabular setting. Also, our work is closely related to Mguni et al. (2018), which focuses on a more restrictive setting where the state transition does not involve the mean-field state. In such a setting, mean-field games are potential games, which is, however, not true in more general settings (Li et al., 2017; Briani and Cardaliaguet, 2018). In comparison, we allow the state transition to depend on the mean-field state. Meanwhile, they propose a fictitious self-play method based on the single-agent actor-critic algorithm and establishes its asymptotic convergence. However, their proof of convergence relies on the assumption that the single-agent actor-critic algorithm converges to the optimal policy, which is unverified therein. In addition, our work is related to Jayakumar and Aditya (2019), where the proposed algorithm is only shown to converge asymptotically to a stationary point of the mean-field game.
Our work also extends the line of work on finding the Nash equilibria of Markov games using MARL. Due to the computational intractability introduced by the large number of agents, such a line of work focuses on finite-agent Markov games (Littman, 1994; 2001; Hu and Wellman, 1998; Bowling, 2001; Lagoudakis and Parr, 2002; Hu and Wellman, 2003; Conitzer and Sandholm, 2007; Perolat et al., 2015; Pérolat et al., 2016b;a; 2018; Wei et al., 2017; Zhang et al., 2018; Zou et al., 2019; Casgrain et al., 2019). See also Shoham et al. (2003; 2007); Busoniu et al. (2008); Li (2018) for detailed surveys. Our work is related to Yang et al. (2018b), where they combine the mean-field approximation of actions (rather than states) and Nash Q-learning (Hu and Wellman, 2003) to study general-sum Markov games with a large number of agents. However, the Nash Q-learning algorithm is only applicable to finite state and action spaces, and its convergence is established under rather strong assumptions. Also, when the number of agents goes to infinity, their approach yields a variant of tabular Q-learning, which is different from our mean-field actor-critic algorithm.
For policy optimization, based on the policy gradient theorem, Sutton et al. (2000); Konda and Tsitsiklis (2000) propose the actor-critic algorithm, which is later generalized to the natural actor-critic algorithm (Peters and Schaal, 2008; Bhatnagar et al., 2009). Most existing results on the convergence of actor-critic algorithms are based on stochastic approximation using ordinary differential equations (Bhatnagar et al., 2009; Castro and Meir, 2010; Konda and Tsitsiklis, 2000; Maei, 2018), which are asymptotic in nature. For policy evaluation, the convergence of primal-dual gradient temporal difference is studied in Liu et al. (2015); Du et al. (2017); Wang et al. (2017b); Yu (2017); Wai et al. (2018). However, this line of work assumes that the feature mapping is bounded, which is not the case in our setting. Thus, the existing convergence results are not applicable to analyzing the critic update in our setting. To handle the unbounded feature mapping, we utilize a truncation argument, which requires more delicate analysis.
Finally, our work extends the line of work that studies model-free RL for LQR. For example, Bradtke (1993); Bradtke et al. (1994) show that policy iteration converges to the optimal policy, Tu and Recht (2017); Dean et al. (2017) study the sample complexity of least-squares temporal-difference for policy evaluation. More recently, Fazel et al. (2018); Malik et al. (2018); Tu and Recht (2018) show that the policy gradient algorithm converges at a linear rate to the optimal policy. See as also Hardt et al. (2016); Dean et al. (2018) for more in this line of work. Our work is also closely related to Yang et al. (2019), where they show that the sequence of policies generated by the natural actorcritic algorithm enjoys a linear rate of convergence to the optimal policy. Compared with this work, when fixing the mean-field state, we use the actor-critic algorithm to study LQR in the presence of drift, which introduces significant difficulties in the analysis. As we show in §3, the drift causes the optimal policy to have an additional intercept, which makes the state- and action-value functions more complicated.
Notations. We denote by ‖M‖∗ the spectral norm, ρ(M) the spectral radius, σmin(M) the minimum singular value, and σmax(M) the maximum singular value of a matrix M . We use ‖α‖2 to represent the `2-norm of a vector α, and (α) j i to denote the sub-vector (αi, αi+1, . . . , αj)
>, where αk is the k-th entry of the vector α. For scalars a1, . . . , an, we denote by poly(a1, . . . , an) the polynomial of a1, . . . , an, and this polynomial may vary from line to line. We use [n] to denote the set {1, 2, . . . , n} for any n ∈ N.

heading: 2 LINEAR-QUADRATIC MEAN-FIELD GAME
text: A linear-quadratic mean-field Na-player game involves Na ∈ N agents, whose state transitions are given by
xit+1 = Ax i t +Bu i t +Axt + d i + ωit, ∀t ≥ 0, i ∈ [Na].
HereA ∈ Rm×m,B ∈ Rm×k, andA ∈ Rm×m are matrices, xit ∈ Rm and uit ∈ Rk are the state and action vectors of agent i, respectively, the vector di ∈ Rm is a drift term, ωit ∈ Rm is an independent random noise term following the Gaussian distribution N (0,Ψω), and xt = 1/Na · ∑Na j=1 x j t is the mean-field state. The agents are coupled through the mean-field state xt. In the linear-quadratic mean-field Na-player game, the cost of agent i ∈ [Na] at time t ≥ 0 is given by
cit = (x i t) >Qxit + (u i t) >Ruit + x > t Qxt,
where Q ∈ Rm×m, R ∈ Rk×k, and Q ∈ Rm×m are matrices, and uit is generated by πi, i.e., the policy of agent i. To measure the performance of agent i following its policy πi under the influence of the other agents, we define the expected total cost of agent i as
J i(π1, π2, . . . , πNa) = lim T→∞ E
( 1
T T∑ t=0 cit
) .
We are interested in finding a Nash equilibrium (π1, π2, . . . , πNa), which is defined by
J i(π1, . . . , πi−1, πi, πi+1, . . . , πNa) ≤ J i(π1, . . . , πi−1, π̃i, πi+1, . . . , πNa), ∀π̃i, i ∈ [Na].
That is, agent i cannot further decrease its expected total cost by unilaterally deviating from its Nash policy.
For the simplicity of discussion, we assume that the drift term di is identical for each agent. We consider taking the infinite-population limit Na →∞, where each agent has an infinitesimal contribution to the dynamics of the system. Thus, the joint policy of all the agents except agent i can be modeled as a mean-field policy π†, and all the agents following such a mean-field policy π† generate the mean-field state Ex†t , where {x † t}t≥0 is generated following the policy π†. By the symmetry of the agents in terms of their state transitions and cost functions, we focus on a fixed agent and drop the superscript i hereafter.
Before we formally present the formulation of linear-quadratic mean-field games, we first introduce the following mean-field LQR (MF-LQR) problem, which aims to find an optimal policy for the fixed agent given the mean-field policy π†. Problem 2.1 (MF-LQR). Given the mean-field policy π†, we consider the following formulation,
xt+1 = Axt +But +AEx†t + d+ ωt, c(xt, ut) = x > t Qxt + u > t Rut + (Ex † t) >Q(Ex†t),
J(π, π†) = lim T→∞ E
[ 1
T T∑ t=0 c(xt, ut)
] ,
where xt ∈ Rm is the state vector, ut ∈ Rk is the action vector generated by the policy π, {x†t}t≥0 is the trajectory generated by the policy π†, ωt ∈ Rm is an independent random noise term following the Gaussian distribution N (0,Ψω), and d ∈ Rm is a drift term. Here the expectation Ex†t is taken across all the agents. We aim to find π∗ such that J(π∗, π†) = infπ∈Π J(π, π†).
Note that a controllable linear system using linear quadratic optimal control is always stable. Further combining the fact that our linear closed-loop dynamics in Problem 2.1 is driven by the Gaussian noise term ωt, we know that the Markov chain of states generated by the policy π† admits a stationary distribution and converges to this stationary distribution. This implies that the mean-field state Ex†t converges to a constant vector µ† as t→∞, which serves as a time-invariant mean-field state. As we consider the ergodic setting, it then suffices to study Problem 2.1 with t sufficiently large. Therefore, the influence of the mean-field policy π† is captured by the mean-field state µ†. By re-formulating Problem 2.1, with slight abuse of notations, we obtain the following drifted-LQR (D-LQR). Problem 2.2 (D-LQR). Given a mean-field state µ ∈ Rm, we consider the following formulation,
xt+1 = Axt +But +Aµ+ d+ ωt, cµ(xt, ut) = x > t Qxt + u > t Rut + µ >Qµ,
Jµ(π) = lim T→∞
E
[ 1
T T∑ t=0 cµ(xt, ut)
] ,
where xt ∈ Rm is the state vector, ut ∈ Rk is the action vector generated by the policy π, ωt ∈ Rm is an independent random noise term following the Gaussian distribution N (0,Ψω), and d ∈ Rm is a drift term. We aim to find an optimal policy π∗µ such that Jµ(π ∗ µ) = infπ∈Π Jµ(π).
Compared with the most studied LQR problem (Lewis et al., 2012), both the state transition and the cost function in Problem 2.2 have drift terms, which act as the mean-field “force” that drives the states away from zero. Such a mean-field “force” introduces additional challenges when solving Problem 2.2 in the model-free setting (see §3.3 for details). On the other hand, the unique optimal policy π∗µ of Problem 2.2 admits a linear form π ∗ µ(xt) = −Kπ∗µxt + bπ∗µ under certain regularity conditions (Anderson and Moore, 2007), where the matrix Kπ∗µ ∈ R k×m and the vector bπ∗µ ∈ R k are the parameters of π∗µ. Motivated by such a linear form of the optimal policy, we define the class of linear-Gaussian policies as
Π = {π(x) = −Kπx+ bπ + σ · η : Kπ ∈ Rk×m, bπ ∈ Rk}, (2.1)
where σ ∈ R and the standard Gaussian noise term η ∈ Rk is included to encourage exploration. To solve Problem 2.2, it suffices to find the optimal policy π∗µ within Π. We define Λ1(µ) = π ∗ µ as the optimal policy under the mean-field state µ.
Assume that all the agents follow the linear policy π(x) = −Kπx + bπ under the mean-field state µ. By plugging ut = π(xt) into the state transition in Problem 2.2, as t → ∞, we know that these agents generate a new mean-field state µnew such that
µnew = (I −A+BKπ)−1(Bbπ +Aµ+ d).
We define Λ2(µ, π) = µnew as such a new mean-field state.
Now, we are ready to present the following linear-quadratic mean-field game (LQ-MFG). Problem 2.3 (LQ-MFG). We consider the following formulation,
xt+1 = Axt +But +Aµ+ d+ ωt, c(xt, ut) = x > t Qxt + u > t Rut + µ >Qµ,
J(π, µ) = lim T→∞
E
[ 1
T T∑ t=0 c(xt, ut)
] ,
where xt ∈ Rm is the state vector, ut ∈ Rk is the action vector generated by the policy π, µ ∈ Rm is the mean-field state, ωt ∈ Rm is an independent random noise term following the Gaussian distribution N (0,Ψω), and d ∈ Rm is a drift term. We aim to find a pair (µ∗, π∗) such that (i) J(π∗, µ∗) = infπ∈Π J(π, µ
∗); (ii) Ex∗t converges to µ∗ as t → ∞, where {x∗t }t≥0 is the Markov chain of states generated by the policy π∗.
The formulation in Problem 2.3 is studied by Lasry and Lions (2007); Bensoussan et al. (2016); Saldi et al. (2018a;b). We propose a more general formulation in Problem C.2 (see §C of the appendix for
details), where an additional interaction term between the state vector xt and the mean-field state µ is incorporated into the cost function. According to our analysis in §C, up to minor modification, the results in the following sections also carry over to Problem C.2. Therefore, for the sake of simplicity, we focus on Problem 2.3 in the sequel.
In Problem 2.3, condition (i) is equivalent to the optimality of the policy π∗ under the mean-field state µ∗, namely, Λ1(µ∗) = π∗. Meanwhile, condition (ii) is equivalent to the invariance of the mean-field state µ∗ given the policy π∗, namely, Λ2(µ∗, π∗) = µ∗. Such equivalence follows from the NCE principle (Huang et al., 2006; 2007), which also motivates the following definition of the Nash equilibrium pair (Saldi et al., 2018a;b). Definition 2.4 (Nash Equilibrium Pair). The pair (µ∗, π∗) ∈ Rm×Π constitutes a Nash equilibrium pair of Problem 2.3 if it satisfies π∗ = Λ1(µ∗) and µ∗ = Λ2(µ∗, π∗). Here µ∗ is called the Nash mean-field state and π∗ is called the Nash policy.
By Definition 2.4, Problem 2.3 aims to find a Nash equilibrium pair (µ∗, π∗).

heading: 3 MEAN-FIELD ACTOR-CRITIC
text: We first characterize the existence and uniqueness of the Nash equilibrium pair of Problem 2.3 under mild regularity conditions, and then propose a mean-field actor-critic algorithm to obtain such a Nash equilibrium. As a building block of the mean-field actor-critic, we propose the natural actor-critic to solve Problem 2.2.

heading: 3.1 EXISTENCE AND UNIQUENESS OF NASH EQUILIBRIUM PAIR
text: We now establish the existence and uniqueness of the Nash equilibrium pair defined in Definition 2.4. We impose the following regularity conditions.
Assumption 3.1. We assume that the following statements hold:
(i) The algebraic Riccati equation X = A>XA+Q−A>XB(B>XB +R)−1B>XA admits a unique symmetric positive definite solution X∗;
(ii) It holds for L0 = L1L3 + L2 that L0 < 1, where L1 = ∥∥[(I −A)Q−1(I −A)> +BR−1B>]−1A∥∥ 2 · ∥∥[K∗Q−1(I −A)> −R−1B>]∥∥ 2 ,
L2 = [ 1− ρ(A−BK∗) ]−1 · ‖A‖2, L3 = [1− ρ(A−BK∗)]−1 · ‖B‖2. Here K∗ = −(B>X∗B +R)−1B>X∗A.
The first assumption is implied by mild regularity conditions on the matrices A, B, Q, and R, which are (1) the positivity of R; (2) the non-negativity of Q = C>C; (3) the observability of (A,C); (4) the stability of (A,B). See De Souza et al. (1986); Lewis et al. (2012) for more details. The second assumption is standard in the literature (Bensoussan et al., 2016; Saldi et al., 2018b), which ensures the stability of the LQ-MFG. In the following proposition, we show that Problem 2.3 admits a unique Nash equilibrium pair.
Proposition 3.2 (Existence and Uniqueness of Nash Equilibrium Pair). Under Assumption 3.1, the operator Λ(·) = Λ2(·,Λ1(·)) is L0-Lipschitz, where L0 is given in Assumption 3.1. Moreover, there exists a unique Nash equilibrium pair (µ∗, π∗) of Problem 2.3.
Proof. See §E.1 for a detailed proof.

heading: 3.2 MEAN-FIELD ACTOR-CRITIC FOR LQ-MFG
text: The NCE principle motivates a fixed-point approach to solve Problem 2.3, which generates a sequence of policies {πs}s≥0 and mean-field states {µs}s≥0 satisfying the following two properties: (i) Given the mean-field state µs, the policy πs is optimal. (ii) The mean-field state becomes µs+1 as t→∞, if all the agents follow πs under the current mean-field state µs. Here (i) requires solving Problem 2.2 given the mean-field state µs, while (ii) requires simulating the agents following the
Algorithm 1 Mean-Field Actor-Critic for solving LQ-MFG. 1: Input:
• Initial mean-field state µ0 and Initial policy π0 with parameters K0 and b0. • Numbers of iterations S, {Ns}s∈[S], {Hs}s∈[S], {T̃s,n, Ts,n}s∈[S],n∈[Ns], {T̃ bs,h, T bs,h}s∈[S],h∈[Hs].
• Stepsizes {γs}s∈[S], {γbs}s∈[S], {γs,n,t}s∈[S],n∈[Ns],t∈[Ts,n], {γbs,h,t}s∈[S],h∈[Hs],t∈[T bs,h]. 2: for s = 0, 1, 2, . . . , S − 1 do 3: Policy Update: Solve for the optimal policy πs+1 with parametersKs+1 and bs+1 of Problem
2.2 via Algorithm 2 with µs, πs, Ns, Hs, {T̃s,n, Ts,n}n∈[Ns], {T̃ bs,h, T bs,h}h∈[Hs], γs, γbs , {γs,n,t}n∈[Ns],t∈[Ts,n], and {γbs,h,t}h∈[Hs],t∈[T bs,h], which gives the estimated mean-field state µ̂Ks+1,bs+1 .
4: Mean-Field State Update: Update the mean-field state via µs+1 ← µ̂Ks+1,bs+1 . 5: end for 6: Output: Pair (πS , µS).
policy πs given the current mean-field µs. Based on such properties, we propose the mean-field actor-critic in Algorithm 1.
Algorithm 1 requires solving Problem 2.2 at each iteration to obtain πs = Λ1(µs) and µs+1 = Λ2(µs, πs). To this end, we introduce the natural actor-critic in §3.3 that solves Problem 2.2.

heading: 3.3 NATURAL ACTOR-CRITIC FOR D-LQR
text: Now we focus on solving Problem 2.2 for a fixed mean-field state µ, we thus drop the subscript µ hereafter. With slight abuse of notations, we write πK,b(x) = −Kx + b + σ · η to emphasize the dependence on K and b, and J(K, b) = J(πK,b) consequently. Now, we propose the natural actor-critic to solve Problem 2.2.
For any policy πK,b ∈ Π, by the state transition in Problem 2.2, we have
xt+1 = (A−BK)xt + (Bb+Aµ+ d) + t, t ∼ N (0,Ψ ), (3.1)
where Ψ = σBB> + Ψω . It is known that if ρ(A − BK) < 1, then the Markov chain {xt}t≥0 induced by (3.1) has a unique stationary distribution N (µK,b,ΦK) (Anderson and Moore, 2007), where the mean-field state µK,b and the covariance ΦK satisfy that
µK,b = (I −A+BK)−1(Bb+Aµ+ d), (3.2) ΦK = (A−BK)ΦK(A−BK)> + Ψ . (3.3)
Meanwhile, the Bellman equation for Problem 2.2 takes the following form
PK = (Q+K >RK) + (A−BK)>PK(A−BK). (3.4)
Then by calculation (see Proposition B.2 in §B.1 of the appendix for details), it holds that the expected total cost J(K, b) is decomposed as
J(K, b) = J1(K) + J2(K, b) + σ 2 · tr(R) + µ>Qµ, (3.5)
where J1(K) and J2(K, b) are defined as J1(K) = tr [ (Q+K>RK)ΦK ] = tr(PKΨ ),
J2(K, b) = ( µK,b b )>( Q+K>RK −K>R −RK R )( µK,b b ) . (3.6)
Here J1(K) is the expected total cost in the most studied LQR problems (Yang et al., 2019; Fazel et al., 2018), where the state transition does not have drift terms. Meanwhile, J2(K, b) corresponds to the expected cost induced by the drift terms. The following two propositions characterize the properties of J2(K, b).
First, we show that J2(K, b) is strongly convex in b.
Proposition 3.3. Given any K, the function J2(K, b) is νK-strongly convex in b. Here νK = σmin(Y > 1,KY1,K + Y > 2,KY2,K), where Y1,K = R
1/2K(I − A + BK)−1B − R1/2 and Y2,K = Q1/2(I − A + BK)−1B. Also, J2(K, b) has ιK-Lipschitz continuous gradient in b, where ιK is upper bounded as ιK ≤ [1− ρ(A−BK)]−2 · (‖B‖22 · ‖K‖22 · ‖R‖2 + ‖B‖22 · ‖Q‖2).
Proof. See §E.4 for a detailed proof.
Second, we show that minb J2(K, b) is independent of K. Proposition 3.4. We define bK = argminb J2(K, b), where J2(K, b) is defined in (3.6). It holds that
bK = [ KQ−1(I −A)> −R−1B> ] · [ (I −A)Q−1(I −A)> +BR−1B> ]−1 · (Aµ+ d). Moreover, J2(K, bK) takes the form of
J2(K, b K) = (Aµ+ d)> [ (I −A)Q−1(I −A)> +BR−1B> ]−1 · (Aµ+ d), which is independent of K.
Proof. See §E.2 for a detailed proof.
Since minb J2(K, b) is independent of K by Proposition 3.4, it holds that the optimal K∗ is the same as argminK J1(K). This motivates us to minimize J(K, b) by first updating K following the gradient direction ∇KJ1(K) to the optimal K∗, then updating b following the gradient direction ∇bJ2(K∗, b). We now design our algorithm based on this idea. We define ΥK , pK,b, and qK,b as
ΥK =
( Q+A>PKA A
>PKB B>PKA R+B >PKB
) = ( Υ11K Υ 12 K
Υ21K Υ 22 K
) ,
pK,b = A >[PK · (Aµ+ d) + fK,b], qK,b = B>[PK · (Aµ+ d) + fK,b], (3.7)
where fK,b = (I − A + BK)−>[(A − BK)>PK(Bb + Aµ + d) −K>Rb]. By calculation (see Proposition B.3 in §B.1 of the appendix for details), the gradients of J1(K) and J2(K, b) take the forms of
∇KJ1(K) = 2(Υ22KK −Υ21K ) · ΦK , ∇bJ2(K, b) = Υ22K (−KµK,b + b) + Υ21KµK,b + qK,b.
Our algorithm follows the natural actor-critic method (Bhatnagar et al., 2009) and actor-critic method (Konda and Tsitsiklis, 2000). Specifically, (i) To obtain the optimal K∗, in the critic update step, we estimate the matrix ΥK by Υ̂K via a policy evaluation algorithm, e.g., Algorithm 3 or Algorithm 4 (see §B.2 and §B.3 of the appendix for details); in the actor update step, we update K via K ← K−γ · (Υ̂22KK− Υ̂21K ), where the term Υ̂22KK− Υ̂21K is the estimated natural gradient. (ii) To obtain the optimal b∗ given K∗, in the critic update step, we estimate ΥK∗ , qK∗,b, and µK∗,b by Υ̂K∗ , q̂K∗,b, and µ̂K∗,b via a policy evaluation algorithm; In the actor update step, we update b via b← b− γ · ∇̂bJ2(K∗, b), where ∇̂bJ2(K∗, b) = Υ̂22K∗(−K∗µ̂K∗,b + b) + Υ̂21K∗ µ̂K∗,b + q̂K∗,b is the estimated gradient. Combining the above procedure, we obtain the natural actor-critic for Problem 2.2, which is stated in Algorithm 2.
One may want to apply gradient method to J(K, b) directly in the joint space of K and b. However, the gradient dominance property of J1(K) in the most studied LQR problem (Yang et al., 2019) no longer holds for J(K, b). Therefore, the convergence of the gradient method to J(K, b) is not guaranteed in our problem.

heading: 4 GLOBAL CONVERGENCE RESULTS
text: The following theorem establishes the rate of convergence of Algorithm 1 to the Nash equilibrium pair (µ∗, π∗) of Problem 2.3.
Algorithm 2 Natural Actor-Critic Algorithm for D-LQR. 1: Input:
• Mean-field state µ and initial policy πK0,b0 . • Numbers of iterations N , H , {T̃n, Tn}n∈[N ], {T̃ bh, T bh}h∈[H]. • Stepsizes γ, γb, {γn,t}n∈[N ],t∈[Tn], {γbh,t}h∈[H],t∈[T bh].
2: for n = 0, 1, 2, . . . , N − 1 do 3: Critic Update: Compute Υ̂Kn via Algorithm 3 with πKn,b0 , µ, T̃n,Tn, {γn,t}t∈[Tn], K0, and b0 as inputs. 4: Actor Update: Update the parameter via
Kn+1 ← Kn − γ · (Υ̂22KnKn − Υ̂ 21 Kn).
5: end for 6: for h = 0, 1, 2, . . . ,H − 1 do 7: Critic Update: Compute µ̂KN ,bh , Υ̂KN , q̂KN ,bh via Algorithm 3 with πKN ,bh , µ, T̃ bh, T bh, {γbh,t}t∈[T bh], K0, and b0. 8: Actor Update: Update the parameter via
bh+1 ← bh − γb · [ Υ̂22KN (−KN µ̂K,bh + bh) + Υ̂ 21 KN µ̂KN ,bh + q̂KN ,bh ] .
9: end for 10: Output: Policy πK,b = πKN ,bH , estimated mean-field state µ̂K,b = µ̂KN ,bH .
Theorem 4.1 (Convergence of Algorithm 1). For a sufficiently small tolerance ε > 0, we set the number of iterations S in Algorithm 1 such that
S > log ( ‖µ0 − µ∗‖2 · ε−1 ) log(1/L0) . (4.1)
For any s ∈ [S], we define εs = min {[ 1− ρ(A−BK∗) ]4(‖B‖2 + ‖A‖2)−4(‖µs‖−22 + ‖d‖−22 ) · σmin(Ψ ) · σmin(R) · ε2,
νK∗ · [ 1− ρ(A−BK∗) ]4 · ‖B‖−22 ·Mb(µs) · ε2, ε} · 2−s−10, (4.2) where νK∗ is defined in Proposition 3.3 and
Mb(µs) = 4 ∥∥∥Q−1(I −A)> · [(I −A)Q−1(I −A)> +BR−1B>]−1 · (Aµs + d)∥∥∥
2 · [ ν−1K∗ + σ −1 min(Ψ ) · σ −1 min(R) ]1/2 . (4.3)
In the s-th policy update step in Line 3 of Algorithm 1, we set the inputs via Theorem B.4 such that Jµs(πs+1) − Jµs(π∗µs) < εs, where the expected total cost Jµs(·) is defined in Problem 2.2, and π∗µs = Λ1(µs) is the optimal policy under the mean-field state µs. Then it holds with probability at least 1− ε5 that
‖µS − µ∗‖2 ≤ ε, ‖KS −K∗‖F ≤ ε, ‖bS − b∗‖2 ≤ (1 + L1) · ε.
Here µ∗ is the Nash mean-field state, KS and bS are parameters of the policy πS , and K∗ and b∗ are parameters of the Nash policy π∗.
Proof Sketch. The proof of the theorem is based on the convergence of the natural actor-critic algorithm 2 and a contraction argument. First, we prove in Theorem B.4 that Algorithm 2 converges linearly to the optimal policy of Problem 2.2. By this, in each iteration of Algorithm 1, we control the error between µs+1 and µ∗s+1 to be εs > 0 with high probability, where µ ∗ s+1 is the mean-field state generated by the optimal policy Λ1(µs); in other words, µ∗s+1 = Λ(µs). Combining the fact from Proposition 3.2 that Λ(·) is a contraction, we deduce that
‖µs+1 − µ∗‖2 ≤ ∥∥Λ(µs)− Λ(µ∗)∥∥2 + ε̃s ≤ L0 · ‖µs − µ∗‖2 + ε̃s
with high probability, where ε̃s > 0 is some error term and is specified in the detailed proof. Moreover, by telescoping sum, and note that the sum ∑S s=1 ε̃s is upper bounded by the desired error ε, we conclude the theorem. See §D.1 for a detailed proof.
We highlight that if the inputs of Algorithm 1 satisfy the conditions stated in Theorem B.4, it holds that Jµs(πs+1) − Jµs(π∗µs) < εs for any s ∈ [S]. See Theorem B.4 in §B.1 of the appendix for details. By Theorem 4.1, Algorithm 1 converges linearly to the unique Nash equilibrium pair (µ∗, π∗) of Problem 2.3. To the best of our knowledge, this theorem is the first successful attempt to establish that reinforcement learning with function approximation finds the Nash equilibrium pairs in mean-field games with theoretical guarantee, which lays the theoretical foundations for applying modern reinforcement learning techniques to general mean-field games.

heading: 5 CONCLUSION
text: For the discrete-time linear-quadratic mean-field games, we provide sufficient conditions for the existence and uniqueness of the Nash equilibrium pair. Moreover, we propose the mean-field actorcritic algorithm with linear function approximation that is shown converges to the Nash equilibrium pair with linear rate of convergence. Our algorithm can be modified to use other parametrized function classes, including deep neural networks, for solving mean-field games. For future research, we aim to extend our algorithm to other variations of mean-field games including risk-sensitive mean-field games (Saldi et al., 2018a; Tembine et al., 2014), robust mean-field games (Bauso et al., 2016), and partially observed mean-field games (Saldi et al., 2019).

heading: A NOTATIONS IN THE APPENDIX
text: In the proof, for convenience, for any invertible matrix M , we denote by M−> = (M−1)> = (M>)−1 and ‖M‖F the Frobenius norm. We also denote by svec(M) the symmetric vectorization of the symmetric matrix M , which is the vectorization of the upper triangular matrix of the symmetric matrix M , with off-diagonal entries scaled by √ 2. We denote by smat(·) the inverse operation. For any matrices G and H , we denote by G ⊗ H the Kronecker product, and G ⊗s H the symmetric Kronecker product, which is defined as a mapping on a vector svec(M) such that (G⊗s H)svec(M) = 1/2 · svec(HMG> +GMH>). For notational simplicity, we write Eπ(·) to emphasize that the expectation is taken following the policy π.

heading: B AUXILIARY ALGORITHMS AND ANALYSIS

heading: B.1 RESULTS IN D-LQR
text: In this section, we provide auxiliary results in analyzing Problem 2.2. First, we introduce the value functions of the Markov decision process (MDP) induced by Problem 2.2. We define the state- and action-value functions VK,b(x) and QK,b(x, u) as follows
VK,b(x) = ∞∑ t=0 { E [ c(xt, ut) |x0 = x ] − J(K, b) } , (B.1)
QK,b(x, u) = c(x, u)− J(K, b) + E [ VK,b(x1) |x0 = x, u0 = u ] , (B.2)
where xt follows the state transition, and ut follows the policy πK,b given xt. In other words, we have ut = −Kxt + b + σηt, where ηt ∼ N (0, I). The following proposition establishes the close forms of these value functions. Proposition B.1. The state-value function VK,b(x) takes the form of
VK,b(x) = x >PKx− tr(PKΦK) + 2f>K,b(x− µK,b)− µ>K,bPKµK,b, (B.3)
and the action-value function QK,b(x, u) takes the form of
QK,b(x, u) = ( x u )> ΥK ( x u ) + 2 ( pK,b qK,b )>( x u ) − tr(PKΦK)− σ2 · tr(R+ PKBB>)
− b>Rb+ 2b>RKµK,b − µ>K,b(Q+K>RK + PK)µK,b + 2f>K,b [ (Aµ+ d)− µK,b ] + (Aµ+ d)>PK(Aµ+ d), (B.4)
where fK,b = (I − A + BK)−>[(A − BK)>PK(Bb + Aµ + d) − K>Rb], and ΥK , pK,b, and qK,b are defined in (3.7).
Proof. See §E.6 for a detailed proof.
By Proposition B.1, we know that VK,b(x) is quadratic in x, while QK,b(x, u) is quadratic in (x>, u>)>. Now, we show that (3.5) holds. Proposition B.2. The expected total cost J(K, b) defined in Problem 2.2 takes the form of
J(K, b) = J1(K) + J2(K, b) + σ 2 · tr(R) + µ>Qµ,
where
J1(K) = tr [ (Q+K>RK)ΦK ] = tr(PKΨ ),
J2(K, b) = ( µK,b b )>( Q+K>RK −K>R −RK R )( µK,b b ) .
Here µK,b is defined in (3.2), ΦK is defined in (3.3), and PK is defined in (3.4).
Proof. See §E.3 for a detailed proof.
The following proposition establishes the gradients of J1(K) and J2(K, b), respectively. Proposition B.3. The gradient of J1(K) and the gradient of J2(K, b) with respect to b take the forms of
∇KJ1(K) = 2(Υ22KK −Υ21K ) · ΦK , ∇bJ2(K, b) = 2 [ Υ22K (−KµK,b + b) + Υ21KµK,b + qK,b ] ,
where ΥK and qK,b are defined in (3.7).
Proof. See §E.5 for a detailed proof.
The following theorem establishes the convergence of Algorithm 2. Theorem B.4 (Convergence of Algorithm 2). Assume that ρ(A − BK0) < 1. Let ε > 0 be a sufficiently small tolerance. We set
γ ≤ [ ‖R‖2 + ‖B‖22 · J(K0, b0) · σ−1min(Ψ ) ]−1 ,
N ≥ C · ‖ΦK∗‖2 · γ−1 · log { 4 [ J(K0, b0)− J(K∗, b∗) ] · ε−1 } ,
Tn ≥ poly ( ‖Kn‖F, ‖b0‖2, ‖µ‖2, J(K0, b0) ) · λ−4Kn · [ 1− ρ(A−BKn) ]−9 · ε−5, T̃n ≥ poly ( ‖Kn‖F, ‖b0‖2, ‖µ‖2, J(K0, b0) ) · λ−2Kn · [ 1− ρ(A−BKn)
]−12 · ε−12, γn,t = γ0 · t−1/2,
γb ≤ min { 1− ρ(A−BKN ), [ 1− ρ(A−BKN ) ]−2 · (‖B‖22 · ‖KN‖22 · ‖R‖2 + ‖B‖22 · ‖Q‖2)}, H ≥ C0 · ν−1KN · (γ b)−1 · log { 4 [ J(KN , b0)− J(KN , bKN ) ] · ε−1 } ,
T bh ≥ poly ( ‖KN‖F, ‖bh‖2, ‖µ‖2, J(KN , b0) ) · λ−4KN · ν −4 KN · [ 1− ρ(A−BKN ) ]−11 · ε−5, T̃ bh ≥ poly ( ‖KN‖F, ‖bh‖2, ‖µ‖2, J(KN , b0) ) · λ−4KN · ν −2 KN · [ 1− ρ(A−BKN )
]−17 · ε−8, γbh,t = γ0 · t−1/2,
where C, C0, and γ0 are positive absolute constants, {Kn}n∈[N ] and {bh}h∈[H] are the sequences generated by Algorithm 2, λKn is specified in Proposition B.6, and νKN is specified in Proposition 3.3. Then it holds with probability at least 1− ε10 that
J(KN , bH)− J(K∗, b∗) < ε, ‖bH − b∗‖2 ≤Mb(µ) · ε1/2, ‖KN −K∗‖F ≤ [ σ−1min(Ψ ) · σ −1 min(R) · ε ]1/2 , ‖µ̂KN ,bH − µK∗,b∗‖2 ≤ ε,
where Mb(µ) is defined in (4.3).
Proof. See §D.2 for a detailed proof.
By Theorem B.4, given any mean-field state µ, Algorithm 2 converges linearly to the optimal policy π∗µ of Problem 2.2.

heading: B.2 PRIMAL-DUAL POLICY EVALUATION ALGORITHM
text: Note that the critic update steps in Algorithm 2 are built upon the estimators of the matrix ΥK and the vector qK,b. We now derive a policy evaluation algorithm to establish the estimators of ΥK and qK,b, which is based on gradient temporal difference algorithm (Sutton et al., 2009a).
We define the feature vector as
ψ(x, u) = ( ϕ(x, u) x− µK,b
u− (−KµK,b + b)
) , (B.5)
where
ϕ(x, u) = svec
[( x− µK,b
u− (−KµK,b + b)
)( x− µK,b
u− (−KµK,b + b)
)>] .
Recall svec(M) gives the symmetric vectorization of the symmetric matrix M . We also define
αK,b =  svec(ΥK) ΥK ( µK,b
−KµK,b + b
) + ( pK,b qK,b ) , (B.6) where ΥK , pK,b, and qK,b are defined in (3.7). To estimate ΥK and qK,b, it suffices to estimate αK,b. Meanwhile, we define
ΘK,b = EπK,b { ψ(x, u) [ ψ(x, u)− ψ(x′, u′) ]>} , (B.7)
where (x′, u′) is the state-action pair after (x, u) following the policy πK,b and the state transition. The following proposition characterizes the connection between ΘK,b and αK,b. Proposition B.5. It holds that(
1 0 EπK,b [ ψ(x, u) ] ΘK,b )( J(K, b) αK,b ) = ( J(K, b) EπK,b [ c(x, u)ψ(x, u) ]) , where ψ(x, u) is defined in (B.5), αK,b is defined in (B.6), and ΘK,b is defined in (B.7).
Proof. See §E.7 for a detailed proof.
By Proposition B.5, to obtain αK,b, it suffices to solve the following linear system in ζ = (ζ1, ζ>2 ) >, Θ̃K,b · ζ = (
J(K, b) EπK,b [ c(x, u)ψ(x, u) ]) , (B.8) where for notational convenience, we define
Θ̃K,b =
( 1 0
EπK,b [ ψ(x, u) ] ΘK,b
) . (B.9)
The following proposition shows that ΘK,b is invertible. Proposition B.6. If ρ(A − BK) < 1, then the matrix ΘK,b is invertible, and ‖ΘK,b‖2 ≤ 4(1 + ‖K‖2F)2 · ‖ΦK‖22. Also, σmin(Θ̃K,b) ≥ λK , where λK only depends on ‖K‖2 and ρ(A−BK).
Proof. See §E.8 for a detailed proof.
By Proposition B.6, ΘK,b is invertible. Therefore, (B.8) admits the unique solution ζK,b = (J(K, b), α>K,b) >.
Now, we present the primal-dual gradient temporal difference algorithm.
Primal-Dual Gradient Method. Instead of solving (B.8) directly, we minimize the following loss function with respect to ζ = ((ζ1)>, (ζ2)>),[
ζ1 − J(K, b) ]2 + ∥∥∥EπK,b[ψ(x, u)]ζ1 + ΘK,bζ2 − EπK,b[c(x, u)ψ(x, u)]∥∥∥2
2 . (B.10)
By Fenchel’s duality, the minimization of (B.10) is equivalent to the following primal-dual min-max problem,
min ζ∈Vζ max ξ∈Vξ
F (ζ, ξ) = { EπK,b [ ψ(x, u) ] ζ1 + ΘK,bζ 2 − EπK,b [ c(x, u)ψ(x, u) ]}> ξ2 (B.11)
+ [ ζ1 − J(K, b) ] · ξ1 − ‖ξ‖22/2,
where we restrict the primal variable ζ in a compact set Vζ and the dual variable ξ in a compact set Vξ, which are specified in Definition B.7. It holds that
∇ζ1F = ξ1 + EπK,b [ ψ(x, u) ]> ξ2, ∇ζ2F = Θ>K,bξ2, ∇ξ1F = ζ1 − J(K, b)− ξ1,
∇ξ2F = EπK,b [ ψ(x, u) ] ζ1 + ΘK,bζ 2 − EπK,b [ c(x, u)ψ(x, u) ] − ξ2. (B.12)
The primal-dual gradient method updates ζ and ξ via
ζ1 ← ζ1 − γ · ∇ζ1F (ζ, ξ), ζ2 ← ζ2 − γ · ∇ζ2F (ζ, ξ) ξ1 ← ξ1 − γ · ∇ξ1F (ζ, ξ), ξ2 ← ξ2 − γ · ∇ξ2F (ζ, ξ). (B.13)
Estimation of Mean-Field State µK,b. To utilize the primal-dual gradient method in (B.13), it remains to evaluate the feature vector ψ(x, u). Note that by (B.5), the evaluation of the feature vector ψ(x, u) requires the mean-field state µK,b. In what follows, we establish the estimator µ̂K,b of the mean-field state µK,b by simulating the MDP following the policy πK,b for T̃ steps, and calculate the estimated feature vector ψ̂(x, u) by
ψ̂(x, u) = ( ϕ̂(x, u) x− µ̂K,b
u− (−Kµ̂K,b + b)
) , (B.14)
where ϕ̂(x, u) takes the form of
ϕ̂(x, u) = svec
[( x− µ̂K,b
u− (−Kµ̂K,b + b)
)( x− µ̂K,b
u− (−Kµ̂K,b + b)
)>] .
We now define the sets Vζ and Vξ in (B.11). Definition B.7. Given K0 and b0 such that ρ(A − BK0) < 1 and J(K0, b0) < ∞, we define the sets Vζ and Vξ as
Vζ = { ζ : 0 ≤ ζ1 ≤ J(K0, b0), ‖ζ2‖2 ≤Mζ,1 +Mζ,2 · (1 + ‖K‖F) · [ 1− ρ(A−BK) ]−1} ,
Vξ = { ξ : |ξ1| ≤ J(K0, b0), ‖ξ2‖2 ≤Mξ · ( 1 + ‖K‖2F )3 · [1− ρ(A−BK)]−1}. Here Mζ,1, Mζ,2, and Mξ are constants independent of K and b, which take the forms of
Mζ,1 = [( ‖Q‖F + ‖R‖F ) + ( ‖A‖2F + ‖B‖2F ) · √ d · J(K0, b0) · σ−1min(Ψω) ] + ( ‖A‖2 + ‖B‖2 ) · J(K0, b0)2 · σ−1min(Ψω) · σ −1 min(Q),
+ [( ‖Q‖2 + ‖R‖2 ) + ( ‖A‖2 + ‖B‖2 )2 · J(K0, b0) · σ−1min(Ψω)] · J(K0, b0) · [ σ−1min(Q) + σ −1 min(R)
] Mζ,2 = ( ‖A‖2 + ‖B‖2 ) · (κQ + κR), Mξ = C · (Mζ,1 +Mζ,2) · J(K0, b0)2 · σ−2min(Q),
where C is a positive absolute constant, and κQ and κR are condition numbers of Q and R, respectively.
We summarize the primal-dual gradient temporal difference algorithm in Algorithm 3. Hereafter, for notational convenience, we denote by ψ̂t the estimated feature vector ψ̂(xt, ut).
We now characterize the rate of convergence of Algorithm 3. Theorem B.8 (Convergence of Algorithm 3). Given K0, b0, K, and b such that ρ(A − BK0) < 1 and J(K, b) ≤ J(K0, b0), we define the sets Vζ and Vξ through Definition B.7. Let γt = γ0t −1/2, where γ0 is a positive absolute constant. Let ρ ∈ (ρ(A − BK), 1). For T̃ ≥ poly0(‖K‖F, ‖b‖2, ‖µ‖2, J(K0, b0)) · (1− ρ)−6 and a sufficiently large T , it holds with probability at least 1− T−4 − T̃−6 that ‖α̂K,b − αK,b‖22 ≤ λ−2K · poly1 ( ‖K‖F, ‖b‖2, ‖µ‖2, J(K0, b0) ) · [ log6 T
T 1/2 · (1− ρ)4 +
log T̃
T̃ 1/4 · (1− ρ)2
] ,
Algorithm 3 Primal-Dual Gradient Temporal Difference Algorithm.
1: Input: Policy πK,b, mean-field state µ, numbers of iteration T̃ and T , stepsizes {γt}t∈[T ], parameters K0 and b0. 2: Define the sets Vζ and Vξ via Definition B.7 with K0 and b0. 3: Initialize the parameters by ζ0 ∈ Vζ and ξ0 ∈ Vξ. 4: Sample x̃0 from the the stationary distribution N (µK,b,ΦK). 5: for t = 0, . . . , T̃ − 1 do 6: Given the mean-field state µ, take action ũt following πK,b and generate the next state x̃t+1. 7: end for 8: Set µ̂K,b ← 1/T̃ · ∑T̃ t=1 x̃t and compute the estimated feature vector ψ̂ via (B.14).
9: Sample x0 from the the stationary distribution N (µK,b,ΦK). 10: for t = 0, . . . , T − 1 do 11: Given the mean-field state µ, take action ut following πK,b, observe the cost ct, and generate the next state xt+1. 12: Set δt+1 ← ζ1t + (ψ̂t − ψ̂t+1)>ζ2t − ct. 13: Update parameters via
ζ1t+1 ← ζ1t − γt+1 · (ξ1t + ψ̂>t ξ2t ), ζ2t+1 ← ζ2t − γt+1 · ψ̂t(ψ̂t − ψ̂t+1)>ξ2t ,
ξ1t+1 ← (1− γt+1) · ξ1t + γt+1 · (ζ1t − ct), ξ2t+1 ← (1− γt+1) · ξ2t + γt+1 · δt+1 · ψ̂t.
14: Project ζt+1 and ξt+1 to Vζ and Vξ, respectively. 15: end for 16: Set α̂K,b ← ( ∑T t=1 γt) −1 · ( ∑T t=1 γt · ζ2t ), and
Υ̂K ← smat(α̂K,b,1), ( p̂K,b q̂K,b ) ← α̂K,b,2 − Υ̂K ( µ̂K,b −Kµ̂K,b + b ) ,
where α̂K,b,1 = (α̂K,b) (k+d+1)(k+d)/2 1 and α̂K,b,2 = (α̂K,b) (k+d+3)(k+d)/2 (k+d+1)(k+d)/2+1.
17: Output: Estimators µ̂K,b, Υ̂K , and q̂K,b.
where λK is defined in Proposition B.6. Same bounds for ‖Υ̂K − ΥK‖2F, ‖p̂K,b − pK,b‖22, and ‖q̂K,b − qK,b‖22 hold. Meanwhile, it holds with probability at least 1− T̃−6 that
‖µ̂K,b − µK,b‖2 ≤ log T̃
T̃ 1/4 · (1− ρ)−2 · poly2
( ‖ΦK‖2, ‖K‖F, ‖b‖2, ‖µ‖2, J(K0, b0) ) .
Proof. See §D.3 for a detailed proof.

heading: B.3 TEMPORAL DIFFERENCE POLICY EVALUATION ALGORITHM
text: Besides the primal-dual gradient temporal difference algorithm, we can also evaluate αK,b by TD(0) method (Sutton and Barto, 2018) in practice, which is presented in Algorithm 4.
Note that in related literature (Bhandari et al., 2018; Korda and La, 2015), non-asymptotic convergence analysis of TD(0) method with linear function approximation is only applied to discounted MDP. As for our ergodic setting, the convergence of TD(0) method is only shown asymptotically (Borkar and Meyn, 2000; Kushner and Yin, 2003) using ordinary differential equation method. Therefore, in the convergence theorem proposed in §3, we only focus on the primal-dual gradient temporal difference method (Algorithm 3) to establish non-asymptotic convergence result.

heading: C GENERAL FORMULATION
text: Compared with Problem 2.3, a more general formulation includes an additional term x>t Pµ in the cost function. For the completeness of this paper, we define this general formulation here. Following from a same argument as in §2, it suffices to study the setting where t is sufficiently large. First, we propose the following general drifted LQR (general D-LQR) problem, which is parallel to Problem 2.2.
Algorithm 4 Temporal Difference Policy Evaluation Algorithm.
1: Input: Policy πK,b, number of iteration T̃ and T , stepsizes {γt}t∈[T ]. 2: Sample x̃0 from the stationary distribution N (µK,b,ΦK). 3: for t = 0, . . . , T̃ − 1 do 4: Take action ũt under the policy πK,b and generate the next state x̃t+1. 5: end for 6: Set µ̂K,b ← 1/T̃ · ∑T̃ t=1 x̃t. 7: Sample x0 from the the stationary distribution N (µK,b,ΦK). 8: for t = 0, . . . , T do 9: Given the mean-field state µ, take action ut following πK,b, observe the cost ct, and generate
the next state xt+1. 10: Set δt+1 ← ζ1t + (ψ̂t − ψ̂t+1)>ζ2t − ct. 11: Update parameters via ζ1t+1 ← (1− γt+1) · ζ1t + γt+1 · ct and ζ2t+1 ← ζ2t − γt+1 · δt+1 · ψ̂t. 12: Project ζt to V ′ζ , where V ′ζ is a compact set. 13: end for 14: Set α̂K,b ← ( ∑T t=1 γt) −1 · ( ∑T t=1 γt · ζ2t ), and
Υ̂K ← smat(α̂K,b,1), ( p̂K,b q̂K,b ) ← α̂K,b,2 − Υ̂K ( µ̂K,b −Kµ̂K,b + b ) ,
where α̂K,b,1 = (α̂K,b) (k+d+1)(k+d)/2 1 and α̂K,b,2 = (α̂K,b) (k+d+3)(k+d)/2 (k+d+1)(k+d)/2+1.
15: Output: Estimators µ̂K,b, Υ̂K , and q̂K,b.
Problem C.1 (General D-LQR). For any given mean-field state µ ∈ Rm, consider the following formulation
xt+1 = Axt +But +Aµ+ d+ ωt, c̃µ(xt, ut) = x > t Qxt + u > t Rut + µ >Qµ+ 2x>t Pµ,
J̃µ(π) = lim T→∞
E
[ 1
T T∑ t=0 c̃µ(xt, ut)
] ,
where xt ∈ Rm is the state vector, ut ∈ Rk is the action vector generated by the policy π, ωt ∈ Rm is an independent random noise term following the Gaussian distribution N (0,Ψω), and d ∈ Rm is a drift term. We aim to find an optimal policy π∗µ such that J̃µ(π ∗ µ) = infπ∈Π J̃µ(π).
In Problem C.1, the unique optimal policy π∗µ(·) still admits a linear form π∗µ(xt) = −Kπ∗µxt + bπ∗µ (Anderson and Moore, 2007), where the matrix Kπ∗µ ∈ R k×m and the vector bπ∗µ ∈ R k are the parameters of the policy π. It then suffices to find the optimal policy in the class Π defined in (2.1).
Parallel to Problem 2.3, we define the general LQ-MFG problem as follows.
Problem C.2 (General LQ-MFG). We consider the following formulation
xt+1 = Axt +But +Aµ+ d+ ωt, c̃(xt, ut) = x > t Qxt + u > t Rut + µ >Qµ+ 2x>t Pµ,
J̃(π, µ) = lim T→∞
E
[ 1
T T∑ t=0 c̃(xt, ut)
] ,
where xt ∈ Rm is the state vector, ut ∈ Rk is the action vector generated by the policy π, µ ∈ Rm is the mean-field state, ωt ∈ Rm is an independent random noise term following the Gaussian distribution N (0,Ψω), and d ∈ Rm is a drift term. We aim to find a pair (µ∗, π∗) such that (i) J̃(π∗, µ∗) = infπ∈Π J̃(π, µ
∗); (ii) Ex∗t converges to µ∗ as t → ∞, where {x∗t }t≥0 is the Markov chain of states generated by the policy π∗.
One can see that Problem C.2 aims to find a Nash equilibrium pair (µ∗, π∗).
Similar to the discussions in §3.2, to solve Problem C.2, one can design an algorithm similar to Algorithm 1, which solves Problem C.1 and obtain the new mean-field state at each iteration. We omit the detailed algorithm here. Now, we focus on solving Problem C.1 in the sequel.
Similar to §3.3, we drop the subscript µ when we focus on Problem C.1 for a fixed µ. We write πK,b(x) = −Kx + b + σ · η to emphasize the dependence on K and b, and J̃(K, b) = J̃(πK,b) consequently. We derive an explicit form of the expected total cost J̃(K, b) in the following proposition.
Proposition C.3. The expected total cost J̃(K, b) in Problem C.1 is decomposed as
J̃(K, b) = J̃1(K) + J̃2(K, b) + σ 2 · tr(R) + µ>Qµ,
where J̃1(K) and J̃2(K, b) take the forms of J̃1(K) = tr [ (Q+K>RK)ΦK ] = tr(PKΨ ),
J̃2(K, b) = ( µK,b b )>( Q+K>RK −K>R −RK R )( µK,b b ) + 2µ>PµK,b.
Here µK,b is given in (3.2), ΦK is given in (3.3), and PK is given in (3.4).
Proof. The proof is similar to the one of Proposition B.2. Thus we omit it here.
Compared with the form of J(K, b) given in (3.5), we see that the only difference is that J̃(K, b) contains an extra term 2µ>PµK,b in J̃2(K, b), which is only a linear term in b (recall that µK,b is linear in b by (3.2)). Thus, J̃2(K, b) is still strongly convex in b, as shown in the proposition below.
Proposition C.4. Given any K, the function J̃2(K, b) is νK-strongly convex in b, here νK = σmin(Y > 1,KY1,K + Y > 2,KY2,K), where Y1,K = R
1/2K(I − A + BK)−1B − R1/2 and Y2,K = Q1/2(I − A + BK)−1B. Also, J̃2(K, b) has ιK-Lipschitz continuous gradient in b, where ιK is upper bounded such that ιK ≤ [1− ρ(A−BK)]−2 · (‖B‖2∗ · ‖K‖2∗ · ‖R‖∗ + ‖B‖2∗ · ‖Q‖∗).
Proof. The proof is similar to the one of Proposition 3.3. Thus we omit it here.
Parallel to Proposition 3.4, we derive a similar proposition in the sequel.
Proposition C.5. Denote by b̃K = argminb J̃2(K, b), then J̃2(K, b̃K) takes the form
J̃2(K, b̃ K) = ( Aµ+ d P>µ )>( S S(I −A)Q−1 Q−1(I −A)>S 3Q−1(I −A)>S(I −A)Q−1 −Q−1 )( Aµ+ d P>µ ) ,
which is independent of K. Here S = [(I − A)Q−1(I − A)> + BR−1B>]−1. And b̃K takes the form
b̃K = [ KQ−1(I −A)> −R−1B> ] · S · [ (Aµ+ d) + (I −A)Q−1P>µ ] −KQ−1P>µ.
Proof. The proof is similar to the one of Proposition 3.4. Thus we omit it here.
Similar to Problem 2.2, we define the state- and action-value functions as
ṼK,b(x) = ∞∑ t=0 { E [ c̃(xt, ut) |x0 = x, ut = −Kxt + b+ σηt ] − J̃(K, b) } ,
Q̃K,b(x, u) = c̃(x, u)− J̃(K, b) + E [ ṼK,b(x ′) |x, u ] ,
where the x′ is the state generated by the state transition after the state-action pair (x, u). A slight modification of Proposition B.1 gives the proposition below.
Proposition C.6. For Problem C.1, the state-value function ṼK,b(x) takes the form
ṼK,b(x) = x >PKx− tr(PKΦK) + 2f̃>K,b(x− µK,b)− (µK,b)>PKµK,b,
and the action-value function Q̃K,b(x, u) takes the form
Q̃K,b(x, u) = ( x u )> ΥK ( x u ) + 2 ( p̃K,b q̃K,b )>( x u ) − tr(PKΦK)− σ2 · tr(R+ PKBB>)− b>Rb
+ 2b>RKµK,b − (µK,b)>(Q+K>RK + PK)µK,b + 2f̃>K,b [ (Aµ+ d)− µK,b ] + (Aµ+ d)>PK(Aµ+ d)− 2µ>PµK,b.
Here the matrix ΥK is given in (3.7), and the vectors p̃K,b, q̃K,b are given as( p̃K,b q̃K,b ) = ( A> [ PK · (Aµ+ d) + f̃K,b ] + Pµ B> [ PK · (Aµ+ d) + f̃K,b ] ) , (C.1) where the vector f̃K,b = (I −A+BK)−>[(A−BK)>PK(Bb+Aµ+ d)−K>Rb+ Pµ].
Proof. The proof is similar to the one of Proposition B.1. Thus we omit it here.
Now we establish the gradients of J̃(K, b) for Problem C.1.
Proposition C.7. The gradient of J̃1(K) and the gradient of J̃2(K, b) w.r.t. b takes the form ∇K J̃1(K) = 2(Υ22KK −Υ21K ) · ΦK , ∇bJ̃2(K, b) = 2 [ Υ22K (−KµK,b + b) + Υ21KµK,b + q̃K,b ] ,
where the matrix ΥK is given in (3.7), and the vector q̃K,b is given in (C.1).
Proof. The proof is similar to the one of Proposition B.3. Thus we omit it here.
Equipped with above results, parallel to the analysis in §3, it is clear that by slight modification of Algorithms 1, 2, and 3, we derive similar actor-critic algorithms to solve both Problem C.2 and Problem C.1, where all the non-asymptotic convergence results hold. We omit the algorithms and the convergence results here.

heading: D PROOFS OF THEOREMS

heading: D.1 PROOF OF THEOREM 4.1
text: We define µ∗s+1 = Λ(µs), which is the mean-field state generated by the optimal policy πK∗(µs),b∗(µs) = Λ1(µs) under the current mean-field state µs. By Proposition 3.4, the optimal K∗(µ) is independent of the mean-field state µ. Therefore, we write K∗ = K∗(µ) hereafter for notational convenience. By (3.2), we know that
µ∗s+1 = (I −A+BK∗)−1 · [ Bb∗(µs) +Aµs + d ] .
We define
µ̃s+1 = (I −A+BKs)−1(Bbs +Aµs + d), which is the mean-field state generated by the policy πs under the current mean-field state µs, where Ks and bs are the parameters of the policy πs. By triangle inequality, we have
‖µs+1 − µ∗‖2 ≤ ‖µs+1 − µ̃s+1‖2︸ ︷︷ ︸ E1 + ‖µ̃s+1 − µ∗s+1‖2︸ ︷︷ ︸ E2 + ‖µ∗s+1 − µ∗‖2︸ ︷︷ ︸ E3 , (D.1)
where µs+1 is generated by Algorithm 1. We upper bound E1, E2, and E3 in the sequel.
Upper Bound of E1. By Theorem B.4, it holds with probability at least 1− ε10 that E1 = ‖µs+1 − µ̃s+1‖2 < εs ≤ ε/8 · 2−s, (D.2)
where εs is given in (4.2).
Upper Bound of E2. By the triangle inequality, we have E2 = ∥∥∥(I −A+BKs)−1(Bbs +Aµs + d)− (I −A+BK∗)−1 · [Bb∗(µs) +Aµs + d]∥∥∥
2 ≤ ∥∥Bb∗(µs) +Aµs + d∥∥2 · ∥∥∥[I −A+BK∗ +B(Ks −K∗)]−1 − (I −A+BK∗)−1∥∥∥2
+ ∥∥(I −A+BKs)−1∥∥2 · ‖B‖2 · ∥∥bs − b∗(µs)∥∥2. (D.3)
By Taylor’s expansion, we have∥∥∥[I −A+BK∗ +B(Ks −K∗)]−1 − (I −A+BK∗)−1∥∥∥ 2
= ∥∥∥(I −A+BK∗)−1[I + (I −A+BK∗)−1B(Ks −K∗)]−1 − (I −A+BK∗)−1∥∥∥
2 ≤ 2 ∥∥(I −A+BK∗)−1B(Ks −K∗)(I −A+BK∗)−1∥∥2. (D.4)
Meanwhile, by Taylor’s expansion, it holds with probability at least 1− ε10 that∥∥(I −A+BKs)−1∥∥2 = ∥∥∥(I −A+BK∗ +B(Ks −K∗))−1∥∥∥2 = ∥∥∥(I −A+BK∗)−1(I + (I −A+BK∗)−1B(Ks −K∗))−1∥∥∥
2 ≤ [ 1− ρ(A−BK∗) ]−1 · (1 + ∥∥(I −A+BK∗)−1B∥∥ 2 · ‖K∗ −Ks‖2 ) ≤ 2 [ 1− ρ(A−BK∗) ]−2 , (D.5)
where the last inequality comes from Theorem B.4. By plugging (D.4) and (D.5) in (D.3), it holds with probability at least 1− ε10 that
E2 ≤ 2 ∥∥Bb∗(µs) +Aµs + d∥∥2 · ∥∥(I −A+BK∗)−1B(Ks −K∗)(I −A+BK∗)−1∥∥2 + ∥∥(I −A+BKs)−1∥∥2 · ‖B‖2 · ∥∥bs − b∗(µs)∥∥2
≤ 2 ∥∥Bb∗(µs) +Aµs + d∥∥2 · [1− ρ(A−BK∗)]−2 · ‖B‖2 · ‖Ks −K∗‖2 (D.6) + 2 [ 1− ρ(A−BK∗)
]−2 · ‖B‖2 · ∥∥bs − b∗(µs)∥∥2. By Proposition 3.4, it holds that∥∥Bb∗(µs) +Aµs + d∥∥2 ≤ L1 · ‖B‖2 · ‖µs‖2 + ‖A‖2 · ‖µs‖2 + ‖d‖2
≤ ( L1 · ‖B‖2 + ‖A‖2 ) · ‖µs‖2 + ‖d‖2, (D.7)
where the scalar L1 is defined in Assumption 3.1. Meanwhile, by Theorem B.4, it holds with probability at least 1− ε10 that
‖Ks −K∗‖F ≤ [ σ−1min(Ψ ) · σ −1 min(R) · εs ]1/2 , ∥∥bs − b∗(µs)∥∥2 ≤Mb(µs) · ε1/2s , (D.8) where Mb(µs) is defined in (4.3). Combining (D.6), (D.7), (D.8), and the choice of εs in (4.2), it holds with probability at least 1− ε10 that
E2 ≤ ε/8 · 2−s. (D.9)
Upper Bound of E3. By Proposition 3.2, we have E3 = ‖µ∗s+1 − µ∗‖2 = ∥∥Λ(µs)− Λ(µ∗)∥∥2 ≤ L0 · ‖µs − µ∗‖2, (D.10)
where L0 = L1L3 + L2 by Assumption 3.1.
By plugging (D.2), (D.9), and (D.10) in (D.1), we know that
‖µs+1 − µ∗‖2 ≤ L0 · ‖µs − µ∗‖2 + ε · 2−s−2, (D.11)
which holds with probability at least 1 − ε10. Following from (D.11) and a union bound argument with S = O(log(1/ε)), it holds with probability at least 1− ε5 that
‖µS − µ∗‖2 ≤ LS0 · ‖µ0 − µ∗‖2 + ε/2,
where we use the fact that L0 < 1 by Assumption 3.1. By the choice of S in (4.1), it further holds with probability at least 1− ε6 that
‖µS − µ∗‖ ≤ ε. (D.12)
By Theorem B.4 and the choice of εs in (4.2), it holds with probability at least 1− ε5 that ‖KS −K∗‖F = ∥∥KS −K∗(µS)∥∥F ≤ [σ−1min(Ψ ) · σ−1min(R) · εS]1/2 ≤ ε. (D.13)
Meanwhile, by the triangle inequality and the choice of εs in (4.2), it holds with probability at least 1− ε5 that
‖bS − b∗‖2 ≤ ∥∥bS − b∗(µS)∥∥2 + ∥∥b∗(µS)− b∗∥∥2
≤Mb(µS) · ε1/2S + L1 · ‖µS − µ ∗‖2 ≤ (1 + L1) · ε, (D.14)
where the second inequality comes from Theorem B.4 and Proposition 3.4, and the last inequality comes from (D.12). By (D.12), (D.13), and (D.14), we conclude the proof of the theorem.

heading: D.2 PROOF OF THEOREM B.4
text: Proof. We first show that J1(KN ) − J1(K∗) < ε/2 with a high probability, then show that J2(KN , bH)− J2(K∗, b∗) < ε/2 with a high probability. Then we have
J(KN , bN )− J(K∗, b∗) = J1(KN ) + J2(KN , bH)− J1(K∗)− J2(K∗, b∗) < ε
with a high probability, which proves Theorem B.4.
Part 1. We show that J1(KN )− J1(K∗) < ε/2 with a high probability. We first bound J1(K1)− J1(K2) for any K1 and K2. By Proposition B.2, J1(K) takes the form of
J1(K) = tr(PKΨ ) = Ey∼N (0,Ψ )(y >PKy). (D.15)
The following lemma calculates y>PK1y − y>PK2y for any K1 and K2. Lemma D.1. Assume that ρ(A − BK1) < 1 and ρ(A − BK2) < 1. For any state vector y, we denote by {yt}t≥0 the sequence generated by the state transition yt+1 = (A − BK2)yt with initial state y0 = y. It holds that
y>PK2y − y>PK1y = ∑ t≥0 DK1,K2(yt),
where
DK1,K2(y) = 2y >(K2 −K1)(Υ22K1K1 −Υ 21 K1)y + y >(K2 −K1)>Υ22K1(K2 −K1)y.
Here ΥK is defined in (3.7).
Proof. See §F.1 for a detailed proof.
The following lemma shows that J1(K) is gradient dominant. Lemma D.2. Let K∗ be the optimal parameter and K be a parameter such that J1(K) < ∞, then it holds that
J1(K)− J1(K∗) ≤ σ−1min(R) · ‖ΦK∗‖2 · tr [ (Υ22KK −Υ21K )>(Υ22KK −Υ21K ) ] , (D.16)
J1(K)− J1(K∗) ≥ σmin(Ψω) · ‖Υ22K ‖−12 · tr [ (Υ22KK −Υ21K )>(Υ22KK −Υ21K ) ] . (D.17)
Proof. See §F.2 for a detailed proof.
Recall that from Algorithm 2, the parameter K is updated via
Kn+1 = Kn − γ · (Υ̂22KnKn − Υ̂ 21 Kn), (D.18)
where Υ̂Kn is the output of Algorithm 3. We upper bound |J1(Kn+1) − J1(K∗)| in the sequel. First, we show that if J1(Kn)− J1(K∗) ≥ ε/2 holds for any n ≤ N , we obtain that
J1(KN ) ≤ J1(KN−1) ≤ · · · ≤ J1(K0), (D.19)
which holds with probability at least 1− ε13. We prove (D.19) by mathematical induction. Suppose that
J1(Kn) ≤ J1(Kn−1) ≤ · · · ≤ J1(K0), (D.20)
which holds for n = 0. In what follows, we define K̃n+1 as
K̃n+1 = Kn − γ · (Υ22KnKn −Υ 21 Kn), (D.21)
where ΥKn is given in (3.7). By (D.21), we have J1(K̃n+1)− J1(Kn) = Ey∼N (0,Ψ ) [ y>(PK̃n+1 − PKn)y ] = −2γ · tr [ ΦK̃n+1 · (Υ 22 KnKn −Υ 21 Kn) >(Υ22KnKn −Υ 21 Kn)
] + γ2 · tr [ ΦK̃n+1 · (Υ 22 KnKn −Υ 21 Kn) >Υ22Kn(Υ 22 KnKn −Υ 21 Kn)
] ≤ −2γ · tr [ ΦK̃n+1 · (Υ 22 KnKn −Υ 21 Kn) >(Υ22KnKn −Υ 21 Kn) ] (D.22)
+ γ2 · ‖Υ22Kn‖2 · tr [ ΦK̃n+1 · (Υ 22 KnKn −Υ 21 Kn) >(Υ22KnKn −Υ 21 Kn) ] ,
where the first equality comes from (D.15), the second equality comes from Lemma D.1, and the last inequality comes from the trace inequality. By the definition of ΥK in (3.7), we obtain that
‖Υ22Kn‖2 ≤ ‖R‖2 + ‖B‖ 2 2 · ‖PKn‖2 ≤ ‖R‖2 + ‖B‖22 · J1(Kn) · σ−1min(Ψ )
≤ ‖R‖2 + ‖B‖22 · J1(K0) · σ−1min(Ψ ), (D.23)
where the second inequality comes from Proposition B.2. By plugging (D.23) and the choice of stepsize γ ≤ [‖R‖2 + ‖B‖22 · J1(K0) · σ−1min(Ψ )]−1 into (D.22), we obtain that
J1(K̃n+1)− J1(Kn) ≤ −γ · tr [ ΦK̃n+1 · (Υ 22 KnKn −Υ 21 Kn) >(Υ22KnKn −Υ 21 Kn) ] ≤ −γ · σmin(Ψ ) · tr [ (Υ22KnKn −Υ 21 Kn) >(Υ22KnKn −Υ 21 Kn)
] ≤ −γ · σmin(Ψ ) · σmin(R) · ‖ΦK∗‖−12 · [ J1(Kn)− J1(K∗) ] < 0, (D.24)
where the last inequality comes from Lemma D.2.
The following lemma upper bounds |J1(K̃n+1)− J1(Kn+1)|. Lemma D.3. Assume that J1(Kn) ≤ J1(K0). It holds with probability at least 1− ε15 that∣∣J1(K̃n+1)− J1(Kn+1)∣∣ ≤ γ · σmin(Ψ ) · σmin(R) · ‖ΦK∗‖−12 · ε/4, where Kn+1 and K̃n+1 are defined in (D.18) and (D.21), respectively.
Proof. See §F.3 for a detailed proof.
Combining (D.24) and Lemma D.3, if J1(Kn) − J1(K∗) ≥ ε/2, it holds with probability at least 1− ε15 that
J1(Kn+1)− J1(Kn) ≤ J1(K̃n+1)− J1(Kn) + ∣∣J1(K̃n+1)− J1(Kn+1)∣∣
≤ −γ · σmin(Ψ ) · σmin(R) · ‖ΦK∗‖−12 · ε/4 < 0. (D.25)
Combining (D.20) and (D.25), it holds with probability at least 1− ε15 that J1(Kn+1) ≤ J1(Kn) ≤ · · · ≤ J1(K0). Finally, following from a union bound argument and the choice of N in Theorem B.4, if J1(Kn)− J1(K
∗) ≥ ε/2 holds for any n ≤ N , we have J1(KN ) ≤ J1(KN−1) ≤ · · · ≤ J1(K0),
which holds with probability at least 1− ε13. Thus, we complete the proof of (D.19). Combining (D.24) and (D.25), for J1(Kn)− J1(K∗) ≥ ε/2, we have
J1(Kn+1)− J1(K∗) ≤ [ 1− γ · σmin(Ψ ) · σmin(R) · ‖ΦK∗‖−12 ] · [ J1(Kn)− J1(K∗) ] ,
which holds with probability at least 1 − ε13. Meanwhile, following from a union bound argument and the choice of N in Theorem B.4, it holds with probability at least 1− ε11 that
J1(KN )− J1(K∗) ≤ ε/2. (D.26)
The following lemma upper bounds ‖KN −K∗‖F. Lemma D.4. For any K, we have
‖K −K∗‖2F ≤ σ−1min(Ψ ) · σ −1 min(R) · [ J1(K)− J1(K∗) ] .
Proof. See §F.4 for a detailed proof.
Combining (D.26) and Lemma D.4, we have ‖KN −K∗‖F ≤ [ σ−1min(Ψ ) · σ −1 min(R) · ε/2 ]1/2 , (D.27) which holds with probability 1− ε11.
Part 2. We show that J2(KN , bH) − J2(K∗, b∗) < ε/2 with high probability. Following from Proposition 3.4, it holds that J2(K∗, b∗) = J2(KN , bKN ). Therefore, it suffices to show that J2(KN , bH)− J2(KN , bKN ) < ε/2. First, we show that if J2(KN , bh)− J2(KN , bKN ) ≥ ε/2 for any h ≤ H , we obtain that J2(KN , bH) ≤ J2(KN , bH−1) ≤ · · · ≤ J2(KN , b1) ≤ J2(KN , b0), (D.28) which holds with probability at least 1− ε13. We prove (D.28) by mathematical induction. Suppose that
J2(KN , bh) ≤ J2(KN , bh−1) ≤ · · · ≤ J2(KN , b0), (D.29) Recall that by Algorithm 2, the parameter b is updated via
bh+1 = bh − γb · ∇̂bJ2(KN , bh). (D.30) Here
∇̂bJ2(KN , bh) = Υ̂22KN (−KN µ̂KN ,bh + bh) + Υ̂ 21 KN µ̂KN ,bh + q̂KN ,bh , (D.31)
where Υ̂KN and q̂KN ,bh are the outputs of Algorithm 3. We define b̃h+1 as
b̃h+1 = bh − γb · ∇bJ2(KN , bh). (D.32) Here
∇bJ2(KN , bh) = Υ22KN (−KNµKN ,bh + bh) + Υ 21 KNµKN ,bh + qKN ,bh , (D.33)
where ΥKN and qKN ,bh are defined in (3.7). We upper bound J2(KN , bh+1)− J2(KN , bKN ) in the sequel. Following from (D.32) and Proposition 3.3, we have
J2(KN , b̃h+1)− J2(KN , bh) ≤ −γb/2 · ∥∥∇bJ2(KN , bh)∥∥22
≤ −νKN · γb · [ J2(KN , bh)− J2(KN , bKN ) ] ≤ −νKN · γb · ε < 0, (D.34)
where νKN is specified in Proposition 3.3. The following lemma upper bounds |J2(KN , bh+1) − J2(KN , b̃h+1)|.
Lemma D.5. Assume that J2(KN , bh) ≤ J2(KN , b0). It holds with probability at least 1− ε15 that∣∣J2(KN , bh+1)− J2(KN , b̃h+1)∣∣ ≤ νKN · γb · ε/2, where bh+1 and b̃h+1 are defined in (D.30) and (D.32), respectively.
Proof. See §F.5 for a detailed proof.
Combining (D.34) and Lemma D.5, we know that if J2(KN , bh)− J2(KN , bKN ) ≥ ε, it holds with probability at least 1− ε15 that
J2(KN , bh+1)− J2(KN , bh) ≤ J2(KN , b̃h+1)− J2(KN , bh) + ∣∣J2(KN , bh+1)− J2(KN , b̃h+1)∣∣
≤ −νKN · γb · ε/2 < 0. (D.35)
Combining (D.29) and (D.35), it holds with probability at least 1− ε15 that J2(KN , bh+1) ≤ J2(KN , bh) ≤ · · · ≤ J2(KN , b0).
Following from a union bound argument and the choice of H in Theorem B.4, if J2(KN , bh) − J2(KN , b
KN ) ≥ ε holds for any h ≤ H , we have J2(KN , bH) ≤ J2(KN , bH−1) ≤ · · · ≤ J2(KN , b0),
which holds with probability at least 1− ε13. Thus, we finish the proof of (D.28). Combining (D.34) and Lemma D.5, for J2(KN , bh)− J2(KN , bKN ) ≥ ε/2, we have
J2(KN , bh+1)− J2(KN , bKN ) ≤ (1− νKN · γb) · [ J2(KN , bh)− J2(KN , bKN ) ] ,
which holds with probability at least 1 − ε13. Meanwhile, following from a union bound argument and the choice of H in Theorem B.4, it holds with probability at least 1− ε11 that
J2(KN , bH)− J2(KN , bKN ) ≤ ε/2. (D.36)
By Proposition 3.3 and (D.36), it holds with probability at least 1− ε11 that
‖bH − bKN ‖2 ≤ (2ε/νK∗)1/2. (D.37) Following from Proposition 3.4, we know that
bKN − b∗ = (KN −K∗)Q−1(I −A)> (D.38) · [ (I −A)Q−1(I −A)> +BR−1B> ]−1 · (Aµ+ d). Combining (D.27), (D.37), and (D.38), it holds with probability 1− ε10 that
‖bH − bKN ‖2 ≤Mb · ε1/2, where
Mb(µ) = 4 ∥∥∥Q−1(I −A)> · [(I −A)Q−1(I −A)> +BR−1B>]−1 · (Aµ+ d)∥∥∥
2 · [ ν−1K∗ + σ −1 min(Ψ ) · σ −1 min(R) ]1/2 .
We finish the proof of the theorem.

heading: D.3 PROOF OF THEOREM B.8
text: Proof. We follow the proof of Theorem 4.2 in Yang et al. (2019), where they only consider LQR without drift terms. Since our proof requires much more delicate analysis, we present it here.
Part 1. We denote by ζ̂ and ξ̂ the primal and dual variables generated by Algorithm 3. We define the primal-dual gap of (B.11) as
gap(ζ̂, ξ̂) = max ξ∈Vξ F (ζ̂, ξ)− min ζ∈Vζ F (ζ, ξ̂). (D.39)
In the sequel, we upper bound ‖α̂K,b − αK,b‖2 using (D.39). We define ζK,b and ξ(ζ) as
ζK,b = ( J(K, b), α>K,b )> , ξ(ζ) = argmax
ξ F (ζ, ξ). (D.40)
Following from (B.12), we know that ξ1(ζ) = ζ1 − J(K, b), ξ2(ζ) = EπK,b [ ψ(x, u) ] ζ1 + ΘK,bζ 2 − EπK,b [ c(x, u)ψ(x, u) ] . (D.41)
The following lemma shows that ζK,b ∈ Vζ and ξ(ζ) ∈ Vξ for any ζ ∈ Vζ . Lemma D.6. Under the assumptions in Theorem B.8, it holds that ζK,b = (J(K, b), α>K,b)> ∈ Vζ . Also, for any ζ ∈ Vζ , the vector ξ(ζ) defined in (D.40) satisfies that ξ(ζ) ∈ Vξ.
Proof. See §F.6 for a detailed proof.
By (B.12), we know that ∇ζF (ζK,b, 0) = 0 and ∇ξF (ζK,b, 0) = 0. Combining Lemma D.6, it holds that (ζK,b, 0) is a saddle point of the function F (ζ, ξ) defined in (B.11).
Following from (D.39), it holds that∥∥∥EπK,b[ψ(x, u)]ζ̂1 + ΘK,bζ̂2 − EπK,b[c(x, u)ψ(x, u)]∥∥∥2 2 + ∣∣ζ̂1 − J(K, b)∣∣2
= F ( ζ̂, ξ(ζ̂) ) = max
ξ∈Vξ F (ζ̂, ξ) = gap(ζ̂, ξ̂) + min ζ∈Vζ F (ζ, ξ̂), (D.42)
where the first equality comes from (D.41), and the second equality comes from the fact that ξ(ζ̂) = argmaxξ∈Vξ F (ζ̂, ξ) by (D.40) and Lemma D.6. We upper bound the RHS of (D.42) and lower bound the LHS of (D.42) in the sequel.
As for the RHS of (D.42), it holds for any ξ ∈ Vξ that
min ζ∈Vζ F (ζ, ξ) ≤ min ζ∈Vζ max ξ∈Vξ F (ζ, ξ) = min ζ∈Vζ
F ( ζ, ξ(ζ) ) = 1
2 min ζ∈Vζ {∥∥∥EπK,b[ψ(x, u)]ζ1 + ΘK,bζ2 − EπK,b[c(x, u)ψ(x, u)]∥∥∥2 2 + ∣∣ζ1 − J(K, b)∣∣2}
= 0, (D.43)
where the first equality comes from the fact that ξ(ζ) = argmaxξ∈Vξ F (ζ, ξ) by (D.40) and Lemma D.6, the second equality comes from (D.41), and the last equality holds by taking ζ = ζK,b ∈ Vζ . Meanwhile, we lower bound the LHS of (D.42) as∥∥∥EπK,b[ψ(x, u)]ζ̂1 + ΘK,bζ̂2 − EπK,b[c(x, u)ψ(x, u)]∥∥∥2 2 + ∣∣ζ̂1 − J(K, b)∣∣2
= ∥∥Θ̃K,b(ζ̂ − ζK,b)∥∥22 ≥ λ2K · ‖ζ̂ − ζK,b‖22 ≥ λ2K · ‖α̂K,b − αK,b‖22, (D.44)
where the first equality comes from the definition of Θ̃K,b in (B.9), and the first inequality comes from Proposition B.6. Here λK is defined in Proposition B.6. Combining (D.42), (D.43), and (D.44), it holds that
‖α̂K,b − αK,b‖22 ≤ λ−2K · gap(ζ̂, ξ̂), (D.45)
which finishes the proof of this part.
Part 2. We now upper bound gap(ζ̂, ξ̂). We denote by z̃t = (x̃>t , ũ>t )> for t ∈ [T̃ ], where x̃t and ũt are generated in Line 6 of Algorithm 3. Following from the state transition in Problem 2.3 and the form of the linear policy, {z̃t}t∈[T̃ ] follows the following transition,
z̃t+1 = Lz̃t + ν + δt, (D.46)
where
ν =
( Aµ+ d
−K(Aµ+ d) + b
) , δt = ( ωt
−Kωt + ση
) , L = ( A B −KA −KB ) .
Note that we have
L = ( A B −KA −KB ) = ( I −K ) (A B) .
Then by the property of spectral radius, it holds that
ρ(L) = ρ ( (A B) ( I −K )) = ρ(A−BK) < 1.
Thus, the Markov chain generated by (D.46) admits a unique stationary distribution N (µz,Σz), where
µz = (I − L)−1ν, Σz = LΣzL> + (
Ψω −ΨωK> −KΨω KΨωK> + σ2I
) . (D.47)
The following lemma characterizes the average
µ̂z = 1/T̃ · T̃∑ t=1 z̃t. (D.48)
Lemma D.7. It holds that
µ̂z ∼ N ( µz + 1
T̃ µT̃ ,
1 T̃ Σ̃T̃
) ,
where ‖µT̃ ‖2 ≤Mµ · (1− ρ) −2 · ‖µz‖2 and ‖Σ̃T̃ ‖F ≤MΣ · (1− ρ) −1 · ‖Σz‖F. Here Mµ and MΣ are positive absolute constants. Moreover, it holds with probability at least 1− T̃−6 that
‖µ̂z − µz‖2 ≤ log T̃
T̃ 1/4 · (1− ρ)−2 · poly
( ‖ΦK‖2, ‖K‖F, ‖b‖2, ‖µ‖2 ) .
Proof. See §F.7 for a detailed proof.
Lemma D.7 gives that
‖µ̂K,b − µK,b‖2 ≤ log T̃
T̃ 1/4 · (1− ρ)−2 · poly
( ‖ΦK‖2, ‖K‖F, ‖b‖2, ‖µ‖2 ) ,
which holds with probability at least 1− T̃−6.
We now apply a truncation argument to show that gap(ζ̂, ξ̂) is upper bounded. We define the event E in the sequel. Following from Lemma D.7, it holds for any z ∼ N (µz,Σz) that
z − µ̂z + 1/T̃ · µT̃ ∼ N (0,Σz + 1/T̃ · Σ̃T̃ ).
By Lemma G.3, there exists a positive absolute constant C0 such that P [∣∣‖z − µ̂z + 1/T̃ · µT̃ ‖22 − tr(Σ̃z)∣∣ > τ] ≤ 2 exp[−C0 ·min(τ2‖Σ̃z‖−2F , τ‖Σ̃z‖−12 )], (D.49)
where we write Σ̃z = Σz + 1/T̃ · Σ̃T̃ for notational convenience. By taking τ = C1 · log T · ‖Σ̃z‖F in (D.49) for a sufficiently large positive absolute constant C1, it holds that
P [∣∣‖z − µ̂z + 1/T̃ · µT̃ ‖22 − tr(Σ̃z)∣∣ > C1 · log T · ‖Σ̃z‖F] ≤ T−6. (D.50)
We define the event Et,1 for any t ∈ [T ] as Et,1 = {∣∣‖zt − µ̂z + 1/T̃ · µT̃ ‖22 − tr(Σ̃z)∣∣ ≤ C1 · log T · ‖Σ̃z‖F}.
Then by (D.50), it holds for any t ∈ [T ] that
P(Et,1) ≥ 1− T−6. (D.51) Also, we define
E1 = ⋂ t∈[T ] Et,1. (D.52)
Following from a union bound argument and (D.51), it holds that
P(E1) ≥ 1− T−5. (D.53)
Also, conditioning on E1, it holds for sufficiently large T̃ that max t∈[T ] ‖zt − µ̂z‖22
≤ C1 · log T · ‖Σ̃z‖F + tr(Σ̃z) + ‖1/T̃ · µT̃ ‖ 2 2 ≤ 2C̃1 · [ 1 +MΣ(1− ρ)−1/T̃ 2 ] · log T · ‖Σz‖2 +Mµ(1− ρ)−2/T̃ 2 · ‖µz‖22
≤ C2 · log T · ( 1 + ‖K‖2F ) · ‖ΦK‖2 · (1− ρ)−1 + C3 · ( ‖b‖22 + ‖µ‖22 ) · (1− ρ)−4 · T̃−2
≤ 2C2 · log T · ( 1 + ‖K‖2F ) · ‖ΦK‖2 · (1− ρ)−1, (D.54)
where C̃1, C2, and C3 are positive absolute constants. Here, the first inequality comes from the definition of E1 in (D.52), the second inequality comes from Lemma D.7, and the third inequality comes from (D.47). Also, we define the following event
E2 = { ‖µ̂z − µz + 1/T̃ · µT̃ ‖2 ≤ C1 } . (D.55)
Then by Lemma D.7, we know that
P(E2) ≥ 1− T̃−6 (D.56)
for T̃ sufficiently large. We define the event E as E = E1 ⋂ E2.
Then following from (D.53), (D.56), and a union bound argument, we know that
P(E) ≥ 1− T−5 − T̃−6.
Now, we define the truncated feature vector ψ̃(x, u) as ψ̃(x, u) = ψ̂(x, u)1E , the truncated cost function c̃(x, u) as c̃(x, u) = c(x, u)1E , and also the truncated objective function F̃ (ζ, ξ) as
F̃ (ζ, ξ) = { E(ψ̃)ζ1 + E [ (ψ̃ − ψ̃′)ψ̃> ] ζ2 − E(c̃ψ̃) }> ξ2 + [ ζ1 − E(c̃) ] · ξ1 − ‖ξ‖22/2, (D.57)
where we write ψ̃ = ψ̃(x, u) and c̃ = c̃(x, u) for notational convenience. Here the expectation is taken following the policy πK,b and the state transition. The following lemma establishes the upper bound of |F (ζ, ξ) − F̃ (ζ, ξ)|, where F (ζ, ξ) and F̃ (ζ, ξ) are defined in (B.11) and (D.57), respectively.
Lemma D.8. It holds with probability at least 1− T̃−6 that∣∣F (ζ, ξ)− F̃ (ζ, ξ)∣∣ ≤ ( 1 2T + log T̃
T̃ 1/4
) · (1− ρ)−2 · poly ( ‖K‖F, ‖b‖2, ‖µ‖2, J(K0, b0) ) .
Proof. See §F.8 for a detailed proof.
Following from (D.39) and Lemma D.8, it holds with probability at least 1− T̃−6 that∣∣gap(ζ̂, ξ̂)− g̃ap(ζ̂, ξ̂)∣∣ ≤ ( 1
2T +
log T̃
T̃ 1/4
) · (1− ρ)−2 · poly ( ‖K‖F, ‖b‖2, ‖µ‖2, J(K0, b0) ) . (D.58)
where we define g̃ap(ζ̂, ξ̂) as
g̃ap(ζ̂, ξ̂) = max ξ∈Vξ F̃ (ζ̂, ξ)− min ζ∈Vζ F̃ (ζ, ξ̂).
Therefore, to upper bound of gap(ζ, ξ), we only need to upper bound g̃ap(ζ, ξ).
Part 3. We upper bound g̃ap(ζ, ξ) in the sequel. We first show that the trajectory generated by the policy πK,b and the state transition in Problem 2.2 is β-mixing. Lemma D.9. Consider a linear system yt+1 = Dyt + ϑ + υt, where {yt}t≥0 ⊂ Rm, the matrix D ∈ Rm×m satisfying ρ(D) < 1, the vector ϑ ∈ Rm, and υt ∼ N (0,Σ) is the Gaussians. We denote by $t the marginal distribution of yt for any t ≥ 0. Meanwhile, assume that the stationary distribution of {yt}t≥0 is a Gaussian distributionN ((I−D)−1ϑ,Σ∞), where Σ∞ is the covariance matrix. We define the β-mixing coefficients for any n ≥ 1 as follows
β(n) = sup t≥0
Ey∼$t [∥∥Pyn(· | y0 = y)− PN ((I−D)−1ϑ,Σ∞)(·)∥∥TV].
Then, for any ρ ∈ (ρ(D), 1), the β-mixing coefficients satisfy that
β(n) ≤ Cρ,D,ϑ · [ tr(Σ∞) +m · (1− ρ)−2 ]1/2 · ρn, where Cρ,D,ϑ is a constant, which only depends on ρ, D, and ϑ. We say that the sequence {yt}t≥0 is β-mixing with parameter ρ.
Proof. See Proposition 3.1 in Tu and Recht (2017) for details.
Recall that by (3.1), the sequence {xt}t≥0 follows
xt+1 = (A−BK)xt + (Bb+Aµ+ d) + t, t ∼ N (0,Ψ ),
where the matrix A−BK satisfies that ρ(A−BK) < 1. Therefore, by Lemma D.9, the sequence {zt}t≥0 is β-mixing with parameter ρ ∈ (ρ(A − BK), 1), where zt = (x>t , u>t )>. The following lemma upper bounds the primal-dual gap for a convex-concave problem.
Lemma D.10. Let X and Y be two compact and convex sets such that ‖x − x′‖2 ≤ M and ‖y − y′‖2 ≤M for any x, x′ ∈ X and y, y′ ∈ Y . We consider solving the following minimax problem
min x∈X max y∈Y
H(x, y) = E ∼$ [ G(x, y; ) ] ,
where the objective function H(x, y) is convex in x and concave in y. In addition, we assume that the distribution$ is β-mixing with β(n) ≤ C ·ρn, where C is a constant. Meanwhile, we assume that it holds almost surely that G(x, y; ) is L̃0-Lipschitz in both x and y, the gradient∇xG(x, y; ) is L̃1-Lipschitz in x for any y ∈ Y , the gradient ∇yG(x, y; ) is L̃1-Lipschitz in y for any x ∈ X , where C , L̃0, L̃1 > 1. Each step of our gradient-based method takes the following forms,
xt+1 = ΓX [ xt − γt+1 · ∇xG(xt, yt; t) ] , yt+1 = ΓY [ yt − γt+1 · ∇yG(xt, yt; t) ] ,
where the operators ΓX and ΓY projects the variables back to X and Y , respectively, and the stepsizes take the form γt = γ0 · t−1/2 for a constant γ0 > 0. Moreover, let x̂ = ( ∑T t=1 γt) −1( ∑T t=1 γtxt) and ŷ = ( ∑T t=1 γt) −1( ∑T t=1 γtyt) be the final output of the gradient method after T iterations, then there exists a positive absolute constant C, such that for any δ ∈ (0, 1), the primal-dual gap to the minimax problem is upper bounded as
max x∈X H(x̂, y)−min y∈Y
H(x, ŷ) ≤ C · (M 2 + L̃20 + L̃0L̃1M) log(1/ρ) · log 2 T + log(1/δ)√ T + C · C L̃0M T ,
which holds with probability at least 1− δ.
Proof. See Theorem 5.4 in Yang et al. (2019) for details.
To use Lemma D.10, we define the function G(ζ, ξ; ψ̃, ψ̃′) as
G(ζ, ξ; ψ̃, ψ̃′) = [ ψ̃ζ1 + (ψ̃ − ψ̃′)ψ̃>ζ2 − c̃ψ̃ ]> ξ2 + (ζ1 − c̃) · ξ1 − 1/2 · ‖ξ‖22,
where ψ̃ = ψ̃(x, u) and ψ̃′ = ψ̃(x′, u′). Note that the gradients of G(ζ, ξ; ψ̃, ψ̃′) take the form ∇ζG(ζ, ξ; ψ̃, ψ̃′) = ( ψ̃>ξ2 + ξ1
ψ̃(ψ̃ − ψ̃′)>ξ2
) , ∇ξG(ζ, ξ; ψ̃, ψ̃′) = ( ζ1 − c̃− ξ1 ψ̃ζ1 + (ψ̃ − ψ̃′)ψ̃>ζ2 − c̃ψ̃ − ξ2 ) .
By Definition B.7 and Lemma D.6, we know that∥∥∇ζG(ζ, ξ; ψ̃, ψ̃′)∥∥2 ≤ poly(‖K‖F, J(K0, b0)) · log2 T · (1− ρ)−2,∥∥∇ξG(ζ, ξ; ψ̃, ψ̃′)∥∥2 ≤ poly(‖K‖F, ‖µ‖2, J(K0, b0)) · log2 T · (1− ρ)−2. (D.59) This gives the Lipschitz constant L̃0 in Lemma D.10 for G(ζ, ξ; ψ̃, ψ̃′). Also, the Hessians of G(ζ, ξ; ψ̃, ψ̃′) take the forms of
∇2ζζG(ζ, ξ; ψ̃, ψ̃′) = 0, ∇2ξξG(ζ, ξ; ψ̃, ψ̃′) = −I,
which follows that ∥∥∇2ζζG(ζ, ξ; ψ̃, ψ̃′)∥∥2 = 0, ∥∥∇2ξξG(ζ, ξ; ψ̃, ψ̃′)∥∥2 = 1. (D.60) This gives the Lipschitz constant L̃1 in Lemma D.10 for ∇ζG(ζ, ξ; ψ̃, ψ̃′) and ∇ξG(ζ, ξ; ψ̃, ψ̃′). Moreover, note that (D.54) provides an upper bound of M , combining (D.59), (D.60) and Lemma D.10, it holds with probability at least 1− T−5 that
g̃ap(ζ̂, ξ̂) ≤ poly
( ‖K‖F, ‖µ‖2, J(K0, b0) ) · log6 T
(1− ρ)4 · √ T
. (D.61)
Combining (D.45), (D.58), and (D.61), we know that ‖α̂K,b − αK,b‖22 ≤ λ−2K · poly1 ( ‖K‖F, ‖b‖2, ‖µ‖2, J(K0, b0) ) · [
log6 T
T 1/2 · (1− ρ)4 +
log T̃
T̃ 1/4 · (1− ρ)2
] .
Same bounds for ‖Υ̂K − ΥK‖2F, ‖p̂K,b − pK,b‖22, and ‖q̂K,b − qK,b‖22 hold. We finish the proof of the theorem.

heading: E PROOFS OF PROPOSITIONS

heading: E.1 PROOF OF PROPOSITION 3.2
text: Proof. We follow a similar proof as in the one of Theorem 1.1 in Sznitman (1991) and Theorem 3.2 in Bensoussan et al. (2016). Note that for any policy πK,b ∈ Π, the parameters K and b uniquely determine the policy. We define the following metric on Π. Definition E.1. For any πK1,b1 , πK2,b2 ∈ Π, we define the following metric,
‖πK1,b1 − πK2,b2‖2 = c1 · ‖K1 −K2‖2 + c2 · ‖b1 − b2‖2,
where c1 and c2 are positive constants.
One can verify that Definition E.1 satisfies the requirement of being a metric. We first evaluate the forms of the operators Λ1(·) and Λ2(·, ·).
Forms of the operators Λ1(·) and Λ2(·, ·). By the definition of Λ1(µ), which gives the optimal policy under the mean-field state µ, it holds that
Λ1(µ) = π ∗ µ,
where π∗µ solves Problem 2.2. This gives the form of Λ1(·). We now turn to Λ2(µ, π), which gives the mean-field state µnew generated by the policy π under the current mean-field state µ. In Problem
2.2, the sequence of states {xt}t≥0 constitutes a Markov chain, which admits a unique stationary distribution. Thus, by the state transition in Problem 2.2 and the form of the linear-Gaussian policy, we have
µnew = (A−BKπ)µnew + (Bbπ +Aµ+ d), (E.1)
where Kπ and bπ are parameters of the policy π. By solving (E.1) for µnew, it holds that
Λ2(µ, π) = µnew = (I −A+BKπ)−1(Bbπ +Aµ+ d).
This gives the form of Λ2(·, ·).
Next, we compute the Lipschitz constants for Λ1(·) and Λ2(·, ·).
Lipschitz constant for Λ1(·). By Proposition 3.4, for any µ1, µ2 ∈ Rm, the optimal K∗ is fixed for Problem 2.2. Therefore, by the form of the optimal bK given in Proposition 3.4, it holds that∥∥Λ1(µ1)− Λ1(µ2)∥∥2 ≤ c2 · ∥∥∥[(I −A)Q−1(I −A)> +BR−1B>]−1A∥∥∥2
· ∥∥∥[K∗Q−1(I −A)> −R−1B>]∥∥∥
2 · ‖µ1 − µ2‖2
= c2L1 · ‖µ1 − µ2‖2, (E.2)
where L1 is defined in Assumption 3.1.
Lipschitz constants for Λ2(·, ·). By Proposition 3.4, for any µ1, µ2 ∈ Rm, the optimal K∗ is fixed for Problem 2.2. Thus, for any π ∈ Π such that π is an optimal policy under some µ ∈ Rm, it holds that ∥∥Λ2(µ1, π)− Λ2(µ2, π)∥∥2 = ∥∥(I −A+BKπ)−1 ·A · (µ1 − µ2)∥∥2
≤ [ 1− ρ(A−BK∗) ]−1 · ‖A‖2 · ‖µ1 − µ2‖2 = L2 · ‖µ1 − µ2‖2, (E.3)
where L2 is defined in Assumption 3.1, and Kπ = K∗ is the parameter of the policy π. Meanwhile, for any mean-field state µ ∈ Rm, and any poicies π1, π2 ∈ Π that are optimal under some mean-field states µ1, µ2, respectively, we have∥∥Λ2(µ, π1)− Λ2(µ, π2)∥∥2 = ∥∥(I −A+BK∗)−1B · (bπ1 − bπ2)∥∥2
≤ [ 1− ρ(A−BK∗) ]−1 · ‖B‖2 · ‖bπ1 − bπ2‖2 = c−12 L3 · ‖π1 − π2‖2, (E.4)
where in the last equality, we use the fact that Kπ1 = Kπ2 = K ∗ by Proposition 3.4. Here L3 is defined in Assumption 3.1, and bπ1 and bπ2 are the parameters of the policies π1 and π2.
Now we show that the operator Λ(·) is a contraction. For any µ1, µ2 ∈ Rm, it holds that∥∥Λ(µ1)− Λ(µ2)∥∥2 =∥∥∥Λ2(µ1,Λ1(µ1))− Λ2(µ2,Λ1(µ2))∥∥∥2 ≤ ∥∥∥Λ2(µ1,Λ1(µ1))− Λ2(µ1,Λ1(µ2))∥∥∥ 2 + ∥∥∥Λ2(µ1,Λ1(µ2))− Λ2(µ2,Λ1(µ2))∥∥∥ 2
≤ c−12 L3· ∥∥Λ1(µ1)− Λ1(µ2)∥∥2 + L2 · ‖µ1 − µ2‖2 ≤ c−12 L3 · c2L1 · ‖µ1 − µ2‖2 + L2 · ‖µ1 − µ2‖2 = (L1L3 + L2) · ‖µ1 − µ2‖2,
where the first inequality comes from triangle inequality, the second inequality comes from (E.3) and (E.4), and the last inequality comes from (E.2). By Assumption 3.1, we know that L0 = L1L3 + L2 < 1, which shows that the operator Λ(·) is a contraction. Moreover, by Banach fixedpoint theorem, we obtain that Λ(·) has a unique fixed point, which gives the unique equilibrium pair of Problem 2.3. We finish the proof of the proposition.

heading: E.2 PROOF OF PROPOSITION 3.4
text: Proof. By the definition of J2(K, b) in (3.6) and the definition of µK,b in (3.2), the problem min b J2(K, b) is equivalent to the following constrained problem,
min µ,b ( µ b )>( Q+K>RK −K>R −RK R )( µ b ) s.t. (I −A+BK)µ− (Bb+Aµ+ d) = 0. (E.5)
Following from the KKT conditions of (E.5), it holds that
2MK ( µ b ) +NKλ = 0, N > K ( µ b ) +Aµ+ d = 0, (E.6)
where
MK = ( Q+K>RK −K>R −RK R ) , NK = ( −(I −A+BK)> B> ) .
By solving (E.6), the minimizer of (E.5) takes the form of( µK,bK bK ) = −M−1K NK(N > KM −1 K NK) −1(Aµ+ d). (E.7)
By substituting (E.7) into the definition of J2(K, b) in (3.6), we have
J2(K, b K) = (Aµ+ d)>(N>KM −1 K NK) −1(Aµ+ d). (E.8) Meanwhile, by calculation, we have
M−1K =
( Q−1 Q−1K>
KQ−1 KQ−1K> +R−1
) .
Therefore, the term N>KM −1 K NK in (E.8) takes the form of
N>KM −1 K NK = (I −A)Q −1(I −A>) +BR−1B>. (E.9) By plugging (E.9) into (E.8), we have
J2(K, b K) = (Aµ+ d)> [ (I −A)Q−1(I −A>) +BR−1B> ]−1 (Aµ+ d).
Also, by plugging (E.9) into (E.7), we have( µK,bK bK ) = ( Q−1(I −A)> KQ−1(I −A)> −R−1B> )[ (I −A)Q−1(I −A)> +BR−1B> ]−1 (Aµ+ d).
We finish the proof of the proposition.

heading: E.3 PROOF OF PROPOSITION B.2
text: Proof. By the definition of the cost function c(x, u) in Problem 2.2 (recall that we drop the subscript µ when we focus on Problem 2.2), we have
Ect = E(x>t Qxt + u>t Rut + µ>Qµ) = E(x>t Qxt + x>t K>RKxt − 2b>RKxt + b>Rb+ σ2η>t Rηt + µ>Qµ) = E [ x>t (Q+K >RK)xt − 2b>RKxt ]
+ b>Rb+ σ2 · tr(R) + µ>Qµ, (E.10) where we write ct = c(xt, ut) for notational convenience. Here in the second line we use ut = πK,b(xt) = −Kxt+ b+σηt. Therefore, combining (E.10) and the definition of J(K, b) in Problem 2.2, we have
J(K, b) = lim T→∞
1
T T∑ t=0 { E [ x>t (Q+K >RK)xt − 2b>RKxt ] + b>Rb+ σ2 · tr(R) + µ>Qµ }
= Ex∼N (µK,b,ΦK) [ x>(Q+K>RK)x− 2b>RKx ] + b>Rb+ σ2 · tr(R) + µ>Qµ
= tr [ (Q+K>RK)ΦK ] + µ>K,b(Q+K
>RK)µK,b − 2b>RKµK,b (E.11) + b>Rb+ σ2 · tr(R) + µ>Qµ.
Now, by iteratively applying (3.3) and (3.4), we have tr [ (Q+K>RK)ΦK ] = tr(PKΨ ), (E.12) where PK is given in (3.4). Combining (E.11) and (E.12), we know that
J(K, b) = J1(K) + J2(K, b) + σ 2 · tr(R) + µ>Qµ,
where J1(K) = tr [ (Q+K>RK)ΦK ] = tr(PKΨ ),
J2(K, b) = ( µK,b b )>( Q+K>RK −K>R −RK R )( µK,b b ) .
We finish the proof of the proposition.

heading: E.4 PROOF OF PROPOSITION 3.3
text: Proof. By calculating the Hessian matrix of J2(K, b), we have
∇2bbJ2(K, b) =B>(I −A+BK)−>(Q+K>RK)(I −A+BK)−1B − [ RK(I −A+BK)−1B +B>(I −A+BK)−>K>R ] +R
= [ R1/2K(I −A+BK)−1B −R1/2 ]>[ R1/2K(I −A+BK)−1B −R1/2 ] +B>(I −A+BK)−>Q(I −A+BK)−1B,
which is a positive definite matrix independent of b. We denote by its minimum singular value as νK . Also, note that ‖∇2bbJ2(K, b)‖2 is upper bounded as∥∥∇2bbJ2(K, b)∥∥2 ≤ [1− ρ(A−BK)]−2 · (‖B‖22 · ‖K‖22 · ‖R‖2 + ‖B‖22 · ‖Q‖2). Therefore, it holds that
ιK ≤ [ 1− ρ(A−BK) ]−2 · (‖B‖22 · ‖K‖22 · ‖R‖2 + ‖B‖22 · ‖Q‖2), where ιK is the maximum singular value of∇2bbJ2(K, b). We finish the proof of the proposition.

heading: E.5 PROOF OF PROPOSITION B.3
text: Proof. Following from Proposition B.2, it holds that
J1(K) = tr(PKΨ ) = Ey∼N (0,Ψ )(y >PKy) = Ey∼N (0,Ψ ) [ fK(y) ] , (E.13)
where fK(y) = y>PKy. By the definition of PK in (3.4), we obtain that ∇KfK(y) = ∇K { y>(Q+K>RK)y + [ (A−BK)y ]> PK [ (A−BK)y ]>} = 2RKyy> +∇K [ fK ( (A−BK)y )] . (E.14)
Also, we have ∇K [ fK ( (A−BK)y )] = ∇KfK ( (A−BK)y ) − 2B>PK(A−BK)yy>. (E.15)
By plugging (E.15) into (E.14), we have ∇KfK(y) = 2 [ (R+B>PKB)K −B>PKA ] yy> +∇KfK ( (A−BK)y ) . (E.16) By iteratively applying (E.16), it holds that
∇KfK(y) = 2 [ (R+B>PKB)K −B>PKA ] · ∞∑ t=0 yty > t , (E.17)
where yt+1 = (A−BK)yt with y0 = y. Now, combining (E.13) and (E.17), it holds that ∇KJ1(K) = 2 [ (R+B>PKB)K −B>PKA ] ΦK = 2(Υ 22 KK −Υ21K ) · ΦK , where ΥK is defined in (3.7). Meanwhile, combining the form of µK,b in (3.2), it holds by calculation that
∇bJ2(K, b) = 2 [ Υ22K (−KµK,b + b) + Υ21KµK,b + qK,b ] ,
where qK,b is defined in (3.7). We finish the proof of the proposition.

heading: E.6 PROOF OF PROPOSITION B.1
text: Proof. From the definition of VK,b(x) in (B.1) and the definition of the cost function c(x, u) in Problem 2.2, it holds that
VK,b(x) = ∞∑ t=0 { E [ x>t (Q+K >RK)xt − 2b>RKxt
+ b>Rb+ σ2η>t Rηt + µ >Qµ |x0 = x
] − J(K, b) } .
Combining (3.1), we know that VK,b(x) is a quadratic function taking the form of VK,b(x) = x>Gx+ r>x+ h, where G, r, and h are functions of K and b. Note that VK,b(x) satisfies that
VK,b(x) = c(x,−Kx+ b)− J(K, b) + E [ VK,b(x ′) |x ] , (E.18)
by substituting the form of c(x,−Kx+b) in Problem 2.2 and J(K, b) in (3.5) into (E.18), we obtain that
x>Gx+ r>x+ h
= x>(Q+K>RK)x− 2b>RKx+ b>Rb+ µ>Qµ (E.19) − [ tr(PKΨ ) + µ>K,b(Q+K >RK)µK,b − 2b>RKµK,b + µ>Qµ+ b>Rb ]
+ [ (A−BK)x+ (Bb+Aµ+ d) ]> G [ (A−BK)x+ (Bb+Aµ+ d) ] + tr(GΨ ) + r> [ (A−BK)x+ (Bb+Aµ+ d) ] + h− σ2 · tr(R).
By comparing the quadratic terms and linear terms on both the LHS and RHS in (E.19), we obtain that
G = PK , r = 2fK,b,
where fK,b = (I −A+BK)−>[(A−BK)>PK(Bb+Aµ+ d)−K>Rb]. Also, by the definition of VK,b(x) in (B.1), we know that E[VK,b(x)] = 0, where the expectation is taken following the stationary distribution generated by the policy πK,b and the state transition. Therefore, we have
h = −2fK,bµK,b − µ>K,bPKµK,b − tr(PKΦK), which shows that
VK,b(x) = x >PKx− tr(PKΦK) + 2f>K,b(x− µK,b)− µ>K,bPKµK,b. (E.20)
For the action-value function QK,b(x, u), by plugging (E.20) into (B.2), we obtain that
QK,b(x, u) = ( x u )> ΥK ( x u ) + 2 ( pK,b qK,b )>( x u ) − tr(PKΦK)− σ2 · tr(R+ PKBB>)
− b>Rb+ 2b>RKµK,b − µ>K,b(Q+K>RK + PK)µK,b + 2f>K,b [ (Aµ+ d)− µK,b ] + (Aµ+ d)>PK(Aµ+ d).
We finish the proof of the proposition.

heading: E.7 PROOF OF PROPOSITION B.5
text: Proof. By Proposition B.1, it holds that QK,b takes the following linear form
QK,b(x, u) = ψ(x, u) >αK,b + βK,b, (E.21)
where βK,b is a scalar independent of x and u. Note that QK,b(x, u) satisfies that QK,b(x, u) = c(x, u)− J(K, b) + EπK,b [ QK,b(x ′, u′) |x, u ] , (E.22) where (x′, u′) is the state-action pair after (x, u) following the policy πK,b and the state transition. Combining (E.21) and (E.22), we obtain that
ψ(x, u)>αK,b = c(x, u)− J(K, b) + EπK,b [ ψ(x′, u′) |x, u ]> αK,b. (E.23)
By left multiplying ψ(x, u) to both sides of (E.23), and taking the expectation, we have EπK,b { ψ(x, u) [ ψ(x, u)− ψ(x′, u′) ]>} · αK,b + EπK,b[ψ(x, u)] · J(K, b) = EπK,b[c(x, u)ψ(x, u)]. Combining the definition of the matrix ΘK,b in (B.7), we have(
1 0 EπK,b [ ψ(x, u) ] ΘK,b )( J(K, b) αK,b ) = ( J(K, b) EπK,b [ c(x, u)ψ(x, u) ]) , which concludes the proof of the proposition.

heading: E.8 PROOF OF PROPOSITION B.6
text: Proof. Invertibility and Upper Bound. We denote by zt = (x>t , u>t )> for any t ≥ 0. Then following from the state transition and the policy πK,b, the transition of {zt}t≥0 takes the form of
zt+1 = Lzt + ν + δt, (E.24)
where L, ν and δ are defined as
L = ( A B −KA −KB ) , ν = ( Aµ+ d −K(Aµ+ d) + b ) , δt = ( ωt −Kωt + σηt ) .
Note that L also takes the form of
L = ( I −K ) (A B) .
Combining the fact that ρ(UV ) = ρ(V U) for any matrices U and V , we know that ρ(L) = ρ(A− BK) < 1, which verifies the stability of (E.24). Following from the stability of (E.24), we know that the Markov chain generated by (E.24) admits a unique stationary distributionN (µz,Σz), where µz and Σz satisfy that
µz = Lµz + ν, Σz = LΣzL > + Ψδ.
where
Ψδ = ( Ψω −ΨωK> −KΨω KΨωK> + σ2I ) .
Also, we know that Σz takes the form of
Σz = Cov [( x u )] = ( ΦK −ΦKK> −KΦK KΦKK> + σ2I ) = ( 0 0 0 σ2I ) + ( I −K ) ΦK ( I −K )> ,
(E.25)
where ΦK is defined in (3.3).
The following lemma establishes the form of ΘK,b. Lemma E.2. The matrix ΘK,b in (B.7) takes the form of
ΘK,b =
( 2(Σz ⊗s Σz)(I − L⊗s L)> 0 0 Σz(I − L)> ) .
Proof. See §F.9 for a detailed proof.
Note that since ρ(L) < 1, both I − L ⊗s L and I − L are positive definite. Therefore, by Lemma E.2, the matrix ΘK,b is invertible. This finishes the proof of the invertibility of ΘK,b. Moreover, by (E.25) and Lemma E.2, we upper bound the spectral norm of ΘK,b as
‖ΘK,b‖2 ≤ 2 max { ‖Σz‖22 · ( 1 + ‖L‖22 ) , ‖Σz‖2 · ( 1 + ‖L‖2 )} ≤ 4 ( 1 + ‖K‖2F )2 · ‖ΦK‖22, which proves the upper bound of ‖ΘK,b‖2.
Minimum singular value. To lower bound σmin(Θ̃K,b), we only need to upper bound σmax(Θ̃ −1 K,b) = ‖Θ̃ −1 K,b‖2. We first calculate Θ̃ −1 K,b. Recall that the matrix Θ̃K,b in (B.8) takes the form of
Θ̃K,b =
( 1 0
EπK,b [ ψ(x, u) ] ΘK,b
) .
By the definition of the feature vector ψ(x, u) in (B.5), the vector σ̃z = EπK,b [ψ(x, u)] takes the form of
σ̃z = EπK,b [ ψ(x, u) ] = ( svec(Σz) 0k+m ) ,
where 0k+m denotes the all-zero column vector with dimension k +m. Also, since ΘK,b is invertible, the matrix Θ̃K,b is also invertible, whose inverse takes the form of
Θ̃−1K,b =
( 1 0
−Θ−1K,b · σ̃z Θ −1 K,b
) .
The following lemma upper bounds the spectral norm of Θ̃−1K,b. Lemma E.3. The spectral norm of the matrix Θ̃−1K,b is upper bounded by a positive constant λ̃K , where λ̃K only depends on ‖K‖2 and ρ(A−BK).
Proof. See §F.10 for a detailed proof.
By Lemma E.3, we know that σmin(Θ̃K,b) is lower bounded by a positive constant λK = 1/λ̃K , which only depends on ‖K‖2 and ρ(A−BK). This concludes the proof of the proposition.

heading: F PROOFS OF LEMMAS

heading: F.1 PROOF OF LEMMA D.1
text: Proof. Following from (3.4), it holds that y>PK2y = ∑ t≥0 y> [ (A−BK2)t ]> (Q+K>2 RK2)(A−BK2)ty. (F.1) Meanwhile, by the state transition yt+1 = (A−BK2)yt, we know that yt = (A−BK2)ty0 = (A−BK2)ty. (F.2) By plugging (F.2) into (F.1), it holds that
y>PK2y = ∑ t≥0 y>t (Q+K > 2 RK2)yt = ∑ t≥0 (y>t Qyt + y > t K > 2 RK2yt). (F.3)
Also, it holds that
y>PK1y = ∑ t≥0 (y>t+1PK1yt+1 − y>t PK1yt) (F.4)
Combining (F.3) and (F.4), we have y>PK2y − y>PK1y = ∑ t≥0 (y>t Qyt + y > t K > 2 RK2yt + y > t+1PK1yt+1 − y>t PK1yt). (F.5) Also, by the state transition yt+1 = (A−BK2)yt, it holds for any t ≥ 0 that y>t Qyt + y > t K > 2 RK2yt + y > t+1PK1yt+1 − y>t PK1yt
= y>t [ Q+ (K2 −K1 +K1)>R(K2 −K1 +K1) ] yt
+ y>t [ A−BK1 −B(K2 −K1) ]> PK1 [ A−BK1 −B(K2 −K1) ] yt − y>t PK1yt
= 2y>t (K2 −K1)> [ (R+B>PK1B)K1 −B>PK1A ] yt
+ y>t (K2 −K1)>(R+B>PK1B)(K2 −K1)yt = 2y>t (K2 −K1)>(Υ22K1K1 −Υ 21 K1)yt + y > t (K2 −K1)>Υ22K1(K2 −K1)yt, (F.6)
where the matrix ΥK1 is defined in (3.7). By plugging (F.6) into (F.5), we have
y>PK2y − y>PK1y = ∑ t≥0 2y>t (K2 −K1)>(Υ22K1K1 −Υ 21 K1)yt + y > t (K2 −K1)>Υ22K1(K2 −K1)yt
= ∑ t≥0 DK1,K2(yt),
where DK1,K2(y) = 2y >(K2 − K1)(Υ22K1K1 − Υ 21 K1 )y + y>(K2 − K1)>Υ22K1(K2 − K1)y. We finish the proof of the lemma.

heading: F.2 PROOF OF LEMMA D.2
text: Proof. We prove (D.16) and (D.17) separately in the sequel.
Proof of (D.16). From the definition of J1(K) in (3.6), we have
J1(K)− J1(K∗) = tr(PKΨ − PK∗Ψ ) = Ey∼N (0,Ψ )(y >PKy − y>PK∗y) = −E [∑ t≥0 DK,K∗(yt) ] , (F.7)
where in the last equality, we apply Lemma D.1 and the expectation is taken following the transition yt+1 = (A−BK∗)yt with initial state y0 ∼ N (0,Ψ ). Here we denote by DK,K∗(y) as
DK,K∗(y) = 2y >(K∗ −K)(Υ22KK −Υ21K )y + y>(K∗ −K)>Υ22K (K∗ −K)y.
Also, we write DK,K∗(y) as
DK,K∗(y) = 2y >(K∗ −K)(Υ22KK −Υ21K )y + y>(K∗ −K)>Υ22K (K∗ −K)y (F.8) = y> [ K∗ −K + (Υ22K )−1(Υ22KK −Υ21K ) ]> Υ22K [ K∗ −K + (Υ22K )−1(Υ22KK −Υ21K ) ] y − y>(Υ22KK −Υ21K )>(Υ22K )−1(Υ22KK −Υ21K )y. Note that the first term on the RHS of (F.8) is positive, due to the fact that it is a quadratic form of a positive definite matrix, we lower bound DK,K∗(y) as
DK,K∗(y) ≥ −y>(Υ22KK −Υ21K )>(Υ22K )−1(Υ22KK −Υ21K )y. (F.9) Combining (F.7) and (F.9), it holds that
J1(K)− J1(K∗) ≤ ∥∥∥∥E(∑
t≥0
yty > t )∥∥∥∥ 2 · tr [ (Υ22KK −Υ21K )>(Υ22K )−1(Υ22KK −Υ21K ) ] = ‖ΦK∗‖2 · tr [ (Υ22KK −Υ21K )>(Υ22K )−1(Υ22KK −Υ21K )
] ≤ ∥∥(Υ22K )−1∥∥2 · ‖ΦK∗‖2 · tr[(Υ22KK −Υ21K )>(Υ22KK −Υ21K )]
≤ σ−1min(R) · ‖ΦK∗‖2 · tr [ (Υ22KK −Υ21K )>(Υ22KK −Υ21K ) ] ,
where the last line comes from the fact that Υ22K = R + B >PKB R. This complete the proof of (D.16).
Proof of (D.17). Note that for any K̃, it holds by the optimality of K∗ that J1(K)− J1(K∗) ≥ J1(K)− J1(K̃) = −E [∑ t≥0 DK,K̃(yt) ] , (F.10)
where the expectation is taken following the transition yt+1 = (A − BK̃)yt with initial state y0 ∼ N (0,Ψ ). By taking K̃ = K − (Υ22K )−1(Υ22KK −Υ21K ) and following from a similar calculation as in (F.8), the function DK,K̃(y) takes the form of
DK,K̃(y) = −y >(Υ22KK −Υ21K )>(Υ22K )−1(Υ22KK −Υ21K )y. (F.11)
Combining (F.10) and (F.11), it holds that J(K)− J(K∗) ≥ tr [ ΦK̃(Υ 22 KK −Υ21K )>(Υ22K )−1(Υ22KK −Υ21K ) ] ≥ σmin(Ψ ) · ‖Υ22K ‖−12 · tr [ (Υ22KK −Υ21K )>(Υ22KK −Υ21K ) ] ,
where we use the fact that ΦK̃ = (A − BK̃)ΦK̃(A − BK̃) > + Ψ Ψ in the last line. This finishes the proof of (D.17).

heading: F.3 PROOF OF LEMMA D.3
text: Proof. By Proposition B.2, we have∣∣J1(K̃n+1)− J1(Kn+1)∣∣ = tr[(PK̃n+1 − PKn+1)Ψ ] ≤ ‖PK̃n+1 − PKn+1‖2 · ‖Ψ ‖F. (F.12) The following lemma upper bounds the term ‖PK̃n+1 − PKn+1‖2.
Lemma F.1. Suppose that the parameters K and K̃ satisfy that ‖K̃ −K‖2 · ( ‖A−BK‖2 + 1 ) · ‖ΦK‖2 ≤ σmin(Ψω)/4 · ‖B‖−12 , (F.13)
then it holds that
‖PK̃ − PK‖2 ≤ 6 · σ −1 min(Ψω) · ‖ΦK‖2 · ‖K‖2 · ‖R‖2 · ‖K̃ −K‖2 (F.14) · ( ‖B‖2 · ‖K‖2) · ‖A−BK‖2 + ‖B‖2 · ‖K‖2 + 1 ) .
Proof. See Lemma 5.7 in Yang et al. (2019) for a detailed proof.
To use Lemma F.1, it suffices to verify that K̃n+1 and Kn+1 satisfy (F.13). Note that from the definitions of Kn+1 and K̃n+1 in (D.18) and (D.21), respectively, we have
‖K̃n+1 −Kn+1‖2 · ( ‖A−BK̃n+1‖2 + 1 ) · ‖ΦK̃n+1‖2
≤ γ · ‖Υ̂Kn −ΥKn‖F · ( 1 + ‖Kn‖2 ) · ( ‖A−BK̃n+1‖2 + 1 ) · ‖ΦK̃n+1‖2. (F.15)
Now, we upper bound the RHS of (F.15). For the term ‖A−BK̃n+1‖2, it holds by the definition of K̃n+1 in (D.21) that
‖A−BK̃n+1‖2 ≤ ‖A−BKn‖2 + γ · ‖B‖2 · ‖Υ22KnKn −Υ 21 Kn‖2 ≤ ‖A−BKn‖2 + γ · ‖B‖2 · ‖ΥKn‖2 · ( 1 + ‖Kn‖2 ) . (F.16)
By the definition of ΥKn in (3.7), we upper bound ‖ΥKn‖2 as ‖ΥKn‖2 ≤ ‖Q‖2 + ‖R‖2 + ( ‖A‖F + ‖B‖F )2 · ‖PKn‖2 ≤ ‖Q‖2 + ‖R‖2 + ( ‖A‖F + ‖B‖F
)2 · J1(K0) · σ−1min(Ψ ), (F.17) where the last line comes from the fact that
J1(K0) ≥ J1(Kn) = tr [ (Q+K>n RKn)ΦKn ] = tr(PKnΨ ) ≥ ‖PKn‖2 · σmin(Ψ ).
As for the term ‖ΦK̃n+1‖2 in (F.15), from the fact that
J1(K0) ≥ J1(K̃n+1) = tr [ (Q+ K̃>n+1RK̃n+1)ΦK̃n+1 ] ≥ ‖ΦK̃n+1‖2 · σmin(Q),
it holds that
‖ΦK̃n+1‖2 ≤ J1(K0) · σ −1 min(Q). (F.18)
Therefore, combining (F.15), (F.16), (F.17), and (F.18), we know that ‖K̃n+1 −Kn+1‖2 · ( ‖A−BK̃n+1‖2 + 1 ) · ‖ΦK̃n+1‖2
≤ poly1 ( ‖Kn‖2 ) · ‖Υ̂Kn −ΥKn‖F. (F.19)
From Theorem B.8, it holds with probability at least 1− T−4n − T̃−6n that
‖Υ̂Kn −ΥKn‖F ≤ poly3
( ‖Kn‖F, ‖µ‖2 ) λKn · (1− ρ)2 · log 3 Tn T 1/4 n
(F.20)
+ poly4
( ‖Kn‖F, ‖b0‖2, ‖µ‖2 ) λKn · log 1/2 T̃n T̃ 1/8 n · (1− ρ) ,
which holds for any ρ ∈ (ρ(A−BKn), 1). Note that from the choice of Tn and T̃n in the statement of Theorem B.4 that
Tn ≥ poly5 ( ‖Kn‖F, ‖b0‖2, ‖µ‖2 ) · λ−4Kn · [ 1− ρ(A−BKn) ]−9 · ε−5, T̃n ≥ poly6 ( ‖Kn‖F, ‖b0‖2, ‖µ‖2 ) · λ−2Kn · [ 1− ρ(A−BKn)
]−12 · ε−12, it holds that
poly3 ( ‖Kn‖F, ‖µ‖2 ) λKn · (1− ρ)2 · log 3 Tn T 1/4 n + poly4 ( ‖Kn‖F, ‖b0‖2, ‖µ‖2 ) λKn · log 1/2 T̃n T̃ 1/8 n · (1− ρ)
≤ min {[ poly1 ( ‖Kn‖2 )]−1 · σmin(Ψω)/4 · ‖B‖−12 , (F.21)[
poly2 ( ‖Kn‖2 )]−1 · ε/8 · γ · σmin(Ψ ) · σmin(R) · ‖ΦK∗‖−12 · ‖Ψ ‖ −1 F } .
Combining (F.19), (F.20), and (F.21), we know that (F.13) holds with probability at least 1 − ε15 for sufficiently small ε > 0. Meanwhile, by (F.16), (F.17), and (F.18), the RHS of (F.14) is upper bounded as
6 · σ−1min(Ψω) · ‖ΦK̃n+1‖2 · ‖K̃n+1‖2 · ‖R‖2 · ‖K̃n+1 −Kn+1‖2 · ( ‖B‖2 · ‖K̃n+1‖2) · ‖A−BK̃n+1‖2 + ‖B‖2 · ‖K̃n+1‖2 + 1 ) ≤ poly2 ( ‖Kn‖2 ) · ‖Υ̂Kn −ΥKn‖F. (F.22)
Now, by Lemma F.1, it holds with probability at least 1− ε15 that
‖PK̃n+1 − PKn+1‖2 ≤ 6 · σ −1 min(Ψω) · ‖ΦK̃n+1‖2 · ‖K̃n+1‖2 · ‖R‖2 · ‖K̃n+1 −Kn+1‖2 · ( ‖B‖2 · ‖K̃n+1‖2) · ‖A−BK̃n+1‖2 + ‖B‖2 · ‖K̃n+1‖2 + 1 ) ≤ poly2 ( ‖Kn‖2 ) · ‖Υ̂Kn −ΥKn‖F
≤ ε/8 · γ · σmin(Ψ ) · σmin(R) · ‖ΦK∗‖−12 · ‖Ψ ‖ −1 F , (F.23)
where the second inequality comes from (F.22), and the last inequality comes from (F.20) and (F.21). Combining (F.12) and (F.23), it holds with probability at least 1− ε15 that∣∣J1(K̃n+1)− J1(Kn+1)∣∣ ≤ γ · σmin(Ψ ) · σmin(R) · ‖ΦK∗‖−12 · ε/4, which concludes the proof of the lemma.

heading: F.4 PROOF OF LEMMA D.4
text: Proof. Note that Υ22K∗K ∗ − Υ21K∗ is the natural gradient of J1 at the minimizer K∗, which implies that
Υ22K∗K ∗ −Υ21K∗ = 0. (F.24)
By Lemma D.1, it holds that
J1(K)− J1(K∗) = tr(PKΨ − PK∗Ψ ) = Ey∼N (0,Ψ )(y >PKy − y>PK∗y) = E {∑ t≥0 [ 2y>t (K −K∗)(Υ22K∗K∗ −Υ21K∗)yt + y>t (K −K∗)>Υ22K∗(K −K∗)yt ]}
= E {∑ t≥0 y>t (K −K∗)>Υ21K∗(K −K∗)yt } , (F.25)
where we use (F.24) in the last line. Here the expectations are taken following the transition yt+1 = (A−BK)yt with initial state y0 ∼ N (0,Ψ ). Also, we have
E {∑ t≥0 y>t (K −K∗)>Υ22K∗(K −K∗)yt }
= tr [ ΦK(K −K∗)>Υ22K∗(K −K∗) ] ≥ ‖ΦK‖2 · ‖Υ22K∗‖2 · tr [ (K −K∗)>(K −K∗)
] ≥ σmin(Ψ ) · σmin(R) · ‖K −K∗‖2F, (F.26)
where we use the fact that ΦK = (A−BK)ΦK(A−BK)+Ψ Ψ and Υ22K∗ = R+B>PK∗B R in the last line. Combining (F.25) and (F.26), we have J1(K)− J1(K∗) ≥ σmin(Ψ ) · σmin(R) · ‖K −K∗‖2F. We conclude the proof of the lemma.

heading: F.5 PROOF OF LEMMA D.5
text: Proof. Following from Proposition 3.3, we have
J2(KN , bh+1)− J2(KN , b̃h+1) ≤ γb · ∇bJ2(KN , b̃h+1)> [ ∇bJ2(KN , bh)− ∇̂bJ2(KN , bh) ] + (γb)2 · νKN /2 ·
∥∥∇bJ2(KN , bh)− ∇̂bJ2(KN , bh)∥∥22, J2(KN , b̃h+1)− J2(KN , bh+1)
≤ −γb · ∇bJ2(KN , b̃h+1)> [ ∇bJ2(KN , bh)− ∇̂bJ2(KN , bh) ] (F.27)
− (γb)2 · ιKN /2 · ∥∥∇bJ2(KN , bh)− ∇̂bJ2(KN , bh)∥∥22,
where νKN and ιKN are defined in Proposition 3.3. Also, following from Proposition B.3, it holds that∥∥∇bJ2(KN , b̃h+1)∥∥2 ≤ poly1(‖KN‖F, ‖bh‖2, ‖µ‖2, J(KN , b0)) · [1− ρ(A−BKN )]−1.
(F.28)
Combining (F.27), (F.28), and the fact that νKN ≤ ιKN ≤ [1− ρ(A− BKN )]−2 · poly2(‖KN‖2), we know that∣∣J2(KN , bh+1)− J2(KN , b̃h+1)∣∣ (F.29)
≤ (γb)2 · poly2 ( ‖KN‖2 ) · ∥∥∇bJ2(KN , bh)− ∇̂bJ2(KN , bh)∥∥22 · [1− ρ(A−BKN )]−2
+ γb · poly1 ( ‖KN‖F, ‖bh‖2, ‖µ‖2, J(KN , b0) ) · ∥∥∇bJ2(KN , bh)− ∇̂bJ2(KN , bh)∥∥2
· [ 1− ρ(A−BKN ) ]−1 .
Note that from the definition of ∇̂bJ2(KN , bh) and ∇bJ2(KN , bh) in (D.31) and (D.33), respectively, it holds by triangle inequality that∥∥∇bJ2(KN , bh)− ∇̂bJ2(KN , bh)∥∥2
≤ ‖Υ̂22KN −Υ 22 KN ‖2 · ‖KN‖2 · ‖µ̂KN ,bh‖2 + ‖Υ 22 KN ‖2 · ‖KN‖2 · ‖µ̂KN ,bh − µKN ,bh‖2
+ ‖Υ̂22KN −Υ 22 KN ‖2 · ‖bh‖2 + ‖Υ̂ 21 KN −Υ 21 KN ‖2 · ‖µ̂KN ,bh‖2 + ‖q̂KN ,bh − qKN ,bh‖2
+ ‖Υ21KN ‖2 · ‖µ̂KN ,bh − µKN ,bh‖2. By Theorem B.8, combining the fact that J2(KN , bh) ≤ J2(KN , b0) and the fact that ‖µKN ,b‖2 ≤ J(KN , b0)/σmin(Q), we know that with probability at least 1− (T bn)−4 − (T̃ bn)−6, it holds for any ρ ∈ (ρ(A−BKN ), 1) that∥∥∇bJ2(KN , bh)− ∇̂bJ2(KN , bh)∥∥2 (F.30)
≤ λ−1KN · poly3 ( ‖KN‖F, ‖bh‖2, ‖µ‖2, J2(KN , b0) ) · [
log3 T bn (T bn) 1/4(1− ρ)2 +
log1/2 T̃ bn
(T̃ bn) 1/8 · (1− ρ)
] .
Following from the choices of γb, T bn, and T̃ b n in the statement of Theorem B.4, it holds that γb · poly1 ( ‖KN‖F, ‖bh‖2, ‖µ‖2, J(KN , b0) ) · λ−1KN · poly3 ( ‖KN‖F, ‖bh‖2, ‖µ‖2, J2(KN , b0) ) · [
log3 T bn (T bn) 1/4(1− ρ)2 +
log1/2 T̃ bn
(T̃ bn) 1/8 · (1− ρ)
] · [ 1− ρ(A−BKN ) ]−1 + [ 1− ρ(A−BKN ) ]−2 · poly3 ( ‖KN‖F, ‖bh‖2, ‖µ‖2, J2(KN , b0) ) · [
log6 T bn (T bn) 1/2(1− ρ)4 +
log T̃ bn
(T̃ bn) 1/4 · (1− ρ)2 ] · (γb)2 · poly2 ( ‖KN‖2 ) · λ−1KN
≤ νKN · γb · ε/2.
Further combining (F.29) and (F.30), it holds with probability at least 1− ε15 that∣∣J2(KN , bh+1)− J2(KN , b̃h+1)∣∣ ≤ νKN · γb · ε/2. We then finish the proof of the lemma.

heading: F.6 PROOF OF LEMMA D.6
text: Proof. We show that ζK,b ∈ Vζ and ξ(ζ) ∈ Vξ for any ζ ∈ Vζ separately.
Part 1. First we show that ζK,b ∈ Vζ . Note that from Definition B.7, we know that ζ1K,b = J(K, b) satisfies that 0 ≤ ζ1K,b ≤ J(K0, b0). It remains to show that ζ2K,b = αK,b satisfies that ‖ζ2K,b‖2 ≤ Mζ . By the definition of αK,b in (B.6), we know that
‖αK,b‖22 ≤ ‖ΥK‖2F + ‖ΥK‖22 · ( ‖µK,b‖22 + ‖µuK,b‖22 ) + ( ‖A‖2 + ‖B‖2
)2 · (‖PK‖2 · ‖Aµ+ d‖2 + ‖fK,b‖2)2 (F.31) where fK,b = (I − A + BK)−>[(A − BK)>PK(Bb + Aµ + d) − K>Rb] and for notational simplicity, we denote by µuK,b = −KµK,b + b. We only need to bound ΥK , µK,b, µuK,b, PK , and fK,b. Note that by Proposition B.2, the expected total cost J(K, b) takes the form of
J(K, b) = tr(PKΨ ) + µ>K,bQµK,b + (µ u K,b) >RµuK,b + σ 2 · tr(R) + µ>Qµ.
Thus, we have
J(K0, b0) ≥ J(K, b) ≥ σmin(Ψω) · tr(PK) ≥ σmin(Ψω) · ‖PK‖2, J(K0, b0) ≥ J(K, b) ≥ µ>K,bQµK,b ≥ σmin(Q) · ‖µK,b‖2, J(K0, b0) ≥ J(K, b) ≥ (µuK,b)>RµuK,b ≥ σmin(R) · ‖µuK,b‖2,
which imply that
‖PK‖2 ≤ J(K0, b0)/σmin(Ψω), ‖µK,b‖2 ≤ J(K0, b0)/σmin(Q), ‖µuK,b‖2 ≤ J(K0, b0)/σmin(R). (F.32)
For ΥK , it holds that
ΥK = ( Q 0 0 R ) + ( A> B> ) PK (A B) ,
which gives
‖ΥK‖F ≤ (‖Q‖F + ‖R‖F) + ( ‖A‖2F + ‖B‖2F ) · ‖PK‖F,
‖ΥK‖2 ≤ (‖Q‖2 + ‖R‖2) + ( ‖A‖2 + ‖B‖2 )2 · ‖PK‖2. Combining (F.32) and the fact that ‖PK‖F ≤ √ m · ‖PK‖2, we know that
‖ΥK‖F ≤ ( ‖Q‖F + ‖R‖F ) + ( ‖A‖2F + ‖B‖2F ) · √ m · J(K0, b0)/σmin(Ψω),
‖ΥK‖2 ≤ ( ‖Q‖2 + ‖R‖2 ) + ( ‖A‖2 + ‖B‖2 )2 · J(K0, b0)/σmin(Ψω). (F.33)
Now, we upper bound the vector fK,b. Note that by algebra, the vector fK,b takes the form of
fK,b = −PKµK,b + (I −A+BK)−T (QµK,b −K>RµuK,b).
Therefore, we upper bound fK,b as
‖fK,b‖2 ≤ J(K0, b0)2 · σ−1min(Ψω) · σ −1 min(Q) +
[ 1− ρ(A−BK) ]−1 · (κQ + κR · ‖K‖F) (F.34) Combining (F.31), (F.32), (F.33), and (F.34), it holds that
‖ζ2K,b‖2 = ‖αK,b‖2 ≤Mζ,1 +Mζ,2 · (1 + ‖K‖F) · [1− ρ(A−BK)]−1.
Therefore, it holds that ζK,b ∈ Vζ .
Part 2. Now we show that for any ζ ∈ Vζ , we have ξ(ζ) ∈ Vξ. Recall that from (D.41), it holds that ξ1(ζ) = ζ1 − J(K, b), ξ2(ζ) = EπK,b [ ψ(x, u) ] ζ1 + ΘK,bζ 2 − EπK,b [ c(x, u)ψ(x, u) ] . (F.35) Then we have ∣∣ξ1(ζ)∣∣ = ∣∣ζ1 − J(K, b)∣∣ ≤ J(K0, b0), (F.36) where we use the fact that since ζ ∈ Vζ , we have 0 ≤ ζ1 ≤ J(K0, b0) by Definition B.7. Also, by (F.35), we have∥∥ξ2(ζ)∥∥ 2 ≤ ∥∥∥EπK,b[ψ(x, u)]ζ1∥∥∥
2︸ ︷︷ ︸ B1 + ‖ΘK,b‖2 · ‖ζ2‖2︸ ︷︷ ︸ B2
+ ∥∥∥EπK,b[c(x, u)ψ(x, u)]∥∥∥
2︸ ︷︷ ︸ B3 . (F.37)
Note that we upper bound B1 as B1 ≤ J(K0, b0) · ∥∥∥EπK,b[ψ(x, u)]∥∥∥
2 . (F.38)
Following from the definition of ψ(x, u) in (B.5), we know that∥∥∥EπK,b[ψ(x, u)]∥∥∥ 2 ≤ ‖Σz‖F, (F.39)
where Σz is defined as
Σz = Cov [( x u )] = ( ΦK −ΦKK> −KΦK KΦKK> + σ2I ) = ( 0 0 0 σ2I ) + ( I −K ) ΦK ( I −K )> .
Combining (F.38) and (F.39), we have
B1 ≤ J(K0, b0) · ‖Σz‖F. (F.40)
By Proposition B.6, we upper bound B2 as B2 ≤ 4(1 + ‖K‖2F)3 · ‖ΦK‖22 · (Mζ,1 +Mζ,2) · [ 1− ρ(A−BK) ]−1 , (F.41)
where we use the fact that ζ ∈ Vζ and Definition B.7. As for the term B3 in (F.37), we utilize the following lemma to provide an upper bound. Lemma F.2. The vector EπK,b [c(x, u)ψ(x, u)] takes the following form,
EπK,b [ c(x, u)ψ(x, u) ] = 2svec[Σzdiag(Q,R)Σz + 〈Σz,diag(Q,R)〉Σz] Σz ( 2QµK,b 2RµuK,b )  + [ µ>K,bQµK,b + (µ u K,b) >RµuK,b + µ >Qµ
](svec(Σz) 0m 0k ) .
Here the matrix Σz takes the form of
Σz = ( ΦK −ΦKK> −KΦK KΦKK> + σ2 · I ) .
Proof. See §F.11 for a detailed proof.
From Lemma F.2 and (F.32), it holds that B3 ≤ 3 [ ‖Q‖F + ‖R‖F + J(K0, b0) · ‖Q‖2/σmin(Q) (F.42) + J(K0, b0) · ‖R‖2/σmin(R) ] · ‖Σz‖22.
Moreover, by the definition of Σz in (E.25), combining the triangle inequality, we have the following bounds for ‖Σz‖F and ‖Σz‖2,
‖Σz‖F ≤ 2(d+ ‖K‖2F) · ‖ΦK‖2, ‖Σz‖2 ≤ 2(1 + ‖K‖2F) · ‖ΦK‖2. (F.43)
Also, we have J(K0, b0) ≥ J(K, b) ≥ tr [ (Q+K>RK)ΦK ] ≥ ‖ΦK‖2 · σmin(Q),
which gives the upper bound for ΦK as follows,
‖ΦK‖2 ≤ J(K0, b0)/σmin(Q). (F.44)
Therefore, combining (F.37), (F.40), (F.41), (F.42), (F.43), and (F.44), we know that∥∥ξ2(ζ)∥∥ 2 ≤ C · (Mζ,1 +Mζ,2) · J(K0, b0)2/σ2min(Q) (F.45)
· ( 1 + ‖K‖2F )3 · [1− ρ(A−BK)]−1. By (F.36) and (F.45), we know that ξ(ζ) ∈ Vξ for any ζ ∈ Vζ . We conclude the proof of the lemma.

heading: F.7 PROOF OF LEMMA D.7
text: Proof. Assume that z̃0 ∼ N (µ†,Σ†). Following from the fact that
z̃t+1 = Lz̃t + ν + δt,
it holds that
z̃t ∼ N ( Ltµ† +
t−1∑ i=0 Li · ν, (L>)tΣ†Lt + t−1∑ i=0 (L>)iΨδL i
) , (F.46)
where
Ψδ = ( Ψω KΨω KΨω KΨωK > + σ2I ) .
From (D.47), we know that µz takes the form of
µz = (I − L)−1ν = ∞∑ j=0 Ljν. (F.47)
Therefore, combining (F.46) and (F.47), we have
E(µ̂z) = µz + 1
T̃ T̃∑ t=1 Ltµ† − 1 T̃ T̃∑ t=1 ∞∑ i=t Liν. (F.48)
We denote by
µT̃ = T̃∑ t=1 Ltµ† − T̃∑ t=1 ∞∑ i=t Liν.
Meanwhile, it holds that∥∥∥∥ T̃∑ t=1 Ltµ† − T̃∑ t=1 ∞∑ i=t Liν ∥∥∥∥ 2
≤ T̃∑ t=1 ρ(L)t · ‖µ†‖2 + T̃∑ t=1 ∞∑ i=t ρ(L)i · ‖ν‖2
≤ [ 1− ρ(L) ]−1 · ‖µ†‖2 + [1− ρ(L)]−2 · ‖ν‖2 ≤Mµ · (1− ρ)−2 · ‖µz‖2, (F.49)
where Mµ is a positive absolute constant.
For the covariance, note that for any random variables X ∼ N (µ1,Σ1) and Y ∼ N (µ2,Σ2), we know that Z = X + Y ∼ N (µ1 + µ2,Σ), where ‖Σ‖F ≤ 2‖Σ1‖F + 2‖Σ2‖F. Combining (F.46), we know that µ̂z ∼ N (Eµ̂z, Σ̃T̃ /T̃ ), where Σ̃T̃ satisfies that
T̃ /2 · ‖Σ̃T̃ ‖F ≤ T̃∑ t=1 ρ(L)2t · ‖Σ†‖F + T̃∑ t=1 t−1∑ i=0 ρ(L)2i · ‖Ψδ‖F
≤ [ 1− ρ(L)2 ]−1 · ‖Σ†‖F + T̃ · [1− ρ(L)2]−1 · ‖Ψδ‖F, which implies that
‖Σ̃T̃ ‖F ≤MΣ · (1− ρ) −1 · ‖Σz‖F, (F.50)
where MΣ is a positive absolute constant. Combining (F.48), (F.49), and (F.50), we obtain that µ̂z ∼ N ( µz + 1
T̃ µT̃ ,
1 T̃ Σ̃T̃
) ,
where ‖µT̃ ‖2 ≤ Mµ · (1− ρ) −2 · ‖µz‖2 and ‖Σ̃T̃ ‖F ≤ MΣ · (1− ρ) −1 · ‖Σz‖F. Moreover, by the Gaussian tail inequality, it holds with probability at least 1− T̃−6 that
‖µ̂z − µz‖2 ≤ log T̃
T̃ 1/4 · (1− ρ)−2 · poly
( ‖ΦK‖2, ‖K‖F, ‖b‖2, ‖µ‖2 ) .
Then we finish the proof of the lemma.

heading: F.8 PROOF OF LEMMA D.8
text: Proof. We continue using the notations given in §D.3. We define F̂ (ζ, ξ) = { E(ψ̂)ζ1 + E [ (ψ̂ − ψ̂′)ψ̂> ] ζ2 − E(cψ̂) }> ξ2 + [ ζ1 − E(c) ] · ξ1 − 1/2 · ‖ξ‖22,
where ψ̂ = ψ̂(x, u) is the estimated feature vector. Here the expectation is only taken over the trajectory generated by the state transition and the policy πK,b, conditioning on the randomness induced when calculating the estimated feature vectors. Thus, the function F̂ (ζ, ξ) is still random, where the randomness comes from the estimated feature vectors. Note that |F (ζ, ξ) − F̃ (ζ, ξ)| ≤ |F (ζ, ξ) − F̂ (ζ, ξ)| + |F̂ (ζ, ξ) − F̃ (ζ, ξ)|. Thus, we only need to upper bound |F (ζ, ξ) − F̂ (ζ, ξ)| and |F̂ (ζ, ξ)− F̃ (ζ, ξ)|.
Part 1. First we upper bound |F (ζ, ξ)− F̂ (ζ, ξ)|. Note that by algebra, we have∣∣F (ζ, ξ)− F̂ (ζ, ξ)∣∣ =
∣∣∣∣{E(ψ − ψ̂)ζ1 + E[(ψ − ψ′)ψ> − (ψ̂ − ψ̂′)ψ̂>]ζ2 − E[c(ψ − ψ̂)]}>ξ2∣∣∣∣ ≤ E ( ‖ψ − ψ̂‖2 ) · [ |ζ1|+ E ( ‖ψ − ψ′‖2 + 2‖ψ̂‖2 ) · ‖ζ2‖2 + E(c) ] · ‖ξ2‖2, (F.51)
where the expectation is only taken over the trajectory generated by the state transition and the policy πK,b. From Lemma D.7, it holds that
P ( ‖µ̂z − µz + 1/T̃ · µT̃ ‖2 ≤ C1 ) ≥ 1− T̃−6. (F.52)
Therefore, combining (F.52), it holds with probability at least 1− T̃−6 that E ( ‖ψ − ψ′‖2 + 2‖ψ̂‖2 ) ≤ poly ( ‖ΦK‖2, ‖K‖F, ‖b‖2, ‖µ‖2, J(K0, b0) ) , (F.53) where the expectation is conditioned on the randomness induced when calculating the estimated feature vectors. Also, we know that
E(c) ≤ poly ( ‖ΦK‖2, ‖K‖F, ‖b‖2, ‖µ‖2, J(K0, b0) ) . (F.54)
Therefore, combining (F.51), (F.53), (F.54), and Definition B.7, it holds with probability at least 1− T̃−6 that∣∣F (ζ, ξ)− F̂ (ζ, ξ)∣∣ ≤ E(‖ψ − ψ̂‖2) · poly(‖ΦK‖2, ‖K‖F, ‖b‖2, ‖µ‖2, J(K0, b0)). (F.55) Following from the definitions of ψ(x, u) in (B.5) and ψ̂(x, u) in (B.14), we upper bound ‖ψ(x, u)− ψ̂(x, u)‖2 for any x and u as
‖ψ(x, u)− ψ̂(x, u)‖22 = ‖µ̂z − µz‖22 + ∥∥z(µ̂z − µz)> + (µ̂z − µz)z>∥∥2F + ‖µzµ>z − µ̂zµ̂>z ‖2F
≤ poly ( ‖ΦK‖2, ‖K‖F, ‖b‖2, ‖µ‖2, J(K0, b0) ) · ‖µ̂z − µz‖22, (F.56)
where µz is defined in (D.47), µ̂z is defined in (D.48), and z = (x>, u>)>. Also, by Lemma D.7, we know that
‖µ̂z − µz‖2 ≤ log T̃
T̃ 1/4 · (1− ρ)−2 · poly
( ‖ΦK‖2, ‖K‖F, ‖b‖2, ‖µ‖2, J(K0, b0) ) , (F.57)
which holds with probability at least 1 − T̃−6. Combining (F.55), (F.56), and (F.57), it holds with probability at least 1− T̃−6 that∣∣F (ζ, ξ)− F̂ (ζ, ξ)∣∣ ≤ log T̃
T̃ 1/4 · (1− ρ)−2 · poly
( ‖K‖F, ‖b‖2, ‖µ‖2, J(K0, b0) ) . (F.58)
Part 2. We now upper bound |F̂ (ζ, ξ)− F̃ (ζ, ξ)| in the sequel. By definitions, we have∣∣F̃ (ζ, ξ)− F̂ (ζ, ξ)∣∣ =
∣∣∣∣{E(ψ̃ − ψ̂)ζ1 + E[(ψ̃ − ψ̃′)ψ̃> − (ψ̂ − ψ̂′)ψ̂>]ζ2 − E(c̃ψ̃ − ĉψ̂)}>ξ2 + E(ĉ− c̃)ξ1∣∣∣∣ ≤ ∣∣∣∣{E(ψ̂)ζ1 + E(ψ̂ψ̂>)ζ2 − E(ĉψ̂)}>ξ2 + E(ĉ)ξ1∣∣∣∣ · 1Ec (F.59)
+ ∣∣∣[E(ψ̂′ψ̂>)ζ2]>ξ2∣∣∣ · 1(E′∩E)c ,
where we define the event E ′ as E ′ = ( ⋂ t∈[T ] {∣∣‖z′t − µz + 1/T̃ · µT̃ ‖22 − tr(Σ̃z)∣∣ ≤ C1 · log T · ‖Σ̃z‖2})⋂ E2, where E2 is defined in (D.55). Combining the fact that P(E2) ≥ 1 − T̃−6 and Lemma G.3, it holds that P(E ′) ≥ 1− T−5 − T̃−6. Following a similar argument as in Part 1, it holds from (F.59) that∣∣F̃ (ζ, ξ)− F̂ (ζ, ξ)∣∣ ≤ ( 1
T +
1
T̃ 1/4
) · poly ( ‖K‖F, ‖b‖2, ‖µ‖2, J(K0, b0) ) (F.60)
for sufficiently large T and T̃ .
Now, combining (F.58) and (F.60), by triangle inequality, it holds with probability at least 1− T̃−6 that ∣∣F (ζ, ξ)− F̃ (ζ, ξ)∣∣ ≤ ( 1
2T +
log T̃
T̃ 1/4
) · (1− ρ)−2 · poly ( ‖K‖F, ‖b‖2, ‖µ‖2, J(K0, b0) ) .
We finish the proof of the lemma.

heading: F.9 PROOF OF LEMMA E.2
text: Proof. Recall that the feature vector ψ(x, u) takes the following form
ψ(x, u) =
( svec [ (z − µz)(z − µz)> ] z − µz ) .
We then have
ψ(x, u)− ψ(x′, u′) = ( svec [ yy> − (Ly + δ)(Ly + δ)> ] y − (Ly + δ) ) , (F.61)
where we denote by y = z − µz , and (x′, u′) is the state-action pair after (x, u) following the state transition and the policy πK,b. Therefore, for any symmetric matrices M , N and any vectors m, n, it holds from (B.7) and (F.61) that(
svec(M) m
)> ΘK,b ( svec(N)
n ) = Ey,δ {( svec(M)
m
)>( svec(yy>)
y
)( svec [ yy> − (Ly + δ)(Ly + δ)> ] y − (Ly + δ) )>( svec(N) n )} = Ey,δ {( 〈M,yy>〉+m>y ) · [ 〈N, yy> − (Ly + δ)(Ly + δ)>〉+ n>(y − Ly − δ)
]} = Ey ( 〈yy>,M〉 · 〈yy> − Lyy>L> −Ψδ, N〉
)︸ ︷︷ ︸ A1 +Ey ( 〈yy>,M〉 · n>(y − Ly) )︸ ︷︷ ︸ A2
(F.62)
+ Ey ( m>y · 〈yy> − Lyy>L> −Ψδ, N〉 )︸ ︷︷ ︸ A3 +Ey [ m>y · n>(y − Ly) ]︸ ︷︷ ︸ A4 ,
where the expectations are taken over y ∼ N (0,Σz) and δ ∼ N (0,Ψδ). We evaluate the terms A1, A2, A3, and A4 in the sequel.
For the terms A2 and A3 in (F.62), by the fact that y = z−µz ∼ N (0,Σz), we know that these two terms vanish. For A4, it holds that
A4 = Ey [ m>y · (y − Ly)>n ] = Ey [ m>yy>(I − L)>n ] = m>Σz(I − L)>n. (F.63)
For A1, by algebra, we have A1 = Ey ( 〈yy>,M〉 · 〈yy> − Lyy>L> −Ψδ, N〉 ) = Ey ( 〈yy>,M〉 · 〈yy> − Lyy>L>, N〉 ) − Ey ( 〈yy>,M〉 · 〈Ψδ, N〉
) = Ey [ y>My · y>(N − L>NL)y ] − 〈Σz,M〉 · 〈Ψδ, N〉
= Eu∼N (0,I) [ u>Σ1/2z MΣ 1/2 z u · u>Σ1/2z (N − L>NL)Σ1/2z u ] − 〈Σz,M〉 · 〈Ψδ, N〉. (F.64)
Now, by applying Lemma G.1 to the first term on the RHS of (F.64), we know that A1 = 2 tr [ Σ1/2z MΣ 1/2 z · Σ1/2z (N − L>NL)Σ1/2z ] + tr(Σ1/2z MΣ 1/2 z ) · tr [ Σ1/2z (N − L>NL)Σ1/2z ] − 〈Σz,M〉 · 〈Ψδ, N〉
= 2 〈 M,Σz(N − L>NL)Σz〉+ 〈Σz,M 〉 · 〈Σz − LΣzL> −Ψδ, N〉
= 2 〈 M,Σz(N − L>NL)Σz 〉 ,
where we use the fact that Σz = LΣzL> + Ψδ in the last equality. By using the property of the operator svec(·) and the definition of the symmetric Kronecker product, we obtain that
A1 = 2svec(M)>svec [ Σz(N − L>NL)Σz ] = 2svec(M)> [ Σz ⊗s Σz − (ΣzL>)⊗s (ΣzL>) ] svec(N)
= 2svec(M)> [ (Σz ⊗s Σz)(I − L⊗s L)> ] svec(N). (F.65)
Combining (F.62), (F.63), and (F.65), we obtain that( svec(M)
m
)> ΘK,b ( svec(N)
n ) = svec(M)> [ 2(Σz ⊗s Σz)(I − L⊗s L)> ] svec(N) +m>Σz(I − L)>n
=
( svec(M)
m
)>( 2(Σz ⊗s Σz)(I − L⊗s L)> 0 0 Σz(I − L)> )( svec(N) n ) .
Thus, the matrix ΘK,b takes the following form,
ΘK,b =
( 2(Σz ⊗s Σz)(I − L⊗s L)> 0 0 Σz(I − L)> ) ,
which concludes the proof of the lemma.

heading: F.10 PROOF OF LEMMA E.3
text: Proof. From the definition of Θ̃K,b in (B.9), it holds that
‖Θ̃−1K,b‖ 2 2 ≤ 1 + ‖Θ−1K,b‖ 2 2 + ‖Θ−1K,bσ̃z‖ 2 2, (F.66)
where σ̃z is defined as
σ̃z = EπK,b [ ψ(x, u) ] = ( svec(Σz) 0k+m ) .
We bound the RHS of (F.66) in the sequel. For the term Θ−1K,bσ̃z , combining Lemma E.2, we have
Θ−1K,bσ̃z =
( 1/2 · (I − L⊗s L)−>(Σz ⊗s Σz)−1 · svec(Σz)
0k+m ) = ( 1/2 · (I − L⊗s L)−>(Σ−1z ⊗s Σ−1z ) · svec(Σz)
0k+m ) = ( 1/2 · (I − L⊗s L)−> · svec(Σ−1z )
0k+m
) , (F.67)
where we use the property of the symmetric Kronecker product in the second and last line. By taking the spectral norm on both sides of (F.67), it holds that
‖Θ−1K,bσ̃z‖2 = 1/2 · ∥∥(I − L⊗s L)−> · svec(Σ−1z )∥∥2
≤ 1/2 · ∥∥(I − L⊗s L)−>∥∥2 · ∥∥svec(Σ−1z )∥∥2
≤ 1/2 · [ 1− ρ2(L) ]−1 · ‖Σ−1z ‖F ≤ 1/2 · √ k +m · [ 1− ρ2(L)
]−1 · ‖Σ−1z ‖2 = 1/2 · √ k +m · [ 1− ρ2(L)
]−1 · σ−1min(Σz), (F.68) where in the third line we use Lemma G.2 to the matrix L⊗sL. Similarly, we upper bound ‖Θ−1K,b‖2 in the sequel
‖Θ−1K,b‖2 ≤ min { 1/2 · [ 1− ρ2(L) ]−1 σ−2min(Σz), [ 1− ρ(L) ]−1 σ−1min(Σz) } . (F.69)
Thus, combining (F.66), (F.68), and (F.69), we obtain that
‖Θ̃−1K,b‖ 2 2 ≤ 1 + 1/2 ·
√ k +m · [ 1− ρ2(L) ]−1 · σ−1min(Σz) + min { 1/2 · [ 1− ρ2(L) ]−1 σ−2min(Σz), [ 1− ρ(L) ]−1 σ−1min(Σz) } . (F.70)
Now it remains to characterize σmin(Σz). For any vectors s ∈ Rm and r ∈ Rk, we have( s r )> Σz ( s r ) = Ex∼N (µK,b,ΦK),u∼πK,b(· | x) {[ s>(x− µK,b) + r>(u+KµK,b − b) ]2} = Ex∼N (µK,b,ΦK),η∼N (0,I) {[ (s−K>r)>(x− µK,b) + σr>η
]2} = Ex∼N (µK,b,ΦK) {[ (s−K>r)>(x− µK,b) ]2} + Eη∼N (0,I) [ (σr>η)2 ] . (F.71)
The first term on the RHS of (F.71) is lower bounded as Ex∼N (µK,b,ΦK) {[ (s−K>r)>(x− µK,b) ]2} = (s−K>r)>ΦK(s−K>r)
≥ ‖s−K>r‖22 · σmin(ΦK) ≥ ‖s−K>r‖22 · σmin(Ψω), (F.72)
where the last inequality comes from the fact that σmin(ΦK) ≥ σmin(Ψω) by (3.3). The second term on the RHS of (F.71) takes the form of
Eη∼N (0,I) [ (σr>η)2 ] = σ2‖r‖22. (F.73)
Therefore, combining (F.71), (F.72), and (F.73), we have( s r )> Σz ( s r ) ≥ ‖s−K>r‖22 · σmin(Ψω) + σ2‖r‖22
≥ σmin(Ψω) · ‖s‖22 + [ σ2 − ‖K‖22 · σmin(Ψω) ] · ‖r‖22.
From this, we know that σmin(Σz) ≥ min { σmin(Ψω), σ 2 − ‖K‖22 · σmin(Ψω) } . (F.74)
Thus, combining (F.70) and (F.74), we know that ‖Θ̃−1K,b‖2 is upper bounded by a constant λ̃K , where λ̃K only depends on ‖K‖2 and ρ(L) = ρ(A−BK). This finishes the proof of the lemma.

heading: F.11 PROOF OF LEMMA F.2
text: Proof. First, note that the cost function c(x, u) takes the following form,
c(x, u) = ψ(x, u)> svec[diag(Q,R)]2QµK,b 2RµuK,b + [µ>K,bQµK,b + (µuK,b)>RµuK,b + µ>Qµ]. For any matrix V and vectors vx, vu, it holds that
EπK,b [ c(x, u)ψ(x, u) ]>(svec(V ) vx vu )
= EπK,b ψ(x, u)> svec[diag(Q,R)]2QµK,b 2RµuK,b ψ(x, u)>(svec(V )vx vu )︸ ︷︷ ︸ D1
(F.75)
+ EπK,b ψ(x, u)>(µ>K,bQµK,b + (µuK,b)>RµuK,b + µ>Qµ) (svec(V ) vx vu )︸ ︷︷ ︸ D2 .
In the sequel, we calculate D1 and D2 respectively.
Calculation of D1. Note that by the definition of ψ(x, u) in (B.5), it holds that
D1 = EπK,b  [ (z − µz)>diag(Q,R)(z − µz) + (z − µz)> ( 2QµK,b 2RµuK,b )]
· [ (z − µz)>V (z − µz) + (z − µz)> ( vx vu )] = EπK,b [ (z − µz)>diag(Q,R)(z − µz) · (z − µz)>V (z − µz) ] (F.76)
+ EπK,b [( 2QµK,b 2RµuK,b )> (z − µz)(z − µz)> ( vx vu )] .
Here z = (x>, u>)> and µz = EπK,b(z). For the first term on the RHS of (F.76), note that z − µz ∼ N (0,Σz). Therefore, by Lemma G.1, we obtain that
EπK,b [ (z − µz)>diag(Q,R)(z − µz) · (z − µz)>V (z − µz) ] = 2 〈 Σzdiag(Q,R)Σz, V 〉 + 〈 Σz,diag(Q,R) 〉 · 〈Σz, V 〉
= svec [ 2Σzdiag(Q,R)Σz + 〈 Σz,diag(Q,R) 〉 · Σz ]> svec(V ). (F.77)
Meanwhile, the second term on the RHS of (F.76) takes the form of
EπK,b [( 2QµK,b 2RµuK,b )> (z − µz)(z − µz)> ( vx vu )] = [ Σz ( 2QµK,b 2RµuK,b )]>( vx vu ) . (F.78)
Combining (F.76), (F.77), and (F.78), we obtain that
D1 = 2svec[Σzdiag(Q,R)Σz + 〈Σz,diag(Q,R)〉Σz] Σz ( 2QµK,b 2RµuK,b ) >(svec(V )vx vu ) . (F.79)
Calculation of D2. By the definition of the feature vector ψ(x, u) in (B.5), we know that
D2 = (µ > K,bQµK,b + (µ u K,b) >RµuK,b + µ >Qµ) (svec(Σz) 0m 0k )>(svec(V ) vx vu ) . (F.80)
Now, combining (F.75), (F.79), and (F.80), it holds that
EπK,b [ c(x, u)ψ(x, u) ] = 2svec[Σzdiag(Q,R)Σz + 〈Σz,diag(Q,R)〉Σz] Σz ( 2QµK,b 2RµuK,b )  + [ µ>K,bQµK,b + (µ u K,b) >RµuK,b + µ >Qµ
](svec(Σz) 0m 0k ) ,
which concludes the proof of the lemma.

heading: G AUXILIARY RESULTS
text: Lemma G.1. Assume that the random variable w ∼ N (0, I), and let U and V be two symmetric matrices, then it holds that
E(w>Uw · w>V w) = 2 tr(UV ) + tr(U) · tr(V ).
Proof. See Magnus et al. (1978) and Magnus (1979) for a detailed proof.
Lemma G.2. Let M , N be commuting symmetric matrices, and let α1, . . . , αn, β1, . . . , βn denote their eigenvalues with v1, . . . , vn a common basis of orthogonal eigenvectors. Then the n(n+ 1)/2 eigenvalues of M ⊗s N are given by (αiβj + αjβi)/2, where 1 ≤ i ≤ j ≤ n.
Proof. See Lemma 2 in Alizadeh et al. (1998) for a detailed proof.
Lemma G.3. For any integer m > 0, let A ∈ Rm×m and η ∼ N (0, Im). Then, there exists some absolute constant C > 0 such that for any t ≥ 0, we have
P [∣∣η>Aη − E(η>Aη)∣∣ > t] ≤ 2 · exp[−C ·min(t2‖A‖−2F , t‖A‖−12 )].
Proof. See Rudelson et al. (2013) for a detailed proof.

Review summaries:
This paper considers the problem of model-free reinforcement learning in discrete-time, linear-quadratic Markovian mean-field games with ergodic costs. The authors begin by establishing the existence and uniqueness of the Nash equilibrium in such a setting, and then proposes an mean-field actor-critic algorithm with linear function approximation. The proposed algorithm is shown to converge linearly with high probability under certain standard conditions.