Introduction

The modeling and optimization of networked systems such as wireless communication networks and traffic networks is a long-standing challenge. Typically analytic models must make numerous assumptions to obtain tractable models as a result of the complexity of the systems, which include many unknown, or unmodeled dynamics. Given the success of Reinforcement Learning (RL) in a wide array of domains such as game play [Silver et al., 2016;Mnih et al., 2015], robotics [Duan et al., 2016], and autonomous driving [Li et al., 2019], it has emerged as a promising tool for tackling the complexity of networked systems. However, when seeking to use RL in the context of the control and optimization of large-scale networked systems, scalability quickly becomes an issue. The goal of this paper is to develop scalable multi-agent RL for networked systems.

Motivated by real-world networked systems like wireless communication, epidemics, and traffic, we consider an RL model of n agents with local interaction structure. Specifically, each agent i has local state s i , local action a i and the agents are associated with an underlying dependence graph G and interact locally, i.e, the distribution of s i (t + 1) only depends on the current states of the local neighborhood of i as well as the local a i (t). Further, each agent is associated with stage reward r i that is a function of s i , a i , and the global stage reward is the average of r i . In this setting, the design goal is to find a decision policy that maximizes the (discounted) global reward. This setting captures a wide range of applications, e.g. epidemics [Mei et al., 2017], social networks [Chakrabarti To illustrate our model and our results, we provide stylized examples of applications in three areas: multi-access wireless communication, epidemics, and traffic signal control in Section 2.2. We conduct numerical experiments to demonstrate the performance of the approach using both synthetic examples and an application to wireless communication in Section 5.

Related Literature. Our problem falls into the category of the "succinctly described" MDPs in Blondel and Tsitsiklis [2000, Section 5.2], where the state/action space is a product space formed by the individual state/action space of multiple agents. As the state/action space is exponentially large, such problems are not scalable in general, even when the problem has structure [Blondel and Tsitsiklis, 2000;Whittle, 1988;Papadimitriou and Tsitsiklis, 1999]. Despite this, there is a large literature on RL/MDPs in multi-agent settings, which we discuss below.

Multi-agent RL dates back to the early work of Littman [1994]; Claus and Boutilier [1998]; Littman [2001]; Hu and Wellman [2003] (see Bu et al. [2008] for a review) and has been actively studied, e.g. Zhang et al. [2018]; Kar et al. [2013]; Macua et al. [2015]; Mathkar and Borkar [2017]; Wai et al. [2018], see a more recent review in Zhang et al. [2021]. Multi-agent RL encompasses a broad range of settings including competitive agents and Markov games. The case most relevant to ours is the cooperative multi-agent RL where typically, the agents can take their own actions but they share a common global state and maximize a global reward [Bu et al., 2008]. This is in contrast to the model we study, in which each agent has its own state and acts upon its own state. Despite the existence of a global state, multi-agent RL still faces scalability issues since the joint-action space is exponentially large. Methods have been proposed to deal with this, including independent learners [Tan, 1993;Claus and Boutilier, 1998;Matignon et al., 2012], where each agent employs a single-agent RL method. While successful in some cases, the independent learner approach can suffer from instability [Matignon et al., 2012]. Alternatively, one can use function approximation schemes to approximate the large Q-table, e.g. linear function approximation [Zhang et al., 2018] or neural networks [Lowe et al., 2017]. Such methods can reduce computation complexity significantly, but it is unclear whether the performance loss caused by the function approximation is small. In contrast, our technique not only reduces computation but also guarantees small performance loss.

Factored MDPs are problems where every agent has its own state and the state transition factorizes in a way similar to our model [Kearns and Koller, 1999;Guestrin et al., 2003;Osband and Van Roy, 2014]. However, they differ from the model we consider in that each agent does not have its own action. Instead, there is a global action affecting every agent. Despite the difference, Factored MDPs still suffer from scalability issues. Similar approaches as in the case of multi-agent RL are used, e.g., Guestrin et al. [2003] proposes a class of "factored" linear function approximators; however, it is unclear whether the loss caused by the approximation is small.

Other Related Work. Our work is also related to weakly coupled MDPs, where every agent has its own state and action but their transition is decoupled [Meuleau et al., 1998]. Additionally, our model shares some similarity with Glauber dynamics in physics [Lokhov et al., 2015;Mezard and Montanari, 2009], though our focus is very different from these works. As we consider the class of localized policies, another related line of work is Partially Observable MDP (POMDP) [Nair et al., 2005;Oliehoek and Amato, 2016;Bertsekas, 2005], though the formulations and results we have are very different from those works.

Finally, this work is related to our earlier work Qu and Li [2019], which assumes the full knowledge of the MDP model (not RL) and imposes strong assumptions on the graph. In contrast, our work here does not need knowledge of the MDP and significantly relaxes the assumptions.

Preliminaries

In this section, we introduce our model, provide a few illustrative examples, and provide important background in RL that underlies our analysis. Throughout this paper, • denotes Euclidean norm and • ∞ denotes infinity norm. Notation t and T are reserved as iteration counters for the inner loop of the algorithm to be introduced later, m and M for the outer loop, and κ is used for counting the hops of neighbors. Notation O(•) hides constants and Õ(•) hides log factors with respect to iteration variables T, M and variable κ. The total variation distance for two distributions π, π over a finite set S is defined as TV(π, π ) = sup E⊂S |π(E) -π (E)|.

Model

We consider a network of n agents that are associated with an underlying undirected graph G = (N , E), where N = {1, . . . , n} is the set of agents and E ⊂ N × N is the set of edges. Each agent i is associated with state s i ∈ S i , a i ∈ A i where S i and A i are finite sets. The global state is denoted as s = (s 1 , . . . , s n ) ∈ S := S 1 × • • • × S n and similarly the global action a = (a 1 , . . . , a n ) ∈ A := A 1 × • • • × A n . At time t, given current state s(t) and action a(t), the next individual state s i (t + 1) is independently generated and is only dependent on neighbors:

where notation N i means the neighborhood of i (including i itself) and notation s N i means the states of the agents in N i . In addition, for integer κ ≥ 0, we use N κ i to denote the κ-hop neighborhood of i, i.e. the nodes whose graph distance to i has length less than or equal to κ. We also let f (κ) = sup i |N κ i |. Each agent is associated with a class of localized policies ζ θ i i parameterized by θ i . The localized policy ζ θ i i (a i |s i ) is a distribution on the local action a i conditioned on the local state s i , and each agent, conditioned on observing s i (t), takes an action a i (t) independently drawn from ζ θ i i (•|s i (t)). We use θ = (θ 1 , . . . , θ n ) to denote the tuple of the localized policies ζ θ i i , and also use ζ θ (a|s) = n i=1 ζ θ i i (a i |s i ) to denote the joint policy, which is a product distribution of the localized policies as each agent acts independently.

Further, each agent is associated with a stage reward function r i (s i , a i ) that depends on the local state and action, and the global stage reward is r(s, a) = 1 n n i=1 r i (s i , a i ). The objective is to find localized policy tuple θ such that the discounted global stage reward is maximized, starting from some initial state distribution π 0 ,

(2)

Remark 1. In the state transition (1) of our model, the distribution of each node's next state is allowed to depend on its neighbors' states s N i (t), but only on its own action a i (t) as opposed to its neighbors' actions a N i (t). This restriction is imposed only for simplicity of exposition. With a simple change of notation, our model, algorithm and analysis can be extended to the more general dependence on a N i (t). Similarly, each agent's reward function r i (s i , a i ) can be generalized to depend on its neighbors' state-action pairs, i.e. r i (s N i , a N i ), and each agent's localized policy can also be generalized to depend on its neighbors' states, i.e.

Remark 2. In the paper, the interaction graph G is undirected, but our model and results can be easily generalized to the directed graph setting without essential changes in the algorithm and the analysis. In detail, to generalize to the directed graph case, the only change needed is to redefine the "neighbors". Specifically, N i needs to be redefined as the "in-neighborhood" of i (including i itself ), i.e. i itself and the set of nodes that have a directed link pointing towards i. In addition, N κ i needs to be redefined as the κ-hop "in-neighborhood" of i, i.e. the nodes whose shortest directed link to i has a length less than or equal to κ (including i itself ).

Examples

In this section, we provide three networked system examples in wireless communication, epidemics, and traffic that feature the local dependence structure we study in this paper. For ease of exposition, we present simple versions of these examples, keeping the essence of the model and highlighting the dependence structure while ignoring some application-specific details.

Wireless Communication. We consider a wireless network with multiple access points [Zocca, 2019], where there is a set of users N = {1, 2, • • • , n}, and a set of network access points Y = {y 1 , y 2 , • • • , y m }. Each user i only has access to a subset Y i ⊆ Y of the access points. We define the interaction graph as the conflict graph, in which two users i and j are neighbors if and only if they share an access point, i.e. the neighbors of user i is

Each user i maintains a queue of packets defined as follows. At time step t, with probability p i , user i receives a new packet with an initial deadline d i . Then, user i can choose to send the earliest packet in its queue to one access point in its available set Y i , or not send anything at all. If an action of sending to y k ∈ Y i is taken, and if no other users send to the same access point at this time, then the earliest packet in user i's queue is transmitted with success probability q k which depends on the access point y k ; however, if another user also chooses to send to y k , then there is a conflict and no transmission occurs. If the packet is successfully transmitted, it will be removed from user i's queue and user i will get a reward of 1. After this, the system moves to the next time step, with all deadlines of the remaining packets decreasing by 1 and packets with deadline 0 being discarded. In this example, the local state s i of user i is a characterization of its queue of packets, and is represented by a d i binary tuple

where for each ∈ {1, . . . , d i }, e ∈ {0, 1} indicates whether user i has a packet with remaining deadline . The action space is A i = {null} ∪ Y i , where null represents the action of not sending. The detailed transition is provided in Table 1, where s i (t + 1) only depends on s N i (t), a N i (t), which fits into the local interaction structure we consider. The local reward is given by r i (s N i (t), a N i (t)) = 1 in the case of the last row of Table 1, and r i (s N i (t), a N i (t)) = 0 in all other cases. The local state space S i , the local action space A i , the local transition probabilities in Table 1 and the local reward function r i (•) form a networked MDP model described in Section 2.1. The above model serves as a basis for more complex multi-access wireless communication models studied in the literature, including those with multiple channels [Block and Van Houdt, 2016], (imperfect) carrier sensing [Kim et al., 2011]. Existing analytical approaches typically require knowledge of modeling details and parameters such as the packet arrival rate [Tassiulas and Ephremides, 1990;Yun et al., 2012], while our RL-based approach does not require such knowledge and learns to improve performance in a model-free manner.

Epidemic Network. We consider an SIS (Susceptible-Infected-Susceptible) epidemic network model [Mei et al., 2017;Ahn, 2014;Azizan Ruhi et al., 2016a], where there is a undirected graph of nodes G = (N , E), and each node has a binary state space S i = {susceptible, infected}, as well all other cases (denote a i (t) = y k ) Flip left most "1" in s i (t) to "0" w.p. q k , then left shift s i (t) and append Bernoulli(p i ).

Table 1: State transition for the wireless communication example. "Left shift" means for a binary tuple, discarding the left most bit. Bernoulli(p i ) means a random variable sampled i.i.d. from the Bernoulli distribution that has probability p i to be 1, and probability 1 -p i to be 0. as a finite action space A i with action a i ∈ A i representing epidemic control measures like different levels of vaccination [Preciado et al., 2013]. The evolution of the states follows a local interaction structure: the probability of a node turning from susceptible to infected depends on the whether its neighboring nodes are infected or not as well as its control action in place [Azizan Ruhi et al., 2016b;Preciado et al., 2013]; the probability of a node turning from infected to susceptible depends on the recovering rate. More precisely, s i (t + 1) only depends on s N i (t) and a i (t), and the state transition is provided by,

where δ i ∈ (0, 1) is a given recovering rate parameter, and β i : A i → (0, 1) is a given transmission rate function and it depends on the control action a i (t), and |{j ∈ N i /{i} : s j (t) = 1}| is the number of neighboring nodes excluding i itself that are infected. The local reward of each node is given by,

which consists of two parts: a positive reward of 1 if the node is free from infection, subtracting a given cost function c i (a i ) on the epidemic control measure a i . In this setup, the expected global reward is a weighted balance between the overall infection level and the epidemic control cost. The above defined local state space S i , local action space A i , local transition probabilities in (3) and local rewards in (4) form a networked MDP model in Section 2.1. We comment that the above SIS model is a basis for more complex models, e.g. those with "exposed" and "recovered" states [Kuznetsov and Piccardi, 1994;Britton, 2010] or more complex control interventions [Morris et al., 2020]. These more complex models can also be captured using the framework above. We reiterate that the approach in our paper does not require knowledge of model specifications and learns in a model-free manner.

Traffic Network. Lastly, we consider a traffic signal control problem adapted from Varaiya [2013]. In this setting, each node i represents a road link and the interaction graph represents the physical connection of the road links. Given road link i, the local state s i = (x i,j ) i→j is a tuple of variables with j ranging from neighboring links i can turn to, and x i,j is the number of vehicles on link i that intend to turn to link j, and can only take values in [S] = {0, 1, . . . , S}. Correspondingly, the local state space is S i = [S] N i /{i} . Similarly, the local action a i = (y i,j ) i→j is the binary traffic signal tuple with y i,j controlling the on-off of turn movement (i → j), and the local action space is A i = {0, 1} {N i }/{i} . At each time, a random amount of vehicles on the queue x i,j will flow into link j when the traffic signal y i,j is on. Meanwhile, link i will receive vehicles from other incoming links, a random fraction of which are then assigned to each of the queues in (x i,j ) i→j . Mathematically,

where [x] S 0 means max(min(x, S), 0), C i,j (t) (and similarly C k,i (t)) is an i.i.d. random variable indicating the random amount of vehicles leaving i for j, and R i,j (t) is an i.i.d. random variable that controls the split of the inflow to link i to the queue x i,j (t). See Varaiya [2013] for the complete details. Given a fixed distribution on the random variables C i,j (t), R i,j (t), (5) provides a complete characterization of the distribution of s i (t + 1) conditioned on s N i (t) and a N i (t), in which the local state at each link s i (t + 1) only depends on its neighbors' current states and current actions, which fits into our local interaction structure (equation (1) and Remark 1). The local reward is a characterization of the congestion level at link i, and one version of the reward is the negative queue length

The local state space S i , the local action space A i , the local transition probabilities (5) and the local rewards (6) form a networked MDP discussed in Section 2.1. We comment that policies for traffic signal control, like the max pressure policy in Varaiya [2013], typically require knowledge of the statistics of the random variables C ij (t) and R ij (t), while our approach learns from data in a model-free manner.

Background in RL

To provide background for the analysis in this paper, we review a few key concepts in RL. First, fixing a localized policy tuple θ = (θ 1 , . . . , θ n ), an important notion is the Q-function, which is defined for policy θ as a "table" of values for each state-action pair (s, a) ∈ S × A and it is the expected infinite horizon discounted reward under policy θ conditioned on the initial state and action being (s, a):

In the last step, we have defined Q θ i (s, a) which is the Q function for the individual reward r i . Both Q θ and Q θ i are exponentially large tables and, therefore, are intractable to compute and store. Additionally, another important concept we use is the policy gradient theorem, which provides a characterization of the gradient of the objective J(θ) and is the basis of many algorithmic results in RL. The policy gradient theorem shows that the gradient of J(θ) depends on Q θ and, therefore, is intractable to compute using the form in Lemma 1.

Lemma 1 (Sutton et al. [2000]). Let π θ be a distribution on the state space given by π θ (s) = (1 -γ) ∞ t=0 γ t π θ t (s), where π θ t is the distribution of s(t) under a fixed policy θ when s(0) is drawn from π 0 . Then, we have,

3 Algorithm Design and Results

In this paper we propose an algorithm, Scalable Actor Critic (SAC), which provably finds an O(ρ κ+1 )-stationary point of the objective J(θ) (i.e. a θ s.t. ∇J(θ) 2 ≤ ε) for some ρ ≤ γ, with complexity scaling in the size of the local state-action space of the largest κ-hop neighborhood. We state our main result formally in Theorem 4 after introducing the details of SAC and the key idea underlying its design.

Key Idea: Exponential Decay of Q-function Leads to Efficient Approximation

Recall that the policy gradient in Lemma 1 is intractable to compute due to the dimension of the Q-function. Our key idea is that exponential decay of the Q function allows efficient approximation of the Q-function via truncation. To illustrate this, we start with the definition of the exponential decay property. Recall that N κ i is the set of κ-hop neighborhood of node i and define N κ -i = N /N κ i , i.e. the set of agents that are outside of i'th κ-hop neighborhood. We write state s as (s N κ i , s N κ -i ), i.e. the states of agents that are in the κ-hop neighborhood of i and outside of the κ-hop neighborhood respectively. Similarly, we write a as (a N κ i , a N κ -i ). The exponential decay property is then defined as follows.

Definition 1. The (c, ρ)-exponential decay property holds if, for any localized policy θ, for any

It may not be immediately clear when the exponential decay property holds. Lemma 2 (a) below highlights that the exponential decay property holds generally with ρ = γ, without any assumption on the transition probabilities except for the factorization structure (1) and the localized policy structure. Further, many MDPs in practice have ergodicity and fast mixing properties, and Lemma 2 (b) shows that when such fast mixing property holds, the (c, ρ)-exponential decay property holds for some ρ < γ depending on the mixing rate. The proof of Lemma 2 is postponed to Appendix A. The condition on mixing rate in Lemma 2 (b) is similar to those used in the literature on the finite time analysis of RL methods, e.g. Zou et al. [2019]. In fact, our condition is weaker than the common mixing rate condition in that we only require the distribution of the local state-action pair (s i (t), a i (t)) to mix, instead of the full state-action pair (s(t), a(t)). We leave it as future work to study such "local" mixing behavior and its relation to the local transition probabilities (1).

Lemma 2. Assume ∀i, r i is upper bounded by r. Then the following holds.

(a) The ( r 1-γ , γ)-exponential decay property holds.

(b) If there exists c > 0 and µ ∈ (0, 1) s.t. under any policy θ, the Markov chain is ergodic and starting from any initial state, TV(π t,i , π ∞,i ) ≤ c µ t , ∀t, where π t,i is the distribution of (s i (t), a i (t)) and π ∞,i is the distribution for (s i , a i ) in stationarity, and recall TV(•, •) is the total variation distance. Then, the ( 2c r 1-γµ , γµ)-exponential decay property holds. Our definition of exponential decay is similar in spirit to the "correlation decay", or "spatial decay" that has been studied in the literature [Gamarnik, 2013;Gamarnik et al., 2014;Bamieh et al., 2002], though these works consider very different settings. For example, Gamarnik [2013] and Gamarnik et al. [2014] study optimization in a graphical model setting (no concept of state and/or time), and show that the effect of cost functions far away on the optimal solution at a particular node shrinks exponentially in their graph distance, under certain weak interaction assumptions. Compared to these works where the optimization problem is static, we focus on an MDP setting which has states that evolve on a time axis. Further, our exponential decay is in terms of the Q functions, as opposed to the optimal solution. That being said, we believe there are deep connections between our results and that in Gamarnik [2013]; Gamarnik et al. [2014], and we leave the investigation of it as future work.

The power of the exponential decay property is that such properties usually lead to scalable and distributed algorithm design, as in Gamarnik [2013]. In our context, the exponential decay property guarantees that the dependence of Q θ i on other agents shrinks quickly as the distance between them grows. This motivates us to consider the following class of truncated Q-functions,

where

) are any non-negative weights satisfying

