Year: 2021

Abstract:
This paper proposes an optimal policy learning algorithm under a novel multi-agent predictive state representation reinforcement learning model. Compared to the state-of-the-art methods, the most striking feature of our approach is to introduce a dynamic interaction graph to the model, which allows us to represent each agent’s predictive state based on its “neighborhood” agents. Methodologically, we develop an online algorithm that simultaneously learns the predictive state representation and agent policies. Theoretically, we propose an upper bound of the L2-norm of the learnt predictive state representation. Empirically, we provide positive numerical results on both a MAMuJoCo robotic learning experiment and a multi-agent particle learning environment to demonstrate the efficacy of the method.

heading: 1 INTRODUCTION
text: Real-world multi-agent systems are considered partially observable since agents do not have a complete perception of the environment states. The predictive state representation (PSR) (Littman et al., 2001) is a representation of the dynamic system state by a vector of predictions of tests conditioned on a history. The tests are sequences of actions and observations, which are true if and only if all the observations occur, given that all the actions taken and the histories are sequences of actions and observations recorded before. It records a sequence of actions and observations and focuses on observable quantities rather than the system’s dynamics. Therefore, compared with learning the hidden-state-based partially observable Markov decision process (POMDP) Kaelbling et al. (1998) models from observation data, learning PSR should be easier and less prone to local minimal problems (Singh et al., 2012). In addition, in the cases where the moment-method is used to learn a model for controlled systems that incorporates actions, PSR is considered as a more general setting to reconstruct the structure of the dynamic system than POMDP (Azizzadenesheli et al., 2016; Hamilton, 2014). Furthermore, (Littman et al., 2001; James & Singh, 2004; Singh et al., 2012) showed that PSR can offer an alternative more compact state representation for POMDP models. Overall, PSR is considered a generalization of POMDP in terms of representation power and the development of learning algorithms.
Due to the above benefits of PSR, it has been used as the basis of some model-based single-agent reinforcement learning (James & Singh, 2004; Boots et al., 2011; Hamilton et al., 2014; Hefny et al., 2018a). Moreover, the predictive states that act as sufficient statistics for the states of multiagent dynamic systems could potentially help multi-agent decision-making in the complicated, partially observable domains. However, there is no method for developing multi-agent predictive state representation (MAPSR) for multi-agent reinforcement learning (MARL). Due to this reason, we focused on building an algorithm for learning the optimal policies by representing the multi-agent dynamic systems with the MAPSR.
An existing algorithm by (Chen et al., 2020) provides a spectral framework for learning MAPSR using a tensor, with its dimensions representing the agents and joint history. A tensor decomposition is implemented to learn the model parameters. This type of algorithm avoids explicit estimation of the latent state and can get an unbiased estimation if an infinite amount of data is available. However, it is not easy to incorporate prior information such as sparsity or structure because each new source of information leads to new moments and a different and more complex set of equations to solve. (Hefny et al., 2015) represent the distributions in PSR via observable statistics and learn PSR by supervised learning, alleviating the above problem. Later, (Hefny et al., 2018a) extend this method to
develop a new one called predictive state controlled model that can represent the system model where the agents can affect the environment. Due to the trend of deep RL, (Hefny et al., 2018b) developed an end-to-end policy learning method by using the states represented by a predictive state controlled model as the belief state to fit the policy function to train RL. Both the PSR and RL parameters are updated by using a combined loss objective. The PSR part maintains a sufficient representation of distributions under the partially observable environment, the deep RL algorithm brings efficiency for learning.
The environment of MARL can be partially observable and always non-stationary (one agent learning while the other agents’ actions can influence the environment) (Bowling & Veloso, 2002). Predictive state controlled model can model the prediction of systems states affected by actions; if we design a state representation model for MARL environment that can predict the future observations of agents upon exerting future actions of other agents and conditioned on past histories, then the non-stationary environment caused by other agents learning can be modeled. Due to the above reason, we seek to develop a MARL-based PSR model and learning algorithm, which we call MAPSRL. It could potentially ameliorate both the problem of partially observable observations and non-stationarity by other agents’ interaction for multi-agent decision making.
The paper’s first goal is to develop a PSR model for MARL environment. In MARL settings, the distribution of observed observations is dependent on its actions and relevant to other agents and histories. So far, many works (Ryu et al., 2020; Liu et al., 2020b; Niu et al., 2021) that take other agents’ actions into consideration have been shown to achieve competitive performance. In order to consider the interactive relationship under the MAPSR, we introduce a dynamic interaction graph to characterize the “neighborhoods” of agents. Such a graph topology, which is often easy to obtain nowadays, enables us to encode the adjacency relationship of agents into our MAPSR model. We build our model under three common graphs, namely static complete graphs, static non-complete graphs, and dynamic graphs, covering as many real-world scenarios as possible. Under the existence of such a dynamic interaction graph, we build predictions of its observations by considering its neighbor’s actions for every agent. We call this primitive interactive predictions. Then an agent’s PSR will be formulated as a linear combination of those primitive interactive predictions. Due to the decentralized formulation, our model allows implementing the learning of the MAPSR on a single agent level while maintaining other agents’ interactions through the encoding of the interactive graph topology.
As set as our second goal, we analytically upper bound theL2-norm difference between the underlying true and learned predictive state representations, providing theoretical guarantees for the performance of our learning algorithm under the MAPSR model in all three graph settings.
Our paper’s third goal is to develop an algorithm that integrates the MAPSR model to the MARL. We develop an online algorithm that simultaneously learns the MAPSR, and the agents’ policies end-toend with a fully differentiable neural network structure. Our algorithm provides a new model-based MARL, and it has the flexibility to be extended to most of the offline MARL algorithms. We focus on the continuous action space, which is more practical and common to real-world applications such as robotic control.
Finally, we test our algorithm through systematic numerical experiments on MAMujoCo robotic learning experiments (de Witt et al., 2020) and multi-agent particle learning environments (Ryu et al., 2020), and compare our proposed method against two baselines as detailed in Section 6. The results demonstrate the efficiency of our method.

heading: 2 RELATED WORKS
text: Partially Observable Environment. Real-world agents often experience situations that the observed signals are aliased and do not fully determine their state in the system. This is particularly true for multiple agents environments where agents have partial observability due to limited communication (Oliehoek & Amato, 2016). In accommodation to the partially observable environment, POMDP has been adopted by (Kaelbling et al., 1998), and the algorithms (Kaelbling et al., 1998; Cassandra, 1998; Thrun, 1999; Pineau et al., 2003; Poupart & Vlassis, 2008; Platt Jr et al., 2010) for determining an optimal policy have shifted to using the probability distributions (belief state) over the state space instead of exact state space. In general, they have high complexity or suffer from local optima. Moreover, the most common POMDP policy learning assumes the agent has access to a priori knowledge of the system. The access to such prior knowledge has a premise that the agent has
considerable domain knowledge (Kaelbling et al., 1998). However, it is expected that the real-world agents learn the system model and thus a planning policy further without knowledge of the domain.
Overview of PSR. (Littman et al., 2001; James & Singh, 2004) introduced the PSR over an expressive and robust framework for modeling dynamical systems and defined PSR as a representation of state by using a vector of predictions of fully observable quantities (tests) conditioned on past events (histories). A predictive model is constructed directly from execution traces in the PSR framework, utilizing minimal prior information about the domain. The PSR paradigm subsumes POMDP as a special case (Littman et al., 2001). PSR is considered much more compact than POMDP (Aberdeen et al., 2007). The spectral learning method has been proved to show success for learning the PSR (Boots et al., 2011; Jiang et al., 2018). There are other classes of dynamical system learning algorithms that are based on likelihood-based optimization or sampling approaches (Frigola et al., 2013), but they are prone to poor local optima. The spectral learning represents the estimated state by sufficient statistics of future observations and estimates the model parameters by method of moments. However, this line of algorithms is hard to incorporate prior information (Hefny et al., 2015). Thus, (Hefny et al., 2015; 2018a) introduce the supervised learning method to learn PSR and proves its convergence. Although many works study PSR in discrete action space (Hsu et al., 2012; Siddiqi et al., 2010; Boots et al., 2011), (Boots et al., 2013) proposes Hilbert space embedding (HSE)-PSR to deal with continuous actions. (Hefny et al., 2018a) uses an approximation of HSE-PSR by Random Fourier transform (RFF) and built a more principled generalization of PSR to deal with high dimensions. However, all of these studies aim for the single agent scenario.
PSR and RL. The predictive states estimated by the PSR are considered as states in a fully observable Markov Decision Process so that the value function is learned on these states. This line of work has been done in the single-agent environment (Boots & Gordon, 2010; Boots et al., 2011; Hamilton et al., 2014; Venkatraman et al., 2017; Hefny et al., 2018b). Especially, (Hefny et al., 2018b) proposes the recurrent predictive state in the RNN network. Moreover, the learning PSR and policy functions are connected with the end-to-end training.
MAPSR and MARL. MAPSR model is formulated by (Chen et al., 2020) and the model parameters are learnt by tensor decomposition. This formulation is aligned with the spectral algorithm, so it also has the same issue as the single-agent case. MARL algorithms have been developed to learn the Markov Game environment (Lyu & Amato, 2020; Son et al., 2019; Zhang et al., 2019; Rashid et al., 2018; Foerster et al., 2018; Lowe et al., 2017). The partially observable states are pervasive in this environment, bringing non-stationarity to the learning algorithm (Bowling & Veloso, 2002) so that the learning needs to consider the non-stationarity. Recent works propose more sophisticated deep MARL algorithms for multi-agent problems under the paradigm of centralized training with decentralized execution (Zhou et al., 2020; Sunehag et al., 2017; Lowe et al., 2017; Foerster et al., 2018; Rashid et al., 2018). Other than their methods, our approach also accounts for the importance of macroscopic measures of underline systems. Also, providing a predictive state as a belief state rather than the partially observable observations can further help the agent learn. Like the single-agent case, we hope our method can further improve these algorithms.

heading: 3 MODEL DESCRIPTION OF MAPSR

