Multi-agent network systems have found applications in various societal infrastructures, such as power systems, tra!c networks, and smart cities [McArthur et al., 2007, Burmeister et al., 1997, Roscia et al., 2013]. One particularly important class of such problems is the cooperative multi-agent network MDP setting, where agents are embedded in a graph, and each agent has its own local state [Qu et al., 2020b].
In network MDPs, the local state transition probabilities and rewards only depend on the states and actions of the agent's Reinforcement Learning (RL) [Sutton, 2018] has emerged as a promising tool for addressing the complex dynamics of these systems [Chen et al., 2024, Nezamoddini and Gholami, 2022, Yan and Xu, 2020]. There are several pioneering works on designing scalable RL algorithms for network systems [Qu et al., 2021, Lin et al., 2021a, Zhang et al., 2023c]. To facilitate scalable control in network control, in [Qu et al., 2021], the authors introduced a key insight, referred to as the exponential decay property of the Q-function. This property suggests that each agent's local Q-function can be well-approximated using only information from its ω-hop neighborhood. We note that a similar property has also been proposed in [Gu et al., 2022] which focuses on reinforcement learning in the mean field multi-agent setting. Leveraging this property, the proposed algorithm concentrates on learning truncated Q-functions and then applying either policy gradient [Qu et al., 2021, Lin et al., 2021a] or policy iteration [Zhang et al., 2023c]. However, although these methods are scalable with respect to the network size, they are limited to the tabular setting, where each agent must store a local Q-table that scales with the state and action spaces of its neighborhood, making it inefficient for large state and action spaces. In fact, due to the inherent complexity in network MDPs, i.e., network size is large and the state and action spaces of each agent are large or continuous, designing e!cient and scalable RL algorithms for such systems remains a long-standing challenge.
There have been several works aimed at addressing scalability in the context of large state and action spaces. A common approach is to use function approximation to find an e!cient representation of the Q-function. For instance, [Stankovic and Stankovic, 2016] explores function approximation to solve network RL problems. However, their setting is simpler than the network MDP considered here, as they assume fully decoupled agent dynamics, whereas we allow an agent's dynamics to depend on the states of neighboring agents. In the broader context of multi-agent learning, function approximation has also been widely studied [Zhang et al., 2018, Dubey andPentland, 2021]. However, these works di"er from ours: [Zhang et al., 2018] focuses on the stochastic game setting, where agents share a common global state, while [Dubey and Pentland, 2021] examines the parallel MDP setting. Additionally, outside the RL domain, there are works on network representation learning [Dong et al., 2020, Li and Pi, 2020, Zhang et al., 2020]. However, it remains unclear whether these techniques can be applied to control and RL in network systems, which presents an interesting open question for future research.
Finding a suitable representation for the Q-function is not a unique problem in network RL. It is also a central challenge in classical centralized or single-agent RL. However, picking the right class of function approximators that can represent the Q function while being sample e!cient to learn is challenging. A natural approach is to use deep neural network (NN) architectures, which have great representational capacities. However, challenges to using deep NNs include sampleine!ciency (deep NNs often require huge amounts of data to train), hyperparameter sensitivity, stability of training (it is known that TD learning with nonlinear function approximation may fail to converge), and it can be di!cult to pick an appropriate architecture for the problem setting at hand. One promising approach arises in the (low-rank) linear MDP setting [Jin et al., 2020a], where the transition kernel of the MDP can be represented as a linear combination of low-rank features. By applying the Bellman equation, the Q-value function can then be represented as a linear combination of these low-rank features. It has been shown in [Jin et al., 2020a] in this setting, e!cient RL can be achieved, with sample complexity depending on the dimension of the feature space rather than the size of the state and action spaces. Moreover, computationally, to realize the theoretical promise of linear MDPs, there has been a line of work [Ren et al., 2022c, Ren et al., 2022a, Zhang et al., 2022] that show that the transition can be e"ectively approximated by a linear decomposition of nonlinear features, with strong empirical performances. Notably, [Ren et al., 2022c] explore the connection between stochastic nonlinear dynamical systems and linear MDPs, showing that under certain noise assumptions, stochastic nonlinear dynamics can be well-approximated by a linear decomposition of finite-dimensional (nonlinear) spectral features through an approach called spectral dynamic embedding. Building on these spectral features, [Ren et al., 2023] developed RL algorithms, with strong theoretical guarantees and empirical performances. Given the existing literature, the following question remains open:
Can we identify an appropriate representation for network MDPs and leverage it for scalability in both the size of the network and state-action space?
Our contribution Building on the existing literature, this paper addresses the critical gap by proposing a spectral dynamic embedding-based representation and developing a multi-agent RL algorithm for network systems that scales e!ciently with both network size and the complexity of state and action spaces, while also providing provable convergence guarantees.
Our approach integrates insights from both network RL and linear MDP/representation-based approaches in centralized RL. Specifically, utilizing the exponential decay property and local nature of the transition dynamics, we show how we can approximate the local Q i -value function linearly via network ω-local spectral features that factorize the ω-hop transition dynamics.
Leveraging this property, we develop a scalable samplee!cient method to learn local Q-functions in continuous network MDPs, followed by policy optimization based on the learned Q-functions.
We provide rigorous sample complexity guarantees for our framework, and to the best of our knowledge, this is the first work to propose a provably e!cient multiagent RL algorithm for network systems that is scalable with respect to both network size and the size of the state and action spaces of individual agents. Finally, we validate our approach with numerical experiments on network thermal control and Kuramoto oscillator synchronization. In both cases, we find that our approach provides benefits over generic neural network function approximations, demonstrating the advantages of our spectral representation-based framework.
Notations For any vectors v 1 , . . . , v n → R d , the notation ↑ n i=1 v i → R nd denotes their tensor product. The inner product of two tensor products is defined as follows. Consider another set of vectors w 1 , . . . , w n → R d . Then, we denote ↓↑ n i=1 v i , ↑ n i=1 w i ↔ := n i=1 ↓v i , w i ↔ . We also use the notation [n] to denote the set {1, . . . , n} for a positive integer n. In addition, when the context is clear, for notational convenience, we may drop the time indices and denote (s(t), a(t), s(t + 1)) as (s, a, s ↑ ).
Network Markov Decision Process (MDP) We consider the network MDP model, where there are n agents 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 ≃ R S and A i ≃ R A are bounded compact sets. At each time t → N, the global state of the network is denoted as s(t) = (s 1 (t), . . . , s n (t)) → S := S 1 ↘ . . . S n . Similarly, the global actuation of the network at each time t is denoted as
We also introduce the following notations related to ω-hop neighborhoods. Let N ω i denote 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 agent's ω-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 ω-hop neighborhood respectively. Similarly, we write a as (a N ω i , a N ω →i ). When ω = 1, for simplicity we denote N i := N 1 i . We assume that the next state of each agent i depends only on the current states and actions of its neighbors, so that the probability transition admits the following factorization
where N i indicates the neighbors of agent i, and s Ni (t) denotes the states of the neighbors of agent i at time t. Further, each agent is associated with a stage reward function r i (s Ni , a Ni ) 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 Ni , a Ni ); for simplicity, in the rest of our paper, we will assume that r i depends only on (s i , a i ), but we note that our analysis carries with minimal changes when r i depends on (s Ni , a Ni ). The objective is to find a (localized) policy tuple ε
) depends only on a ω ε -hop neighborhood, such that the discounted global stage reward is maximized, starting from some initial state distribution µ 0 ,
Next, we give the Kuramoto oscillator synchronization problem as an example of continuous state-action network MDPs. This example will be used in our simulations in Section 5 later. For space reasons, we defer another example, that of thermal control of multi-zone buildings, to Appendix 7.2.
Example 1 (Kuramoto oscillator synchronization). The Kuramoto model [Acebrón et al., 2005, Dorfler andBullo, 2012] is a well-known model of nonlinear coupled oscillators, and has been widely applied in various fields, ranging from synchronization of neurons in the brain [Cumin and Unsworth, 2007], to synchronization of frequency of the alternating current (AC) generators or oscillators [Filatrella et al., 2008]. Concretely, we consider here a Kuramoto system with n agents, with an underlying graph G = (N , E), where N = {1, . . . , n} is the set of agents and E ↗ N ↘ N is the set of edges. The state of each agent i is its phase ϖ i → [⇒ε, ε], and the action of each agent is a scalar a i → A i ≃ R in a bounded subset of R. The dynamics of each agent is influenced only by the states of its neighbors as well as its own action, satisfying the following form in discrete time [Mozafari et al., 2012]:
Above, ϱ i denotes the natural frequency of agent i, dt is the discretization time-step, K ij denotes the coupling strength between agents i and j, a i (t) is the action of agent i at time t, and ς i (t) ⇑ N (0, φ 2 ) is a noise term faced by agent i at time t. We note that this fits into the localized transition considered in network MDPs. For the reward, we consider frequency synchronization to a fixed target ϱ target . In this case, the local reward of each agent can be described as r i (ϖ Ni , a i ) = ⇒ πi ⇒ϱ target .
To provide context for what follows, we review a few key concepts in RL. First, fixing a localized policy tuple ε = (ε 1 , . . . , ε n ), the Q-function for this policy ε is:
In the last step, we defined the local Q-functions Q ε i (s, a) which represent the Q functions for the individual reward r i . Correspondingly, we can also define the local value function
We note that the global Q(s, a) function can be obtained by averaging n local Q i (s, a) functions. This plays an important role due to the following result known as the policy gradient theorem, which states that the policy gradient can be computed with knowledge of the Q(s, a) function.
Lemma 0 ([Sutton et al., 1999]
A natural approach to learning the Q(s, a) function in the networked case is for each agent to learn its local Q i (s, a) function and share it across the network to form a global average. However, this poses a significant challenge when (i) the network size n is large, and (ii) the individual state and action spaces S i and A i are continuous. Even if S i and A i are finite, representing
n entries, which grows exponentially with n. This challenge worsens with continuous spaces, which have infinite cardinality.
To address this, we first explore the exponential decay property from prior work, which improves scalability with network size. We then present our main contribution: integrating the exponential decay property with spectral representations from single-agent RL to derive scalable local Q i -value function representations for continuous state-action network MDPs. We begin by discussing the exponential decay property.
The exponential decay property [Qu et al., 2020b, Qu et al., 2020a, Lin et al., 2021b] is defined as follows. Definition 1. Given any c > 0 and 0 < ↼ < 1, the (c, ↼)-exponential decay property holds for a policy ε if given any natural number ω, for any i → N , s
As an immediate corollary, it follows that
We defer discussion about when the exponential decay property holds to Appendix 7.3. The power of the exponential decay property is that it immediately guarantees that the dependence of Q ε i on other agents shrinks quickly as the distance between them grows, such that the true local Q i (s, a)-functions can be approximated by truncated Qi (s N ω i , a N ω i )-functions up to an error that decays exponentially with ω. The truncated Qi function is significantly easier to represent in the finite state-action setting since each agent only needs to keep track of (|S i | ↘ |A i |) ω entries. However, continuous state and action space problems still pose a significant challenge. To overcome this, we will use the idea of spectral representations from linear MDPs and show how this can be adapted to the networked setting to yield truncated functions.
3 Spectral representations for truncated approximations of local Q i -value functions To recap, the key question we face is this: how can we derive scalable local Q i -value function representations in network problems with continuous state-action spaces, and integrate them into a scalable control framework? This forms the main contribution of our work. In this section, we tackle this question by demonstrating that the spectral representation of local transition kernels provides an e"ective representation for the local Q i -value functions (see Lemma 3 below). We first motivate our analysis by reviewing representation learning in centralized RL via spectral decompositions [Jin et al., 2020b, Ren et al., 2022b]. From such works, we know that if the global P (s ↑ | s, a) admits a linear decomposition in terms of some spectral features ↽(s, a) and µ(s ↑ ), then the Q(s, a)-value function can be linearly represented in terms of the spectral features ↽(s, a). In the case of representing local Q i -functions, this property can be stated as follows. Lemma 1 (Representing local Q i -value functions via spectral decomposition of P (Linear MDP in [Jin et al., 2020b])). Suppose the probability transition P (s ↑ | s, a) of the next state s ↑ given the current (s, a) pair can be linearly decomposed as P (s ↑ | s, a) = ↽(s, a) ≃ µ(s ↑ ) for some features ↽(s, a) → R D and µ(s ↑ ) → R D , which we also refer to as spectral representations. Then, the local Q i -value function admits the linear representation
where ↽i (s, a) := [r i (s i , a i ), ↽(s, a)],
Remark 1. We note that Lemma 1 requires the assumption of the existence of a linear decomposition of the transition kernel. One significant such example of where a linear decomposition of the transition kernel is possible was discussed in [Ren et al., 2023]. The authors in [Ren et al., 2023] showed that for a wide class of stochastic control setting with Gaussian noise (or more generally, noise which take the form of a positive-definite kernel), the transition kernel admits an exact (but infinite-dimensional) linear decomposition. By considering a finite-dimensional truncation of these infinite-dimensional features, it can be shown that the linear decomposition holds approximately with finite-dimensional features, and rigorous approximation error bounds can be shown for these finite-dimensional truncations, as shown later in Lemma 5 of our paper as well as [Ren et al., 2023].
More generally, a linear decomposition of the kernel also exists when the transition displays a particular latent variable structure [Ren et al., 2022a].
The benefit of the spectral decomposition property is that the Q i -value functions can be represented by a (D+ 1)-dimensional representation ↽i (s, a) comprising the spectral representation ↽(s, a) → R D and local reward r i (s i , a i ) → R; as demonstrated in [Jin et al., 2020a], under appropriate normalization conditions on the norm of the features, the sample complexity of using RL using such features will only depend polynomially on the feature dimension D, rather than the number of the states and actions.. However, applying this result directly in the networked case poses significant challenges, since the required feature dimension D + 1 may be high. To see why this is the case, recall that the probability transition in our networked setting admits the following factorization:
Suppose for the sake of discussion that each agent's transition probability satisfies the following Property
1, which states that the local transition has an exact d-dimensional spectral/linear decomposition. As discussed earlier in Remark 1, such a property can be expected to hold (at least approximately) for many problems, in particular stochastic dynamics under Gaussian noise. Property 1. For any i → [n] and any state-action-next state tuple (s, a, s ↑ ), there exist features ↽i (s Ni , a Ni ) → R d and μi (s
Given the factorization of the dynamics, this implies that
We note that the above expression of P (s ↑ | s, a) as an inner product between two tensor products follows from the definition of the tensor product introduced earlier in our notations; in particular, note that ↽(s, a) and μ(s ↑ ) are both in the space R d n . Intuitively, to interpret why the tensor product arises, we note that in the absence of any structure and correlations between the agents, the features factorizing the transition of the entire network are essentially the outer/tensor product of the local transition-factorization features of the n individual agents. Agnostically, Property 1 means that representing the global network dynamics may require using the d n -dimensional features ↽(s, a) := n i=1 ↽i (s Ni , a Ni ), which even for small d is undesirable due to an exponential dependence on the network size n.
While the exponential decay property suggests that the Qi -function can be approximated by considering a ω-hop neighborhood of agent i, it is unclear how we can combine this with the spectral decomposition property to derive scalable representations for the local Q i -value functions.
To resolve this, we combine insights from both the exponential decay and spectral decomposition property, which intuitively, suggests that what matters in determining Q ε i (s, a) is the probability transition dynamics within a ω-hop neighborhood of agent i. In fact, due to the local factorization property of the dynamics, the evolution of ω-hop neighborhood only depends on the ω + 1-hop neighborhood, which, when Property 1 holds, admits the following spectral decomposition.
Property 2 (Network ω-local spectral features). For any i → [n] and any state-action-next state tuple (s, a, s ↑ ), there exist some positive integer d i,ω and features
Property 2 is a statement that a linear decomposition of the ω-hop neighborhood transition kernel for any agent i is possible. As we see in Lemma 3 later, the features ↽ i,ω (s
that arise in this decomposition can represent the local Q i -value function with error decaying exponentially in ω. While Property 2 can be viewed as a standalone property independent of Property 1, the following lemma shows that when Property 1 is true, Property 2 automatically holds, with ↽ i,ω and µ i,ω given by appropriate tensor products of the original ↽i and μi from the factorization of the local dynamics. We defer the proof to Appendix 7.4. Lemma 2. Property 2 holds whenever Property 1 holds, by setting
i μj (s j (t + 1)). Remark 2. While the tensor product representation in Lemma 2 can be used to give a factorization of the ωtransition dynamics in Property 2, for specific problems, there may exist problem-specific alternative ↽ i,ω and µ i,ω features that may be lower-dimensional and thus more tractable to use.
Property 2 presents us with a path towards scalable representation of the local Q i via factorization of the local ω-hop neighborhood dynamics and approximating Q i (s, a) by network local representations. We first formalize this in the case when the spectral decomposition is exact and error-free. When this holds, we have the following lemma which shows how Q ε i (s, a) can be approximated by network local representations. We defer the proof to Appendix 7.5.1. Lemma 3 (Local Q i approximation via network ω-local spectral features). Suppose the (c, ↼)-exponential decay property holds. Suppose Property 2 also holds. Then, for any (s, a) pair, agent i, and natural number ω, there exists an approximation Qε i which depends linearly on network ω-local spectral features ↽i,ω (s
where
with the definitions ↽i,ω (s
where
Approximation. In general, it may be impossible to find ↽ i,ω and µ i,ω that can exactly factorize the transition kernel. However, in both the unknown-model and the known-model cases, there exist ways to approximate the kernel. In the model-free case, we may leverage representation-learning techniques to approximate the spectrum of the ω-hop transition kernel, such as the spectral decomposition in [Ren et al., 2022b] which seeks to approximate the SVD of the kernel. In the model-based case, in the case when the local transition evolves according to a known dynamics function subject to Gaussian noise, we may approximate the kernel by random or Nystrom features [Ren et al., 2023].
We provide below a unified analysis for the error bound of approximating
Lemma 4. For any distribution ⇀ o over the space
, the following holds for some approximation error ς P > 0:
where
Then, by setting ↽i,ω (s
Proof. Suppose (1) holds. Then, define w
The approximation error in the bound above relies on the condition in (1) to hold. In the case when the local transition evolves according to a known dynamics function subject to a positive-definite kernel noise (e.g. Gaussian noise), we may approximate the ω-hop transition kernel with random features such that (1) holds with high probability. For clarity of exposition, we focus on the approximation error of random features for Gaussian kernels [Rahimi and Recht, 2007] 1 . In this case, our error bound is shown in the following result, whose proof we defer to Appendix 7.5.2. Lemma 5. Fix any i → [n]. Suppose the local dynamics take the form
In addition, define
We note that our result easily generalizes to any positivedefinite transition kernel noise (e.g. Laplacian, Cauchy, Matérn, etc; see Table 1 in [Dai et al., 2014] for more examples)
for some δ > 0. Then, with probability at least 1 ⇒ δ, the condition in (1) holds for every i → [n] and any distribution ⇀ o over S ↘ A, with μi,ε(s
where p ς (s ↑
) is a
Gaussian distribution with standard deviation φ ς .
The key takeaway from the above result is that under Gaussian noise and known reward and dynamics function, there exists finite-dimensional features that can, with high probability, approximately factorize the local ω-transition kernels, satisfying the condition in (1) with high probability. Moreover, from this result, we note that the required number of features to achieve this is
, which only depends on the dimension of states in the largest ω-hop neighborhood. We note that the tunable ⇁ in Lemma 5 allows greater flexibility and may be tuned to improve empirical performance [Ren et al., 2023].
As suggested in Lemma 3, ↽i,ω serves as a good representation for the local Q i -functions. Based upon this observation, this section focuses on how the local Q ifunction and subsequently a good localized policy can be learned. The algorithm contains three major steps: feature generation, policy evaluation and policy gradient.
The first step is feature generation (Lines 1 through 3), where we generate the appropriate features ↽i,ω . This comprises the local reward function as well as the spectral features ↽i,ω (s N ω i ,a N ω i ) coming from the factorization of the local ω-hop dynamics. In the case of known dynamics and Gaussian noise, we know from Lemma 5 that ↽i,ω (s
) can be derived by random features which factorize the local ω-hop dynamics with high probability. In this case, we note that our spectral features are scalable with respect to both the network size and the continuous state-action space, since the required number of features only depend on the dimensions of the ω-hop neighborhoods.
The second step is policy evaluation, where we use the feature ↽i,ω and apply LSTD to find a set of weights
we first sample M s samples from the stationary distribution of ε (k) (Line 5), and then perform a LSTD update for each agent i → [n] to learn the appropriate weights for the local Q i -functions (Line 6).
Finally, the last step is updating policy using policy gradient (Lines 6 to 7). For each agent i → [n], with the learned { Qj } j↘N ωε +ω i , we perform a gradient step to update the local policy weights ϖ i , and update to the new policy. We note that this update is scalable since from the perspective of each agent i, it only requires knowledge of the local Qj for agents j in a (ω ε + ω)hop neighborhood of agent i. In practice, the ω-hop spectral representation we introduce can be combined flexibly with any actor in a distributed cooperative actor-critic framework that requires knowledge of the local Q i -functions.
We have the following result on the policy evaluation error with our network ω-local features. We defer the details of the proof (including preliminary results required for the proof) to Appendix 7.6. Lemma 6 (Policy Evaluation Error). Suppose equation (1) in Lemma 4 holds. Suppose the sample size
. Then, with probability at least
where
is defined as in equation ( 4).
From the above result, we note that the policy evaluation error comprises three components, with one being the statistical error due to using finite samples, which decays with the square root of the sample size M s , and the truncation error from considering a truncated ωhop neighborhood (this decays exponentially in ω), as well as the approximation error of the spectral features in approximating the ω-hop transition (ς P ).
. Suppose with probability at least 1 ⇒ δ, for all i → [n], the following holds for some features ↽i,ω → R m and μi,ω → R m :
,where θi,ε(s
) satisfies the condition in (1).
Policy evaluation and update
i using least square temporal di"erence (LSTD) as follows:
where ⇀i,ε(j) := θi,ε(s
Update approximate
Policy gradient for control
for some ς P > 0 and distribution ⇀ o over S ↘ A, where
where L ↑ is the Lipschitz continuity parameter of
From the above result, we see that our algorithm can achieve convergence to an approximate stationary point of the global objective J as the number of rounds K increases, up to an error term depending on ς J , which depends on the policy evaluation error ς Q from Lemma 6. As we observed before, the policy evaluation error comprises a statistical error, a truncation error decaying exponentially as ω increases, and a feature approximation error term ς P . Consequently, the convergence error to an approximation stationary point also depends on these three terms.
We consider a stochastic linear dynamical system modeling the thermal control of a 50-zone building. We assume that the network is connected, with each agent having 2 neighbors. The dynamics of each agent is only a"ected by its neighbors, and subject to Gaussian noise. We also assume access to the model dynamics and reward function. In this problem, to implement our algorithm, we utilize random features that factorize the ω-hop Gaussian transition (cf. Lemma 5), and perform least squares, followed by normalized gradient descent. The controller is parameterized to be linear. More details on our experimental setup can be found in Appendix 7.9.
Since the dynamics are assumed to be linear, we have access to the cost of the optimal controller, making this a good way to benchmark the performance of our algorithm. The performance of our algorithm is shown in Figure 1 below. As we can see, as ω ε increases, our algorithm is indeed able to approximate the performance of the optimal controller. Moreover, the speed at which it converges is faster than Qi approximations that leverage a generic two-hidden layer neural network (NN) to represent the (truncated) local Qi value functions; we note that in both cases, the algorithms utilize the same learning rate for the policy gradient step, and have access to the rewards and dynamics function.
Due to the nonlinearity in this problem, we adopt the Soft-Actor-Critic (SAC) framework for this problem, and compare the performance of a generic deep NN critic with our network spectral local-ω critic. In this problem, we consider the more realistic and di!cult setting where the dynamics is unknown. In this problem, the network has 40 agents in total, and the network graph is connected, with each agent having 2 neighbors. We set the goal for the agents to synchronize to a target frequency of 0.2.
In both the generic SAC and our spectral SAC implementation, the local critic for Q i -value function considers a ω-hop neighborhood, i.e. approximate Q i by Qi (s
) is a two-hidden layer neural network. However, for our approach (spectral + SAC), we add a feature step that regularizes the feature ↽i (s N ω i , a N ω i ) towards factorizing the local dynamics, i.e. minimizing the objective in Condition 1 in Lemma 4. We defer more details on the problem setup as well as experimental details to Appendix 7.9.
In Figure 2, we compare the performance of our approach (Spectral + SAC) with a generic SAC with two-hidden layer NN critic. As we can see, our approach leads to significantly higher rewards. Moreover, we observe that our approach leads to qualitatively better synchronization behavior when starting from the same initial condition, as indicated in Figure 3. Finally, we note that in the model-based setting, our algorithm (utilizing random features) achieves a performance comparable to that of generic NN approaches. The model-based results are deferred to Appendix 7.9. (a) SAC Controller (b) Spectral+ SAC Controller Figure 3: Synchronization of frequency ( θ) under SAC and Spectral + SAC controller, for 800 time steps with time interval dt = 0.01.
In this work, utilizing local spectral representations, we provide the first provably e!cient algorithm for scalable network control in continuous state-action spaces. We validate our results numerically, where we find that utilizing ω-local spectral features can achieve e"ective control on a thermal network control problem as well as a Kuramoto nonlinear coupled oscillator control problem. Moreover, in both cases, we demonstrate that our approach has benefits over generic neural network approximations for local Q i -value functions.
Collectively, our theoretical and empirical results demonstrate the validity and importance of a representation-based viewpoint in achieving more e"ective and scalable control in continuous state-action network MDPs. by NSF ECCS 2401391, NSF IIS 2403240, as well as Dolby. games with decoupled dynamics: Price of anarchy and sample complexity. In 2023 62nd IEEE Conference on Decision and Control (CDC), pages 8100-8107. IEEE. [Zhang et al., 2023b] Zhang, R. C., Li, W., and Li, N. (2023b). On the Optimal Control of Network LQR with Spatially-Exponential Decaying Structure. In 2023 American Control Conference (ACC), pages 1775-1780. ISSN: 2378-5861. [Zhang et al., 2022] Zhang, T., Ren, T., Yang, M., Gonzalez, J., Schuurmans, D., and Dai, B. (2022). Making linear mdps practical via contrastive representation learning. In International Conference on Machine Learning, pages 26447-26466. PMLR. [Zhang et al., 2016] Zhang, X., Shi, W., Li, X., Yan, B., Malkawi, A., and Li, N. (2016). Decentralized temperature control via hvac systems in energy e!cient buildings: An approximate solution procedure. In 2016 IEEE Global Conference on Signal and Information Processing (GlobalSIP), pages 936-940. IEEE. [Zhang et al., 2023c] Zhang, Y., Qu, G., Xu, P., Lin, Y., Chen, Z., and Wierman, A. (2023c). Global convergence of localized policy iteration in networked multi-agent reinforcement learning. Proceedings of the ACM on Measurement and Analysis of Computing Systems, 7(1):1-51. [Zocca, 2019] Zocca, A. (2019). Temporal starvation in multi-channel csma networks: an analytical framework. ACM SIGMETRICS Performance Evaluation Review, 46(3):52-53. Checklist 1. For all models and algorithms presented, check if you include: (a) A clear description of the mathematical setting, assumptions, algorithm, and/or model. Yes (b) An analysis of the properties and complexity (time, space, sample size) of any algorithm. Yes (c) (Optional) Anonymized source code, with specification of all dependencies, including external libraries. Yes 2. For any theoretical claim, check if you include: (a) Statements of the full set of assumptions of all theoretical results. Yes (b) Complete proofs of all theoretical results. Yes (c) Clear explanations of any assumptions. Yes 3. For all figures and tables that present empirical results, check if you include: (a) The code, data, and instructions needed to reproduce the main experimental results (either in the supplemental material or as a URL). Yes (b) All the training details (e.g., data splits, hyperparameters, how they were chosen). Yes (c) A clear definition of the specific measure or statistics and error bars (e.g., with respect to the random seed after running experiments multiple times). Yes (d) A description of the computing infrastructure used. (e.g., type of GPUs, internal cluster, or cloud provider). Yes 4. If you are using existing assets (e.g., code, data, models) or curating/releasing new assets, check if you include: (a) Citations of the creator If your work uses existing assets. Yes (b) The license information of the assets, if applicable. Not applicable (c) New assets either in the supplemental material or as a URL, if applicable. Yes (d) Information about consent from data providers/curators. Not applicable (e) Discussion of sensible content if applicable, e.g., personally identifiable information or offensive content. Not applicable 5. If you used crowdsourcing or conducted research with human subjects, check if you include: (a) The full text of instructions given to participants and screenshots. Not Applicable (b) Descriptions of potential participant risks, with links to Institutional Review Board (IRB) approvals if applicable. Not applicable (c) The estimated hourly wage paid to participants and the total amount spent on participant compensation. Not applicable 7.1 More Related Works Multi-agent reinforcement learning (MARL) for Markov games The study of MARL dates back to the early work of [Littman , 1994, Claus and Boutilier, 1998, Hu and Wellman, 2003]. One classical setting considered is the stochastic game setting [Shapley, 1953], where agents can take their own actions but they share a common global state and maximize a global reward. There are many recent works that studies reducing the computational and sample complexity in this scenario, e.g. [Zhang et al., 2018, Leonardos et al., 2022, Song et al., 2021, Zhang et al., 2024]. However, one restriction of their approach is that it assumes access to the global state s, as opposed to our setting where each agent only observes a local state s N ω i , thus their approach won't extend easily to the network setting considered in our paper.
Distributed Control for Network Linear Quadratic Regulator (LQR) Apart from reinforcement learning, there is also another line of work tackling the network system from the control perspective. Distributed control is a classical topic that is widely discussed in literature (c.f. [Saeks, 1979, Bamieh et al., 2002, Rotkowitz and Lall, 2005, Rantzer, 2011]) Notably, some recent works [Zhang et al., 2023b, Shin et al., 2023, Olsson et al., 2024] also leverage a similar exponential decaying property to address the near-global-optimality of distributed controllers. However, these works are primarily focused on linear dynamical systems, as opposed to our work which is targeted to potentially nonlinear dynamics via representation learning Other multi-agent dynamical models There are also other settings in multi-agent/network systems such as the weakly coupled MDP [Meuleau et al., 1998, Zhang et al., 2023a], where agents' transition dynamics are fully decoupled and the only coupling is through the reward function; factored MDP [Kearns and Koller, 1999, Guestrin et al., 2003, Osband and Van Roy, 2014], where there is a global action a"ecting every agent's individual local state. (See the 'Related Literature' in [Qu et al., 2021] for a more detailed summary and comparison.)
Linear MDPs There has been a sequence of work on sample-e!cient RL via the linear MDPs approach [Jin et al., 2020a]. In a linear MDP, the transition kernel of the MDP can be represented as a linear combination of low-rank features. By applying the Bellman equation, the Q-value function can then be represented as a linear combination of these low-rank features. It has been shown in [Jin et al., 2020a] in this setting, sample-e!cient RL can be achieved, with sample complexity depending on the dimension of the feature space rather than the size of the state and action spaces. Computationally, to realize the theoretical promise of linear MDPs, there has been a line of works with strong empirical performances that represents the Q-value function using finite-dimensional features that factorize the transition kernel, where the features are learnt in di"erent ways, e.g. a latent variable approach [Ren et al., 2022a], a noise-contrastive approach [Zhang et al., 2022] building on noise contrastive estimation [Gutmann andHyvärinen, 2010, Ma andCollins, 2018], and a di"usion-inspired approach [Shribak et al., 2024]. A key question in the study of linear MDPs is the circumstances under which a lowrank linear factorization of possibly nonlinear features that represent the transition kernel indeed exists. Towards answering this question, the works in [Ren et al., 2022c, Ren et al., 2023] show that for a wide range of problems in stochastic nonlinear control, specifically stochastic nonlinear control problems where the transition noise takes the form of a positive-definite kernel (e.g. Gaussian noise), there exists infinite-dimensional spectral features which exactly factorizes the transition kernel. However, to enable tractable control, finite-dimensional features are required. Tn [Ren et al., 2023], the authors propose a finite-dimensional approximation of the infinite-dimensional spectral features via random features [Rahimi and Recht, 2007] and Nystrom features [Williams and Seeger, 2000], characterize the approximation error of these finite-dimensional truncation approaches, and provide end-to-end theoretical guarantees for a actor-critic framework building on these finite-dimensional features, with strong empirical performance.
Example 2 (Thermal control in buildings). The problem of thermal control of multiple zones in a building can also be cast as a network MDP. Consider a multi-zone building with a Heating Ventilation and Air Conditioning (HVAC) system. Each zone is equipped with a sensor that can measure the local temperatures and can adjust the supply air flow rate of its associated HVAC system. For simplicity, we consider a discrete-time linear thermal dynamics model based on [Zhang et al., 2016, Li et al., 2021], where for any i → [n],
where x i (t) denotes the temperature of zone i at time t, a i (t) denotes the control input of zone i that is related with the air flow rate of the HVAC system, ϖ o (t) denotes the outdoor temperature, ε i represents a constant heat from external sources to zone i, w i (t) represents random disturbances, ! is the time resolution, v i is the thermal capacitance of zone i, ▷ i represents the thermal resistance of the windows and walls between zone i and the outside environment, and ▷ ij represents the thermal resistance of the walls between zone i and j. Again, we note that the transition dynamics of each agent depends only on its neighbors (and itself ). At each zone i, there is a desired temperature ϖ → i set by the users. The local reward function is composed of the (negative) deviation from the desired temperature and the control cost, i.e.
where ⇁ i > 0 is a trade-o" parameter.
It may not be immediately clear when the exponential decay property holds. The following lemma (cf. Appendix A in [Qu et al., 2020b]) highlights that for a local policy where each agent's actions depend only only on its and its neighbors' states (i.e.
), the exponential decay property holds generally, with ↼ = ϑ. We defer the proof to our appendix.
Lemma 7. Suppose ⇔i → [n], agent i adopts a localized policy, i.e.
Suppose also that the local rewards are bounded such that 0 ↭ r i ↭ r. Then the r 1↓ϖ , ϑ -exponential decay property holds.
We note that under some mixing time assumptions on the MDP [Qu et al., 2020b], the exponential decay property may in fact hold for ↼ < ϑ depending on the system parameters, making it applicable to problems with large discount factors or even in the average-reward setting [Qu et al., 2020a].
We proceed now to prove Lemma 7, which shows that the exponential decay property holds for localized policies and bounded rewards. We note that this was first shown in [Qu et al., 2020b], and we provide the proof here for completeness.
Proof. Consider any i, and choose any natural number ω. For an arbitrary (s, a) = (s N ω i , s N ω →i , a N ω i , a N ω →i ), consider any state-action pair (s ↑ , a ↑ ) that di"ers with (s, a) only outside the
For any natural number t, let p t,i denote the distribution of s i (t), a i (t) conditional on s(0) = s, a(0) = a, and let p ↑ t,i denote the distribution of s i (t), a i (t) conditional on s(0) = s ↑ , a(0) = a ↑ . Then,
Above, (i) is a direct application of the definition of p t,i and p ↑ t,i . Meanwhile, (ii) utilizes the fact that for any 0 ↭ t ↭ ω, we have p t,i ⇐ p ↑ t,i . This is because of (a) localized policy, such that a i (t) depends only on s Ni (t ⇒ 1), and (b) factorized localized dynamics, that s N j i (t) depends only only on s N j+1 i (t) for any natural number j; hence an iterative argument shows that for any t, p t,i and p ↑ t,i both only depend on s N t i (0) and a N t i (0). Thus, since (s, a) and (s ↑ , a ↑ ) share identical s N ω i (0) and a N ω i (0), it follows that p t,i ⇐ p ↑ t,i for t ↭ ω. Finally, (iii) uses the fact that bounded reward assumption, i.e. 0 ↭ r i ↭ r. The proof then concludes by rerunning the argument on
Next, we state and prove the following elementary technical result, which bounds any two truncated Q(or V )-functions with di"erent weights.
Lemma 8. Suppose the (c, ↼)-exponential decay property holds. Then, for any two di"erent weights
where
Similarly, for any two di"erent weights w i (s
we have
where
Proof. Compare both truncated Q-functions to a Q-function evaluated at any specific state action pair where the states and actions of the agents in N ω i are s N ω i and a N ω i respectively. The desired result then follows by Definition 1. A similar argument works for the V -function.
Lemma 2. Property 2 holds whenever Property 1 holds, by setting
Proof. To see that, observe that
Above, (iv) follows from Property 1, while (v) uses the definition of the inner product of two tensor products. Thus, when Property 1 holds, Property 2 holds by setting
7.5 Approximation error of spectral features
We recall and prove this result, which bounds the approximation error of using the truncated spectral features to approximate the local Q i -function, in the case when there is no approximation error in the spectral features in representing the ω-hop transition. Lemma 3 (Local Q i approximation via network ω-local spectral features). Suppose the (c, ↼)-exponential decay property holds. Suppose Property 2 also holds. Then, for any (s, a) pair, agent i, and natural number ω, there exists an approximation Qε i which depends linearly on network ω-local spectral features ↽i,ω (s
where
with the definitions ↽i,ω (s
where
Proof. For notational convenience, we omit the t and (t + 1) in the parentheses of the state and action notations, and instead use a superscript + to denote (t + 1), e.g. s + to denote s(t + 1). Observe that
Scalable spectral representations for multiagent reinforcement learning in network MDPs
In (vi) above, we used the notation
and in (vii), we recall that we defined
Since the (c, ↼)-exponential decay property holds, applying Lemma 8, we have
Thus our desired result holds by setting
We first state the following result on uniform convergence of random Fourier features, adapted from Claim 1 in [Rahimi and Recht, 2007]. Lemma 9 (Uniform convergence of Fourier features). Let M be a compact subset of R d with diameter diam(M).
Let k be a positive definite shift-invariant kernel k(x, y) = k(x ⇒ y). Define the mapping z, where
where ϱ 1 , . . . , ϱ D → R d are D iid samples from p, where p is the Fourier transform of k, i.e. p(ϱ) =
1 2ε
e ↓j⇀ ↔ ↽ k(δ)dδ, and b 1 , . . . , b D are D are iid samples from Unif(0, 2ε). We assume that k is suitably scaled such that p is a probability distribution.
Then, for the mapping z defined above, we have
,
is the second moment of the Fourier transform of k. Further, sup
Then, for a positive integer m, define the m-dimensional features ↽i,ω (s
In addition, define
for some δ > 0. Then, with probability at least 1 ⇒ δ, the condition in (1) holds for every i → [n] and any distribution ⇀ o over S ↘ A, with μi,ε(s
where p ς (s ↑
) is a Gaussian distribution with standard deviation φ ς .
Proof. Observe that P s +
In addition, in this proof we denote
where
it follows that with probability at least 1 ⇒ δ,
When ϱ := 0, this simplifies to
. However, we allow a general 0 ↭ ϱ < 1 because it gives greater flexibility both theoretically and empirically.
Here we overload notation to use fi,ε(S
and {ϱ ϱ }'s are drawn iid from N (0, φ ↓2 I d ), {b ϱ }'s are drawn iid from Unif(0, 2ε). It follows then for any
where we recall that ↽i,ω (x) := g ς ⇁ d ↽i,ω (x), μi,ω (y) := p ς (y)μ i,ω (y). The proof then follows by rescaling ς P := ς P /n, and taking a union bound over all i → [n].
For simplicity, we assume throughout the analysis that we are solving the LSTD step in the policy evaluation exactly, i.e. we take the number of least square solves, T , to infinite. Moreover, we drop the i subscript in the notation of ↽i,ω , and use ⇀ to denote ⇀ ε (k) . At round k, the algorithm output of the policy parameter w i of agent i is given by
For notational convenience, when the context is clear, we drop the k-superscript indicating the current round k, and denote
We define an intermediate variable wi as follows:
and further define
(s, a). From Lemma 3 and Lemma 4, we have that there exists
such that with probability at least 1 ⇒ δ, for every i → [n],
Lemma 10 (Bellman Error). On the event that condition (1) in Lemma 4 holds, for every i → [n], we have
From Bellman equation we have that
On both side multiply by ↽i,ω (s
) and take expectation over s, a ⇑ ⇀, we have that
Further, given that
on the event that condition (1) in Lemma 4 holds.
where for the final inequality we used (5).
Lemma 11 (Statistical Error).
Proof. Again, for notational simplicity we drop the k-superscript. We first bound the di"erences
From the Matrix Bernstein inequality (see Lemma 15 in Appendix 7.8) we have that when
Thus with probability 1 ⇒ 2δ
which completes the proof.
Combining the above statement we can get the following Lemma for policy evaluation error, which is a restatement of our result in Lemma 6.
where denoting ◁i,ω := ↽i,ω (s
↙ ◁i,ω ↙, where
Proof. Suppose the condition (1) in Lemma 4 holds. Consider any i → [n] and k → [K]. From Lemma 10 and 11 we have that with probability at least 1 ⇒ 2δ,
The desired result then follows by rescaling δ := δ/(Kn) and taking an union bound over all i → [n] and k → [K].
We show now that our algorithm can find an approximate stationary point of the averaged discounted cumulative reward function J(ε (ϑ) ). For notational convenience, for a given set of policy parameters ϖ = (ϖ 1 , . . . , ϖ n ), we define
where we recall that r(s, a) := 1 n n i=1 r i (s i , a i ). From Lemma 0, we have that
We first provide the following result, which shows that assuming Lipschitz continuity of the gradient of the objective function J as well as the gradients of log ε ϑ , there exists the following bound on the following weighted sum of the squared gradient norms.
Lemma 13. Suppose that ⇓J(ϖ) is L ↑ -Lipschitz continuous. Suppose that for each i → [n], k) . Then,
Proof. By the Lipschitz continuity of ⇓J(ϖ), we have
By rearranging and using a telescoping sum, we obtain
Recall that J(ϖ (K+1) ) ⇒ J(ϖ (1) ) ↭ r/(1 ⇒ ϑ). Hence, by the given assumption on the bound on the derivative term ⇓ ϑi log ε ϑi i , it follows that ↙⇓J(•)↙ ↭ Lε r 1↓ϖ . The desired bound then follows by plugging this in as well as using the triangle inequality to decompose
As we can tell from the above result, the crux to bounding the average stationarity gap after K rounds of optimization is the di"erence between the true gradient ⇓J(ϖ (k) ) and the learned gradient ĝ(k) used in the update.
In this next result, we bound this error, assuming that the truncated local Qi -functions are learned up to some error.
Lemma 14. For any optimization round k → [K], let ⇀(k) denote the empirical distribution of the samples used during round k, i.e. {s(j), a(j)} j↘ [Ms] where (s(j), a(j)) ⇑ ⇀ ε (k) . Suppose that for each 0 → [n], the learnt Qϱ -value function satisfies the following error bound:
where the i-th component of the approximate gradient
is defined in Line 9 of Algorithm 1.
Proof. For notational convenience, in the proof, we fix the optimization round k → [K], and hence, denote ĝi := ĝ(k) i , ϖ := ϖ (k) and Qϱ := Q(k) ϱ unless otherwise specified. Moreover, we also denote Q ϑ := Q ε ϑ for simplicity. From Lemma 0, for any agent i → [n], we have that
To bound the di"erence between ĝi and ⇓ ϑi J(ϖ), we define the following intermediate terms.
We define the terms
Then, we decompose the error as
We proceed now to bound the three error terms in (8).
Error term E J,h . We can bound the term E J,h as follows. For any 0 → [n] and positive integer ω, we define
where we let w((s N ω ϱ ) ↑ , (a N ω ϱ ) ↑ ) denote the uniform weight over the space S N ω →ϱ ↘ A N ω →ϱ . From Lemma 8, we know that
We then have
To bound E J,h,1 , utilizing the bound in (9) as well as the bound ⇓ ϑi log ε ϑi i (a i | s ω Ni ) ↭ L i,ε in the statement of the lemma, we have that
Meanwhile, observe that by definition, for any 0
Thus E J,h,2 = 0. This implies then that
Error term E h,g . To bound E h,g , we may use standard concentration inequalities. Observe that
.
Since
it follows that E E h,g(j) = 0. Moreover, using the fact that for any 0 → [n], ϖ and (s, a) pair, 0 ↭ Q ϑ ϱ (s, a) ↭ r 1↓ϖ , and the bound
in the assumption, we have ↙E h,g (j)↙ ↭ 2rLε 1↓ϖ . Using the i.i.d. assumption between the samples j → [M s ], we may apply Bernstein's concentration inequality for vectors (see Lemma 15) to find that for any δ > 0, with probability at least 1 ⇒ δ,
where d ϑ is the dimension of ϖ i .
Error term E g,ĝ . Observe that
Above, (ix) follows from the bound in (6), as well as the bound
Combining the bounds for E J,h , E h,g and E g,ĝ , we find that with probability at least 1 ⇒ δ,
The final result then follows by applying a union bound over k → [K].
We are now ready to state our main convergence result.
Theorem 2 (Restatement of Theorem 1). Suppose the sample size M s ↫ log 2dω (↽/Kn)
. Suppose with probability at least 1 ⇒ δ, for all i → [n], the following holds for some features ↽i,ω and μi,ω :
Proof. Fix a δ > 0. Suppose the condition in (1) holds with probability at least 1 ⇒ δ for all i → [n] for a distribution ⇀ o over S ↘ A. In other words, with probability at least 1 ⇒ δ, for all i → [n], the following holds:
for some ς P > 0. Then, by Lemma 6, it follows that with probability at least 1 ⇒ 3δ, for every i → [n] and optimization round k → [K], we have
Note that by Lemma 13, with probability at least 1 ⇒ δ, for every optimization round k → [K], we have
Thus, by picking η = O(1/ ↖ K), using union bound, with probability at least 1 ⇒ 4δ, we have
Lemma 15 (Matrix Bernstein). Suppose {M k } n k=1 are i.i.d random matrices where M k → R d1⇐d2 and that
then for a given δ → (0, 1) and n ↫ log
Now we can apply the matrix Bernstein inequality (Theorem 6.1.1 in [Tropp et al., 2015]) and get that Pr
which completes the proof.
All code for this project is available as a zip folder with the supplementary material.
Problem setup details. In the simulations, we consider a discrete-time linear thermal dynamics model adapted from [Zhang et al., 2016, Li et al., 2021], where for any i → [n],
where x i (t) denotes the temperature of zone i at time t, a i (t) denotes the control input of zone i that is related with the air flow rate of the HVAC system, ϖ o (t) denotes the outdoor temperature, ε i represents a constant heat from external sources to zone i, w i (t) represents random disturbances, ! is the time resolution, v i is the thermal capacitance of zone i, ▷ i represents the thermal resistance of the windows and walls between zone i and the outside environment, ▷ ij represents the thermal resistance of the walls between zone i and j, and ⇁ i and 1 i denote scaling factors on the input and noise respectively. The local reward is defined as
where ϖ → i is the target temperature and ↼ i is a trade-o" parameter. The parameters in the dynamics and rewards are set as follows. For simplicity, we center the temperatures at 0, and hence set the target ϖ → i to be 0. We set ↼ i = 3. We set the following parameters for the dynamics: ! = 20, v i = 200, ▷ ij = 1, ▷ i = 1 2 , ⇁ i = 1 7 , 1 i = vi ! , ϖ 0 = 0. We also assume w i (t) to be drawn iid from N (0, 1). We set the discount factor in the problem to be 0.75, and (when collecting data) set the horizon length of each episode to be 20.
In this problem, there are n = 50 agents, and the agents have circular connectivity and has two neighbors each, such that agent 1 is connected to agents N and agent 2, agent 2 is connected to agents 1 and 3, so on and so forth.
Experimental details. We assume knowledge of the dynamics and rewards. For policy truncation parameter ω ε = 0, 1, 2, 3, we use ω = 0, 1, 2, 2 respectively 4 as the evaluation ω parameter. We now explain the simulation setup for our implementation of Algorithm 1 with random features, as well as the benchmark algorithm using a two-hidden layer NN.
1. (Spectral embedding generation step). For Algorithm 1 with random features, for each agent i, we use random feature dimension of m = 30, 50, 800, 800 for each of the four experiments (ω ε = 0, 1, 2, 3) to represent the function T ε (s
)↓ri(si,ai) ϖ
. For the NN implementation, we used a two-hidden layer NN with 128 neurons to represent the function T ε (s N ω+1 i , a N ω+1 i ).
2. (Policy evaluation step) We used M s = 100, 200, 500, 1000 episodes respectively for each of the four experiments (ω ε = 0, 1, 2, 3) to perform the policy evaluation. For the random features implementation, we used the least squares method in Algorithm 1 to compute the new weights for the local value functions. For the NN implementation, we ran batch gradient descent, and used a target network with update rate of 0.005.
For both implementations, we normalize the policy gradient, and run gradient descent with η = 0.2.
Problem setup details. We recall the setup described earlier in the paper. We consider here a Kuramoto system with n agents, with an underlying graph G = (N , E), where N = {1, . . . , n} is the set of agents and E ↗ N ↘ N is the set of edges. The state of each agent i is its phase ϖ i → [⇒ε, ε], and the action of each agent is 4
We found in practice that using κ = 3 for κω = 3 performed less well in this specific example.
Scalable spectral representations for multiagent reinforcement learning in network MDPs a scalar a i → A i ≃ R in a bounded subset of R. The dynamics of each agent is influenced only by the states of its neighbors as well as its own action, satisfying the following form in discrete time [Mozafari et al., 2012]:
ϱ i (t)+a i (t)+ j↘Ni K ij sin(ϖ j ⇒ϖ i ) := θi(t)
+ς i (t).
Above, ϱ i denotes the natural frequency of agent i, dt is the discretization time-step, K ij denotes the coupling strength between agents i and j, a i (t) is the action of agent i at time t, and ς i (t) ⇑ N (0, φ 2 ) is a noise term faced by agent i at time t. We note that this fits into the localized transition considered in network MDPs. For the reward, we consider frequency synchronization to a fixed target ϱ target . In this case, the local reward of each agent can be described as r i (ϖ Ni , a i ) = ⇒ πi ⇒ ϱ target .
The parameters in the dynamics and rewards are set as follows. We set the target ϱ target to be 0.75. We set the action space as [⇒3, 3]. For agents i and j that are connected, we sample K ij uniformly at random from [0.2, 1.2].
For the natural frequency ϱ i 's, we sample them iid uniformly at random from [0, 1.5]. For the noise, we sample ς i (t) ⇑ N (0, 0.01 2 ). The time resolution is dt = 0.01.
We also assume w i (t) to be drawn iid from N (0, 1). We set the discount factor in the problem to be 0.99, and set the horizon length to be 800 steps.
Connectivity. In this problem, there are n = 20 agents, and the agents have circular connectivity and has two neighbors each, such that agent 1 is connected to agents N and agent 2, agent 2 is connected to agents 1 and 3, so on and so forth.
Experimental details (model-free). In this case, we do not assume access to the dynamics function. We now explain the simulation setup for our implementation of Algorithm 1 with spectral features, as well as the benchmark algorithm using a two-hidden layer NN.
1. (Policy evaluation step) For both the spectral feature and NN implementation, the features are the last layer of a two-hidden layer neural network with hidden dimension 256. At each iteration, for each agent i, we draw a batch (of 128 transitions) from the replay bu"er and we run 1 step of gradient descent on the least square bellman error, and used a target network with update rate of 0.005.
For each agent i, the policy is parameterized to be a 3-hidden layer NN which outputs the mean and standard deviation of the agent's action, and the input is s Ni , i.e. the states of the neighborhood of agent i. We update the policy parameters {ϖ i } n i=1 by taking one gradient descent step on the following objective:
where D is a set of data from the replay bu"er, Z(s) is a normalization constant. Above, we assume the temperature parameter 2 to be 1. This objective is identical to the implementation in Soft-Actor-Critic (SAC) [Haarnoja et al., 2018] but for factored policies, as well as using the learned Qi -value functions to approximate the value function.
3. (Feature step). For the spectral features, for each agent i, we add an additional feature step, which seeks to regularize the features such that they approximate the top left eigenfunctions of the probability transition P (s ↑
), by taking a gradient descent step on the following objective to update agent i's feature
1. Fix a positive integer L.
2. For each i → [L], sample a random function ϱ i (s ↑ ) → R, e.g. ϱ i (s ↑ ) = cos(⇁ ≃ i s ↑ + 1 i ), where ⇁ i ⇑ N (0, I S ) and 1 i ⇑ Unif([0, 2ε]).
3. For each i → [L], learn a ↽ i (s, a) that approximates P ϱ i (s, a) := s ↑ P (s ↑ | s, a)ϱ i (s ↑ )ds ↑ as follows: We note that the final objective is equivalent to solving the L ↽ i 's jointly which is single-agent analogue of the objective in ( 12