With the definition of the truncated Q-function, our key insight is the following Lemma 3, which says when the exponential decay property holds, the truncated Q-function (9) approximates the full Q-function with high accuracy and can be used to approximate the policy gradient. The proof of Lemma 3 is postponed to Appendix B.

Lemma 3. Under the (c, ρ)-exponential decay property, the following holds:

(a) Any truncated Q-function in the form of (9) satisfies,

(b) Given i, define the following truncated policy gradient,

where Qθ j can be any truncated Q-function in the form of (9). Then, if

Algorithm 1: SAC: Scalable Actor Critic Input: θ i (0); parameter κ; T , length of each episode; step size parameters h, t 0 , η.

i to be the all zero vector.

Each agent i calculates approximated gradient,

12 Each agent i conducts gradient step

The power of this lemma is that the truncated Q function has a much smaller dimension than the true Q function, and is thus scalable to compute and store. However, despite the reduction in dimension, the error resulting from the approximation is small. In the next section, we use this idea to design a scalable algorithm.

Algorithm Design: Scalable Actor Critic (SAC)

The good properties of the truncated Q-function open many possibilities for algorithm design. For instance, one can first obtain the truncated Q-function in some way (which could be much easier than directly computing the full Q-function) and then do a policy gradient step using the Lemma 3. In this subsection, we propose one particular approach using the actor critic framework. Our approach, Scalable Actor Critic (SAC), uses temporal difference (TD) learning to obtain the truncated Q-function and then uses policy gradient for policy improvement. The pseudocode of the proposed algorithm is given in Algorithm 1.