heading: 3.1 REVIEW OF SINGLE AGENT PSR
text: Predictive State Representation (PSR). A prediction of a state is defined as the conditional probability of seeing a test’s observations in sequence given that actions of the test are taken in sequence from a history (Littman et al., 2001). Givn a finite observation space O (o ∈ O) and action space A (a ∈ A). A test of length k at time t, is defined as a sequence of action-observation pairs that starts at time t and ends at time t+ k − 1, {(ol, al)}t+k−1l=t = {ot, at, ot+1, at+1, · · · , ol, al, · · · , ot+k−1, at+k−1}. A history, at time t is a sequence of action-observation pairs that start from the beginning of time and ends at time t− 1, {(ol, al)}t−1l=1 = {o1, a1, · · · , ol, al, · · · , ot−1, at−1}. Hilbert Space Embedding Predictive State Controlled Model. In this work, we are interested in extending PSR to decision makings of controlled systems with continuous actions. So we use the model introduced by (Hefny et al., 2018a). In this model, predictive state Qt satisfies Qtψ a t = E[ψot |ψat ;ψht ] and extended predictive state Pt satisfies Ptξat = E[ξot |ξat ;ψht ] (i.e., Qt and Pt are conditional linear expectation operators which maps to the conditional expectation of future observations), where ψot := φO({ol}t+k−1l=t ) and ψat := φA({al} t+k−1 l=t ) are feature maps by kernels kO, kA over future observation and action features. The extended predictive state compared
to predictive state adds one more pair of {(at+k, ot+k)} to the prediction. The ξat and ξot are the corresponding extended feature maps, which satisfy ξot = ψ o t ⊗ φot and ξat = ψat ⊗ φat . Here φot := φo(ot) and φ a t := φa(at) are the shorthands for one time feature map by ko and ka. We use ⊗ to denote the transposed Khatri–Rao product for two matrices with the same number of rows, and each row of the resultant matrix is the vectorzied outer product of the corresponding row vectors in the two matrices. Also, we use the ψht := ψ
h({(ol, al)}t−1l=1) to define a set of features extracted from previous observations and actions (typically from a fixed length window ending at t− 1).

heading: 3.2 MAPSR FORMULATION
text: For n-agent systems, we denote joint observations as {ot,i}ni=1 and joint actions as {at,i}ni=1 at any time t. The joint test is then a collection of tests for all agents {(ol,1, al,1), · · · , (ol,n, al,n)}t+k−1l=t and the joint history is a collection of histories for all agents {(ol,1, al,1), · · · , (ol,n, al,n)}t−1l=1 . A prediction of multi-agent state is thus the conditional probability of a joint test given a joint history.
Our Formulation of MAPSR. We bring the HSE-PSR model to the multi-agent settings. For any agent i in an n-agent system, given kernels {kOi}ni=1,{kAi}ni=1,{koi}ni=1,{kai}ni=1, we use ψot,i, ψat,i, ξot,i, ξ a t,i, φ o t,i, φ a t,i and ψ h t,i to denote its corresponding feature maps, same as the previous definitions for single agent. The superscripts o, a, h are local observations, actions, and histories for agent i. Moreover, we use the Qt,i and Pt,i to denote its predictive state and extended predictive state.

heading: 3.3 THE ESTIMATION FOR MAPSR MODEL COMPONENTS
text: As stated by (Littman et al., 2001), a complete PSR model can build a recursive rule to update itself. In other words, if we are given Qt, we can get Qt+1 via the Bayes’ rules. We use a supervised learning method proposed by (Hefny et al., 2018a) with three types of models to update the predictive states: For any agent i, Pt,i = Wi(Qt,i), Qt+1,i = Fi(Pt,i, ot,i, at,i), and ot,i = Zi(Qt,i, at,i). Typically, Wi, Zi are learnable linear maps, and Fi is non-linear but known in advance. Since the model components, Qt,i and Pt,i are operators, so they have to be learned as well. The Wi, Zi can be learned from regression if we know the estimation of Qt,i and Pt,i, see Appendix A.2.2 for details of estimation of Wi, Zi. We focus on the details of estimation of Qt,i, Pt,i to section 4.1.1.

heading: 4 DYNAMIC INTERACTION GRAPH FOR MAPSR MODEL
text: Many works have considered graph representation of the multi-agent network (Liu et al., 2020b; Ryu et al., 2020). In general, the relationship between agents is characterized as an undirected graph. We introduce a dynamic interaction graph to represent the MAPSR by considering the interaction between agents.
Definition 1. Let G = (V,E), including the set V of nodes and set E of the edges. Each node represents the agent entry, and the edge represents the relationship between the two adjacent agents.
Here we suppose the graph structure is given, that means the number of nodes (n), the number of edges (m), the edge weights, and the maximum number of degrees (k) are available to us. We think that this kind of presupposition is very reasonable because in the real world, for example, there are multiple robots; we can quickly get the geographic position of the robots through sensors, then calculate the structure of the graph formed by them.

heading: 4.1 STATIC COMPLETE GRAPH
text: We starts with a static complete graph Gc, where the relationship between nodes are invariant to time change. A complete graph has m = n(n − 1)/2 edges, where m = |E|. For each agent, we then represent its PSR by considering other agents’ interactions by
Qt,i := g ( {Qt,i,j}nj=1 ) = ∑ j Qt,i,j , (1)
same for Pt,i. To consider the interactive behavior between agents, we introduce two additional notations Qt,i,j and Pt,i,j . Let {(i, j)}ni,j=1 represents a pair of agents on a n multi-agent system. The Qt,i,j is a primitive predictive state of i’s observation (ψot,i) by intervening agent j’s action (ψ a t,j) and observing agent i’s observation history and agent j’s action history (ψht,i,j). And Pt,i,j is the extended counterpart. If i == j, then it becomes an exact single agent scenario. Similarly, we use the same approach to represent Pt,i. Each agent’s PSR Qt,i and extended PSR Pt,i are modeled by fully considering all other available agents.
Based on equation (1), in practice, estimating Qt,i and Pt,i denoted as Q̂t,i and P̂t,i requires us to get Q̂t,i,j and P̂t,i,j at first.
4.1.1 ESTIMATION OF Qt,i,j , Pt,i,j
To estimate the Pt and Qt for single agent, (Hefny et al., 2018a) use the supervised learning method. They show that Q̂t = Mψot |ψat ;ψht = Cψotψat |ψht (Cψat ψat |φht + λI)
−1, where MA|B;c is a linear operator that satisfies E[A|B = b;C = c] = MA|B;cb, and CXY := E[φ(X) ⊗ φ(Y )] is the uncentered covariance operator, and CXY |z is covariance of X and Y given Z = z. They estimate Cψotψat |ψht , Cψat ψat |φht by sampling data, then get the estimation Q̂t. The same procedure is used to estimate P̂t by replacing the features ψot and ψ a t with their extend counterparts ξ o t and ξ a t .
Similar to the single-agent case above, the representation of Qt,i,j can be achieved as
Q̂t,i,j = Mψot,i|ψat,j ;ψht,i,j = Cψot,iψat,j |ψht,i,j ( Cψat,jψat,j |φht,i,j + λI )−1 (2)
To estimate Cψot,iψat,j |ψht,i,j and Cψat,jψat,j |ψht,i,j , we learn two linear maps Ti,j and Ui,j such that Ti,j ( ψht,i,j ) ≈ Cψot,iψat,j |ψht,i,j and Ui,j ( ψht,i,j ) ≈ Cψat,jψat,j |ψht,i,j . The training examples for Ti,j and Ui,j consist of pairs (ψht,i,j , ψ o t,i ⊗ ψat,j) and (ψht,i,j , ψat,j ⊗ ψat,j). The learning of P̂t,i,j can be done in a similar way. More details can be found in Appendix A.2.1
After calculation of Q̂t,i,j and P̂t,i,j , under the static complete graph setting where the interaction is considered, defined in equation (1), we can get the estimate of Q̂t,i and P̂t,i.
4.1.2 THEORETICAL GUARANTEE OF ESTIMATION OF Qi , Pi UNDER THE STATIC COMPLETE GRAPH
Theoretically, we show that the difference between Qi and its estimator Q̂i is bounded with high probability. Theorem 1. Let πΘ be a data collection policy and H is the range of πΘ on joint histories. If Equation 1 and 2 used, then for all h ∈ H and any ∈ (0, 1), such thatN > t2Aj log(2dAj / )
v ( Cψa
j ) where N is the number of time points we collect sample, then
∥∥∥Q̂i −Qi∥∥∥ is bounded as below with probability at least 1− 3 ,
∥∥∥Q̂i −Qi∥∥∥ ≤ n∆ (3) where
∆ = √√√√√ u ( C ψo i |ψh i,j
) v ( C ψa j |ψh i,j )3 · ‖∆1‖ 2+2u ( C ψa j |ψh i,j ) ‖∆1‖+λ v ( C ψa j |ψh i,j ) (1−γ)+λ + ∥∥∥∥Cψo i ψa j |ψh i,j ∥∥∥∥‖∆1‖+‖∆2‖∥∥∥∥Cψa j |ψh i,j ∥∥∥∥+‖∆2‖‖∆1‖ v ( C ψa j |ψh i,j )2 (1−γ)2+λ .
Here ∆1, ∆2 are two other relevant bounds, we provide them in Appendix E.1. u(·), v(·) denote the largest, smallest eigenvalue of a matrix. And γ = t2Aj log(2dAj / )
v ( Cψa
j
) N < 1 is a constant that depends
on the magnitude of the norm of ψaj (we assume ∥∥ψaj ∥∥ ≤ tAj ), the dimension of ψaj (dAj ), the
(uncentered) covariances (Cψaj := E [ ψajψ aT j ] ) and the sample size (N ).
Theorem 1 says we need at least N samples for the bound in equation (3) to be valid. We give the proofs in Appendix E. It is not hard to obtain a bound for Pi using the same approach. We omit that.

heading: 4.2 EXTENSION TO STATIC NON-COMPLETE GRAPH
text: In large-scale multi-agent systems, the number of agents is large, and not all agents need to interact with each other. A static non-completed graph can perfectly represent such a situation. For example, in a given static non-complete graph Gs, we know its maximum number of degree k, and we use the binary n× n matrix with each entry as Ii,j to indicate the interaction between two agents. Then the MAPSR for each agent will be
Qt,i := g ( {Qt,i,j}nj=1 ) = ∑ j Ii,jQt,i,j . (4)
Lemma 1. Under the same environment depicted in Theorem 1 and given a Gs with k maximum number of the degree to represent agents, then the bound in Theorem 1 can be rewritten as∥∥∥Q̂i −Qi∥∥∥ ≤ k∆. The conclusion of Lemma 1 is evident since we replace the n with the k neighbors, the total error bound is also decreased approximately as kn .