Overall structure. The overall structure of SAC is a for-loop from line 1 to line 13. Inside the outer loop, there is an inner loop (line 4 through line 9) that uses temporal difference learning to get the truncated Q-function, which is followed by a policy gradient step that does policy improvement.

The Critic: TD-inner loop. Line 4 through line 9 is the policy evaluation inner loop that obtains the truncated Q function, where line 7 and 8 are the temporal difference update. We note that steps 7 and 8 use the same update equation as TD learning, except that it "pretends" (s N κ i , a N κ i ) is the true state-action pair while the true state-action pair should be (s, a). As will be shown in the theoretic analysis, such a TD update implicitly gives an estimate of a truncated Q function.

The Actor: policy gradient. Line 10 through line 12 define the actor actions. Here, each agent calculates an estimate of the truncated gradient based on (11), and then conducts a gradient step.

Communication. To implement our training algorithm, each agent needs to communicate with other agents in its κ-hop neighborhood in line 7 and line 11; after training is done, each agent implements its localized policy that does not need communication. This communication requirement is weaker than the "centralized training with decentralized execution" paradigm in the multi-agent RL literature [Lowe et al., 2017], where in the training phase, global communication is used. We also comment that when κ = 0, our algorithm does not need communication and is effectively the same as the independent learner approach in the literature [Tan, 1993;Lowe et al., 2017], as each agent simply runs a single-agent actor critic method based on its local state and local action. When κ > 1, our algorithm requires communication with agents beyond the direct 1-hop neighbors, which may be unrealistic for some applications. An interesting future direction is to reduce the communication requirements, e.g. potentially using consensus schemes like in Zhang et al. [2018], and also techniques that only communicate quantized bits as opposed to real numbers [Magnússon et al., 2020].

Discussion. Our algorithm serves as an initial concrete demonstration of how to make use of the truncated Q-functions to develop a scalable RL method for networked systems. There are many extensions and other approaches that could be pursued, either within the actor critic framework or beyond. One immediate extension is to do a warm start, i.e., initialize Q0

i as the final estimate QT i in the previous outer-loop. Additionally, one can use the TD-λ variant of TD learning, incorporate variance reduction schemes like the advantage function (Advantage Actor Critic), or incorporate function approximation. Further, beyond the actor critic framework, another direction is to develop Q-learning/SARSA type algorithms based on the truncated Q-functions. An appealing aspect of Q-learning/SARSA algorithms is that they may exhibit better convergence properties, but the challenge is that, unlike the actor critic framework, it is not straightforward to enforce the policy to be local in Q-learning/SARSA algorithms. These are interesting topics for future work.

Remark 3 (Model-based vs model-free). We note that by using an actor critic framework, the proposed approach is model-free, meaning it does not explicitly estimate the transition probabilities and the reward function. This is in contrast to model-based RL, which explicitly estimates the transition probabilities (or parameters that determine the transition probabilities, like the δ i , β i (•) parameter in the epidemic example in Section 2.2). On one hand, it is known that model-based RL can be more sample efficient than model-free RL in certain circumstances [Tu and Recht, 2019].

Additionally, for specific applications, model-based control design may come with properties like robustness [Varaiya, 2013]. On the other hand, model-free RL offers more flexibility since it does not impose assumptions on the model class. The comparison and tradeoff between model-based and model-free approaches is an open research question and is beyond the scope of this paper. We refer the reader to Tu and Recht [2019]; Qu et al. [2020b] for more details.

Approximation Bound

In this section, we state and discuss the formal approximation guarantee for SAC. Before stating the theorem, we first state the assumptions we use. The first assumption is standard in the RL literature and bounds the reward and state/action space size.

Assumption 1 (Bounded reward and state/action space size). The reward is upper bounded as 0 ≤ r i (s i , a i ) ≤ r, ∀i, s i , a i . The individual state and action space size are upper bounded as

Assumption 2 (Exponential decay). The (c, ρ)-exponential decay property holds for some ρ ≤ γ.

Note that under Assumption 1, Assumption 2 automatically holds with ρ = γ, cf. Lemma 2 (a). However, we state the exponential decay property as an assumption to account for the more general case that ρ could be strictly less than γ, cf. Lemma 2 (b).

Our third assumption can be interpreted as an ergodicity condition which ensures that the state-action pairs are sufficiently visited.

Assumption 3 (Sufficient local exploration). There exists positive integer τ and σ ∈ (0, 1) s.t. under any fixed policy θ and any initial state-action (s, a)

Assumption 3 requires that every state action pair in the κ-hop neighborhood must be visited with some positive probability after some time. This type of assumption is common for finite time convergence results in RL. For example, in Srikant and Ying [2019]; Li et al. [2020], it is assumed that every state-action pair is visited with positive probability in the stationary distribution and the state-action distribution converges to the stationary distribution with some rate. This implies our assumption which is weaker in the sense that we only require local state-action pair (s N κ i , a N κ i ) to be visited as opposed to the full state-action pair (s, a). Having said that, we note that by making Assumption 3, we do not consider the exploration-exploitation tradeoff, which is a challenging issue even in single-agent RL. One potential way to relax Assumption 3 is to use Upper Confidence Bound (UCB) bonuses to encourage exploration, which has been proposed in single-agent RL [Jin et al., 2018]. We leave the study of the exploration-exploitation tradeoff in the multi-agent networked setting as future work.

Finally, we assume boundedness and Lipschitz continuity of the gradients, which is standard in the RL literature.

Assumption 4 (Bounded and Lipschitz continuous gradient). For any i, a i , s i and θ i , we assume

With these assumptions in hand, we are ready to state our convergence result.

Theorem 4. Under Assumption 1, 2, 3 and 4, for any δ ∈ (0, 1), M ≥ 3, suppose the critic step size

where

The proof of Theorem 4 is deferred to Section 4. To interpret the result, note that the first term in (13) converges to 0 in the order of Õ(foot_0 √ M ) and the second term, which we denote as ε κ , is the bias caused by the truncation of the Q-function and it scales in the order of O(ρ κ+1 ). As such, our method SAC will eventually find an O(ρ κ+1 )-approximation of a stationary point of the objective function J(θ), which could be very close to a true stationary point even for small κ as ε κ decays exponentially in κ.

In terms of complexity, (13) gives that, to reach a O(ε κ )-approximate stationary point, the number of outer-loop iterations required is

), which scales polynomially with the parameters of the problem. We emphasize that it does not scale exponentially with n. Further, since the left hand side of ( 12) decays to 0 as T increases in the order of Õ( 1 √ T ) and the right hand side of ( 12) is in the same order as O(ε κ ), the inner-loop length required is

). This iteration complexity for the inner loop can potentially be further reduced if we do a warm start for the inner-loop, as the Q-estimate from the previous outer-loop should be already a good estimate for the current outer-loop. We leave the finite time analysis of the warm start variant as future work.

In the complexity bound, a key parameter is σ, which we recall is defined in Assumption 3 and it roughly means the probability that a state-action pair in a κ-hop neighborhood is visited. Suppose we interpret σ to scale with σ ∼ interpreted as follows. In multi-agent RL, the typical way to handle the curse of dimensionality is to use function approximation of the Q function. For example, Zhang et al. [2018] uses linear function to approximate the Q-function. Alternatively, it has also been popular to use neural networks to do the approximation [Lowe et al., 2017]. However, in these approximation methods, the resulting steady state error in the actor critic framework depends on the approximation error, and it is generally unclear how to choose the function approximator that is both computationally tractable and also accurate in representing the true Q function. In the context of these works, our truncated Q-functions can be viewed as a specific way of function approximation that exploits the local interaction structure in our problem, and our method is not only computationally tractable but also has a small approximation error.

Proof of Main Result

In this section, we provide the proof of the main result Theorem 4 with some auxiliary derivations postponed to the appendix. As our algorithm is an actor critic algorithm, the proof is divided into two parts: firstly, we provide an analysis of the critic, i.e. TD learning that estimates the truncated Q-function; secondly, we analyze the actor and finish the proof of Theorem 4.

Analysis of the critic. The first part of the analysis concerns the critic, and we show that the critic inner loop converges to an estimate with steady-state error exponentially small in κ. Specifically, recall that within iteration m the inner loop update is

i is initialized to be all zero, and α t = h t+t 0 is the step size. We note that when implementing (14) within outer loop iteration m, trajectory (s(t), a(t)) is generated by the agents taking a fixed policy θ(m). Let Q θ(m) i ∈ R S×A be the true Q-function for reward r i under this fixed policy θ(m) as defined in (7). Given the above notation, we prove the following theorem on the critic, which bounds the error between the approximation QT i generated by ( 14) and the true Q θ(m) i . Theorem 5. Assume Assumption 1, 2, 3 are true and suppose t 0 , h satisfies, h ≥ 1 σ max(2, 1 1-√ γ ) and t 0 ≥ max(2h, 4σh, τ ). Then, inside outer loop iteration m, for each i ∈ N , with probability at least 1 -δ, we have the following error bound,

where

The proof of Theorem 5 is postponed to Section 4.1. The result in Theorem 5 is an upper bound of the infinity-norm error between the truncated Q-function QT i obtained by TD learning and the true Q-function. This error bound can be further decomposed into two parts -a transient part that converges to zero in the order of Õ( 1 √ T + 1 T ), and a steady state error that is exponentially small in κ. In proving Theorem 5, we use two key techniques. The first is using the exponential decay property (cf. Definition 1, Lemma 2) to show that in the "steady state", the error of the truncated Q-function is bounded by O(ρ κ+1 ). This is possible due to our results in Lemma 3 which shows the class of truncated Q-functions are good approximations of the full Q-function. Our second proof technique is that we develop novel finite time analysis tools for TD learning to obtain the finite time error bound. Our proof uses a novel recursive decomposition of the error. Compared to existing work on finite time analysis on TD learning in Srikant and Ying [2019], our result does not need the ergodicity assumption (the weaker Assumption 3 instead), and obtains error bounds in the infinity norm as opposed to the (weighted) Euclidean norm in Srikant and Ying [2019]. This finite time proof technique could be of independent interest.

Analysis of the actor. We view the actor step as a biased stochastic gradient step, with the bias characterized by our result in the critic (Theorem 5). Under this viewpoint, we finish the analysis of the actor and the proof of the main result in Section 4.2.

Analysis of the Critic: Proof of Theorem 5

Since Theorem 5 is entirely about a particular outer-loop iteration m, inside which the policy is fixed to be θ(m), to simplify notation we drop the dependence on m and θ(m) throughout this section. Particularly, we refer to Q θ(m) i as Q * i , which is the true Q-function for reward r i under policy θ(m) (cf. ( 7)). We also introduce short-hand notation z = (s, a) ∈ Z = S × A to represent a particular state action pair (s, a) ∈ S × A. Similarly, we define z i = (s i , a i ) ∈ Z i = S i × A i , and