heading: 4.3 DYNAMIC GRAPH
text: Real-world multi-agents can also formulate a time-dependent dynamic graph Gd, rather than a static graph. A dynamic graph has its structure dynamically changing with time. In other words, the edges can be inserted or deleted across time. The dynamic graph brings more challenges to the representation as the interaction relationship among agents changes constantly.
Braha & Bar-Yam (2009), Ma et al. (2017) and Zhao et al. (2010) consider a dynamic graph as a set of ordered static graphs. For each time point, we are given a static graph such that Gd = {Gd1, Gd2, · · · , Gdt}, and a time-dependent given binary matrix with It,i,j indicating the interaction between two agents at each time. Then we have the agent-wise PSR as
Qt,i := g ( {Qt,i,j}nj=1 ) = ∑ j It,i,jQt,i,j . (5)
Compared to (4), the coefficient It,i,j is time-dependent, which brings a challenge to our theoretical bound. We consider a dynamic graph experiences a trajectory path, assuming every node has a chance at least p to interact by connection with another node at any time point. For example, for a node i, if we take the union set of the nodes interacted with i over the path, then the union set could form a static complete graph; in other words, if i connects j, then we can obtain a valid sample to estimate Qi,j as the complete static graph does, if not, then we skip to the next time point. For the complete static graph, we need the trajectory to run at least N time points to collect enough data to estimate our conditional operator accurately Qi,j . Furthermore, The total number of time points needed by i until the N th interaction with j follows a negative binomial distribution NB(N, p). On average, we need Np time points before we see i, j completely connecting N times. Now we consider node i could interact with every node in a set of nodes (J : |J | = n− 1) for N number of time points. We assume the interaction between two nodes (i, j ∈ J) does not affect their interaction with other nodes. Thus we have a set of independent negative binomial random variables {Jl}n−1l=1 ∼ NB(N, p) to characterize the interaction of i with j ∈ J . So we are interested in the expectation of the maximum of J1, . . . , Jn−1, a statistics that tells us the expected maximum number of time points of collection of measurements needed for the node i to be able to connect with every node j ∈ J for at least N number of time points. We denote it as J{1,...,n−1} and we have
E{J{1,...,n−1}} = E{max(J1, · · · , Jn−1)} = ∑ N≥0 ( qN +NpqN−1 + · · ·+ ( N N − 1 ) pN−1q )n−1 . (6)
Lemma 2. Under the same environment depicted in 1 and given a dynamic graph Gd with every node has a chance at least p to interact with another node in a one-time point, let N be the number of time points we collect data of measurements in order to get the bound in Theorem 1, if Equation 5 and 2 used, then we need at least
N ′ = ( N − 1
2
) +K (q, n,N)− γ
log1/q(1/q) + F [K (q, n,N)] +O(1), (7)
total number of time points, where q = 1 − p, K (q, n,N) := log1/q(n − 1) + (N − 1) log1/q [ log1/q(n− 1) ] + (N − 1) log1/q p − log1/q(N − 1)!, F is a periodic C∞-function of
period 1 and mean value 0 whose Fourier-coefficients are given by F̂ (k) = − 1 log( 1q ) Γ(− 2kπi log( 1q ) ) for k ∈ Z \ {0}. Then ∥∥∥Q̂i −Qi∥∥∥ achieves the same bound as in Theorem 1 with probability at least

heading: 1− 3 . In other words,
text: ∥∥∥Q̂i −Qi∥∥∥ ≤ n∆. (7) is an asymptotic expansion of the right-hand of (6).
We give the proof in Appendix E.3.
Lemma 2 gives the worst-case bound for our estimation under the dynamic graph. The result tells us that if we need 1 more sample of measurement for our algorithm to converge on the complete static graph, we need roughly log1/q log1/q(n− 1) more samples on the dynamic graph. So as long as we allow enough learning time, the algorithm can converge with high probability.
Complexity with Increasing Agents: For the complete static graph, we need to evaluate Q̂i,j and P̂i,j every time point, which requiresO(n2) operations and space. Overall, with an increased number of agents, our MAPSR has a polynomial O(n2) scaled complexity, which is feasible for learning in a large number of agents environment and is more efficient compared to the centralized MAPSR (Chen et al., 2020) theoretically, which is combinatorially sample complex, the analysis is shown in Appendix A.4. For a non-complete static graph with k maximum number of degrees for k n, which is more common in the real world, because a very far-away robot will not likely affect the targeted robot, the operation will be significantly decreased to O(k2).
After we get the Q̂t,i and P̂t,i, we build the update rule by usingWi, Fi, Zi, as described in section 3.3.

heading: 5 DECISION-MAKING FRAMEWORK WITH MAPSR MODEL
text: Previous work on single-agent proposed an end-to-end training algorithm (Hefny et al., 2018a) for PSR model and policy learning. Here we design an algorithm for multi-agent settings and incorporate our MARPSR model containing the interactive graph component.
We propose an online learning algorithm to learn the MAPSR and agent policies simultaneously. Our algorithm is shown in Appendix B.2. We call the algorithm MAPSRL. The algorithm can be divided into two parts. Firstly, we have a state tracking component, which lets the agents explore the environment and learn an initial predictive state and its parameters Wi and Zi. The agents first enter the initialization phase (Line 1-4), where we initialize MAPSR by executing an exploration policy to collect data. The agents explore the environment and learn the MAPSR parameters under the given interactive graph. Secondly, we also have a policy component, in which, by starting from the initialized MAPSR and a random policy, the agents experience a T length trajectory under the current policy function that maps predictive states {Qt,i}ni=1 to actions {at,i}ni=1 (line 5-18). The agents then update to the next predictive state as introduced in subsection 3.3. Agents also save trajectories (actions, observations, predictive states, next predictive states, and rewards) over the path. Finally, a MARL algorithm conducts the policy learning component. We use a diagram in Appendix B.3 to illustrate this process. We update the MARL and MAPSR model parameter Θ = {ΘMAPSR,ΘMARL} (line 19) by minimizing the following objective function.
L(Θ) = α1l1(ΘMARL) + α2l2 (ΘMAPSR,MARL)
= α1l1(ΘMARL) + α2 1
n T∑ t=0 n∑ i=1 Ep(τi|Θ) [∥∥Zi (Fi(Wi(g (Qt−1,i,j)))⊗ ψat,i)− ψot,i∥∥22] .(8) Here the l1(ΘMARL) is the loss for MARL, and l2(ΘMAPSR,MARL) is the MSE between prediction and actual observation. And p(τi|Θ) is the distribution over trajectories induced by the policy and MAPSR. ΘMAPSR = {Wi, Zi}ni=1 denotes MAPSR’s parameters, and ΘMARL denotes the MARL part. Usually, in a continuous environment, each agent’s policy will be parameterized by its actor network that outputs the mean and diagonal covariance of a Gaussian distribution over the continuous action space. We could also give a parameter sharing actor network that maps individual PSR to parameters of a Gaussian distribution over the individual action space if agents are homogeneous.
Now we devise our MAPSR based MARL algorithm in detail. Our algorithm follows the actorcritic framework. We first develop an algorithm under partially observable environments, where each agent has its own critic and actor independently, the gradient of the policy is written as ∇θiJ(θi) = E[∇θi log πθi(ai|Qi)(ri + γV πi(Qt+1,i) − V πi(Qi))], and the independent critic is updated by minimizing the loss L = E[(V (Qi)− yi)2]. Here the predictive states Qi estimated by the MAPSR are considered as states to fit the value and policy functions. We call it MAPSRL-1. As we know, independent actor-critic (IAC) (Foerster et al., 2018) would not work well since the environment is not stationary under multi-agent setting, the MADDPG by (Lowe et al., 2017) solves a non-stationary environment by considering other agents actions; however, every agent needs a separate critic that has the global information. Instead, we use a centralized critic like as in COMA (Foerster et al., 2018) and LICA (Zhou et al., 2020) to save space. Our centralized critic minimizes the loss L = E[(V π(Q1, · · · , Qn) − y)2]. Unlike MADDPG using the deterministic policy, the gradient is taken on the policy parameter space, which is reflected in LICA (Zhou et al., 2020). So our policy gradient is written as ∇θiJ(θi) = E[∇θi log πθi(ai|Qi)∇θiV π(Q πθ1 1 , · · · , Q πθn n )], Also, to solve the non-stationary environments when each agent is learning, it uses the joint predictive states as input. The benefits are 1. PSR (Q) is a function of action (a), bringing agents’ action information into the policy gradient, 2. Also, our interactive MAPSR brings the effect of other agents’ actions
to the policy gradient, so it does not learn a deterministic policy like MADDPG rather a policy that affects the predictions of future observations considering other agents’ actions 3. Without the need to input other agents’ actions, which saves the input space since the predictive state is a conditional operator that considers other agents’ actions. We call it MAPSRL-2. We put our integration details in Appendix B.4.
6 EXPERIMENTS (a) 2-Agents Swimmer
(b) 3-Agents Hopper
Environments. We evaluate the performance of our MAPSRL on a collection of MARL tasks under some OpenAI Gym MAMujoco environments (de Witt et al., 2020), such as multi-agent swimmer, hopper, and ant. Each robotic agent is represented as a body graph, where vertices (joints) are connected by adjacent edges (body segments) as shown in Appendix Figure 5. Each agent controls its joints based on the local information observed. All tasks are learned under the partially observable environments by manually hiding some observations for each agent. The goals of the multi-agent systems are aligned with their corresponding single-agent ones. However, different from the single-agent system, the agents in the multi-agent system need to collaborate to reach their goals. For the interactive graph, for simplicity, we considered a complete static graph, where we assume every agent is connected with all other agents in the graph. We defer the details of the robotic agents to Appendix C.1 and the details of the experimental setup to Appendix C.2.
Baselines and Evaluation. We run 50 iterations for each experiment and collect M = 100 trajectories in every iteration with a maximum of 1000 steps in every trajectory. After each iteration, we compute the average return R = 1 M ∑n i=1 ∑M b=1 ∑Tb t=1 r t i,b on a batch of M tra-
jectories, where, Tb is the length of the bth trajectory. We repeat this process using ten different random seeds and report the average and a standard deviation. To verify the effectiveness of interactive graph, we introduce a baseline called Independent PSR learning (IPSR); in this model, we do not consider the graph, so we formulate n independent single PSR without considering their interactions, which means there are no Qi,j any more but only Qi. The architecture of MAPSRL-1 and MAPSRL-2 remains the same. To verify the advantage of PSR, we introduce another baseline (MARL) where we take out the MAPSR entirely, so it matches with the MARL run on a partially observable environment. For MAPSRL-1, its MARL baseline is IAC (Foerster et al., 2018), and for MAPSRL-2, is MADDPG (Lowe et al., 2017).
Results And Discussion Figure 1 illustrates the empirical average return vs. the number of interactions with the environment measured in time steps. Our MAPSRL methods consistently outperform IPSR and get the highest rewards under partially observable environments, which justifies the representation power of the interactive graph to assist agents in learning in the non-stationary environment with limited observation when the MARL algorithm has a defect (IAC). Moreover, it can also boost the performance of the existing good MARL algorithm (MADDPG). To further verify the effect of learning the PSR part, we also plotted the predicted trajectory to verify the MAPSR’s performance for predicting the observations in Figure 2. We plotted the predicted observations vs. actual observations in iterations 1, and 40, respectively, for MAPSRL-2. We plotted a row × columns figure, with each row representing the observation feature and each column representing each agent. As we can see from the figures, the first iteration of the learning does not predict the actual observation very well; it has some mismatches. However, as learning progresses, the predictions
get increasingly more accurate. Note that the actual trajectory is changing according to the current policy, and the current policy is optimized based on the further accurate learning of MAPSR.
To enrich our algorithm environment, we also test our algorithm (MAPSRL-2) in multi-agent particle environments, using the benchmark by (Lowe et al., 2017), please check Appendix D.2 for experiment details. We test our algorithm for large n cases. In this environment, the agents can have cooperative goals such that all agents must maximize a shared return and conflicting competitive goals. We set up the environments where agents can only perform physical actions but not communication; however, to achieve the goals, agents need explicit communication about others’ locations to achieve the best reward. These partially observable environments give us the motivation to test our method. We report the rewards in Figure 3. We see that MAPSRL outperforms IPSR and baseline, in terms of the convergence speed and final attained rewards, with a different number of agents. The predictive states convey the information that can help communication between agents in limited communication and observation environments.

heading: 7 CONCLUSION
text: We propose a MAPSR model, extending ideas from single-agent predictive state representations to a multi-agent scenario, during the process, we intro-
duce the dynamic interactive graph to model agents’ interactions. Furthermore, we provide the theoretical guarantees of the MAPSR model. Finally, a learning algorithm that supports gradientbased deep MARL methods is developed. Our method provides a model-based MARL framework under a partially observable environment. The experiments proved that our model assumption is valid by observing the highest return while reducing the observations’ prediction error over trajectories.

heading: ETHICS STATEMENT
text: The authors of this paper claimed that we have read the Code of Ethics, adhered to it, and explicitly acknowledged this during the submission process. We claimed that we did not perform any of the following studies that involve human subjects, practices to data set releases, potentially harmful insights, methodologies and applications, potential conflicts of interest and sponsorship, discrimination/bias/fairness concerns, privacy and security issues, legal compliance, and research integrity issues.

heading: REPRODUCIBILITY STATEMENT
text: We have an appendix section after the main paper to include the model supplements, algorithm details, method diagram, experiment settings, and proofs for our main theorems. We provide a link which directs to our codebase in the supplementary material.

heading: APPENDICES TO REINFORCEMENT LEARNING UNDER A MULTI-AGENT PREDICTIVE STATE REPRESENTATION MODEL: METHOD AND THEORY

heading: A METHOD DETAILS

heading: A.1 MATHEMATICAL MODEL OF PSR
text: Here, we give the mathematical model for PSR. As in section 3.1, we let τt := {(ol, al)}t+k−1l=t as a test (a sequence of action-observation pairs) at time t, then we have a subset τot := {ol}t+k−1l=t for observations, and τat := {al}t+k−1l=t for actions. And we let ht := {(ol, al)} t−1 l=1 as the history at time t. A test is executed at time t if we intervene to select the sequence of actions specified by the test. It is said to succeed at time t if it is executed and the sequence of observations in the test matches the observations by the system (Boots et al., 2013).
The prediction of a length-k test τt is defined as the probability distribution that the test succeed at time t, given history ht: Pr(τot |τat , ht) = Pr(τot ,τ a t |ht)
Pr(τat |ht)
Let T = {τz}, z ∈ N is a finite set of core tests. And we use To = {τo1 , · · · , τo|T|} and T a = {τa1 , · · · , τa|T|} to indicate its observation and action parts.

heading: A.1.1 LINEAR PSR
text: A linear PSR is the prediction vector Qt := [ Pr(τo1 |τa1 , ht), . . .Pr(τo|T||τ a |T|, ht) ]> (9)
that contains the probabilities of success of the tests in T, if and only if for any test τ Pr(τo|τa, ht) = fτ (Qt), (10) Where fτ : [0, 1]|T| 7→ [0, 1] is a linear projection function. For simplicity of notation, we use the same notations for QtPt,W, F , which only follows the meaning defined in this section.
It means that knowing the probabilities for the tests in T is sufficient for computing the probabilities for all other tests in the system. The prediction vector is a sufficient statistic for the system at time t, so we call it as state for the PSR at time t. Therefore, PSR can represent state by using a vector of predictions of fully observable quantities (tests) conditioned on past events (histories). The prediction function fτ is linear and has one-to-one relationship to a test τ such that fτ (Qt) = f>τ Qt ∀t. The linear PSR can still represent systems with nonlinear dynamics.
To maintain predictions in T, we need to update the state Qt. To do that, we predict the success of any core test τz prepended by a new action a and observation o at time t+ k− 1, which we call oaτz . Based on Bayes rule, we have
Pr(τoz,t+1|τaz,t+1, ht+1 = (ht, o, a))) = Pr(oτoz,t+1|aτaz,t+1, ht) Pr(o|a, ht) = f>aoτzQt f>aoQt (11)
faoτz , fao ∈ R|T| are linear operators such that (∀τz ∈ T,∀a ∈ A,∀o ∈ O). Then let Faoτ be the matrix with its columns as faoτz for all τz ∈ T. Then the updated state in PSR is obtained by
Qt+1 = F>aoτQt f>aoQt
(12)
Given the initial prediction vectorQ1, the PSR can update with equation 12. This recursive application of Bayes rule to a belief state is called the Bayes filter. Now we have seen the extended prediction vector, so we define extended predictive state as:
Pt := [ Pr(aτo1 |aτa1 , ht), . . .Pr(aτo|T||aτ a |T|, ht) ]> (13)
Clearly, we can see there exists a linear extension operator W such that Pt = W (Qt). And given the new observation, we can have a Bayes filter F such that Qt+1 = F (Pt, o, a). We also introduced these operators in section 3.3.

heading: A.1.2 REPRESENTATION OF STATE AS CONDITIONAL EXPECTATION OF SUFFICIENT STATISTICS
text: Instead of learning the distribution, here we recover the idea of representation of states as a conditional expectation of sufficient statistics and using the supervised learning method to learn (Hefny et al., 2015). Let history hot is a sequence of observation that starts from the beginning of time and ends at time t− 1. We define the belief state bt = Pr(st|hot ), where st is the current state of the world. bt+1 = Pr(st+1|hot+1) is the next time belief state. We call bt as "belief state", which represents the probability distribution over state space, also represents the knowledge and uncertainty about the true state of the system. In dynamic system, the task of getting the updated bt+1 with given bt and new observation ot is called filtering. The task of estimating the Pr(st+1|hot ) with given current bt without incorporating any new observation is called one-step prediction.
Instead maintaining a belief bt over states, spectral algorithms try to recover observable operators that can be used to perform filtering and prediction directly, by maintaining the expected value of a sufficient statistic of future observations.
Let a test τot is a length k sequence of observation defined before. As the recursive Bayes rule holds, we can get the next time prediction vector Pr(τot+1|hot+1) using the new observation ot+k and the current prediction Pr(τot |hot ), we also define an extended prediction vector Pr(τot+1|hot ). These prediction vectors characterize the state of the system and they can be estimated by observable quantities. Given Pr(τot+1|hot ), filtering becomes the task of getting the updated prediction vector Pr(τot+1|hot+1), "conditioning" on ot. One-step prediction becomes getting the Pr(τot+1|hot ), "marginalizing" over ot.
Therefore, the spectral algorithms avoid explicitly estimating the latent state or the initial, transition, or observation distributions. We let Qt = E [νt|hot ], where νt a vector of features that determines the distribution of future observations Pr(τot |hot ). For simplicity of notation, we use the same notations for Qt, Pt,W, F . Let Pt = E [ϕt|hot ], where ϕt is a vector of features that determines the distribution of observations Pr(τot+1|hot ). We call Qt is the transformed predictive state. So Qt+1 = E [ νt+1|hot+1 ] , is the updated predictive state, and Pt is the extended predictive state. Let hot be the feature vector of history h o t .
In Hidden Markov Models and Kalman filters, the extended state Pt is linearly related to the predictive stateQt. Which is Pt = WQt. Estimation of theW can be done using linear regression with samples νt and ϕt, however, due to the overlap between observation windows, the noise terms on νt and ϕt are correlated, which will cause biased estimate. The instrumental regression (Pearl et al., 2000; Stock et al., 2012) is employed. hot is a instrumental variable that do not overlap with sequence {l}l=t+k−1l=t and {l} l=t+k l=t . The correlation corr(h o t , e(ϕt)) = 0 and corr(h o t , e(νt)) = 0, where e is a measure of error. By taking the conditional expectation of Pt = WQt given hot , we have E[Pt|hot ] = E[WQt|hot ] (14)
E[E[ϕt|hot ]|hot ] = WE[E[νt|hot ]|hot ] E[ϕt|hot ] = WE[Qt|hot ]
Based on the above relationship, we first estimate the E[ϕt|hot ] and E[Qt|hot ] by sample hot , νt, and ϕt, we then use the estimates to compute W . So if we start with Q1, we can compute P1 = WQ1, and get the Qt+1 = F (P1, o1), where F is the Bayes filter to update the state.

heading: A.2 LEARNING FOR MAPSR IN DETAILS
text: Here we introduce the details about learning MAPSR, which supplements the section 3.3 and 4.1.1.
A.2.1 LEARNING FOR Q̂t,i,j AND P̂t,i,j
To calculate Q̂t,i,j , we introduce two sets of linear operators Ti,j and Ui,j , such that Ti,j(ψht,i,j) ≈ Cψot,iψat,j |ψht,i,j and Ui,j(ψ h t,ij) ≈ Cψat,jψat,j |ψht,i,j . We estimate them by using two ridge regressions:
arg min Ti,j T∑ t=1 L(Ti,j(ψht,i,j)), ψot,i ⊗ ψat,j) +R(Ti,j) (15)
arg min Ui,j T∑ t=1 L(Ui,j(ψht,i,j), ψat,j ⊗ ψat,j) +R(Ui,j) (16)
L represents the ridge regression loss, and R represents the regularizer. After learning Ti,j and Ui,j , we can get the estimate of Cψot,iψat,j |ψht,i,j and Cψat,jψat,j |ψht,i,j . Then we get the Q̂t,i,j by equation 2.
Similarly, we use extended features to obtain P̂t,i,j .