A vector v ∈ R Z means a vector of dimension |Z| that is indexed by z ∈ Z, with its z'th entry denoted by v(z). For example, the full Q-function Q * i will be treated as a vector in R Z with its z'th entry denoted by

, and its z N κ i 'th entry will be denoted by v(z N κ i ). For example, the truncated Q-functions Qt i will be treated as vectors in R Z N κ i . Following a similar convention, a matrix A ∈ R Z×Z will be a |Z|-by-|Z| matrix indexed by (z, z ) ∈ Z × Z, with its (z, z )'th entry denoted by A(z, z ).

Theorem 5 essentially says that the critic iterate Qt i in (14) will become a good estimate of Q * i as t increases. We note that the full Q-function Q * i must satisfy the Bellman equation [Bertsekas and Tsitsiklis, 1996],

where TD : R Z → R Z is the standard Bellman operator for reward r i and P ∈ R Z×Z is the transition probability matrix from z(t) to z(t + 1) under policy θ(m). Note in (15), without causing any confusion, r i is interpreted as a vector in R Z although r i only depends on z i .

Our proof is divided into 3 steps. In Step 1, we rewrite (14a) and (14b) in a linear update form (cf. ( 17)), study its behavior (Lemma 6), and decompose the error into a recursive form (Lemma 7). In Step 2, we bound the noise sequences in the error decomposition (Lemma 10 and Lemma 11). Finally, in Step 3, we use the recursive error decomposition and the bound on the noise sequences to prove Theorem 5.

Step 1: error decomposition. Define e z N κ i to be the indicator vector in R Z N κ i , i.e. the z N κ i 'th entry of e z N κ i is 1 and other entries are zero. Then, the critic update equation ( 14) can be written as,

with Q0 i being the all zero vector in R

Qt-1 i , we can make the following definition (where we omit the dependence on i in notation A z,z , b z ),

and rewrite ( 16) in a linear form

We define the following simplifying notations,

. Let F t be the σ-algebra generated by z(0), . . . , z(t). Then, clearly A t-1 is F t -measurable and b t-1 is F t-1 measurable. As a result, Qt i is F t -measurable. Let τ > 0 to be the integer in Assumption 3. Let

τ measurable random vectors (matrices). With these notations, (17) can be rewritten as,

where in the last step, we have defined sequences t-1 , φ t-1 ∈ R Z N κ i , which are noise sequences that we will later control in Step 2. For now, we focus on the term Āt-1

Qt-1 i + bt-1 , and show the following Lemma 6. The proof of Lemma 6 is similar to the analysis of fixed points for TD learning with linear function approximation, see e.g. Van Roy and Tsitsiklis [1995]. We postpone the detailed proof of Lemma 6 to Appendix C.1. Lemma 6. ∀t, there exists diagonal matrix

i and operator g t-1 : R

Āt-1

Qt-

where D t-1 satisfies D t-1 σI (recall σ > 0 is from Assumption 3) with its z N κ i 'th diagonal entry being dt-1 (z