heading: A.2.2 EXTENSION, FILTERING AND PREDICTION FUNCTIONS
text: Filtering. To obtain Q̂t+1,i from Q̂t,i, P̂t,i, we use filtering. We denote Fi as the filtering function. We describe the filtering process as below. From {oi, ai}ni=1, we obtain the embedding {φot,i, φat,i}ni=1. We then compute the observation covariance
Cot,i,ot,i|ht,i,at,i = Mφot,i⊗φot,i|φat,i;ψht,iφ a t,i. (17)
We then multiply the extended state by inverse observation covariance to change predicting φot,i into conditioning on φot,i. Mψot+1,i|ψat+1,i,φot,i,φat,i;ψht,i = Mψot+1,i⊗φot,i|ψat+1,i,φat,i;ψht,i ×φot,i (Cot,i,ot,i|ht,i,at,i + λI) −1. (18)
× here is to denote n−mode (matrix or vector) product, ×φot,i means multiplying the tensor by a matrix (or vector) in mode φot,i.
We condition on φot,i and φ a t,i to obtain shifted state.
Qt+1,i := Mψot+1,i|ψat+1,i;φot,i,φat,i,ψht,i = Mψot+1,i|ψat+1,i,φot,i,φat,i;ψht,i ×φot,i φ o t,i ×φat,i φ a t,i. (19)
Based on the updating rule, Qt+1,i = Fi(Pt,i, ot,i, at,i), and Pt,i = Wi(Qt,i), we write the filtering equation. Qt+1,i = Fi(Pt,i, ot,i, at,i)
:= Mψot+1,i⊗φot,i|ψat+1,i,φat,i;ψht,i ×φot,i (Mφot,i⊗φot,i|φat,i;ψht,iφ a t,i + λI) −1 ×φot,i φ o t,i ×φat,i φ a t,i,
where Pt,i := Mψot+1,i⊗φot,i|ψat+1,i,φat,i;ψht,i . Fi usually is known because it is obtained through the above calculation using known quantities, however, Wi, Zi must have to be learned by using regressions.
Extension. The Wi can be learned by kernel regression if we know the Q̂t,i and P̂t,i. Previous work (Hefny et al., 2018a) demonstrated the kernel regression model for learning single-agent PSR, here we extend to the MAPSR. We set the model parameter Wi.
We optimized a ridge regression problem for Wi.
arg min Wi T∑ t=1 L(Wi(Q̂t,i), P̂t,i) +R(Wi)
Prediction. We can also get the prediction about the next one time observation ot by the regression function such as: ôt,i := E(ot,i|Qt,i, at,i) = Zi(Qt,i ⊗ ψat,i) We solve the prediction regression function Zi by another ridge regression:
arg min Zi T∑ t=1 L(Zi(Q̂t,i ⊗ ψat,i), ψot,i) +R(Zi)

heading: A.3 TENSOR DECOMPOSITION OF MAPSR
text: Existed Formulation. (Chen et al., 2020) use a n + 1 multi-dimensional tensor called system dynamics tensor D ∈ R|T1|×···×|Tn|×|H| to represent the system dynamics of the MAPSR, with |T1| representing the cardinality of the tests for agent 1, n denoting the number of agents, and the |H| being the cardinality of the joint history. Each element of the tensor is a probability of a joint test given joint histories. Given the system dynamics tensor, finding the latent predictive state can be
transferred into finding a minimal linearly independent set from the system dynamics tensor, and it can be solved by spectral method such as tensor decomposition.
As mentioned earlier, this formulation can not satisfy our needs. In Proposition 1 of Appendix A.4, we also give the sample complexity analysis such that the sample size needed to formulate D scales exponentially with the length of tests and number of agents.
Here we give a summary of tensor decomposition of their method. Given a system dynamic tensor D such that:
D ≈ [λ;D1, · · · , Dn, F ] = R∑ r=1 λrD 1 r ◦ · · · ◦Dnr Fr (20)
Here ◦ is the outer product. D1, · · · , Dn, F are matrices. The factor matrices D1, Dn, F consist of the vectors, i.e., D1 = [D1:1D 1 :2 · · ·D1:R] ∈ R|T1|×R. A colon is used to indicate all elements of a mode, thus, theRth column ofD1 is denoted byD1:R. For any i1 ∈ {1, . . . , |T1|}, in ∈ {1, . . . , |Tn|}, and k ∈ {1, . . . , |H|}, after the decomposition, we get
Di1,··· ,in,k = R∑ r=1 λrD 1 i1r · · ·D n inrFkr, (21)
where {i1, · · · , in, k} are index corresponding to the specific dimension of theD. The last dimension of the tensorD is compressed in a matrix F , and its row vector xk = [xk(1) . . . xk(R)] ∈ R1×R, k ∈ {1, . . . , |H|} is a summary of joint history and can be considered as a compressed version of the system predictive state vector p(Q|hk), the joint history hk ∈ H(k ∈ [1, |H|]) at time step s = |hk|, where Q is the core joint test set. The whole fibers listed in the set Q form a basis of the space spanned by the mode-(n+1) fibers of tensor D. Thus, by constructing the vector m = (λ ∗D1i1: ∗ · · ·D n in:
)T , where ∗ is Hadamard product. D1i1: denotes the i1-th row vector of D
1. Then we could rewrite the previous equation as
Di1,··· ,in,k = R∑ r=1 m(r)xk(r) = xk(λ ∗D1i1: ∗ · · ·D n in:) T = xkm (22)
m(r) is a scalar that m(r) = λrD1i1r · · ·D n inr . And xk is the system state vector and m is the prediction parameter, both of them are obtained by the tensor decomposition.

heading: A.4 ANALYSIS OF SAMPLE COMPLEXITY FOR FORMULATING THE SYSTEM DYNAMIC TENSOR
text: The paper (Chen et al., 2020) does not analyze the sample complexity to construct D, which is a n + 1 multi-dimensional tensor D ∈ R|T1|×···×|Tn|×|H|. We give this analysis. The Di1,··· ,in,k := Pr(t1i1 , · · · , tnin |hk) is an element of that tensor D such that t1i1 is the i1-th test of agent 1, similarly, the tnin is the in-th test of agent n, and hk is the joint history.
Proposition 1. In a n-agents system, assume every agent has the same observation and action space |O|, |A|, for a length-k test, to formulate a complete system dynamics tensorD defined in equation 21. Assume each entry of the tensor needs S samples to give a sufficient estimation using Monte-Carlo roll-out method, then the total sample size is at least (|O||A|)knS, if the agents are homogeneous, in other words, they are permutation invariant such that identity does not matter, then the total sample size is (|O||A|k+n−1 n ) S.
Proof. At one time step, for any agent, it has |O||A| different combinations for the joint test, then for a length of k tests, it follows that (|O||A|)k number of different choices. Then the tensor D would need (|O||A|)kn elements to cover all the possible length k tests. So the total sample size is (|O||A|)knS. If the agents are homogeneous, then the ordering does not matter, for each agent, we have (|O||A|)k number of different choices for length k test, so the total choices for n agents are(|O||A|k+n−1
n
) , then we need (|O||A|k+n−1 n ) S samples.
The sample complexity is exponentially scaled with the number of agents and length of tests.

heading: B ALGORITHM AND INTEGRATING WITH MARL METHOD
text: We first introduce some backgrounds on two multi-agent frameworks: Multi-agent Markov Decision Process (MMDP) and Multi-agent Partially observable Markov Decision Process (MPOMDP) since the algorithms in our paper are developed based on the framework of MPOMDP.

heading: B.1 MULTI-AGENT MDP AND MULTI-AGENT POMDP MODEL
text: B.1.1 MMDP MMDP model is a tuple ( S,N, {Oi}i∈[n]{Ai}i∈[n], T,R ) , where S and N are finite sets of states and agents. Ai is a finite set of actions available to agent i; T : S ×A1 × · · · ×An × S 7→ [0, 1] is a transition function; and R : S 7→ R is the reward function. Each agent i obtains reward as function of the state and agent’s action ri : S ×Ai 7→ R, and receives a private observation from the state by the observation channel oi : S 7→ Oi. The state has distribution d : S 7→ [0, 1]. Each agent aims to maximize its own total expected return ri = ∑T t=0 γ
trti where γ is a discounted factor and T is the time horizon. In a shared reward situation, there is a team reward function r : S ×Ai · · · ×An 7→ R, agents aim to maximize one shared total expected return r = ∑T t=0 γ trt.

heading: B.1.2 MPOMDP
text: A MPOMDP model is a tuple ( S,N, {Oi}i∈[n]{Ai}i∈[n], {Ωi}i∈[n], T,R ) , where (S,Ai, Ti, Oi,Ωi, R) describe a single-agent POMDP. Oi is the set of observations the agent i can make. Ωi : S × Ai × Oi 7→ [0, 1] is the agent’s observation channel function, which specifies probabilities of observations given agent’s actions and resulting states. (S,Ai, Ti, Ri) describes a single agent MDP; and each agent i obtaines reward as function of the state and agetn’s action ri : S ×Ai 7→ R. Each agent aims to maximize its own total expected return ri = ∑T t=0 γ
trti where γ is a discounted factor and T is the time horizon. In POMDP, an agent’s belief about the sate is represented as probability distribution over S. The agent has prior belief b0,i The agent’s current belief, bt,i over S, is continuously revised based on new observations and expected results of performed actions. The belief update takes into account changes in initial belief, bt−1,i, due to action at,i, executed at time t − 1, and the new observation, ot,i. The new time belief state can be obtained from basic probability theory as follows: bi(st) = βΩi(ot,i, st, at−1,i) ∑ st−1∈S bt−1,i(st−1)T (st, at,i, st−1), where β is the normalizing factor.

heading: B.2 ALGORITHM: MAPSRL
text: We give the details of MAPSRL in Algorithm 1.

heading: B.3 MAPSRL IN DIAGRAM
text: We also use a diagram (Fig 4) to depict the algorithm. The algorithm runs k iterations; each iteration first uses the policy to roll out data and uses the regressions to obtain the PSR parameters. Then it executes the policy phase by using the PSR as the input of policy to generate action and using the current PSR parameters to update the predictive state. At the end of the policy phase, it updates both the PSR and policy parameters. The policy parameters and the linear relationship captured by Wi and Zi are parameterized by the neural network, section C.3 has the network architectures. The loss is a composite loss that includes loss from actor-critic and the loss from the predictive state representation. The next iteration will re-learn the PSR parameters using the newly generated data based on the current policy obtained from the previous iteration; then, it does a soft update to update the PSR parameters with ones obtained at the previous iteration.
B.4 INTEGRATING TWO COMMON MARL ALGORITHMS INTO MAPSR MODEL
Here we provide a brief intuitive introduction about how we connect existed MARL algorithms with MAPSR. Please look at Fig 4 for demonstration. ψot,i, ψ a t,j , and ψ h t,i,j are embedding of (ot:t+k−1,i, at:t+k−1,j), and (o1:t−1,i, a1:t−1,j). We also have embedding for extended part, labeled as ξ. Qt,i,j and Pt,i,j are estimated by tensor regression using the embedding vectors. G represents the given
Algorithm 1 MAPSRL 1: Input: Learning rate η, a graph G, a static complete graph Gc or static non-compete graph Gs
or dynamic graph Gt 2: Initialize MARL Policy ΘMARL randomly 3: for k = 1, 2, 3, · · · iterations do 4: State Tracking Phase 5: Sample b = 1, 2, 3, · · ·M batch of initial trajectories: {(obt , abt)}Mb=1 from existed policy obtained from previous iteration k − 1: {πk−1i }ni=1 6: Let πθMARL = {πk−1i }ni=1 if available or the initial policy; 7: Given G, calculate Q̂i and P̂i, and obtain the initial Wi, Zi, Fi: 8: (1).Regression Qt,i,j = Ti,j ◦ Ui,j(ht,i,j) to get the Q̂t,i,j and P̂t,i,j 9: (2).Given the graph G, using equation 1, 4, or 5 to obtain Q̂t,i, P̂t,i
10: (3).Given the Q̂t,i, P̂t,i, to obtain Wi, Zi, Fi 11: 12: Policy Phase 13: Initialize MAPSR parameters ΘMAPSR = {Q1,i,Wi, Zi}ni=1 from state tracking phase and
previous iteration by a soft-update: Wi = βWi + (1− β)W k−1i , Zi = βZi + (1− β)Z k−1 i , Fi = Fi
14: for b = 1, 2, 3, · · ·M batch of trajectories from {πk−1i }ni=1 do 15: Reset episode: ab0,i, o b 0,i 16: for t = 0, 1, 2, · · ·T roll-in in each trajectory do 17: for Each agent i do 18: Get observation obt,i and reward r b t,i 19: Extension P bt,i = Wi(Q b t,i) 20: Filtering Qbt+1,i = Fi(Q b t,i, a b t,i, o b t,i,Wi) 21: Execute abt+1,i ∼ π k−1 i (Q b t+1,i) 22: Predict ôbt,i = Zi(Q b t,i, a b t,i) 23: Collect obt,i, ô b t,i, a b t,i, r b t,i, Q b t,i, Q b t+1,i 24: end for 25: end for 26: end for 27: Update Θ using D = {{{obt,i, ôbt,i, abt,i, rbt,i, Qbt,i}ni=1}Tt=1}Mb=1:
Θ← Update(Θk−1, D, η) as in Equation (8), get W ki ,Zki , and πki , the linear operator Wi, Zi, and nonlinear operator πi are all parameterized by neural network.
28: end for 29: Output: Return Θ = (ΘMAPSR,ΘMARL)
graph, the estimation ofQt,i and Pt,i are based on graphG and equation 1, equation 4 and equation 5 for reference. Please also go to subsection 3.2 and section 4.1.1 for a verbal description. The policy network uses the predictive state as input to return the action. The agent takes the action to get the observation. The filter Fi takes predictive state, action, and observation as inputs to get the next predictive state. We shows a centralized critic and gives the description in paragraph B.4.2. The loss is composed into two parts, LθMAPSR = ‖Zi{Fi [Wi(Qt−1,i), at−1,i, ot−1,i] , at,i} − ot,i‖ 2 2, and Lθcritic , which we gives detailed explanation in equation 29. Please also go to section 5 and Algorithm B.2 for more details about the framework.
B.4.1 MAPSRL-1
It is based on the IAC (Foerster et al., 2018) which directly applies the single-agent policy gradient to have each agent learn independently, with the idea behind independent Q-learning (Tan, 1993), with actor-critic in place of Q-learning.
IAC trains an actor-critic pair for each agent, resulting in actors πi(ai|oi) and critics Vi(oi, ai). ∇θiJ(π) = Eπ [ ∇θi log πi(ai|oi) ( ri + γV πi i (ot+1,i, at+1,i)− V πi i (oi, ai) )] (23)
r i J ( i ) = E[r i log i; i (ai jQi )( r i + V i i (Qt +1 ;i ) V i i (Qi ))] (24) If the agents are homogeneous, we can share the critic and actor network. We have the critic loss as
L = E[(V i (Qi ) yi )2] yi = r i + bV i (Qt +1 ;i ) (25) Here bV is the target value function. B.4.2 MAPSRL-2 The algorithm borrows the idea from MADDPG (Lowe et al., 2017) and LICA (Zhou et al., 2020) using the gradient of the PSR-value function with respect to policy to direct the policy gradient update. MADDPG (Lowe et al., 2017) is an extension of deep deterministic actor-critic policy gradient (DDPG) (Lillicrap et al., 2015) to multi-agent setting such that let each agent’s own critic is augmented with extra information about the actions of other agents, while their individual actor maintains a local state or observation. The gradient of each agent is:
r i J ( i ) = E[r i i (ai joi )r a i V a
i (o1; ; on ; a1; ; an )jai = (oi )] (26) The action-value functionVi is updated as
L i = E[(V a i (o1; ; on ; a1; ; an ) y) 2] y = r i + bVi
a0
(o01; ; o 0 n ; a 0 1; ; a 0 n ) (27)
Wherea0 is the set of target policies with delayed parameters. Unlike MADDPG using the deterministic policy, the gradient is taken on the policy parameter space, which is re ected in LICA Zhou et al. (2020).
r i J ( i ) = E[r i log i (ai jQi )r i V (Q 1 1 ; ; Q n n )] (28)
Our method also uses the centralized critic to save space. Also, in order to solve the non-stationary environments when each agent is learning, it uses the joint predictive states as input without considering other agents’ actions, which saves the input space since the predictive state is a conditional operator that considers other agents’ actions.
The critic loss is defined as below: Lθcritic = E[(V π(Q1, · · · , Qn)− y)2]
y = ∑ i ri + γV̂ π(Q πθ1 t+1,1, · · · , Q πθn t+1,n) (29)
Even though the original MADDPG has the current actions as input to the value function so that the environment is stationary, we take out this operation because the predictive state already considers the effect actions can bring.

heading: C ENVIRONMENT AND EXPERIMENT

heading: C.1 MAMUJOCO ENVIRONMENT SETUP
text: Here, we give the details of setting up our multi-agent environment - MAMuJoCo. As introduced in section 6, the many agents are constructed by separating a existed single agent into parts, and each agent will only control a part of the whole agent (Figure 5).
Partially Observable space: MAMuJoCo is a simulated robotic environment, the partially observable property is achieved by only allowing partial information to the agents. For all environments, only the angles of the agent’s joints are visible to the network; the velocities are hidden.
Action: Each agent’s action space in MAMuJoCo is given by the joint action space overall motors controllable by that agent.
Observation: For each agent i, observations are constructed by inferring which body segments and joints are observable by an agent i. Each agent can always observe all joints within its sub-graph. A configurable parameter k ≥ 0 determines the maximum graph distance to the agent’s subgraph at which joints are observable. For example, k = 0 means agents can only observe their own joints and body parts, while k = 1 means it can observe its adjacent joints, which has 1 graph distance to the agent. The agent observation is then given by a fixed order concatenation of each observable graph element’s representation vector. Depending on the environment and configuration, representation vectors may include attributes such as position, velocity, and external body forces. In addition to joint and body segment-specific observation categories, agents can also be configured to observe the robot’s central torso’s position and velocity attributes.

heading: C.2 EXPERIMENT SETUP
text: We select three experiments from MAMuJoCo and give a detailed description of the experiments’ setup. In all environments, the agent has the goal to maximize the velocity of the first coordinate for the team. We use k to denote the maximum observation distances to the subgraph. We use ∆d to denote the first coordinate position difference between a time difference ∆t. Finally, we use R to denote the reward function. Table 1 has the configuration details for these parameters.

heading: C.3 NEURAL NETWORK ARCHITECTURE
text: We implement all algorithms using deep neural networks as function approximators. We ensure that all policy, value, and action-value functions have the same neural network architecture among all algorithms to the extent each algorithm allows for a fair comparison.
For our experiments with continuous action spaces, a Gaussian distribution with a diagonal covariance matrix is used. The policy network maps from the input feature to a Gaussian distribution vector µ. Moreover, µ = [mean, std], where mean is a vector specifies the action means, and std vector specifies the standard deviation.
In the implementation of all actor-critic method, all the actor-network is parameterized by a multilayer perceptron (MLP) with two hidden layers of size 400 and 300 respectively and ReLU activation, which takes in the individual agent’s predictive state and outputs the mean and covariance of a Gaussian policy. The critic network is also an MLP with two hidden layers with 400 and 300 units, respectively. For MAPSRL-1, the critic network is used to approximate per-agent utilities, which receives each agent’s predictive state as input. For IAC, same as MAPSRL-1, it receives agent local observation and individual action as input.
In MAPSRL-2, there is a shared critic network that approximates all agents utilities, which receives all agents’ predictive states as input. In MADDPG, the critic receives the global state and the joint action of all agents as input. The global state consists of the complete state information from the original OpenAI Gym environment. Each decentralized actor (i.e., policy) network takes in each agent’s observation and outputs the agent’s action vector.

heading: C.4 MODEL PARAMETERS
text: The hyper parameters for the the MAMuJoCo environments are in Table 2.

heading: D SUPPLEMENT EXPERIMENTS

heading: D.1 MAMUJOCO SUPPLEMENT EXPERIMENTAL RESULTS
text: We plotted the predictive observations compared to actual observations in Figure 2 at the beginning of the learning process (iteration 1) and end of the learning (iteration 40). We also show the results of iteration ten and iteration 20 in Figure 6. By comparing with iterations 1 and iterations 40 in Figure 2, we see that the iteration 40 has the smallest difference between predictive observation and true observation, and the difference gets increased as the iteration goes to the earlier stage of the learning process. So the predictive accuracy is improved incrementally with the learning progresses.