From Lemma 6, we can see the first two terms in (18) can be written as (I -

, which roughly speaking drives the iterate to Q * ,t-1 i , the fixed point of operator g t-1 . Though depending on t, Q * ,t-1 i is a good approximation of the true Q-function Q * i regardless of t, as shown in Lemma 6. As such, moving towards Q * ,t-1 i at each time step should eventually produce a good estimate of the full Q-function Q * i . In the following, we turn the above intuition into an error decomposition. In details, we unroll (18) and get,

where we have used the following short-hand notations

Since every diagonal entry of D is lower bounded by σ almost surely (Lemma 6), we have every entry of B k,t is upper bounded by β k,t and every entry of Bk,t is upper bounded by βk,t almost surely. With these short-hand notations, we state the error decomposition below, which is a consequence of ( 20) and the property of operator g k (•) and its fixed point Q * ,k i in Lemma 6. The proof of Lemma 7 is postponed to Appendix C.2.

The following recursion holds almost surely,

where b k,t (z

From Lemma 7, it is clear that to bound the error ξ t , we need to bound t-1 k=τ α k Bk,t k ∞ and t-1 k=τ α k Bk,t φ k ∞ , which is the focus of the next step.

Step 2: bound the k and the φ k -sequence. The goal of this step is to bound

We first show some simple properties of Qt i , t and φ t in Lemma 8. Since every entry of α k Bk,t is upper bounded by β k,t , we also show some properties of β k,t , βk,t in Lemma 9. The proofs of the two lemmas are postponed to Appendix C.3.

Lemma 8. We have almost surely, (a)

Lemma 9. If α t = h t+t 0 , where t 0 ≥ h > 2 σ and t 0 ≥ 4σh, and t 0 ≥ τ , then β k,t , βk,t satisfies (a)

We now bound t-1 k=τ α k Bk,t k ∞ . Clearly, t-1 is F t -measurable, and satisfies

where the last equality is due to the definition of Āt-1 and bt-1 and the fact Qt-τ i is F t-τ -measurable. Equation ( 21) shows that t-1 is a "shifted" martingale difference sequence (it is not a standard martingale difference sequence which would require E t-1 |F t-1 = 0). Therefore, t-1 k=τ α k Bk,t k ∞ can be controlled by Azuma-Hoeffding type inequalities, as shown by Lemma 10. We comment that Bk,t is also random and Bk,t k is no longer a martingale difference sequence. As a result, to prove Lemma 10 requires more than the direct application of the Azuma-Hoeffding bound. For more details, see the full proof of Lemma 10 in Appendix C.4.

Lemma 10. We have with probability 1 -δ,

Finally we bound sequence

, which is small due to the step size selection. The proof of Lemma 11 can also be found in Appendix C.4.

Lemma 11. We have almost surely,

Step 3: bound the critic error.

We are now ready to use the error decomposition in Lemma 7 as well as the bounds on the k , φ k -sequences in Lemma 10 and Lemma 11 to bound the error of the critic. Recall that Theorem 5 states with probability 1 -δ,

where C 0 = 2cρ κ+1 1-γ , and

To prove (22), we start by applying Lemma 10 to t ≤ T with δ replaced by δ/T . Then, using a union bound, we get with probability 1 -δ, for any t ≤ T ,

Combining this with Lemma 7 and using Lemma 11, we get with probability 1 -δ, for all τ ≤ t ≤ T ,

We now condition on ( 23) is true and use induction to show (22). Eq. ( 22) is true for t = τ , as

Then, assume ( 22) is true for up to k ≤ t -1, we have by ( 23),

We use the following auxiliary Lemma, whose proof is provided in Section C.5.

Lemma 12. Recall α k = h k+t 0 , and b k,t (z

and α 0 ≤ 1 2 , then, for any z N κ i , and any 0 < ω ≤ 1,

With Lemma 12, and using the bound on βτ-1,t in Lemma 9 (a), we have

To finish the induction, it suffices to show F t ≤ Ca √ t+t 0 and F t ≤ C a t+t 0 . To see this, note

1-γ , one can check our selection of C a and C a satisfies the above three inequalities, and so the induction is finished and the proof of Theorem 5 is concluded.

Analysis of the Actor and Proof of Main Result (Theorem 4)

With the error of the critic bounded in Theorem 5, the second part of the proof focuses on the actor, i.e. the policy gradient step. Recall that at iteration m, the policy gradient step is given by θ

where Qm,T i is the final estimate of the Q-function for r i at the end of the critic loop in iteration m, where we have added an additional superscript m to Qm,T i to indicate its dependence on m; {s(t), a(t)} T t=0 is the state-action trajectory with s(0) drawn from π 0 (the initial state distribution defined in the objective function J(θ), cf. ( 2)) and the agents taking policy θ(m). Our goal is to show that ĝi (m) is approximately the right gradient direction, ∇ θ i J(θ(m)), which by Lemma 1 can be written as,

where π θ(m) t

is the distribution of s(t) under fixed policy θ(m) when the initial state is drawn from π 0 ; Q θ(m) is the true Q function for the global reward r under policy θ(m), cf. (7).

To bound the difference between ĝi (m) and the true gradient ∇ θ i J(θ(m)), we define the following additional sequences,

where

is the true Q function for r i under policy θ(m). We also use notation h(m), g(m), ĝ(m) to denote the respective h i (m), g i (m), ĝi (m) stacked into a larger vector. The following result is an immediate consequence of Assumption 1 and Assumption 4, whose proof is postponed to Appendix D.1.

Lemma 13. We have almost surely, ∀m ≤ M , max( ĝ(m) , g(m) , h(m) , ∇J(θ(m))

With these definitions, we decompose the error between the gradient estimator ĝ(m) and the true gradient ∇J(θ(m)) into the following three terms,

In the following, we will provide bounds on e 1 (m), e 2 (m), e 3 (m), and then combine these bounds to prove our main result Theorem 4.

Bounds on e 1 (m). Notice that the difference between ĝi (m) and g i (m) is that the critic estimate Qm,T j is replaced with the true Q-function Q θ(m) j

. By Theorem 5, we have Qm,T with high probability when T is large enough, based on which we can bound e 1 (m) , which is formally provided in Lemma 14. The proof of Lemma 14 is postponed to Appendix D.2. Lemma 14. When T is large enough s.t.

with ¯ = 4 r 1-γ + 2r, then we have with probability at least

Bounds on e 2 (m). Let G m be the σ-algebra generated by the trajectories in the first m outer-loop iterations. Then, θ(m) is G m-1 measurable, and so is h i (m). Further, by the way that the trajectory

), e 2 (m) is a martingale difference sequence w.r.t. G m , and we have the following bound in Lemma 15 which is a direct consequence of Azuma-Hoeffding bound. The proof of Lemma 15 is postponed to Section D.3. Lemma 15. With probability at least 1 -δ/2, we have

Bounds on e 3 (m). The term e 3 (m) is the error between h(m) and the true gradient ∇J(θ(m)), where h(m) is similar to the truncated policy gradient in Lemma 3 and therefore should be close to the true gradient by Lemma 3. We provide Lemma 16 below to bound e 3 (m) . Due to technical issues, h(m) is not exactly the same as the truncated policy gradient defined in Lemma 3, but a variant of it instead. Therefore the conclusion of Lemma 3 can not be directly used in the proof of Lemma 16. Nevertheless, the proof of Lemma 16 follows essentially the same arguments as in Lemma 3 and is postponed to Appendix D.4.

, we have almost surely, e 3 (m) ≤ 2 Lc (1-γ) ρ κ+1 .

Combining the bounds and proof of Theorem 4. With the above bounds on e 1 (m), e 2 (m) and e 3 (m), we are now ready to prove the main result Theorem 4. Since ∇J(θ) is L Lipschitz continuous, we have

Using the decomposition of ĝ(m) in ( 28), we get, Plugging the above bounds on ĝ(m) 2 and ∇J(θ(m)), ĝ(m) into (29), we have,

where ε m,0 = ∇J(θ(m)), e 2 (m) , ε m,1 = ∇J(θ(m)) ( e 1 (m) + e 3 (m) ), ε m,2 = 2L ( e 1 (m) 2 + e 2 (m) 2 + e 3 (m) 2 ). Doing a telescope sum for (30), we get

where we have used

We now apply our bounds on e 1 (m), e 2 (m), e 3 (m). By Lemma 15, we have with probability

By Lemma 14 and Lemma 16, we have with probability 1

By Lemma 13, we have almost surely, max( e 1 (m) , e 2 (m) , e 3 (m) ) ≤ 2 rL (1-γ) 2 , and hence,

Using a union bound, we have with probability 1 -δ, all three events (33), ( 34) and ( 35) hold, which when combined with (32) implies

). Further we use the bound J(θ(M )) ≤ r 1-γ and J(θ(0)) ≥ 0 almost surely. Combining these results, we get with probability 1 -δ,

This concludes the proof of Theorem 4.

Numerical Studies

In this section, we first conduct numerical studies in Section 5.1 using a synthetic example where the optimal solution is known in closed form to verify our theoretic results. Then, in Section 5.2, we demonstrate our approach using the wireless communication example introduced in Section 2.2.

Synthetic Experiments

We first study a synthetic example where the interaction graph is a line of n nodes N = {1, 2, . . . , n} with the left most node labeled as 1 and the right most as n. Each node has a binary local state space and local action space S i = A i = {0, 1}. The left most node (node 1) has a reward 1 whenever s 1 = 1, and all other reward values of node 1 and the rewards of all other nodes are 0. Further, s 1 (t + 1) = 1 with probability 1 when the second node has state s 2 (t) = 1, and in all other cases s 1 (t + 1) = 0 with probability 1. For the i'th node with 2 ≤ i ≤ n -1,

all other cases.

For the last node i = n, s n (t + 1) = 1 w.p. 1 when a n (t) = 1, and s n (t + 1) = 0 w.p. 1 when a n (t) = 0. The initial states of all nodes are s i (0) = 1.

In this example, the rewards transitions are designed in a way so that the optimal policy can be determined explicitly. Only agent 1 has a non-zero reward when its state is s 1 = 1, and it will stay in state s 1 = 1 only when s 2 = 1, which happens with high probability when a 2 = 1 and s 3 = 1, and so on so forth. Therefore, it is clear that the optimal policy is for all nodes to always take action a i (t) = 1 regardless of the local state, in which case the states of all agents will stay as s i (t) = 1, and the resulting optimal global discounted reward is

In our experiments, we set n = 8, γ = 0.7. We use the softmax policy for the localized policies, which is a standard policy parameterization that encompasses all stochastic policies from S i to A i [Sutton et al., 1998]. We run the SAC algorithm with κ = 0 up to κ = 7. We plot the global discounted reward throughout the training process in Figure 1, and we also plot the optimality gap in Figure 2 computed as the difference of the optimal discounted reward ( 1 n 1 1-γ ) and the discounted reward achieved by the algorithm with the respective κ values. Figure 1 shows that when κ increases, the performance of the algorithm also increases, though the improvement becomes small when κ > 1. This is also confirmed by Figure 2, which shows the optimality gap is decaying exponentially when κ increases. This is consistent with our theoretic result in Theorem 4. We note that the exponential decay appears to stop in Figure 2 when κ > 4, and this is due to the increasing complexity in training for large κ as discussed in Section 3.3, and may also be due to the fact that any further potential improvement may be too small (in the order of 1 × 10 -3 ) to be noticed because of the noise.

Multi-Access Wireless Communication

We next study the multi-access wireless communication example discussed in Section 2.2. We consider a grid of users in Figure 3, where each user has access points on the corners of the block it is in. In the experiments, we set the grid size as 6 × 6, deadline as d i = 2, and all parameters p i (packet arrival probability for user i) and q k (success transmission probability for access point y k ) are generated uniformly random from [0, 1]. We set γ = 0.7 and the initial state to be uniformly random, and run the SAC algorithm with κ = 0, 1, 2 to learn a localized softmax policy, starting from an initial policy where the action is chosen uniformly random. We compare the proposed method with a benchmark based on the localized ALOHA protocol [Roberts, 1975], where each user has a certain probability of sending the earliest packet and otherwise not sending at all. When it sends, it sends the packet to a random access point in its available set, with probability proportional to the success transmission probability of this access point and inverse proportional to the number of users that share this access point. The results are shown in Figure 4. It shows that the proposed algorithm can outperform the ALOHA based benchmark, despite the proposed algorithm does not have access to the transmission probability q k which the benchmark has access to. It also shows that the SAC with κ = 1, 2 outperforms κ = 0, which as we mentioned in Section 3.2 corresponds to the independent learner approach in the literature [Tan, 1993;Lowe et al., 2017]. SAC with κ = 2 outperforms κ = 1, but the improvement is small, which is consistent with the results in the synthetic experiments in Section 5.1. We also study a case with the same 6 × 6 grid of access points, but the 36 users are assigned randomly to the square blocks in the grid. We plot the results in Figure 5. Similar phenomenons can be observed in Figure 5 as in Figure 4, with SAC outperforms the benchmark and SAC with κ = 1, 2 outperforms κ = 0, the independent learner approach.

Concluding Remarks and Extensions

This paper proposes a SAC algorithm that provably finds a close-to-stationary point of J(θ) in time that scales with the local state-action space size of the largest κ-hop neighborhood, which can be much smaller than the full state-action space size when the graph is sparse. This represents the first scalable RL method for localized control of multi-agent networked systems with such a provable guarantee. There are many possible extensions and open questions, which we discuss below.

Average Reward. One extension is to consider average reward instead of discounted reward, i.e., to consider

Under appropriate ergodicity assumptions, the average reward above is equivalent to the reward under the stationary distribution. Average reward is common in applications where the performance is measured in stationarity, e.g. throughput or waiting time in communication and queueing networks. Despite the importance in applications, average reward RL is known to be more challenging even in single-agent settings, see e.g. Tsitsiklis andVan Roy [1999, 2002]. For example, the Q function needs to be defined in a different way,

and because of the lack of a discounting factor γ, the associated Bellman operator no longer has γ as the natural contraction factor. In ongoing work summarized in Qu et al. [2020a], we have begun to study an average reward multi-agent RL setting with the same local interaction structure (1) as in this paper. While similar exponential decay properties on the Q-functions can be defined, in the average reward setting, Qu et al. [2020a] shows that the exponential decay only holds under a form of bounded interaction strength assumption. Under this assumption, Qu et al. [2020a] proposes a variant of the SAC algorithm that achieves similar guarantees as the algorithm in this paper. However, the bounded interaction strength assumption in Qu et al. [2020a] may be restrictive, and searching for weaker assumptions that guarantee the exponential decay (or weaker forms of decay) is an interesting open future direction. Time-varying dependence structure. In this paper, the interaction graph that defines the dependence structure (1) is a fixed graph, while in many real world applications the graph is time-varying. In other words, if the interaction graph at time t is G t , then (1) becomes,

where N i (G t ) means the set of neighbors of i in graph G t . In another of our extensions of this work [Lin et al., 2021], we study a stochastic graph setting where, at each time step, the graph G t is sampled from a fixed graph distribution in which each link is present with probability exponentially decreasing in a predefined distance measure between the two nodes. Lin et al. [2021] shows that a weaker form of exponential decay holds in this setting, termed µ-decay, and a similar SAC algorithm can be used. Beyond Lin et al. [2021], an open question is how to handle the case where G t can be arbitrarily chosen, rather than stochastically sampled. An issue in this setting is that the transition kernel may not be time-homogeneous anymore, which means that is a challenging open problem.

Other future directions. While SAC is an actor critic based algorithm, the idea underlying SAC, including the exponential decay property and the truncated Q-function (9), is a contribution in its own right, and can potentially lead to other scalable RL methods for networked systems. For example, within the actor critic framework, one idea is to change the critic to variants like TD-λ, change the actor to include the advantage function, or use the simultaneous update version of the actor critic algorithm (as opposed to the inner-loop structure used in this paper). Beyond the actor critic framework, one can develop policy-iteration type algorithms (e.g. Q-learning/SARSA type methods) on our truncated Q-functions. Further, our proof technique in showing the finite time convergence of TD learning can be of independent interest. In fact, we have already considered preliminary applications of our technique in the context of Q-learning and TD learning with state aggregation [Lin et al., 2021], and we expect other applications to emerge in the coming years. Additionally, the setting we consider here does not provide an investigation of the tradeoff between exploration and exploitation. Adding consideration of this tradeoff is an interesting direction for future work. Finally, other future directions include studying the landscape of our policy optimization problem, and studying the robustness of the trained policies.

A The Exponential Decay Property

Our main results depend on the (c, ρ)-exponential decay of the Q-function (cf. Definition 1), which means that for any i, any s

In Section 3.1, we have pointed out in Lemma 2 that the (c, ρ)-exponential decay property always holds with ρ ≤ γ, assuming the rewards r i are upper bounded. We now provide the proof of Lemma 2.

Proof of Lemma 2. We first prove part (a). For notational simplicity, denote s = (s

) and a = (a

). Let π t,i be the distribution of (s i (t), a i (t)) conditioned on (s(0), a(0)) = (s, a) under policy θ, and let π t,i be the distribution of (s i (t), a i (t)) conditioned on (s(0), a(0)) = (s , a ) under policy θ. Then, we must have π t,i = π t,i for all t ≤ κ. The reason is that, due to the local dependence structure (1) and the localized policy structure, π t,i only depends on (s N t i , a N t i ) (the initial state-action of agent i'th t-hop neighborhood) which is the same as (s N t i , s N t i ) when t ≤ κ per the way the initial state (s, a), (s , a ) are chosen. With these definitions, we expand the definition of

where TV(π t,i , π t,i ) is the total variation distance between π t,i and π t,i which is upper bounded by 1. The above inequality shows that the ( r 1-γ , γ)-exponential decay property holds and concludes the proof of Lemma 2 (a).

The proof of part (b) is almost identical to that of part(a). The only change is that in step (37), we use TV(π t,i , π t,i ) ≤ 2c µ t .

B Proof of Lemma 3

We first show part (a), the truncated Q function is a good approximation of the true Q function.

To see that, we have for any (s, a) ∈ S × A, by ( 9) and (10),

where in the last step, we have used the (c, ρ)-exponential decay property, cf. Definition 1. Next, we show part (b). Recall by the policy gradient theorem (Lemma 1),

where we have used

) by the localized policy structure. With the above equation, we can compute ĥi (θ) -∇ θ i J(θ),

We claim that E 2 = 0. To see this, consider for any j ∈ N κ -i ,

= s,a 1 ,...,a i-1 ,a i+1 ,...,an

where in the last equality, we have used Qθ j (s N κ j , a N κ j ) does not depend on a i as i ∈ N κ j ; and

where in the last step, we have used (38) and the upper bound ∇ θ i log ζ θ i i (a i |s i ) ≤ L i . This concludes the proof of Lemma 3.

C Proof of Auxiliary Results for the Analysis of the Critic C.1 Proof of Lemma 6

The goal of Lemma 6 is to understand Āt-1

Qt-

Recall that P is the transition matrix from z(t-1) to z(t). Given any distribution d on the stateaction space Z, we define Ãz = E z ∼P (•|z) A z,z and Ād = E z∼d Ãz , bd = E z∼d b z . Then, Āt-1 and bt-1 can be rewritten as Āt-1 = Ād t-1 , bt-1 = bd t-1 . In what follows, we provide characterizations of Ād , bd for general distributions d.

Firstly, Ãz is given by,

where P (•|z) is understood as the z'th row of P and is treated as a row vector. Also, we have defined Φ ∈ R Z×Z N κ i to be a matrix with each row indexed by z ∈ Z and each column indexed by

Further, the z'th row of Φ is the indicator vector e z N κ i , in other words Φ(z, z

Then, Ād , bd are given by,

where diag(d) ∈ R Z×Z is a diagonal matrix with the z'th diagonal entry being d(z); in the last equation, r i is understood as a vector over the entire state-action space Z, though it only depends on z i . With the above characterizations, we show the following property of Ād and bd in Lemma 17. Lemma 17 (with d set as d t-1 ) will directly lead to the results in Lemma 6, with D t-1 being the D in Lemma 17 and it satisfies D t-1 σI due to Assumption 3; g t-1 and Q * ,t-1 i being g d t-1 and Qd t-1 i in Lemma 17.

Lemma 17. Given distribution d on state-action pair z whose marginalization onto z N κ i is non-zero for every z N κ i , we have Ād Qi + bd can be written as

where

i is a diagonal matrix, with the z N κ i 'th entry being the marginalized distribution of z N κ i under distribution d; g d (•) is given by g d ( Qi ) = Π d TDΦ Qi , where

i is a diagonal matrix, and the z N κ i 'th diagonal entry is the marginal probability of z N κ i under d, which is non-zero by the assumption of the lemma. Therefore, Φ diag(d)Φ ∈ R

i is invertable and matrix Π d = (Φ diag(d)Φ) -1 Φ diag(d) is well defined. Further, the z N κ i 'th row of Π d is in fact the conditional distribution of the full state z given z N κ i . So, Π d must be a stochastic matrix and is non-expansive in infinity norm.

By the definition of Ād and bd , we have,

where TD is the Bellman operator for reward r i defined in (15), and operator g d is given by

Notice that Φ is non-expansive in • ∞ norm since each row of Φ has precisely one entry being 1 and all others are zero. Also since Π d is non-expansive in • ∞ norm and TD is a γ-contraction in • ∞ norm, we have g d = Π d TDΦ is a γ contraction in • ∞ norm. As a result, g d has a unique fixed point Qd i . Finally, we show (42), which bounds the distance between Φ Qd i and Q * i , where Q * i is the true Q-function for reward r i and it is the unique fixed point of TD operator (15). We have,

where the equality follows from the fact that Qd i is the fixed point of Π d TDΦ, Q * i is the fixed point of TD; the last inequality is due to ΦΠ d TD is a γ contration in infinity norm. Therefore,

Next, recall that the z N κ i 's row of Π d is the distribution of the state-action pair z conditioned on its N κ i coordinates being fixed to be z N κ i . We denote this conditional distribution of the states outside of

And therefore,

Further, we have

where the last inequality is due to the exponential decay property (cf. Definition 1 and Assumption 2). Therefore,

Combining the above with (43), we get the desired result

C.2 Proof of Lemma 7

Recall that the z N κ i 'th diagonal entry of B k,t is b k,t (z N κ i ), and we define similarly the z N κ i 'th diagonal entry of Bk,t to be bk,t (z N κ i ). Using these notations, equation ( 20) can be written as,

Notice that by definition, bτ-1,t (z

where in the thrid inequality, we have used that g k is γ-contraction in infinity norm with fixed point Q * ,k i , and in the last inequality, we have used the property of Q * ,k i in Lemma 6. Combining the above with (44), we have

C.3 Proof of Lemma 8 and Lemma 9

In this section, we provide proofs of the two auxiliary lemmas, Lemma 8 and Lemma 9. We start with the proof of Lemma 8.

Proof of Lemma 8. First, notice that A

Part (a) can be proved by induction. Part (a) is true for t = 0 as Q0

or in other words,

And for other entries of Qt i , it stays the same as Qt-1

which finishes the induction and the proof of part (a).

For part (b), notice that t = (A t -Āt ) Qt+1-τ i + b t -bt . Therefore, it is easy to check that by part (a), t ∞ ≤ 4 r 1-γ + 2r = ¯ .

For part (c), notice that, for any k

Therefore, by triangle inequality,

As a consequence,

Proof of Lemma 9. Notice that log(1 -x) ≤ -x for all x < 1. Then,

(

which leads to the bound on β k,t and βk,t .

For part (b),

where we have used (k + 1 + t 0 ) 2σh ≤ 2(k + t 0 ) 2σh , which is true when t 0 ≥ 4σh. Then,

,

where in the last inequality we have used 2σh -1 > σh.

For part (c), notice that for k

where we have used (k + 1 + t 0 ) σh ≤ 2(k + t 0 ) σh , and σh -1 > 1 2 σh.

C.4 Proof of Lemma 10 and Lemma 11

Since

and Lemma 9 (c) directly leads to the bound in Lemma 11. So, in this section, we focus on the proof of Lemma 10. We start by stating a variant of the Azuma-Hoeffding bound that handles our "shifted" Martingale difference sequence.

Lemma 18. Let X t be a F t -adapted stochastic process, satisfying EX t |F t-τ = 0. Further, |X t | ≤ Xt almost surely. Then with probability 1 -δ, we have,

Proof of Lemma 18. Let be an integer between 0 and τ -1. For each , define process Y k = X τ k+ , scalar Ȳ k = Xkτ+ , and define Filtration F k = F τ k+ . Then, Y k is F k -adapted, and satisfies

Therefore, applying Azuma-Hoeffding bound on Y k , we have

i.e. with probability at least 1 -δ τ , k:kτ + ≤t

Using the union bound for = 0, . . . , τ -1, we get that with probability at least 1 -δ,

where the last inequality is due to Cauchy-Schwarz.

Recall that Lemma 10 is an upper bound on

with d (z

Proof of Lemma 19. Let p k be a scalar sequence defined as follows. Set p τ = 0, and

, and to prove Lemma 19 we need to bound

We claim that for all k ≥ k 0 + 1, p k and pk have the same sign, and |p k | ≤ |p k |. This is obviously true for k = k 0 + 1. Suppose it is true for for k -1. Without loss of generality, suppose both p k-1 and pk-1 are non-negative. Since k -1 > k 0 and by the definition of k 0 , we must have

|. These imply pk ≥ p k > 0. The case where both p k-1 and pk-1 are negative are similar. This finishes the induction, and as a result,

By the definition of k 0 , we have

, where in the last step, we have used the upper bound on k 0 ∞ in Lemma 8 (b). As a result,

With the above preparations, we are now ready to prove Lemma 10.

Proof of Lemma 10. Fix z N κ i and τ ≤ k 0 ≤ t -1. As have been shown in (21), k (z

As a result, we can use the Azuma-Hoeffding bound in Lemma 18 to get with probability 1 -δ,

By a union bound on τ ≤ k 0 ≤ t -1, we get with probability 1 -δ,

C.5 Proof of Lemma 12

Throughout the proof, we fix z N κ i and prove the desired upper bound. For notational simplicity, we drop the dependence on z N κ i and write b k,t and dk instead, and we will use the property dk ≥ σ. Define the sequence

We use induction to show that e t ≤ 1 √ γ(t+t 0 ) ω . The statement is clearly true for t = τ + 1, as e τ +1 = b τ,τ +1 1 (τ +t 0 ) ω = α τ dτ 1 (τ +t 0 ) ω ≤ 1 √ γ(τ +1+t 0 ) ω (last step needs α τ ≤ 1 2 , (1 + 1 t 0 ) ω ≤ 2 √ γ , implied by t 0 ≥ 1, ω ≤ 1). Let the statement be true for t -1. Then, notice that,

where the inequality is based on induction assumption. Then, plug in α t-1 = h t-1+t 0 and use dt-1 ≥ σ, we have,

Now using the inequality that for any x > -1, (1 + x) ≤ e x , we have, 1 -σh t -1 + t 0

(1 -√ γ) 1 + 1 t -1 + t 0

where we have used that Qm,T The upper bounds for g(m) , h(m) and ∇J(θ(m)) can be obtained in an almost identical way and their proof is therefore omitted.

D.2 Proof of Lemma 14

Let G m be the σ-algebra generated by the trajectories in the first m outer-loop iterations. Then, Theorem 5 implies that, fixing each m ≤ M and i ∈ N , conditioned on G m-1 , following event happens with probability at least 1 -δ:

where

with ¯ = 4 r 1-γ + 2r. We can take expectation and average out G m-1 , and apply union bound over 0 ≤ m ≤ M -1 and i ∈ N , getting with probability at least 1 -

As a result, sup 0≤m≤M -1 ĝ(m) -g(m) ≤ 4cLρ κ+1

(1 -γ) 3 , which is true conditioned on event (47) is true that happens with probability at least 1 -δ 2 .

D.3 Proof of Lemma 15

By The rest of the proof is essentially the same as Lemma 3. For completeness we provide a proof below. Combining the above two equations, we have,

Clearly, the second term satisfies E 2 ≤ L i r (1-γ) 2 γ T +1 . For E 1 , we have as defined in (9). We claim E 4 is zero. To see this, consider for any j ∈ N κ -i and any t,

where in the last equality, we have used Qθ(m) j (s N κ j , a N κ j ) does not depend on a i as i ∈ N κ j ; and

(a i |s i ) = ∇ θ i 1 = 0. For E 3 , by the exponential decay property, the truncated Q function has a small error, cf. ( 38 and as a result,

Therefore,