heading: D.2 MULTI-AGENT PARTICLE ENVIRONMENT
text: We also test MAPSR into another environment, multi-agent particle environment (Lowe et al., 2017). The agents are displaced into a 2-dimensional coordinate. This environment does not assume that all agents have identical action and observation spaces. We run experiments using a different number of agents on two environments, the predator-prey, and cooperative-push. We use the same configuration
in (Liu et al., 2020a). The network has the same design as MAMuJoCo in section C.3, except we use MLP with two hidden layers with the same 128 units respectively, so we do not repeat the description. In Table 3, we report the hyperparameters for multi-agent particle environments.

heading: E DETAILED PROOFS
text: We give the proof of Theorem 1 in Appendix E.1, and we first introduce the following Lemmas to prepare the proof.
Definition 2. Let X1 . . . Xk be independent random variables of dimensionality dX1 . . . dXk such that ‖Xk‖ < txk . Let {(x1j . . . , xkj )}Nj=1 be the N i.i.d samples from distribution of X1, . . . , Xk,
Lemma 3 ((Tropp, 2015)). Let aj be a finite sequence of independent random, Hermitian matrices with dimension d. Assume that 0 ≤ v(aj) and u(aj) ≤ L for each j. Let S = ∑ j aj , then for any η ∈ [0, 1], it follows that
Pr(v(S) ≤ (1− η)v(E[S])) ≤ d [ e−η (1− η)1−η ]v(E[S])/L ≤ 2de−ηv(E[S])/L (30) Corollary 1. Let X be a random variable, for any ∈ (0, 1) such that N > t
2 x log(2dX/ ) v(CX)
the following holds with probability at least 1−
v(ĈX) > t2x log(2dX/ )
v(CX)N

heading: In other words, if N large enough, then ĈX and CX will be close enough.
text: Proof. Define Sj = 1/NxjxTj . Then it follows that u(Sj) ≤ L = t2x/N and define := 2dXe −σNv(CX)/t2x , which implies that σ = t 2 x log(2dX/ ) v(CX)N . Then it follows from Matrix Chernoff Inequality in Lemma 3 that Pr(v(ĈX) ≤ (1− σ)v(CX)) ≤
Lemma 4 ((Tropp, 2015)). A finite sequence {aj} of independent, random matrices with common dimensions a× b, and assume that E[aj ] = 0 and ‖aj‖ ≤ L for each j, let S = ∑ j aj as a random
matrix. Let V ar(S) be the variance statistics such that V ar(S) = max{ ∥∥E[SST ]∥∥ ,∥∥E[STS]∥∥}, then Pr(‖S‖ > c) ≤ (a+ b)e −c2/2 V ar(S)+Lc/3 (31) Corollary 2. With at least probability 1− that∥∥∥ĈY X − CY X∥∥∥ ≤√2 log(dY + dX)/ V ar N + 2 log((dY + dX)/ )L
3N where L = tytx + ‖CY X‖ ≤ 2tytx and V ar = max{t2y ‖CX‖ , t2x ‖CY ‖}+ ‖CY X‖ 2 ≤ 2t2yt2x
Proof. Let X,Y be two random variables, and let a finite sequence {aj} of independent random matrices to satisfy aj = yjxTj − CY X . So the aj will have dimensions dX × dY . Let random matrix S = ∑ j aj . It follows that E[aj ] = 0 and ‖aj‖ = ∥∥yjxTj − CY X∥∥ ≤ ‖yk‖ ‖xk‖ + ‖CY X‖ ≤
tytx + ‖CY X‖∥∥E[SST ]∥∥ = ∥∥∥∥∥∥ ∑ i,j ( E[yixTi xjyTj ]− CY XCXY )∥∥∥∥∥∥ =
∥∥∥∥∥∥ ∑ i ( E[‖xi‖2 yiyTi ]− CY XCXY ) + ∑ i 6=j ( E[yixTi ]E[xjyTj ]− CY XCXY )∥∥∥∥∥∥ =
∥∥∥∥∥∑ i (E[‖xi‖2 yiyTi ]− CY XCXY ) ∥∥∥∥∥ ≤ N(t2x ‖CY ‖+ ‖CY X‖ 2 )
Similarly, ∥∥E[SST ]∥∥ ≤ N(t2y ‖CX‖ + ‖CY X‖2). By applying lemma 4, we have = Pr(‖S‖ ≥
Nc) ≤ (dX + dY )e
( −Nc2/2 V ar+Lc/3 ) and therefore, it implies that
c ≤ log((dX + dY )/ )L 3N +
√ (log(dX + dY / ))2L2
9N2 +
2 log((dX + dY )/ )V ar
N
≤ 2 log((dX + dY )/ )L 3N +
√ 2 log((dX + dY )/ )V ar
N
Corollary 3. For random variable X with dimensionality dX and ‖X‖ ≤ tx, with probability 1− , it follows that ∥∥∥C−1/2X (ĈX)− CX∥∥∥ ≤ 2tx √ 2 log(2dX/ )
N +
2 log(2dX/ )L
3N
where L = t 2 x√
v(CX) + tx
Proof. The proof is similarly to the the proof of corollary 2, define aj = ∑−1/2 X xjx T j − C 1/2 X ,
S = ∑ aj then it follows that E[aj ] = 0 and ‖aj‖ ≤ t 2 x√
v(CX) + tx∥∥E[STS]∥∥ = ∥∥E[SST ]∥∥ ≤ N(t2x + ‖CX‖2) ≤ 2Nt2x (32)
Applying lemma 4 to get
= Pr(‖S‖ ≥ Nc) ≤ 2dXe −Nc2/2
2t2x+Lc/3 (33) it follows that
c ≤ 2 log(2dX/ )L 3N + 2tx
√ log(2dX/ )
N (34)
Lemma 5. For two random variables X,Y , let ĈY X = CY X + ∆Y X , and ĈX = CX + ∆X where E[∆Y X ] and E[∆X ] are not necessarily zero and ĈX is symmetric positive semidefinite. Define A = CY XC −1 X and Â = ĈY X(ĈX + λ)
−1. Then it follows that:∥∥∥Â−A∥∥∥ ≤√u (CY ) v(CX) √v(CX) ∥∥∥C−1/2X ∥∥∥∆X + λ v(ĈX) + λ + ‖∆Y X‖ v(ĈX) + λ

heading: Proof.
text: Â−A = CY X((CX + ∆X + λI)−1 − C−1X ) + ∆Y X(CX + ∆X + λI) −1 = M1 +M2
It follows that ‖M2‖ ≤
∆Y X
v(ĈX) + λ
For M1, by using facts U−1 − V −1 = U−1(V − U)V −1 and CY X = C1/2Y PC 1/2 X , where P is a correlation matrix with ‖P‖ ≤ 1, M1 = −CY XC−1X (∆X + λI)(CX + ∆X + λI) −1
= −C1/2Y PC −1/2 X (∆X + λI)(CX + ∆X + λI) −1
‖M1‖ ≤ √ u (CY ) ∥∥∥C−1/2X ∆X∥∥∥+ λ ∥∥∥C−1/2X ∥∥∥ v(ĈX) + λ
=
√ u(CY )
v(CX)
√ v(CX) ∥∥∥C−1/2X ∆X∥∥∥+ λ v(ĈX) + λ
Corollary 4. Let {(xk, yk)}Nk=1 be i.i.d samples from two random variables X,Y with dimensions dX and dY and (uncentered) covariances CX and CY . Assume ‖X‖ ≤ tx and ‖Y ‖ ≤ ty. Define A = CY XC −1 X and Â = ĈY X(ĈX + λ) −1. For any ∈ (0, 1) such that N > t 2 x log(2dX/ ) v(CX)
the following holds with probability at least 1− 3 ∥∥∥Â−A∥∥∥ ≤√u(CY )
v(CX)
( √ v(CX)α+ λ
v(CX)(1− γ) + λ
) +
β
v(CX)(1− γ) + λ
where
α = 2tx
√ 2 log(2dX/ )
N +
2 log(2dX/ )
3N ( c2x√ v(CX) + tx)
β = 2tytx
√ log(dY + dX)/
N +
3tytx log((dY + dX)/ )
3N
γ = t2x log(2dX/ )
v(CX)N
Proof. It follows by applying Corollaries 1,2,3 to Lemma 5. By union bound, each condition has probability 1− , so the total events are bounded
Pr(bounds satisfied) := 1−Pr( 3⋃ i=1 Ai) ≥ 1− 3∑ i=1 Pr(Ai) = 1− 3
Lemma 6. For two random variables X,Y , let ĈY X = CY X + ∆Y X , and ĈX = CX + ∆X where E[∆Y X ] and E[∆X ] are not necessarily zero and ĈX is symmetric but not positive semidefinite. Define A = CY XC−1X and Â = ĈY XĈX(Ĉ 2 X + λ)
−1. Then it follows that:∥∥∥Â−A∥∥∥ ≤√ u(CY ) v(CX)3 ‖∆x‖2 + 2u(CX) ‖∆X‖+ λ v(ĈX) + λ + ‖CY X‖ ‖∆X‖+ ‖∆Y X‖ ‖CX‖+ ‖∆Y X‖ ‖∆X‖ v(ĈX)2 + λ

heading: Proof.
text: Â−A = (CY X + ∆yx)(CX + ∆X)((CX + ∆X)2 + λI)−1 − CY XCXC−2X
= CY XCX(((CX + ∆X) 2 + λI)−1 − C−2X ) + (CY X∆X + ∆Y XCX + ∆Y X∆X)((CX + ∆X) 2 + λI)
= M1 +M2
For M1, by using facts U−1 − V −1 = U−1(V − U)V −1 and CY X = C1/2Y PC 1/2 X , where P is a correlation matrix with ‖P‖ ≤ 1, it follows that M1 = −C1/2Y PC −3/2 X (∆ 2 X + CX∆X + ∆XCX + λI)((CX + ∆X) 2 + λI)−1
Therefore,
‖M1‖ ≤
√ u(CY )
v(CX)3 ‖∆x‖2 + 2u(CX) ‖∆X‖+ λ
v(ĈX) + λ
‖M2‖ ≤ ‖CY X‖ ‖∆X‖+ ‖CY X‖ ‖CX‖+ ‖CY X‖ ‖∆X‖
v(ĈX)2 + λ
Corollary 5. Let {(xk, yk)}Nk=1 be i.i.d samples from two random variables X,Y with dimensions dX and dY and (uncentered) covariances CX and CY . The E[∆Y X ] and E[∆X ] is not necessarily zero and Cx is symmetric but not necessarily positive semidefinite. Assume ‖X‖ ≤ tx and ‖Y ‖ ≤ ty . DefineA = CY XC−1X and Â = ĈY XĈX(Ĉ 2 X+λ) −1. For any ∈ (0, 1) such thatN > t 2 x log(2dX/ ) v(CX)
the following holds with probability at least 1− 3 ∥∥∥Â−A∥∥∥ ≤√ u(CY ) v(CX)3 ‖∆x‖2 + 2u(CX) ‖∆X‖+ λ v(CX)(1− γ) + λ + ‖CY X‖ ‖∆X‖+ ‖∆Y X‖ ‖CX‖+ ‖∆Y X‖ ‖∆X‖ v(CX)2(1− γ)2 + λ
where
γ = t2x log(2dX/ )
v(CX)N
Proof. It follows by applying Corollaries 1,2,3 to Lemma 6. Also by union bound, so the total events are bounded 1− 3
Theorem 2. Grabner 1997 (3.3) (Grabner & Prodinger, 1997) Consider there are n independent copies X1, . . . Xn i.i.d negative binomial random variables, with parameters defined as NB(b, p), and our goal is to calculate the expectation of the maximum of these N random variables En = E{max(X1, . . . , Xn)} then we have following asymptotic solution En = log 1
q (n) + (b− 1) log 1 q log 1 q (n) + (b− 1) log 1 q p+ (b− 1)− log 1 q (b− 1)! + 1 2 +
γ
log 1 q (1/q)
+ F (log 1 q (n) + (b− 1) log 1 q log 1 q (n) + (b− 1) log 1 q p− log 1 q (b− 1)!) + o(1) (35) (where F is a periodic C∞ − function of period 1 and mean value 0 whose Fourier-coefficients are given by F̂ (k) = − 1
log( 1q ) Γ(− 2kπi log( 1q ) ) for k ∈ Z \ {0}, and q = 1− p)
We omit the proof, interested readers could go to (Grabner & Prodinger, 1997) for details.

heading: E.1 PROOF OF THEOREM 1
text: We first prove the bound for ∥∥∥Q̂i,j −Qi,j∥∥∥.
Proposition 2. Let πΘ be a data collection policy and H is the range of πΘ on joint histories. If Equation 2 used, then for all h ∈ H and any ∈ (0, 1), ∥∥∥Q̂i,j(ψh)−Qi,j(ψh)∥∥∥ is bounded as below with probability at least 1− 3 .
∥∥∥Q̂i,j −Qi,j∥∥∥ ≤ √√√√√√ u ( Cψoi |ψhi,j ) v ( Cψaj |ψhi,j )3 ‖∆1‖ 2 + 2u ( Cψaj |ψhi,j ) ‖∆1‖+ λ v ( Cψaj |ψhi,j ) (1− γ) + λ
+ ∥∥∥Cψoi ψaj |ψhi,j∥∥∥ ‖∆1‖+ ‖∆2‖ ∥∥∥Cψaj |ψhi,j∥∥∥+ ‖∆2‖ ‖∆1‖ v ( Cψaj |ψhi,j )2 (1− γ)2 + λ
where ∆1 follows the bound 38 and ∆2 follows the bound 36, and γ = t2Aj log(2dAj / ) v ( Cψa
j
) N
Proof. Let Ti,j is the tensor such that Cψoi ψaj |ψhi,j = Ti,j ×h ψ h i,j , and Ui,j is the tensor such that Cψaj |ψhi,j = Ui,j ×h ψ h i,j and for simplicity without loss meaning, we use Cψaj |ψhi,j to denote Cψaj ψaj |ψhi,j . Then we have∥∥∥Ĉψoi ψaj |ψhi,j − Cψoi ψaj |ψhi,j∥∥∥ ≤ ∥∥∥T̂i,j − Ti,j∥∥∥∥∥ψhi,j∥∥∥∥∥Ĉψaj |ψhi,j − Cψaj |ψhi,j∥∥∥ ≤ ∥∥∥Ûi,j − Ui,j∥∥∥∥∥ψhi,j∥∥ We finish the above proof by proofing the ‖Ti,j − Ti,j‖ and ‖Ui,j − Ui,j‖ are bounded by using Corollary 4.
∥∥∥T̂i,j − Ti,j∥∥∥∥∥ψhi,j∥∥ ≤ th √√√√√u ( Cψoi ,ψaj ) v ( Cψhi,j )  √ v ( Cψhi,j ) α+ λ v ( Cψhi,j ) (1− γ) + λ + β v ( Cψhi,j ) (1− γ) + λ ,
(36)
where
α = 2th
√ 2 log(2dh/ )
N +
2 log(2dh/ )
3N
 t2h√ v ( Cψhi,j ) + th  ,
β = 2tOitAj th
√ log(dOidAj + dh)/
N +
4tOitAj th log((dOidAj + dh)/ )
3N ,
γ = t2h log(2dh/ ) v ( Cψhi,j ) N ; (37)
and
∥∥∥Ûi,j − Ui,j∥∥∥ ∥∥ψhi,j∥∥ ≤ th √√√√√ u ( Cψaj ) v ( Cψhi,j )  √ v ( Cψhi,j ) α+ λ v ( Cψhi,j ) (1− γ) + λ + β v ( Cψhi,j ) (1− γ) + λ ,
(38)
where
α = 2th
√ 2 log(2dh/ )
N +
2 log(2dh/ )
3N
 t2h√ v ( Cψhi,j ) + th  ,
β = 2tAj th
√ log(dAj + dh)/
N +
4tAj th log((dAj + dh)/ )
3N ,
γ = t2h log(2dh/ ) v ( Cψhi,j ) N . (39)
Then using the equation 2 and corollary 5 to obtain the bound for Qi,j
∥∥∥Q̂i,j −Qi,j∥∥∥ ≤ √√√√√√ u ( Cψoj |ψhi,j ) v ( Cψaj |ψhi,j )3 ‖∆1‖ 2 + 2u ( Cψaj |ψhi,j ) ‖∆1‖+ λ v ( Cψaj |ψhi,j ) (1− γ) + λ
+ ∥∥∥Cψoi ψaj |ψhi,j∥∥∥ ‖∆1‖+ ‖∆2‖ ∥∥∥Cψaj |ψhi,j∥∥∥+ ‖∆2‖ ‖∆1‖ v ( Cψaj |ψhi,j )2 (1− γ)2 + λ ,
where ∆1 follows the bound 38 and ∆2 follows the bound 36, and γ = t2Aj log(2dAj / ) v ( Cψa
j
) N
Now we start to prove our theorem 1.
Proof. The equation 1 says: Qt,i := g ( {Qt,i,j}nj=1 ) = ∑ j Qt,i,j , here we assume the agents are homogeneous, in other words, each pair of {(oi, ai)}ni=1 coming from the same spaces O,A. They are permutation invariant and their identities do not matter. Thus, the bound of Qi,j is invariant to agents. Under the assumption of static fully complete graph, for any agent Qi =
∑n j=1Qi,j , thus∥∥∥Q̂i −Qi∥∥∥ = ∥∥∥∑j Q̂i,j −∑j Qi,j∥∥∥
≤ n ∥∥∥Q̂i,j −Qi,j∥∥∥ (40)

heading: E.2 PROOF OF LEMMA 1
text: Here we prove Lemma 1.
Proof. Equation 4 says Qt,i := g ( {Qt,i,j}nj=1 ) = ∑ j Ii,jQt,i,j , where Ii,j is an indicator function to denote if two agents are connected. Under the assumption that the static non-complete graph has maximum of number of degrees k and the agents are homogeneous, we have
∑ j Ii,j <= k. Thus,∥∥∥Q̂i −Qi∥∥∥ = ∥∥∥∑j Ii,jQ̂i,j −∑j Ii,jQi,j∥∥∥
≤ k ∥∥∥Q̂i,j −Qi,j∥∥∥ (41)
E.3 INFORMAL PROOF OF LEMMA 2
Here we make some intuitions for the proof of Lemma 2. As we already give the proof sketch in our main paper. Lemma 2 is a direct application of Theorem 2. Here our random variable J1, · · · , Jn ∼ NB(r, p), where {(Ji)}ni=1 represents the number of time points node i needs before it meets node j number of r times.
For the complete static graph, we need at least N sample for the bound in equation 3 to be valid; in other words, we need the trajectory to run at least N time points to collect enough data to estimate our conditional operator accurately Qi,j .
For the dynamic graph, each time t, the two nodes are randomly connected with probability p, if it connects, then we can obtain a valid sample to estimate Qi,j ; if not, then we skip to the next
time step. For the dynamic graph node i, if we take the union set of the nodes i connected over the trajectory path, then the union set could form a static complete graph. We say the two graphs are equivalent. The number of time points needed by i until the N th connection with j for each of the pairs (i, j) follows the same distribution J ∼ NB(N, p). Then we are interested in the expectation of the maximum of J1, . . . , Jn−1, we denote it as J{1,...,n−1}.
E{J{1,...,n−1}} means on average, how many time points (N ′) we need for all nodes other than i at least meets N times with i. Obviously, this N ′ >= N since p ∈ [0, 1]. The calculation of this expectation is solved by (Grabner & Prodinger, 1997). And we also put their result in Theorem 2.

Review summaries:
The authors present a framework and method in which predictive state representations for multiple agents simultatneously acting and interacting within an environment. This is presented in a general way, where predictive states are Hilbert space operators which when applied to sequences of observations and actions appropriately predict the predictions of these state representations. The key advance in this work is to apply this existing PSR framework to networks of agents, with 3 types of agent network: static complete graph (all agents affect all others experience); static non-complete graph (only some agents affect one another); and dynamic non-complete graph (agents affect one another in a time varying way). A number of theoretical results are presented, including PAC bounds for the approximators in the framework.  The authors then present two closely related methods to learn policies alongside these multi agent PSRs in an online way. The first MAPSRL-1 is akin to the independent actor critic (Foerster et al., 2017) and thus may suffer from the apparent non-stationarity of the environment from the perspective of any one agent. MAPSRL-2 addresses this by incorporating the PSR information from other agents into the policy gradient update.  The paper presents a series of experiments based on environments encoded in the  OpenAI Gym MAMujoco system. These are environments presented in previous papers and there are a broad selection of these.  

The paper introduces a predictive state representation (PSR)-based MARL framework. The framework uses a graph representation to model the interactions between agents. Performance bounds are given for the learned predictive state representation.  With the MAPSR, individual policies are trained by replacing the partial observation o with PSR Q. The experiments results show the advantage of the proposed method compared with the baselines experimented in the paper.