
Efficient PAC proofs
In this section we provide proofs for the Efficient-PAC theorems presented in Section . We also note that in this work we use the median-of-means estimate for the sample-mean (Definition REF ), instead of using the more standard empirical-mean estimate. It was recently demonstrated that the former has improved convergence results over the latter in the absence of information about the underlying distribution [1]}, [2]}, [3]}.
Definitions
First, we define the sample set as a collection of previously agent-collected samples. We note that only the learner has access to the full sample set, which is used for value iteration.
[]mydefsampleset

For a state-action \((s,a)\) , the sample set \(u(s,a)\)  is a set of up to \(k\)  tuple samples \(x_{l}\triangleq (s,a,\widetilde{R}_{n(l),t(l)},s_{n(l),t(l)+1})\) , where \(l\)  is the sample's index, \(n(l)\)  is the collecting agent and \(t(l)\)  is the time at which the sample was collected. \(\widetilde{R}_{n(l),t(l)}\)  is the corresponding noisy reward, and \(s_{n(l),t(l)+1}\)  is the value of the next state.

For simplicity of representation, we define the saturated \(Q\) -function to be the restriction of some \(Q\) -function to the set \([0,Q_{max}]\)  for each state-action, for some given \(Q\)  and \(Q_{max}\) .
[]mydefQsat

Given a function \(Q:S\times A \rightarrow \mathbb {R} \)  and a positive value \(Q_{max}\) , the saturated \(Q\) -function is,
\(Q_{\operatornamewithlimits{sat}}(s,a)=\max \lbrace 0,\min \lbrace Q_{\max },Q(s,a))\rbrace \rbrace \) 

The next operator \(F\)  is used to compute the approximate Bellman operator in Algorithm , and is essentially an empiric estimate of the Bellman operator with an additive Upper Confidence Bound (UCB) term (as stated before, we use the median-of-means estimate for the sample mean corresponding to \(R(s,a,s^{\prime })+\gamma Q(s^{\prime },\pi (s^{\prime }))\)  instead of the empirical average). For given integers \(k,k_{m}\)  such that \(k/k_{m}\)  is some power of 2, given that the number of samples in the sample set \(u(s,a)\)  is \(2^{p}k_{m}\)  with some integer \(p\) , we divide the group arbitrarily to \(k_{m}\)  groups with \(2^{p}\)  samples each, calculate the sample average resulting from this group as the operator \(G\) , and then define \(F\)  to be the median of the various \(G\)  operators, with an additive UCB term.
[]mydefFoperator

Let \(k_{m} \ge 1\)  be an integer parameter, \(\pi \)  some fixed policy, \(\epsilon _{b}\)  a real value and \(u(s,a)\)  be a sample set with \(|u(s,a)| \triangleq 2^{p}k_{m}\)  samples for some integer \(p \ge 0\) . The operator \(F^{\pi }(Q,u(s,a))\)  is defined as
\({\begin{array}{c}F^{\pi }(Q,u(s,a)) \triangleq \frac{\epsilon _{b}}{\sqrt{|u(s,a)|}}+\mathrm {median}\lbrace G^{\pi }(Q,u(s,a),1),...,G^{\pi }(Q,u(s,a),k_{m})\rbrace ~,\end{array}}\) 
where
\({\begin{array}{c}G^{\pi }(Q,u(s,a),j) \triangleq \frac{k_{m}}{|u(s,a)|}\sum ^{j\frac{|u(s,a)|}{k_{m}}}_{l=1+(j-1)\frac{|u(s,a)|}{k_{m}}}\left(\widetilde{R}_{n(l),t(l)}+\gamma Q(s_{n(l),t(l)+1},\pi (s_{n(l),t(l)+1}))\right)~,\end{array}}\) 
and \((\widetilde{R}_{n(l),t(l)},s_{n(l),t(l)+1})\)  is the \(l\) -th noisy sample in \(u(s,a)\)  collected by agent \(n(l)\)  at time \(t(l)\) . We will use \(\widetilde{F}\)  to denote \(\widetilde{F}^{\pi ^{Q}}\) , where \(\pi ^{Q}\)  is the greedy policy over \(Q\) .

The approximate Bellman operator \(\widetilde{B}\)  from Algorithm  is defined as a mere saturation.
[]mydefBoperator

For a state-action \((s,a)\) , the approximate Bellman operator \(\widetilde{B}^{\pi }\)  for policy \({\pi }\)  is
\({\begin{array}{c}\widetilde{B}^{\pi }Q(s,a) \triangleq \max \lbrace 0,\min \lbrace Q_{\max },F^{\pi }(Q,u(s,a))\rbrace \rbrace \end{array}}\) 
We will use \(\widetilde{B}\)  to denote \(\widetilde{B}^{\pi ^{Q}}\) . When \(|u(s,a)|=0\)  we set \(\widetilde{B}^{\pi }Q(s,a) \triangleq Q_{\max }\) .

The variance of the Bellman operator \(\sigma \)  is defined to an upper bound on the possible values of the variance for all terms that can be encountered by the algorithm.
[]mydefsig 
The variance of the Bellman operator \(\sigma \)  is defined to be the minimal constant satisfying
\(\sum _{s^{^{\prime }}}p(s^{\prime }|s,a)\left(R(s,a,s^{\prime })+\gamma Q(s^{^{\prime }},\pi ^{Q}(s^{^{\prime }}))-B^{\pi ^{Q}}Q(s,a)\right)^{2} \le \sigma ^{2} \quad \forall (s,a,\pi ^{Q},Q),\) 
where \(Q\)  can be any \(Q\) -function produced by the learner during the run of the algorithm, and \(\pi ^{Q}\)  is the corresponding greedy policy.

An efficient exploration algorithm is defined as an algorithm for which the computational complexity, the space complexity and the sample complexity are bounded. Intuitively, the sample complexity can be thought of as the number of samples that need to be collected before the algorithm achieves the exploratory goal, and is commonly defined as the number of sub-optimal steps in which the algorithm is `far' from its goal. We define a measure for the sample complexity termed the Total Cost of Exploration (TCE) which is shown to dominate the number of sub-optimal steps measure [4]}. This measure takes into account how sub-optimal the steps are by summing the difference between the value function under the policy produced by the algorithm, and the value of the optimal policy. It can be easily shown that if the TCE of an algorithm is bounded by some term \(L\) , then the number of sub-optimal steps measure is bounded by \(L/\epsilon \) .
[]mydefTCE

Let \(\pi \)  be a possibly non-stationary and history dependent policy, \((s_{1},s_{2},...)\)  a random path generated by \(\pi \) , \(\epsilon \)  a positive constant, \(T\)  the (possibly infinite) set of time steps for which \(V^{\pi }(s_{t})<V^{*}(s_{t})-\epsilon \) , and define
\(\begin{aligned}\epsilon _{e}(t) &=V^{*}(s_{t})-V^{\pi }(s_{t})-\epsilon ,\quad \forall t \in T\\\epsilon _{e}(t) &=0,\quad \forall t \notin T\end{aligned}\) 
The Total Cost of Exploration (TCE) for parameter \(\epsilon \)  is defined as the un-discounted infinite sum
\(TCE_{\epsilon } \triangleq \sum _{t=0}^{\infty }\epsilon _{e}(t)\) . To avoid encumbering the notation, we use \(TCE\)  instead of \(TCE_{\epsilon }\)  when the context is clear.

A single-agent exploration algorithm is defined to be Probably Approximately Correct (PAC) if its sample complexity, computational complexity and space complexity are all bounded by some polynomial.
[]mydefPAC

A single agent algorithm is said to be efficient PAC-MDP, if for any \(\varepsilon >0\)  and \(0<\delta <1\) , its sample complexity, its per-time step computational complexity, and its space complexity, are less than some polynomial in the relevant quantities \((S,A,1/\varepsilon ,1/\delta ,1/1-\gamma )\) , with probability of at least \(1-\delta \) .

We will define a multi-agent algorithm to be efficient PAC-MDP, if Definition REF  holds with regards to a sample complexity that is defined as the sum of sample complexities for all agents. By using the TCE measure, the TCE of a multi agent system is thus TCE\(=\sum ^{N}_{i=1}\) TCE\((i)\)  where \(TCE(i)\)  is the TCE for agent \(i\)  for \(\epsilon \) .
[]mydefMAPAC

A multi-agent algorithm is said to be efficient PAC-MDP, if for any \(\varepsilon >0\)  and \(0<\delta <1\) , its sample complexity TCE\(=\sum ^{N}_{i=1}\) TCE\((i)\) , per-time step computational complexity, and space complexity, are all less than some polynomial in the relevant quantities \((S,A,1/\varepsilon ,1/\delta ,1/1-\gamma )\) , with probability of at least \(1-\delta \) . The sample complexity is the sum of sample complexities of all agents.

A flow chart for the dependencies between the various lemmas and theorems we state here can be found in Figure REF .
<FIGURE>
Space and computational complexity bounds
Since Algorithm  is a sub-algorithm for Algorithm , in the following theorems we only state results for the latter.
[]mytheospacecomplexity

The space complexity of Algorithm  is \(O\left((k+N)SA\right)\) .


Each agent in Algorithm  only needs access to the Q-function, and the estimated Q-function is calculated as a weighted average of \(d_{i}+1\)  samples (from the learner and from the other agents). Thus, only one value per state-action has to be saved, leading to a term of \(O\left(NSA\right)\)  for all \(N\)  agents. The learner keeps a copy of the \(Q\) -function, and for each state-action saves up to \(O(k)\)  samples. This results in a term of \(O\left(kSA\right)\) . Overall, the space complexity is \(O\left((k+N)SA\right)\) 

In order to bound the computational complexity of the algorithm, we first quote Lemma 13.3 from [4]} which bounds the number of value iterations needed in order to converge to a fixed point given some confidence level.
[]mylemmaBconv

Let \(\hat{B}\)  be a \(\gamma \) -contraction with fixed point \(\hat{Q}\) , and \(Q\)  the output of \((1/(1-\gamma ))\ln Q_{max}/\epsilon \) 
iterations of value iteration using \(\hat{B}\) . Then if \(Q,Q_{0} \in [0,Q_{max}]\)  where \(Q_{0}\)  is the initial Q-function, then
\(|Q(s,a)-\hat{B}Q(s,a)|\le \epsilon \quad \forall (s,a)\) 

We now show that the operator \(\widetilde{B}\)  is indeed a \(\gamma \) -contraction.
[]mylemmacontraction

\(\widetilde{B}\)  is a \(\gamma \) -contraction under the maximum norm.


given some real parameter \(x\) , note that if we define the saturation function
\(f_{sat}(x)= \max \lbrace f_{min},\min \lbrace f_{max},x\rbrace \rbrace \) 
for some \(f_{max} \ge f_{min}\) , then given some \(x_{0}\ge 0\) , we have that \(f_{sat}(x+x_{0}) \le f_{sat}(x)+x_{0}\) . It also holds that \(f_{sat}\)  is non-decreasing. Now define \(||Q_{1}-Q_{2}||_{\infty }=\epsilon \) . To prove the lemma, we must show that \(||\widetilde{B}Q_{1}-\widetilde{B}Q_{2}||_{\infty }\le \gamma \epsilon \) .
First, note that any function \(G(Q,u(s,a,),j)\)  (see Definition REF ) for \(1 \le j\le k_{m}\)  is a simple average of terms containing maximization over the Q-function, and so under the definition of \(\epsilon \)  above
\(G(Q_{1},u(s,a),j) \le G(Q_{2},u(s,a),j) + \gamma \epsilon \) 
For all state-actions and values of \(j\) . Since the operator F from Definition REF  contains a mere median function over the \(G\)  functions, the following also holds
\(F(Q_{1},u(s,a)) \le F(Q_{2},u(s,a)) + \gamma \epsilon \) 
And now for the approximate Bellman operator we have
\(\widetilde{B}Q_{1}(s,a) = f_{sat}(F(Q_{1},u(s,a))) \le f_{sat}(F(Q_{2},u(s,a)))+\gamma \epsilon =\widetilde{B}Q_{2}(s,a)+\gamma \epsilon \) 
Similarly we can show that \(\widetilde{B}Q_{2}(s,a) \le \widetilde{B}Q_{1}(s,a)+\gamma \epsilon \) , which proves the theorem.

We use the previous results to prove a bound over the computational complexity of the algorithm.
[]mytheocomputationalcomplexity

The per step computational complexity of Algorithm  is
\(O\left(NSAd_{max}+\frac{kSA^{2}}{1-\gamma }\ln \frac{Q_{max}}{\epsilon _{a}}\right)\) 
where \(d_{max}=\max _{i}d_{i}\) .


Lines 3-12 in Algorithm  are performed once for each of the \(N\)  agents. For each agent, we first calculate a weighted average over \(d_{i}\)  samples, and this is done \(SA\)  times - once per every state-action. We then use a greedy policy that requires a search over \(A\)  actions. This results in Algorithm  having a computational complexity bounded by \(O\left(NSAd_{max}\right)\)  where \(d_{max}\triangleq \max _{i}d_{i}\) .
Next, we have that lines 5-12 in Algorithm  are also performed \(N\)  times, once for each sample collected by an agent. Assuming that nullifying an array of samples can be performed in some constant time independent of the number of elements, we have a contribution of \(O(N)\)  to the computational complexity from lines 6-12. Since the operator \(\widetilde{B}\)  is a \(\gamma \) -contraction and since the \(Q\) -function satisfies \(Q \in [0,Q_{max}]\)  at all times by the definition of \(\widetilde{B}\) , we can use Lemma REF  and deduce that the value iteration in lines 15-17 will be performed at most \(\frac{1}{1-\gamma }\ln \frac{Q_{max}}{\epsilon _{a}}\)  times for each state-action. By definitions REF ,REF , A single execution of \(\widetilde{B}\)  for a given state-action requires \(k_{m}\)  calculations of an average at the cost of \(\frac{k}{k_{m}}A\)  operations each, and is thus performed \(kA\)  times.
Thus, the overall computational complexity is
\(O\left(NSAd_{max}+\frac{kSA^{2}}{1-\gamma }\ln \frac{Q_{max}}{\epsilon _{a}}\right)~.\) 

Sample complexity bounds
We now introduce various lemmas that will help us prove a sample complexity bound in Theorem REF , which will be used to prove Theorem .
First, for completeness, we quote two concentration inequalities that we use in theorems REF  and REF . The first one is Cantelli's inequality [6]}.
[]mylemmaCantelli

Let \(X\)  be a real valued random variable with a finite variance \(\sigma ^{2}\) . Then for any \(\lambda > 0 \)
\(P(X-EX \ge \lambda ) \le \frac{\sigma ^{2}}{\sigma ^{2}+\lambda ^{2}},\quad ;\quad P(X-EX \le -\lambda ) \le \frac{\sigma ^{2}}{\sigma ^{2}+\lambda ^{2}},\) 

The second one is McDiarmid's inequality [7]}.
[]mylemmaMcdiarmid

Let \(X_{1},...,X_{n} \in \mathbb {M}\)  be independent random variables, and \(f\)  be a mapping \(f:\mathbb {M}^{n} \rightarrow \mathbb {R}\) . If, for all \(i \in \lbrace 1,...,n\rbrace \) , and for all \(x_{1},...,x_{n},x_{i^{^{\prime }}} \in \mathbb {M}\) , the function \(f\)  satisfies
\(|f(x_{1},...,x_{i-1},x_{i},x_{i+1},...,x_{n})-f(x_{1},...,x_{i-1},x_{i^{^{\prime }}},x_{i+1},...,x_{n})| \le c_{i}~,\) 
for some positive constants \(c_{1},...,c_{n}\) , then \(\forall t \ge 0\)
\(P(f(x_{1},...,x_{n})-Ef(x_{1},...,x_{n}) \ge t) \le \exp \left(\frac{-2t^{2}}{\sum _{i=1}^{n}c_{i}^{2}}\right)\) 
and
\(P(f(x_{1},...,x_{n})-Ef(x_{1},...,x_{n}) \le -t) \le \exp \left(\frac{-2t^{2}}{\sum _{i=1}^{n}c_{i}^{2}}\right)\) 

In lemma REF  we show that the function \(G\)  from Definition REF  is close to the true Bellman operator with some probability.
[]mylemmaGtoB

Let \(\sigma \)  be defined as in Definition REF , \(\epsilon _{b}^{2}=4k_{m}(\sigma ^{2}+\sigma ^{2}_{R})\) , and \(k_{m}={5.6\ln \frac{8N{1+\log _{2}\frac{k}{k_{m}}}(SA)^{2}}{\delta }}\) . Under assumption  given a fixed \(Q\) , a fixed state action \((s,a)\)  and a fixed \(j \in \lbrace 1,...,k_{m}\rbrace \)  we have
\(P\left(G^{\pi }(Q,u(s,a),j)-B^{\pi }Q(s,a) \le -\frac{\epsilon _{b}}{\sqrt{|u(s,a)|}}\right) \le \frac{1}{5}\) 
and
\({\begin{array}{c}P\left(G^{\pi }(Q,u(s,a),j)-B^{\pi }Q(s,a) \ge \frac{\epsilon _{b}}{\sqrt{|u(s,a)|}}\right) \le \frac{1}{5}\end{array}}\) 


Using the definition of \(G\)  (REF ) we have
\({\begin{array}{c}G^{\pi }(Q,u(s,a),j) = \frac{k_{m}}{|u(s,a)|}\sum ^{j\frac{|u(s,a)|}{k_{m}}}_{l=1+(j-1)\frac{|u(s,a)|}{k_{m}}}\left(R_{l}+n_{R,l}+\gamma Q(s^{^{\prime }}_{l},\pi (s^{^{\prime }}_{l}))\right)\end{array}}\) 
Where for convenience we denote \(R_{l}\)  to be the actual reward of sample \(l\)  in this group (which is independent of the collecting agent), \(n_{R,l}\)  to be the reward noise of sample \(l\)  collected by agent \(n(l)\) , and \(s^{^{\prime }}_{l}\)  to be the next state value for sample \(l\) . Since the reward-noise is zero-mean, by taking the expectation operator we have
\({\begin{array}{c}EG^{\pi }(Q,u(s,a),j) = B^{\pi }Q(s,a)\end{array}}\) 
Now, let us denote by \(X_{1},...,X_{\frac{|u(s,a)|}{k_{m}}}\)  the random variables such that \(X_{l}=R_{l}+\gamma Q(s^{^{\prime }}_{l},\pi (s^{^{\prime }}_{l}))\) , and by \(Y_{1},...,Y_{\frac{|u(s,a)|}{k_{m}}}\)  the random variables such that \(Y_{l}=n_{R,l}\) . Since we assume a given state action \((s,a)\)  and index \(j\) , a given \(Q\) -function and use the independence assumptions of assumption , we have that the random variables \(X_{1},...,X_{\frac{|u(s,a)|}{k_{m}}},Y_{1},...,Y_{\frac{|u(s,a)|}{k_{m}}}\)  are all mutually independent. To understand this, note that since we fixed the values of \((s,a),j,Q\)  - the randomness defining these random varibales depends only on the next-state values \(s^{\prime }_{l}\)  and on the noise values \(n_{R,l}\) . Both of these terms are mutually independent given (s,a) from the markov property.
By Definition REF  and assumption  defining \(\sigma \)  and \(\sigma _{R}\) , it is also true that \(VarX_{l} \le \sigma ^{2}\)  and \(VarY_{l} \le \sigma ^{2}_{R}\)  for any sample \(l\) , and from the independence property \(VarG \le \frac{k_{m}}{|u(s,a)|}(\sigma ^{2}+\sigma _{R}^{2})\) .
Now using Cantelli's inequality  REF
\(\begin{aligned}P\left(G^{\pi }(Q,u(s,a),j)-B^{\pi }Q(s,a) \le -\frac{\epsilon _{b}}{\sqrt{|u(s,a)|}}\right) &= \\ P\left(G^{\pi }(Q,u(s,a),j)-EG^{\pi }(Q,u(s,a),j) \le -\frac{\epsilon _{b}}{\sqrt{|u(s,a)|}}\right) &\le \\ \frac{VarG}{VarG+\epsilon _{b}^{2}/|u(s,a)|} \le \frac{\frac{k_{m}}{|u(s,a)|}(\sigma ^{2}+\sigma _{R}^{2})}{\frac{k_{m}}{|u(s,a)|}(\sigma ^{2}+\sigma _{R}^{2})+4\frac{k_{m}}{|u(s,a)|}(\sigma ^{2}+\sigma _{R}^{2})}&=\frac{1}{5}\end{aligned}\) 
We prove the other inequality by the same argument.

[]mylemmaFlemma

Let \(\sigma \)  be defined as in Definition REF , \(\epsilon _{b}\) , and \(k_{m}\)  as in Lemma REF . Under assumption , given a fixed \(Q\)  and a fixed state action \((s,a)\)  we have
\(P\left(F^{\pi ^{*}}(Q,u(s,a))-B^{\pi ^{*}}Q(s,a) \le 0 \right) \le \frac{\delta }{8N{1+\log _{2}\frac{k}{k_{m}}}(SA)^{2}}\) 
and
\({\begin{array}{c}P\left(F^{\pi ^{Q}}(Q,u(s,a))-B^{\pi ^{Q}}Q(s,a) \ge 2\frac{\epsilon _{b}}{\sqrt{|u(s,a)|}}\right) \le \frac{\delta }{8N{1+\log _{2}\frac{k}{k_{m}}}(SA)^{2}}\end{array}}\) 


Define \(Z_{1},...,Z_{k_{m}}\)  as random variables such that \(Z_{j}\)  represents the joint distribution of all the samples in the \(j\) -th group used by the function \(G^{\pi }(Q,u(s,a),j)\) . We now define the following two counter functions. \(f^{+}\)  counts the number of groups in which \(G^{\pi }(Q,u(s,a),j)-B^{\pi }Q(s,a) \le -\frac{\epsilon _{b}}{\sqrt{|u(s,a)|}}\)
\(f^{+}(Z_{1},...,Z_{k_{m}})=\sum ^{k_{m}}_{j=1}\mathbf {1}\lbrace G^{\pi }(Q,u(s,a),j)-B^{\pi }Q(s,a) \le -\frac{\epsilon _{b}}{\sqrt{|u(s,a)|}}\rbrace \) 
Where \(\mathbf {1}\lbrace \dot{\rbrace }\)  is the indicator function, and \(f^{-}\)  counts the number of groups in which \(G^{\pi }(Q,u(s,a),j)-B^{\pi }Q(s,a) \ge \frac{\epsilon _{b}}{\sqrt{|u(s,a)|}}\)
\(f^{-}(Z_{1},...,Z_{k_{m}})=\sum ^{k_{m}}_{j=1}\mathbf {1}\lbrace G^{\pi }(Q,u(s,a),j)-B^{\pi }Q(s,a) \ge \frac{\epsilon _{b}}{\sqrt{|u(s,a)|}}\rbrace \) 
By lemma REF  We have that
\(\begin{aligned}Ef^{+} \le \frac{k_{m}}{5}\\Ef^{-} \le \frac{k_{m}}{5}\end{aligned}\) 
Since \(f^{+},f^{-}\)  count events, it is obvious that a change in the value of one random variable \(Z_{i}\)  results in a bounded change of these functions.
\(\begin{aligned}|f^{+}(Z_{1},...,Z_{i-1},Z_{i},Z_{i+1},...,Z_{k_{m}})-f^{+}(Z_{1},...,Z_{i-1},Z_{i^{^{\prime }}},Z_{i+1},...,Z_{k_{m}})| \le 1 \\|f^{-}(Z_{1},...,Z_{i-1},Z_{i},Z_{i+1},...,Z_{k_{m}})-f^{-}(Z_{1},...,Z_{i-1},Z_{i^{^{\prime }}},Z_{i+1},...,Z_{k_{m}})| \le 1\end{aligned}\) 
By using the same arguments as in lemma REF  the variables \(Z_{1},...,Z_{k_{m}}\)  are independent given \((s,a)\)  and the \(Q\) -function, so we can use Mcdiarmid's inequality REF .
\(\begin{aligned}&P\left(\mathrm {median}\lbrace G^{\pi ^{*}}(Q,u(s,a),1),...,G^{\pi ^{*}}(Q,u(s,a),k_{m})\rbrace -B^{\pi ^{*}}Q(s,a) \le -\frac{\epsilon _{b}}{\sqrt{|u(s,a)|}}\right) \le \\&P(f^{+} \ge \frac{k_{m}}{2}) = P(f^{+} -\frac{k_{m}}{5}\ge \frac{3k_{m}}{10}) \le P(f^{+} -Ef^{+}\ge \frac{3k_{m}}{10}) \le \\ &\exp ^{\frac{-2\left(\frac{3k_{m}}{10}\right)^{2}}{k_{m}}} \le \frac{\delta }{8N{1+\log _{2}\frac{k}{k_{m}}}(SA)^{2}}\end{aligned}\) 
Where the last inequality stems from our definition of \(k_{m}\) . By the same argument
\(\begin{aligned}&P\left(\mathrm {median}\lbrace G^{\pi ^{Q}}(Q,u(s,a),1),...,G^{\pi ^{Q}}(Q,u(s,a),k_{m})\rbrace -B^{\pi ^{Q}}Q(s,a) \ge \frac{\epsilon _{b}}{\sqrt{|u(s,a)|}}\right) \le \\&P(f^{-} \ge \frac{k_{m}}{2}) = P(f^{-} -\frac{k_{m}}{5}\ge \frac{3k_{m}}{10}) \le P(f^{-} -Ef^{-}\ge \frac{3k_{m}}{10}) \le \\ &\exp ^{\frac{-2\left(\frac{3k_{m}}{10}\right)^{2}}{k_{m}}} \le \frac{\delta }{8N{1+\log _{2}\frac{k}{k_{m}}}(SA)^{2}}\end{aligned}\) 
The rest follows from the definition of \(F\)  as a median with an additive exploration term (REF ).

We have seen that the operator \(F\)  is close to the Bellman operator with some probability for the real \(Q\) -function the learner possesses. Now we shall show that for each agent, the same holds for the estimated Q-function it holds. This means that we have to show that the above relation holds for \(\hat{Q}^{i}_{sat}\)  and for all agents.
We will use the following assumption over the noise terms.
[]myassumpConcentrationBound

For all agents \(i \in \lbrace 1,...,N\rbrace \) , estimated \(Q\) -functions \(\hat{Q}^{i}\)  and state actions \((s,a)\)  during the run of the algorithm, the noise terms \(n^{i}\)  in (REF ) satisfy
\( \begin{aligned}P(|n^{i}(s,a)-En^{i}(s,a)|\le n_{c}^{i}) &\ge 1-\delta _{L}\\ |En^{i}(s,a)| &\le \mu _{c}^{i}\end{aligned}\) 
for some \(\delta _{L},\mu _{c}^{i},n_{c}^{i}\)  given \((s,a)\) .Note that this relation holds for a noise term given the state-action \((s,a)\) , and does not include the probability of encountering \((s,a)\)  itself

For convenience, we define the following events
\(\begin{aligned}I_{1} &\triangleq \lbrace F^{\pi ^{*}}(Q,u(s,a))-B^{\pi ^{*}}Q(s,a) \le 0\rbrace  ~,\\J_{1}&\triangleq \left\lbrace F^{\pi ^{Q}}(Q,u(s,a))-B^{\pi ^{Q}}Q(s,a) \ge 2\frac{\epsilon _{b}}{\sqrt{|u(s,a)|}}\right\rbrace  ~,\\I_{2} &\triangleq \lbrace F^{\pi ^{*}}(\hat{Q}^{i}_{sat},u(s,a))-B^{\pi ^{*}}\hat{Q}^{i}_{sat}(s,a) \le -2\gamma (n_{c}^{i}+\mu _{c}^{i})\rbrace  ~,\\J_{2} &\triangleq \left\lbrace F^{\pi ^{\hat{Q}^{i}_{sat}}}(\hat{Q}^{i}_{sat},u(s,a))-B^{\pi ^{\hat{Q}^{i}_{sat}}}\hat{Q}^{i}_{sat}(s,a) \ge 2\frac{\epsilon _{b}}{\sqrt{|u(s,a)|}}+2\gamma (n_{c}^{i}+\mu _{c}^{i})\right\rbrace ~,\\L_{2} &\triangleq \lbrace |\hat{Q}^{i}_{sat}(s,a)-\widetilde{B}\hat{Q}^{i}_{sat}(s,a)| > \epsilon _{a}+(1+\gamma ) (n_{c}^{i}+\mu _{c}^{i})\rbrace ~.\end{aligned}\) 
[]mylemmalemma2

Let \(\sigma \)  be defined as in Definition REF ,\(\epsilon _{b}\) , and \(k_{m}\)  as in Lemma REF . Under assumptions  and REF , given a fixed \(Q\) , a fixed state action \((s,a)\)  and a fixed agent \(i\)
\(\begin{aligned}P(I_{2}) &\le \delta _{L}SA+\frac{\delta }{8N{1+\log _{2}\frac{k}{k_{m}}}(SA)^{2}} \\P(J_{2}) &\le \delta _{L}SA+\frac{\delta }{8N{1+\log _{2}\frac{k}{k_{m}}}(SA)^{2}}\\P(L_{2}) &\le \delta _{L}SA\end{aligned}\) 


Let us denote
\(W_{\hat{Q}^{i}} \triangleq \left\lbrace \forall (\overline{s},\overline{a}) :|n^{i}(\overline{s},\overline{a})-En^{i}(\overline{s},\overline{a})|\le n_{c}^{i} \right\rbrace \) 
As the event for which (REF ) holds for all state actions corresponding to \(\hat{Q}^{i}\)  stemming from the fixed \(Q\) . Since there are exactly \(SA\)  noise terms in this event, we have from the Union Bound that \(P(W_{\hat{Q}^{i}}) \ge 1- \delta _{L}SA\) . We shall show that
\(\begin{aligned}P(I_{2}|W_{\hat{Q}^{i}})&\le \frac{\delta }{8N{1+\log _{2}\frac{k}{k_{m}}}(SA)^{2}} \\P(J_{2}|W_{\hat{Q}^{i}})&\le \frac{\delta }{8N{1+\log _{2}\frac{k}{k_{m}}}(SA)^{2}}\\P(L_{2}|W_{\hat{Q}^{i}})&=0\end{aligned}\) 
We then have
\(\begin{aligned}P(I_{2}) &=P(I_{2}|W_{\hat{Q}^{i}})P(W_{\hat{Q}^{i}})+P(I_{2}|\overline{W_{\hat{Q}^{i}}})P(\overline{W_{\hat{Q}^{i}}}) \le \frac{\delta }{8N{1+\log _{2}\frac{k}{k_{m}}}(SA)^{2}}\cdot 1+1\cdot \delta _{L}SA\end{aligned}\) 
and similarly for \(J_{2}\) , and
\(\begin{aligned}P(L_{2})=P(L_{2}|W_{\hat{Q}^{i}})P(W_{\hat{Q}^{i}})+P(L_{2}|\overline{W_{\hat{Q}^{i}}})P(\overline{W_{\hat{Q}^{i}}}) \le 0\cdot 1+1\cdot \delta _{L}SA = \delta _{L}SA~.\end{aligned}\) 
Now, to prove the first probability bound above, we assume that the event \(W_{\hat{Q}^{i}}\)  is true and that \(\overline{I}_{1}\)  is also true, where for event \(X\) , \(\overline{X}\)  is its complement. We then have
\({\begin{array}{c}\forall (\overline{s},\overline{a}): |n^{i}(\overline{s},\overline{a})-En^{i}(\overline{s},\overline{a})|\le n_{c}^{i}\end{array}}\) 
implying that
\(\begin{aligned}En^{i}(\overline{s},\overline{a})-n_{c}^{i}&\le n^{i}(\overline{s},\overline{a}) \le n_{c}^{i}+En^{i}(\overline{s},\overline{a}) \\-\mu _{c}^{i}-n_{c}^{i}&\le n^{i}(\overline{s},\overline{a}) \le n_{c}^{i}+\mu _{c}^{i}\end{aligned}\) 
By using Definition REF  for the saturated \(Q\) -function, we know that for each \((\overline{s},\overline{a})\)
\(\begin{aligned}\hat{Q}^{i}_{sat}(\overline{s},\overline{a}) &=\max \lbrace 0,\min \lbrace Q_{\max },\hat{Q}^{i}(\overline{s},\overline{a})\rbrace \rbrace  \\&\ge \max \lbrace 0,\min \lbrace Q_{\max },Q(\overline{s},\overline{a})-n_{c}^{i}-\mu _{c}^{i}\rbrace \rbrace  \\&\ge Q(\overline{s},\overline{a}) -n_{c}^{i}-\mu _{c}^{i}\end{aligned}\) 
By using the definition of the Bellman operator
\(\begin{aligned}B^{\pi ^{*}}\hat{Q}^{i}_{sat}(s,a) &= \sum _{s^{\prime }}p(s^{\prime }|s,a)\left[R(s,\pi ^{*}(s),s^{\prime })+\gamma \hat{Q}^{i}_{sat}(s^{\prime },\pi ^{*}(s^{\prime }))\right] \\&\le \sum _{s^{\prime }}p(s^{\prime }|s,a)\left[R(s,\pi ^{*}(s),s^{\prime })+\gamma Q(s^{\prime },\pi ^{*}(s^{\prime })) + \gamma n_{c}^{i}+\gamma \mu _{c}^{i}\right] \\&= B^{\pi ^{*}}Q(s,a) + \gamma n_{c}^{i}+\gamma \mu _{c}^{i}\end{aligned}\) 
where we used the assumption that the noise is bounded for all state-actions. By using the same trick for the operator \(G^{\pi ^{*}}\)  we have that for all \(j\)
\(\begin{aligned}G^{\pi ^{*}}(\hat{Q}^{i}_{sat},u(s,a),j)&= \frac{k_{m}}{|u(s,a)|}\sum _{l}(\widetilde{R}_{n(l),t(l)}+\gamma \hat{Q}^{i}_{sat}(s_{l}^{^{\prime }},\pi ^{*}(s_{l}^{^{\prime }})))\\&\ge \frac{k_{m}}{|u(s,a)|}\sum _{l}(\widetilde{R}_{n(l),t(l)}+\gamma Q(s_{l}^{^{\prime }},\pi ^{*}(s_{l}^{^{\prime }}))-\gamma n_{c}^{i}-\gamma \mu _{c}^{i}) \\ &=G^{\pi ^{*}}(Q,u(s,a),j) -\gamma n_{c}^{i}-\gamma \mu _{c}^{i}\end{aligned}\) 
We then have that
\(\begin{aligned}F^{\pi ^{*}}(\hat{Q}^{i}_{sat},u(s,a)) &= \frac{\epsilon _{b}}{\sqrt{|u(s,a)|}}+\overline{\mathrm {med}} \\&\ge \frac{\epsilon _{b}}{\sqrt{|u(s,a)|}}+\overline{\mathrm {med}} - \gamma n_{c}^{i}-\gamma \mu _{c}^{i}\\&=F^{\pi ^{*}}(Q,u(s,a)) - \gamma n_{c}^{i}-\gamma \mu _{c}^{i}\end{aligned}\) 
where \(\overline{\mathrm {med}}=\mathrm {median}\lbrace G^{\pi ^{*}}(\hat{Q}^{i}_{sat},u(s,a),1),...,G^{\pi ^{*}}(\hat{Q}^{i}_{sat},u(s,a),k_{m})\rbrace \) .
Therefore
\(\begin{aligned}F^{\pi ^{*}}(\hat{Q}^{i}_{sat},u(s,a))-B^{\pi ^{*}}\hat{Q}^{i}_{sat}(s,a) &\ge F^{\pi ^{*}}(Q,u(s,a))-B^{\pi ^{*}}Q(s,a)\\&\quad -2\gamma (n_{c}^{i}+\mu _{c}^{i}) \\&\ge -2\gamma (n_{c}^{i}+\mu _{c}^{i})\end{aligned}\) 
Since \(\overline{I}_{1} \subseteq \overline{I}_{2}\)  given \(W_{\hat{Q}^{i}}\)
\(P(I_{2}|W_{\hat{Q}^{i}}) \le P(I_{1}|W_{\hat{Q}^{i}}) \le \frac{\delta }{8N{1+\log _{2}\frac{k}{k_{m}}}(SA)^{2}}\) 
Note that from assumption , the noise terms \(n^{i}\)  are independent of \(Q\)  and of the samples collected so far, which means that the events \(I_{1}\)  and \(W\)  are independent so that \(P(I_{1}|W)=P(I_{1})\) . We use the same principle to prove the second inequality. Note that in that case, since we use the Bellman operator and the \(F\)  operator with the greedy policy, we use the fact that \(\max _{a^{\prime }}\hat{Q}^{i}_{sat}(s^{\prime },a^{\prime }) \ge \max _{a^{\prime }}Q(s^{\prime },a^{\prime }) - n_{c}^{i}-\mu _{c}^{i}\) . For the third inequality, note that we know that the Bellman iterations for the \(Q\) -function \(Q\)  have converged with probability 1,
\({\begin{array}{c}|Q(s,a)-\widetilde{B}Q(s,a)| \le \epsilon _{a}\end{array}}\) 
and using the definition of \(\widetilde{B}\)
\(\begin{aligned}\widetilde{B}\hat{Q}^{i}_{sat}(s,a)&=\min \lbrace Q_{\max },\max \lbrace 0,\widetilde{F}(\hat{Q}^{i}_{sat},u(s,a))\rbrace \rbrace  \\&\le \min \lbrace Q_{\max },\max \lbrace 0,\widetilde{F}(Q,u(s,a))+\gamma (n_{c}^{i}+\mu _{c}^{i})\rbrace \rbrace  \\&\le \widetilde{B}Q(s,a) + \gamma (n_{c}^{i}+\mu _{c}^{i})\end{aligned}\) 
and
\(\begin{aligned}\widetilde{B}\hat{Q}^{i}_{sat}(s,a)&=\min \lbrace Q_{\max },\max \lbrace 0,\widetilde{F}(\hat{Q}^{i}_{sat},u(s,a))\rbrace \rbrace  \\&\ge \min \lbrace Q_{\max },\max \lbrace 0,\widetilde{F}(Q,u(s,a))-\gamma (n_{c}^{i}+\mu _{c}^{i})\rbrace \rbrace  \\&\ge \widetilde{B}Q(s,a) - \gamma (n_{c}^{i}+\mu _{c}^{i})\end{aligned}\) 
Therefore, in that case
\(\begin{aligned}\hat{Q}^{i}_{sat}(s,a)-\widetilde{B}\hat{Q}^{i}_{sat}(s,a) &\le Q(s,a)-\widetilde{B}Q(s,a) +(1+\gamma )(n_{c}^{i}+\mu _{c}^{i}) \\ &\le \epsilon _{a}+(1+\gamma )(n_{c}^{i}+\mu _{c}^{i})\end{aligned}\) 
And by the same arguments
\(\begin{aligned}\hat{Q}^{i}_{sat}(s,a)-\widetilde{B}\hat{Q}^{i}_{sat}(s,a) &\ge Q(s,a)-\widetilde{B}Q(s,a) -(1+\gamma )(n_{c}^{i}+\mu _{c}^{i}) \\&\ge -\epsilon _{a}-(1+\gamma )(n_{c}^{i}+\mu _{c}^{i})\end{aligned}\) 
Therefore, \(P(L_{2}|W_{\hat{Q}^{i}})=0\)  as desired.
We have proven a probability bound for a given \(Q\) -function, state-action \((s,a)\)  and an agent (which received the above mentioned version of the \(Q\) -function). Next, we use the union bound to show that this holds for all \(Q\) -functions and agents during the run of the algorithm.
[]mylemmalemma3

Let \(\sigma \)  be defined as in Definition REF , \(\epsilon _{b}^{2}=4k_{m}(\sigma ^{2}+\sigma ^{2}_{R})\) , and \(k_{m}={5.6\ln \frac{8N{1+\log _{2}\frac{k}{k_{m}}}(SA)^{2}}{\delta }}\) . Assume that assumptions  and REF  hold for all \(Q\) -functions, state-actions and agents that can be encountered during the run of the algorithm. Then the events \(\overline{I}_{2},\overline{J}_{2},\overline{L}_{2}\)  defined before Lemma REF  occur for all the noisy versions that can be encountered during the run of Algorithm  \(\hat{Q}^{i}\) , for all state-actions and agents simultaneously with probability larger than \(1-3N{1+\log _{2}\frac{k}{k_{m}}}(SA)^{3}\delta _{L}-\frac{\delta }{4}\) .


We use a union bound over the number of agents, the number of different \(Q\) -functions the learner produces that can be encountered during the run of the algorithm, and the number of state-actions. Note that a particular \(Q\) -function \(Q\)  is determined by the samples that have been used to calculate it (via Value Iteration) up to that point. As seen in Algorithm , there are at most \({1+\log _{2}\frac{k}{k_{m}}}\)  distinct sample sets with \(|u(s,a)|>0\)  for a given state-action, meaning that there are at most \({1+\log _{2}\frac{k}{k_{m}}}SA\)  distinct \(Q\) -functions encountered in practice.
Therefore, there are at most \(N(SA)^{2}{1+\log _{2}\frac{k}{k_{m}}}\)  possibilities to include in the union bound sum, for each one of the three events above. This lemma is then proved by using the union bound over the results of Lemma REF .

In lemma REF , we show that the events \(\overline{I}_{2},\overline{J}_{2}\)  bounding the difference between \(F^{\pi ^{*}}\)  to \(B^{\pi ^{*}}\)  and \(F^{\pi ^{\hat{Q}^{i}_{sat}}}\)  to \(B^{\pi ^{\hat{Q}^{i}_{sat}}}\)  hold not only for the function \(F\) , but also for the estimated \(Q\)  functions \(\hat{Q}^{i}_{sat}\)  used by the different agents.
[]mylemmalemma4

Let \(\sigma \)  be defined as in Definition REF , \(\epsilon _{b}^{2}=4k_{m}(\sigma ^{2}+\sigma ^{2}_{R})\) , and \(k_{m}={5.6\ln \frac{8N{1+\log _{2}\frac{k}{k_{m}}}(SA)^{2}}{\delta }}\) . Assume that assumptions  and REF  hold for all \(Q\) -functions, state-actions and agents that can be encountered during the run of the algorithm. Then
\(\begin{aligned}\hat{Q}^{i}_{sat}(s,a)-B^{\pi ^{*}}\hat{Q}^{i}_{sat}(s,a) &> -\epsilon _{a}-(n_{c}^{i}+\mu _{c}^{i})(1+3\gamma ) \\\hat{Q}^{i}_{sat}(s,a)-B^{\pi ^{Q}}\hat{Q}^{i}_{sat}(s,a) &< \epsilon _{a}+\frac{2\epsilon _{b}}{\sqrt{|u(s,a)|}}+(n_{c}^{i}+\mu _{c}^{i})(1+3\gamma )\end{aligned}\) 
for all \(Q\) -functions, state-actions \((s,a)\)  and for all agents simultaneously with probability larger than \(1-3N{1+\log _{2}\frac{k}{k_{m}}}(SA)^{3}\delta _{L}-\frac{\delta }{4}\)  .


Until collected samples are used for value iteration, we have that for all agents and state-actions \(\hat{Q}^{i}_{sat}=Q_{\max }\) . Since \(B^{\pi ^{*}}\hat{Q}^{i}_{sat}(s,a) \le Q_{\max }\) . we have that \(\hat{Q}^{i}_{sat}(s,a)-B^{\pi ^{*}}\hat{Q}^{i}_{sat}(s,a) \ge 0 > -\epsilon _{a}-(n_{c}^{i}+\mu _{c}^{i})(1+3\gamma )\)  as wanted (as we shall see in Theorem REF , we won't be needing the second inequality in this lemma for the case of such a \(Q\) -function).
Otherwise, we have that under the probability mentioned above, for all \(Q\) -functions, state-actions and agents that can be encountered during the run of the algorithm \((s,a,\hat{Q}^{i})\)  (see Lemma REF ),
\(\begin{aligned}B^{\pi ^{*}}\hat{Q}^{i}_{sat}(s,a) &= \min \lbrace Q_{\max },\max \lbrace 0,B^{\pi ^{*}}\hat{Q}^{i}_{sat}(s,a)\rbrace \rbrace  \\&\le \min \left\lbrace Q_{\max },\max \lbrace 0,F^{\pi ^{*}}(\hat{Q}^{i}_{sat},u(s,a)))+2\gamma (n_{c}^{i}+\mu _{c}^{i})\rbrace \right\rbrace  \\&\le \widetilde{B}^{\pi ^{*}}\hat{Q}^{i}_{sat}(s,a)+2\gamma (n_{c}^{i}+\mu _{c}^{i}) \\&\le \widetilde{B}\hat{Q}^{i}_{sat}(s,a)+2\gamma (n_{c}^{i}+\mu _{c}^{i}) \\&\le \hat{Q}^{i}_{sat}(s,a)+\epsilon _{a}+(1+3\gamma )(n_{c}^{i}+\mu _{c}^{i})\end{aligned}\) 
The first equality stems from the fact that we use the true Bellman operator on a \(Q\) -function that is bounded in \([0,Q_{max}]\) , and the result of applying the operator to this function is also bounded. The inequality that follows stems from \(\overline{I}_{2}\) . The second and third inequalities are due to the definition of the approximate Bellman operator \(\widetilde{B}\)  (Definition REF ), and the last inequality stems from \(\overline{L}_{2}\) . Similarly
\(\begin{aligned}B^{\pi ^{\hat{Q}^{i}_{sat}}}\hat{Q}^{i}_{sat}(s,a) &= \min \lbrace Q_{\max },\max \lbrace 0,B^{\pi ^{\hat{Q}^{i}_{sat}}}\hat{Q}^{i}_{sat}(s,a)\rbrace \rbrace  \\&\ge \min \lbrace Q_{\max },\max \lbrace 0,F^{\pi ^{\hat{Q}^{i}_{sat}}}(\hat{Q}^{i}_{sat},u(s,a)))-2\gamma (n_{c}^{i}+\mu _{c}^{i})-\widetilde{\epsilon }_b\rbrace \rbrace  \\&\ge \widetilde{B}^{\pi ^{\hat{Q}^{i}_{sat}}}\hat{Q}^{i}_{sat}(s,a)-2\gamma (n_{c}^{i}+\mu _{c}^{i})-\widetilde{\epsilon }_b \\&= \widetilde{B}\hat{Q}^{i}_{sat}(s,a)-2\gamma (n_{c}^{i}+\mu _{c}^{i})-\widetilde{\epsilon }_b \\&\ge \hat{Q}^{i}_{sat}(s,a)-\epsilon _{a}-(1+3\gamma )(n_{c}^{i}+\mu _{c}^{i}) -\widetilde{\epsilon }_b\end{aligned}\) 
where \(\widetilde{\epsilon }_b=2\epsilon _{b}/\sqrt{|u(s,a)|}\) .

Before we continue, we quote two lemmas from [4]} that will be useful ahead. Lemma REF  is a basic result regarding Bernoulli random variables.
[]mylemmaBern

Let \(Y_{i}\)  for \(i \in \lbrace 1,...,n\rbrace \)  be independent Bernoulli random variables such that \(P(Y_{i}=1) \ge p_{i} \quad \forall i\) . for some \(m,\delta > 0\)  , if \(\frac{2}{m}\ln \frac{1}{\delta }<1\)  and
\(\sum _{i=1}^{n} \ge \frac{m}{1-\sqrt{\frac{2}{m}\ln \frac{1}{\delta }}}\) 
then
\(P\left(\sum _{i=1}^{n}Y_{i} \ge m \right) \ge 1- \delta \) 

Lemma REF  is a basic result regarding RL in MDPs, and bounds the difference between the optimal value function and a value function resulting from a given greedy policy, given some MDP in which the Bellman error is bounded.
[]mylemmaBellmanDiff

Given some \(Q\) -function \(Q\)  and non-negative constants \(\epsilon _{*},\epsilon _{\pi ^{Q}},\epsilon _{1},...,\epsilon _{n}\)  such that \(\epsilon _{\pi ^{Q}}\le \epsilon _{i} \quad \forall i\) , let \(Q(s,a)-B^{\pi ^{*}}Q(s,a) \ge \epsilon _{*}\)  for all \((s,a)\) , and let \(X_{1},...,X_{n}\)  be sets of state-actions such that \(Q(s,a)-B^{\pi ^{Q}}Q(s,a) \le \epsilon _{i}\)  for all states-actions \((s,a)\)  in \(X_{i}\) , where \(\pi ^{Q}\)  is the greedy policy over \(Q\) . Also, let \(Q(s,a)-B^{\pi ^{Q}}Q(s,a) \le \epsilon _{\pi ^{Q}}\)  for all \((s,a) \notin \cup _{i=1}^{n}X_{i}\) . Let \(T_{H}={\frac{1}{1-\gamma }\ln \frac{Q_{\max }}{\epsilon _{s}}}\)  and define \(H=\lbrace 1,2,4,...,2^{i}\rbrace \)  where i is the largest integer such that \(2^{i} \le T_{H}\) . Define \(p_{h,i}(s)\)  for an integer \(h \ge 0\)  to be the probability of encountering exactly \(h\)  state-actions \((s,a)\)  for which \((s,a) \in X_{i}\)  when starting from state \(s\)  and following \(\pi ^{Q}\)  for a totall of \(\min \lbrace T,T_{H}\rbrace \)  for some \(T\) . Finally, let \(p^{e}_{h,i}(s)=\sum _{m=h}^{2h-1}p_{m,i}(s)\) . Then
\(V^{*}(s)-V^{\pi ^{Q}}(s) \le \frac{\epsilon _{*}+\epsilon _{\pi ^{Q}}}{1-\gamma } +\epsilon _{s}+\epsilon _{e}\) 
Where \(\epsilon _{e} \triangleq 2\sum _{i=1}^{n}\sum _{h\in H}hp^{e}_{h,i}(s)(\epsilon _{i}-\epsilon _{\pi ^{Q}})+\gamma ^{T}Q_{max}\) .

We will use lemma REF  together with REF  and REF  to prove the main theorem REF . Lemma REF  bounds the term \(\epsilon _{e}\)  in lemma REF  for our algorithm.
[]mylemmalemma5

Let \((s_{1,i},s_{2,i},...,)\)  for \(i \in \lbrace 1,...,N\rbrace \)  be the random paths generated on some execution of Algorithm . Let \(\tau (t)\)  be the number of steps from step t to the next step step \(t^{\prime }\)  for which the policy changes. Let \(T_{H}={\frac{1}{1-\gamma }\ln \frac{Q_{\max }}{\epsilon _{s}}}\)  and define \(H=\lbrace 1,2,4,...,2^{i}\rbrace \)  where i is the largest integer such that \(2^{i} \le T_{H}\) . Let \(K_{a}=\lbrace k_{m},2k_{m},2^{2}k_{m},...,k\rbrace \) . Let \(k_{a}^{-}\)  be the largest value in \(K_{a}\)  that is strictly smaller than \(k_{a}\) , or 0 if such a value does not exist. Let \(X_{k_{a}}(t)\)  be the set of state-actions \((s,a)\)  at step t for which \(k_{a}^{-}=|u(s,a)|\) . Define \(p_{h,k_{a}}(s_{t,i})\)  for \(k_{a} \in K_{a}\)  to be the following conditional probability: Given \(\hat{Q}^{i}_{sat}\)  at step t, exactly \(h\)  state-actions in \(X_{k_{a}}(t)\)  are encountered by agent \(i\)  during the next \(\min \lbrace T_{H},\tau (t)\rbrace \)  steps. Let \(p_{h,k_{a}}^{e}(s_{t,i}) \triangleq \sum _{m=h}^{2h-1}p_{m,k_{a}}(s_{t,i})\) . If \(\frac{2T_{H}}{NSA}\ln \frac{2|K_{a}|}{\delta }<1\)  and \(SA\ge 2\) , with probability \(1-\delta /2\)
\({\begin{array}{c}\sum _{i=1}^{N}\sum _{t=0}^{\infty }\sum _{h \in H}hp_{h,k_{a}}^{e}(s_{t,i}) < \frac{SA(k_{a}+N)(1+\log _{2}T_{H})T_{H}}{1-\sqrt{\frac{2T_{H}}{SA(k_{a}+N)}\ln \frac{2|K_{a}|}{\delta }}}\end{array}}\) 
for all \(k_{a} \in K_{a}\)  simultaneously.


Let us fix some \(k_{a}\)  and denote by \(Y^{e}_{h,k_{a}}(s_{t,i})\)  the Bernoulli random variables corresponding to the probabilities \(p_{h,k_{a}}^{e}(s_{t,i})\) . Since we condition these variables and are only interested in their outcome given this condition, we have from the Markov property that variables \(Y^{e}_{h,k_{a}}(s_{t,i})\)  at least \(T_{H}\)  time steps apart are independent, and variables for different agents are independent as well. Define \(T^{H}_{j}\)  for \(j \in \lbrace 0,1,...,T_{H}-1\rbrace \)  to be the infinite set of time steps such \(T^{H}_{j} = \lbrace j,j+T_{H},j+2T_{H},...\rbrace \) . From Algorithm , \(k_{a}\)  samples will be stored in the temporary sample set \(u^{tmp}(s,a)\)  before a state action \((s,a)\)  with \(|u(s,a)|=k_{a}^{-}\)  advances to have \(k_{a}\)  samples. And since \(N\)  agents are exploring in parallel, at most \(k_{a}+N\)  samples will be collected for such a state-action before it progresses to have \(|u(s,a)|=k_{a}\)  (this is the worst case where all the agents visit the same \((s,a)\)  at the last time step before its sample set is updated).
Let us assume that there exists a \(j \in \lbrace 0,1,...,T_{H}-1\rbrace \)  and \(h \in H\)  such that
\(\sum _{i=1}^{N}\sum _{t\in T^{H}_{j}}p_{h,k_{a}}^{e}(s_{t,i}) \ge \frac{SA(k_{a}+N)}{h\left(1-\sqrt{\frac{2h}{SA(k_{a}+N)}\ln \frac{2|K_{a}|}{\delta }}\right)}\) 
However, since this is a sum of probabilities of independent Bernoulli variables, we have from lemma REF  that
\(P\left(\sum _{i=1}^{N}\sum _{t\in T^{H}_{j}}Y_{h,k_{a}}^{e}(s_{t,i}) \ge \frac{SA(k_{a}+N)}{h}\right) \ge 1-\frac{\delta }{2|K_{a}|}\) 
This event is a contradiction, since it means that more than \(SA(k_{a}+N)\)  samples are collected for state-actions with \(k_{a}^{-}\)  samples. Therefore, we have that
\(\sum _{i=1}^{N}\sum _{t\in T^{H}_{j}}p_{h,k_{a}}^{e}(s_{t,i}) < \frac{SA(k_{a}+N)}{h\left(1-\sqrt{\frac{2h}{SA(k_{a}+N)}\ln \frac{2|K_{a}|}{\delta }}\right)}\) 
For all \(j \in \lbrace 0,1,...,T_{H}-1\rbrace \)  and \(h \in H\)  simultaneously, with a probability larger than \(1-\frac{\delta }{2|K_{a}|}\) . By summing over the values of \(h\)  and \(j\)  we have
\(\sum _{i=1}^{N}\sum _{t=0}^{\infty }\sum _{h \in H}hp_{h,k_{a}}^{e}(s_{t,i}) < \frac{SA(k_{a}+N)(1+\log _{2}T_{H})T_{H}}{1-\sqrt{\frac{2T_{H}}{SA(k_{a}+N)}\ln \frac{2|K_{a}|}{\delta }}}\) 
with a probability larger than \(1-\frac{\delta }{2|K_{a}|}\) . Since this is true for a fixed value of \(k_{a}\) , we now use the Union Bound to conclude that the previous equation holds for all \(k_{a} \in K_{a}\)  simultaneously with probability larger than \(1-\frac{\delta }{2}\) .
Note that Since we have defined the settings of Algorithm  so that a learner only sends a new \(Q\) -function to the agents when its own \(Q\) -function is updated, the agents have a noisy but constant version of the \(Q\) -function \(\hat{Q}^{i}_{sat}\)  for \(\tau (t)\)  steps - which is crucial for the proof. It also agrees with the determination principle for multi-agent systems defined in [9]}.

We are now ready to present the main sample complexity bound. As stated before, this theorem combines the lemmas REF ,REF ,REF  to prove PAC results.
[]mytheosamplecomplexity

Under assumptions ,  and REF , assume further that \(\delta _{L}\le \frac{\delta }{12N{1+\log _{2}\frac{k}{k_{m}}}(SA)^{3}}\) . Then with probability at least \(1-\delta \) , for all \(t\)  and \(i\) ,
\({\begin{array}{c}V^{\widetilde{\pi }_{i}}(s_{t,i}) \ge V^{*}(s_{t,i}) - \frac{2\epsilon _{a}+2(1+3\gamma )(n_{c}^{i}+\mu _{c}^{i})}{1-\gamma } -3\epsilon _{s} -\epsilon _{e}^{i}(t)\end{array}}\) 
where
\({\begin{array}{c}TCE = \sum _{i=1}^{N}\sum _{t=0}^{\infty }\epsilon _{e}^{i}(t) = \widetilde{O}\left(\left(N\left(Q_{\max }+\sqrt{\sigma ^{2}+\sigma ^{2}_{R}}\right)+\frac{\sigma ^{2}+\sigma ^{2}_{R}}{\epsilon _{s}(1-\gamma )}\right)SA/(1-\gamma )\right)~,\end{array}}\) 
and \(\widetilde{O}\)  stands for a big-O up to logarithmic terms.


In short, the non-stationary policy of each one of the agents can be broken up into fixed-policy segments, in which we follow an estimated \(Q\) -function greedily. Given that the Bellman error for each segment is acceptably small, lemma REF  states that the greedy policy for that segment has a bounded error with respect to the optimal policy . We use Lemma REF  to show that with high probability, the Bellman errors of all \(Q\) -functions for all agents during the run of the algorithm are bounded, and then combine it with Lemma REF  which bounds the number of times we can encounter state-actions for which we haven't collected enough samples yet. Combining these three theorems, we get to a sample complexity bound for the algorithm as a whole.
Corresponding to lemma REF , the groups we divide the state-actions to for each agent are \(X_{k_{a}}\)  as defined in Lemma REF . We now have to show that state-actions in each group have the same bound on their Bellman error (while different groups may have different bounds). According to Lemma REF  we have that with probability larger than \(1-\frac{\delta }{2}\) , for all state-actions during the run of the algorithm
\({\begin{array}{c}\hat{Q}^{i}_{sat}(s,a)-B^{\pi ^{*}}\hat{Q}^{i}_{sat}(s,a) > -\epsilon _{a}-(n_{c}^{i}+\mu _{c}^{i})(1+3\gamma )\end{array}}\) 
For the required upper bound on the Bellman error, we will divide the state-actions to groups.

For the first group \(X_{k_{m}}\)  note that none of the previous theorems holds, since these are state-actions with no samples in their sample set. Therefore, we simply use the fact that the \(Q\) -functions are all bounded by our definition:
\(\forall (s,a) \in X_{k_{m}}: \hat{Q}^{i}_{sat}(s,a)-B^{\pi ^{\hat{Q}^{i}_{sat}}}\hat{Q}^{i}_{sat}(s,a) \le Q_{\max }\) 

For all \((s,a) \in X_{k_{a}},k_{a}>0\) , we have again from Lemma REF  that
\(\begin{aligned}\hat{Q}^{i}_{sat}(s,a)-B^{\pi ^{Q}}\hat{Q}^{i}_{sat}(s,a) < \epsilon _{a}+(n_{c}^{i}+\mu _{c}^{i})(1+3\gamma )+2\sqrt{\frac{8k_{m}(\sigma ^{2}+\sigma ^{2}_{R})}{k_{a}}}\end{aligned}\) 
Where we substitute \(\epsilon _{b}^{2}=4k_{m}(\sigma ^{2}+\sigma _{R}^{2})\) , and use the fact that the fact that \(k_{a}^{-}\ge \frac{k_{a}}{2}\) .

Finally, for state-actions with \(k\)  samples in the sample set, we can use the assumption \(k \ge {\epsilon _{b}^{2}/((1-\gamma )^2\epsilon _{s}^2)}\)  and have
\({\begin{array}{c}\hat{Q}^{i}_{sat}(s,a)-B^{\pi ^{Q}}\hat{Q}^{i}_{sat}(s,a) < \epsilon _{a}+(n_{c}^{i}+\mu _{c}^{i})(1+3\gamma )+2(1-\gamma )\epsilon _{s}\end{array}}\) 

Now, we know that even though the policy \(\widetilde{\pi }_{i}\)  is non-stationary, it is comprised of stationary segments. Using the definitions from Lemma REF , starting from step t for agent \(i\) , \(\widetilde{\pi }_{i}\)  is equal to \(\pi ^{\hat{Q}^{i}_{sat}}\)  for at least \(\tau (t)\)  steps. Combining the results from Lemma REF  with the bounds we have listed above, we have that with a probability larger than \(1-\delta \) , for all time steps and all agents:
\({\begin{array}{c}V^{\widetilde{\pi }_{i}}(s_{t,i}) \ge V^{*}(s_{t,i}) - \frac{2\epsilon _{a}+2(1+3\gamma )(n_{c}^{i}+\mu _{c}^{i})}{1-\gamma } -3\epsilon _{s} -\epsilon _{e}^{i}(t)\end{array}}\) 
with
\(\begin{aligned}\epsilon _{e}^{i}(t) = &\gamma ^{\tau (t)}Q_{\max }+2\sum _{h \in H}\left(hp^{e}_{h,k_{m}}(s_{t,i})\right)Q_{\max }+\\&\sum _{k_{a} \in \lbrace K_{a}-k_{m}\rbrace }2\sum _{h \in H}\left(hp^{e}_{h,k_{a}}(s_{t,i})\right)2\sqrt{\frac{8k_{m}(\sigma ^{2}+\sigma ^{2}_{R})}{k_{a}}}\end{aligned}\) 
We now calculate a bound over the \(TCE\) . With a probability larger than \(1-\delta \)  we have
\(\begin{aligned}\sum _{i=1}^{N}\sum _{t=0}^{\infty }\epsilon _{e}^{i}(t) = &\sum _{i=1}^{N}\sum _{t=0}^{\infty }\left(\gamma ^{\tau (t)}Q_{\max }+2\sum _{h \in H}\left(hp^{e}_{h,k_{m}}(s_{t,i})\right)Q_{\max }\right) \\ &+\sum _{i=1}^{N}\sum _{t=0}^{\infty }\sum _{k_{a} \in \lbrace K_{a}-k_{m}\rbrace }2\sum _{h \in H}\left(hp^{e}_{h,k_{a}}(s_{t,i})\right)2\sqrt{\frac{8k_{m}(\sigma ^{2}+\sigma ^{2}_{R})}{k_{a}}} \\&\le \frac{NSAQ_{max}{1+\log _{2}\frac{k}{k_{m}}}}{1-\gamma }+\frac{2Q_{max}SA(k_{m}+N)(1+\log _{2}T_{H})T_{H}}{1-\sqrt{\frac{2T_{H}}{SA(k_{m}+N)}\ln \frac{2|K_{a}|}{\delta }}}\\&+2\sum _{k_{a} \in \lbrace K_{a}-k_{m}\rbrace }\frac{2\sqrt{\frac{8k_{m}(\sigma ^{2}+\sigma ^{2}_{R})}{k_{a}}}SA(k_{a}+N)(1+\log _{2}T_{H})T_{H}}{1-\sqrt{\frac{2T_{H}}{SA(k_{a}+N)}\ln \frac{2|K_{a}|}{\delta }}}\end{aligned}\) 
Here we used the fact there are at most \({1+\log _{2}\frac{k}{k_{m}}}SA\)  policy changes for the first term, and substituted the result from lemma REF  for the second and third term. Continuing
\(\begin{aligned}\sum _{i=1}^{N}\sum _{t=0}^{\infty }\epsilon _{e}^{i}(t) <&\frac{Q_{max}SA(2k_{m}+N{3+\log _{2}\frac{k}{k_{m}}})(1+\log _{2}T_{H})T_{H}}{1-\sqrt{\frac{2T_{H}}{SA(k_{m}+N)}\ln \frac{2|K_{a}|}{\delta }}}\\&+\frac{2\sum _{k_{a} \in \lbrace K_{a}-k_{m}\rbrace }2\sqrt{\frac{8k_{m}(\sigma ^{2}+\sigma ^{2}_{R})}{k_{a}}}SA(k_{a}+N)(1+\log _{2}T_{H})T_{H}}{1-\sqrt{\frac{2T_{H}}{SA(k_{m}+N)}\ln \frac{2|K_{a}|}{\delta }}} \\&=\frac{Q_{max}SA(2k_{m}+N{3+\log _{2}\frac{k}{k_{m}}})(1+\log _{2}T_{H})T_{H}}{1-\sqrt{\frac{2T_{H}}{SA(k_{m}+N)}\ln \frac{2|K_{a}|}{\delta }}}\\&+\frac{2\sum _{k_{a} \in \lbrace K_{a}-k_{m}\rbrace }2\sqrt{8k_{m}(\sigma ^{2}+\sigma ^{2}_{R})}SA(\sqrt{k_{a}}+\frac{N}{\sqrt{k_{a}}})(1+\log _{2}T_{H})T_{H}}{1-\sqrt{\frac{2T_{H}}{SA(k_{m}+N)}\ln \frac{2|K_{a}|}{\delta }}} \\&\le \frac{Q_{max}SA(2k_{m}+N{3+\log _{2}\frac{k}{k_{m}}})(1+\log _{2}T_{H})T_{H}}{1-\sqrt{\frac{2T_{H}}{SA(k_{m}+N)}\ln \frac{2|K_{a}|}{\delta }}}\\&+\frac{\sqrt{8k_{m}(\sigma ^{2}+\sigma ^{2}_{R})}SA(1+\log _{2}T_{H})T_{H}}{1-\sqrt{\frac{2T_{H}}{SA(k_{m}+N)}\ln \frac{2|K_{a}|}{\delta }}}\\&\times \left(4\sqrt{k}\left(1+\sum _{l=0}^{\infty }\left(\frac{1}{2^{l}}+\frac{1}{2^{l}\sqrt{2}}\right)\right)+4N\left(\sum _{l=1}^{\infty }\left(\frac{1}{2^{l}}+\frac{1}{2^{l}\sqrt{2}}\right)\right)\right)\end{aligned}\) 
Here we used the fact that the number of possible values for \(k_{a}\)  is bounded. By bounding the sums in the last equation we get to a simpler expression
\(\begin{aligned}\sum _{i=1}^{N}\sum _{t=0}^{\infty }\epsilon _{e}^{i}(t) <&\frac{Q_{max}SA(2k_{m}+N{3+\log _{2}\frac{k}{k_{m}}})(1+\log _{2}T_{H})T_{H}}{1-\sqrt{\frac{2T_{H}}{SA(k_{m}+N)}\ln \frac{2|K_{a}|}{\delta }}}\\&+\frac{\sqrt{8k_{m}(\sigma ^{2}+\sigma ^{2}_{R})}SA(1+\log _{2}T_{H})T_{H}}{1-\sqrt{\frac{2T_{H}}{SA(k_{m}+N)}\ln \frac{2|K_{a}|}{\delta }}}(18\sqrt{k}+10N) \\&=\frac{SA(1+\log _{2}T_{H})T_{H}}{1-\sqrt{\frac{2T_{H}}{SA(k_{m}+N)}\ln \frac{2|K_{a}|}{\delta }}}\\&\times \left(Q_{max}(2k_{m}+N{3+\log _{2}\frac{k}{k_{m}}})+(18\sqrt{k}+10N)\sqrt{8k_{m}(\sigma ^{2}+\sigma ^{2}_{R})}\right)\end{aligned}\) 
We now assume that for a state-action space of realistic size, the denominator is approximately 1 and write
\(\begin{aligned}\sum _{i=1}^{N}\sum _{t=0}^{\infty }\epsilon _{e}^{i}(t)&< \left(SA(1+\log _{2}T_{H}){\frac{1}{1-\gamma }\ln \frac{Q_{max}}{\epsilon _{s}}}\right)\\&\times \left(Q_{max}(2k_{m}+N{3+\log _{2}\frac{k}{k_{m}}})+(18\sqrt{k}+10N)\sqrt{8k_{m}(\sigma ^{2}+\sigma ^{2}_{R})}\right)\end{aligned}\) 
In order to get to a big \(O\)  notation, we write the various terms in the previous equation explicitly and ignore logarithmic terms. Note that \(k_{m}\)  is assume to be a logarithm term, and we can assume that \(k = C{\epsilon _{b}^{2}/((1-\gamma )^2\epsilon _{s}^2)}\)  for some constant \(C\)  (this does not violate assumption ). Therefore
\(\begin{aligned}\sum _{i=1}^{N}\sum _{t=0}^{\infty }\epsilon _{e}^{i}(t) = \widetilde{O}\left(\left(N\left(Q_{\max }+\sqrt{\sigma ^{2}+\sigma ^{2}_{R}}\right)+\frac{\sigma ^{2}+\sigma ^{2}_{R}}{\epsilon _{s}(1-\gamma )}\right)SA/(1-\gamma )\right)\end{aligned}\) 

Using Theorem REF , the following theorem abandons the concentration bound assumption REF  over the learner-agent noise terms, and instead assumes the noise to be sub-Gaussian with a bounded mean and bounded parameter. We note that we can actually prove such a theorem for the more general case of noise with a bounded mean and variance. However, in such a case we replace a logarithmic term with a non-logarithmic one which worsens the bound. This is to be expected, as we have less information regarding the tails of a general distribution.
*

From the definition of the noise in (REF ) and from assumption , we have that \(b^{i}(s,a)\)  is subgaussian with
\({\begin{array}{c}P(|b^{i}(\overline{s},\overline{a})|\ge n_{c}^{i}) \le 2e^{-\frac{n^{2}_{c}}{2\sigma ^{2}_{c}}}\end{array}}\) 
Setting
\({\begin{array}{c}2e^{-\frac{n^{2}_{c}}{2\sigma ^{2}_{c}}} = \frac{\delta }{12N{1+\log _{2}\frac{k}{k_{m}}}(SA)^{3}}\end{array}}\) 
results in \(n_{c}^{i}=\sigma ^{i}_{c}\sqrt{2\ln \left(\frac{24N{1+\log _{2}\frac{k}{k_{m}}}(SA)^{3}}{\delta }\right)}\) . Now given the event that \(|b^{i}(\overline{s},\overline{a})|\le n_{c}^{i}\) , we know that \(|b^{i}(\overline{s},\overline{a})+m^{i}(\overline{s},\overline{a})|\le n_{c}^{i}+\Delta Q^{i}_{c}\) . Therefore
\(\begin{aligned}P\left\lbrace |b^{i}(\overline{s},\overline{a})+m^{i}(\overline{s},\overline{a})|\ge n_{c}^{i}+\Delta Q^{i}_{c}\right\rbrace  \le P(|b^{i}(\overline{s},\overline{a})|\ge n_{c}^{i}) \le \frac{\delta }{12N{1+\log _{2}\frac{k}{k_{m}}}(SA)^{3}}\end{aligned}\) 
so that the total equivalent noise \(n=b+m\)  is of mean 0 and has the right probability bound. Substituting this back in the bound for Theorem REF  concludes our proof.

Compared to Theorem REF , Theorem  requires information about nothing but the first and second moments of the noise, enables a bounded noise term \(m^{i}(s,a)\)  to be present, and suggests that the sample complexity bound is proportional to a bound on the the sub-Gaussian parameter of the noise, and on the bound on \(m^{i}(s,a)\) .

Estimation methods proofs and properties
In this section, we provide proofs and further information about the optimal linear estimation theorems from Section . The proofs all consist of finding the parameter of the equivalent sub-Gaussian noise resulting from the weighted sum (REF ), and minimizing it with regard to the estimation weights. For simplicity and without loss of generalization, we will present the proofs for agent \(N\)  as the receiving agent, and denote \(\sigma ^{N}_{L} \triangleq \sigma _{L},\sigma ^{j,N}_{A} \triangleq \sigma _{A,j}, \sigma _{c}^{N} \triangleq \sigma _{c}\) .
General additive noise model
*

We know that a linear sum of independent sub-Gaussian random variables is a sub-Gaussian with the sum of their collective parameters. Furthermore, since the weights all sum up to 1, by (REF ),(REF ) (REF ) we can write
\(\begin{aligned}\hat{Q}^{N}(s,a) &= Q(s,a)+n^{N}(s,a)\\n^{N}(s,a) &= w_{NN}n^{N}_{L}(s,a) + \sum _{j=1}^{d_{N}}w_{jN}n^{j,N}_{A}(s,a)\end{aligned}\) 
with \(En^{N}(s,a)=0\)  and parameter
\(\begin{aligned}\sigma _{c}^{2}(\mathbf {w}) = \sum _{j=1}^{d_{N}}w_{jN}^{2}\left(\sigma _{L}^{2}+\sigma _{A,j}^{2}\right)+\left(1-\sum _{j=1}^{d_{N}}w_{jN}\right)^{2}\sigma _{L}^{2}~.\end{aligned}\) 
This means the sample complexity bound of Theorem REF  holds. Since the term \(\sigma _{c}(\mathbf {w})\)  is the only one depending on the weights in the above-mentioned bound, minimizing this term over the weights results in optimizing the sample complexity bound.
Calculating the derivative, we get that
\(\frac{\partial \sigma _{c}^{2}(\mathbf {w})}{\partial \mathbf {w}} = A\mathbf {w}-\sigma _{L}~,\) 
meaning that the original function can be written in the following form
\(\sigma _{c}^{2}(\mathbf {w}) = \frac{1}{2}\left(\mathbf {w}-A^{-1}\sigma _{L}\right)^{T}A\left(\mathbf {w}-A^{-1}\sigma _{L}\right)+C\) 
where C is some constant. We can use the Matrix determinant lemma from [10]} (where \(X\)  is a matrix and \(\mathbf {u},\mathbf {v}\)  are some vectors) \(\mathrm {det}\left(X+\mathbf {u}\mathbf {v}^{T}\right)= \left(1+\mathbf {v}^{T}X^{-1}\mathbf {u}\right)\mathrm {det}(X)\)  to show that the principle minors of A (which are all of the same form) are all positive, meaning that A is positive definite by Sylvester's criterion. This means that \(\sigma _{c}(\mathbf {w})\)  has a single global minimum at
\(\mathbf {w}^{*}=A^{-1}\sigma _{L}~.\) 
Now, by substituting \(\sigma _{c}^{2}(\mathbf {0})=\sigma _{L}^{2}=\frac{1}{2}\left(A^{-1}\sigma _{L}\right)^{T}A\left(A^{-1}\sigma _{L}\right)+C\)  we can easily see that
\(\sigma _{c}^{*} \triangleq \sigma _{c}(\mathbf {w}^{*}) = \sigma _{L}\sqrt{w_{NN}^{*}}~.\) 
Properties of the optimal solution
We list some properties exhibited by the general optimal solution.

We have that \(\sigma _{c}(\mathbf {w}) \ge 0 \)  by definition, and also that \(\sigma _{c}(\mathbf {w}^{*}) \le \sigma _{c}(\mathbf {0})=\sigma _{L}^{2}\) , and therefore \(0 \le w_{NN}^{*} \le 1\) . This also means that the optimal solution ensures that we get a lower sample complexity compared to the case without communication, where the variance is \(\sigma _{L}^{2}\) .

The optimal solution results in \(\sigma _{c}^{*}\)  smaller than the case of a uniform average.

We can furthermore show that
\(\forall j \in \lbrace 1,...,d_{N}\rbrace : w_{jN}^{*} \in [0,1]\) 
which is intuitive.

From the mathematical form of \(\sigma _{c}(\mathbf {w})\) , it is easy to deduce that the larger the variance \(\sigma _{A,j}^{2}\)  is, the smaller the corresponding optimal weight \(w_{jN}^{*}\)  gets.

In the case where \(\sigma _{A,j}=\sigma _{A,l}\)  for some \(j\)  and \(l\) , we have from the symmetry property of \(\sigma _{c}(\mathbf {w})\)  that \(w_{jN}^{*}=w_{lN}^{*}\) .

Let us now also consider two extreme cases and see whether our solution corresponds to intuition.

In the case where \(\sigma _{A,j} \gg \sigma _{L}\)  for some \(j \in \lbrace 1,...,d_{N}\rbrace \) , it is evident from both intuition and the explicit solution that \(w_{jN}^{*} \ll 1\) .

In the case where \(\sigma _{L} \gg \sigma _{A,j}\)  for all \(j \in \lbrace 1,...,d_{N}\rbrace \) , we have that
\(\sigma _{c}(\mathbf {w}) = \sum _{j=1}^{d_{N}}w_{jN}^{2}\sigma _{L}^{2}+\left(1-\sum _{j=1}^{d_{N}}w_{jN}\right)^{2}\sigma _{L}^{2}\) 
which is a case where all the weights are identical. By intuition and by substituting in the explicit solution, we have that \(w_{jN}^{*}=1/(d_{N}+1)\)  for all \(j \in \lbrace 1,...,d_{N}\rbrace \) .


Special examples
We exemplify the properties introduced previously over simple cases.

First, let us consider the case of three fully-connected agents. By calculating the optimal solution we can learn about the effects difference noise variances have on the solution. We can see in Figure REF  that, as expected, the self optimal weight \(w^{*}_{33}\)  (proportional to \(\sigma _{c}^{*}\)  itself) decreases as the inter-agent noise parameters \(\frac{\sigma _{A,1}^{2}}{\sigma _{L}^{2}},\frac{\sigma _{A,2}^{2}}{\sigma _{L}^{2}}\)  decrease, meaning that it is beneficial to not only use the less noisy \(Q\) -function received from the learner. The optimal weight for the case of no inter-agent noise is \(w^{*}_{33}=\frac{1}{3}\) , but the larger the noise between the agents is, the better it is to rely more on the less-noisy \(Q\) -function. We can also see that the larger the relation \(\frac{\sigma _{A,1}^{2}}{\sigma _{A,2}^{2}}\)  gets, the smaller \(\frac{w_{13}^{*}}{w_{23}^{*}}\)  becomes, as explained earlier.

In order to describe the dependence of the optimal weights on the number of agents, we examine a case in which there are 2 groups of agents of size \(N_{A}\) , the first one contains agents whose noise parameters are all \(\sigma _{A,1}^{2}\) , and the second one with noise parameters all equal to \(\sigma _{A,2}^{2}\)  (all the agents are fully connected). This is an extension of the previous case in which \(N_{A}=1\) . Besides these two groups, agent number \(N\)  on which we are focusing is also present, so there are overall \(N=2N_{A}+1\)  agents.
In Figure REF A, we can see that the more agents there are - the smaller \(w_{33}^{*}\)  becomes, meaning that it is beneficial to rely on more agents. Since we keep \(\sigma _{A,2}=\sigma _{L}\)  constant in this case, we see that as the ratio \(\frac{\sigma _{A,1}^{2}}{\sigma _{A,2}^{2}}\)  increases, \(w_{33}^{*}\)  increases as well as explained before. In Figure REF B we can see that the ratio \(\frac{w^{*}_{13}}{w^{*}_{23}}\)  remains constant as long as \(\frac{\sigma _{A,1}^{2}}{\sigma _{A,2}^{2}}\)  does not change, even for a varying number of agents.

<FIGURE><FIGURE>
Same parameters additive noise model
We present a special case of equal noise parameters. Theorem REF  is a direct consequence of Theorem REF .
*

By substituting \(\forall i,j: \sigma ^{j}_{L}=\sigma _{L},\sigma ^{j,i}_{A}=\sigma _{A}\)  for the optimal weights of Theorem REF , we get
\({\begin{array}{c}w_{ii}^{*} = \frac{\sigma ^{2}_{L}+\sigma ^{2}_{A}}{(d+1)\sigma ^{2}_{L}+\sigma ^{2}_{A}}=\frac{1}{1+\frac{d}{1+\sigma ^{2}_{A}/\sigma ^{2}_{L}}} \quad ; \quad w^{*}_{ji} =\frac{1-w_{ii}^{*}}{d} \quad \forall j \in \lbrace 1,...,d_{i}\rbrace \end{array}}\) 
where weights of \(Q\) -functions with the same noise parameter all have the same weight from symmetry.
For the second part of the proof, by substituting \(\forall j \in \lbrace 1,...,d_{i}\rbrace :w_{ji}=\frac{1}{d+1}\)  and \(w_{ii}=\frac{1}{d+1}\)  in \(\sigma _{c}(\mathbf {w})\)  from the proof of Theorem REF , we get
\({\begin{array}{c}\widetilde{\sigma _{c}}^{2} = \frac{\sigma _{L}^{2}}{d+1}+\frac{d}{(d+1)^{2}}\sigma _{A}^{2}\end{array}}\) 
And we can see that using a uniform average is preferable to using the less noisy \(Q\) -function from the learner alone, if and only if \(\widetilde{\sigma _{c}}^{2} \le \sigma ^{2}_{L}\)  which leads to \(d+1 \ge \frac{\sigma ^{2}_{A}}{\sigma ^{2}_{L}}\) .

Quantization noise model
We now prove Theorem REF .
*

Without the loss of generality, we will again prove the theorem for agent \(N\) .
First part:
We partition the noise in the weighted \(Q\) -function \(\hat{Q}^{N}\)  to two terms:
\(\begin{aligned}\hat{Q}^{N}(s,a) &= Q(s,a)+b^{N}(s,a)+m^{N}(s,a)\\b^{N}(s,a) &\triangleq w_{NN}b_{L}^{N}(s,a) + \sum _{j = 1}^{d_{N}}w_{jN}(b_{L}^{j}+b^{j,N}_{A}(s,a))\\m^{N}(s,a) &\triangleq w_{NN}m_{L}^{N}(s,a) + \sum _{j = 1}^{d_{N}}w_{jN}(m_{L}^{j}+m^{j,N}_{A}(s,a))\end{aligned}\) 
We know that a linear sum of independent sub-Gaussian random variables is a sub-Gaussian with the sum of their collective parameters. Therefore, \(b^{N}(s,a)\)  is a mean 0 sub-Gaussian with some parameter \(\sigma _{c}^{2}\)  as in Theorem REF . Regarding the quantization noise, we have
\({\begin{array}{c}|m^{N}(s,a)| \le \Delta Q^{N}|w_{NN}|+ \sum _{j = 1}^{d_{N}}\left(\Delta Q^{j}+\Delta Q^{j,N}\right)|w_{jN}| \triangleq \Delta Q_{c}\end{array}}\) 
Therefore, the conditions of Theorem REF  hold.

Second Part:
In the case of equal noise properties, due to the symmetry of the problem, we know that the optimal weights will be such that the \(Q\) -function that agent \(N\)  receives from the learner has some weight \(w \triangleq w_{NN}\) , and the rest have weights of identical values \(\frac{1-w}{d}\) . Therefore, \(b^{N}(s,a)\)  is a mean 0 sub-Gaussian with parameter \(\sigma _{c}^{2} \le w^{2}\sigma _{L}^{2}+\frac{(1-w)^{2}}{d}(\sigma _{L}^{2}+\sigma _{A}^{2})=\sigma _{L}^{2}\left(w^{2}+2\frac{(1-w)^{2}}{d}\right)\) . And regarding the quantization noise:
\({\begin{array}{c}|m^{N}(s,a)| \le \Delta Q \left(|w|+\sum _{j = 1}^{d}|\frac{1-w}{d}|(1+1)\right)=\Delta Q \left(|w|+2|1-w|\right)\end{array}}\) 
In order to find the optimal weight \(w^{*}\)  that minimizes the sample complexity bound, define
\({\begin{array}{c}g(w) \triangleq \sigma _{L}f\sqrt{w^{2}+2\frac{(1-w)^{2}}{d}} +\Delta Q \left(|w|+2|1-w|\right)\end{array}}\) 
Where the function \(f\)  is defined in Theorem REF . From the sample complexity bound
\({\begin{array}{c}V^{\widetilde{\pi }_{i}}(s_{t,i}) \ge V^{*}(s_{t,i}) - \frac{2\epsilon _{a}+2(1+3\gamma )g(w)}{1-\gamma }-3\epsilon _{s} -\epsilon _{e}^{i}(t)\end{array}}\) 
It is evident that the function \(g(w)\)  is the only term depending on the weight \(w\) , such that minimizing \(g(w)\)  over \(w\)  will lead to maximization of the bound.
Viewing \(g(w)\) , we can learn a few properties regarding the optimal solution. Notice that the first term is a parabola-like function with a minimum value at \(w=\frac{2}{d+2}\)  as we have seen in Theorem REF , and the second term has a few cases of discontinuous change.

For \(w<\frac{2}{d+2}\)  we have a sum of two decreasing function, therefore \(g(w)\)  itself decreases.

For \(w>1\) , we have a sum of two increasing functions - and therefore \(g(w)\)  is increasing.

Therefore, there is a single minimal value of \(g(w)\)  in the interval \(w \in \left[\frac{2}{d+2},1\right]\) .

By calculating the derivative of \(f\)  in \([0,1]\)  and finding the value of \(w\)  for which it vanishes, we have that there is a single solution in this interval, satisfying
\({\begin{array}{c}w=\frac{2}{d+2}\left(1+\frac{d}{\sqrt{(d+2)(2\sigma _{L}^{2}f^{2}/\Delta Q)^{2}-2d}}\right)\end{array}}\) 
Under the condition that \(\frac{f\sigma _{L}}{\Delta Q} \ge \frac{d}{\sqrt{d+2}}\) . Otherwise, there is no minimal value in this interval. Furthermore, we can see that the expression for the solution is decreasing with \(f\) , and that substituting \(f\sigma _{L}= \Delta Q\)  results in the minimal value being \(w=1\) .
Therefore, we conclude that:

For \(f\sigma _{L}\le \Delta Q\)  the optimal solution is \(w^{*}=1\)  (either the minimal value is outside of \([0,1]\)  and thus contradicts our initial assumption, or there is no minimal value at all and \(g(w)\)  is decreasing in \([0,1]\) ).

For \(f\sigma _{L}\ge \Delta Q\) , the minimal value of \(g(w)\)  is achieved at \(w\)  from (REF ), which is indeed inside the interval \(\left[\frac{2}{2+d},1\right]\) , and gets closer to 0 as \(d\)  grows larger.


Computational results
In this section we give more details about the computational results shown in Figure REF .
In this simulation, 4 agents are randomly initiated in an unknown warp-around \(5\times 5\)  grid world, and get a reward of 1 whenever they land on the top-right corner. All other grid locations contain a reward of 0. When an agent reaches the top-right corner, it is re-initiated at a random location in the grid world. \(Q\) -functions are noisy with an additive Gaussian noise, under the conditions of Theorem REF , and we assume no reward noise. The agents communicate via a fully-connected communication network \(\Gamma \) . The results are shown for 150 parallel experiments, 10 episodes and 50 steps per episode, where at the beginning of an episode - each agent is re-initiated at a random location in the grid world. We calculate the average accumulated reward and standard deviation for an episode per agent, where the average is calculated over parallel experiments and over all agents (continuous lines). We also use a simulated agent that can move greedily in the grid world using the exact non-noisy \(Q\) -function (dashed lines). This is done to separate the effects of the noise from those of exploration. In Figure REF , the learner-agent communication noise variance is fixed at \(\sigma ^{2}_{L}=0.1\) , and each one of the figures A,B,C represents a different agent-agent noise variance \(\sigma ^{2}_{A}\) , in order to show the tradeoff between the \(\sigma ^{2}_{A}/\sigma ^{2}_{L}\)  ratio and the estimation weighting scheme. The red lines represent a weight \(w=1\)  for the \(Q\) -function received from the learner, meaning there is no averaging. The blue lines represent a uniform averaging \(w=0.25\) , and the black lines represent the optimal weighting scheme suggested by Theorem REF . We note that for an algorithm to be able to perform these estimation method, prior knowledge over the noise parameters is needed. Given such knowledge, each agent can decide on its optimal weights before the exploration begins and use them throughout the whole loop. The green line represents the more realistic case where the algorithm has no information about the parameters of the environmental noise, and estimates the variance on the go in order to use it at the optimal weighting scheme of Theorem REF . We use an estimation scheme similar to that suggested by [11]}, such that for agent \(i\)
\(\begin{aligned}\hat{w}_{ii}(t) &= \frac{(\hat{\sigma }^{i}_{LA}(t-1))^{2}}{(\hat{\sigma }^{i}_{LA}(t-1))^{2}+(\hat{\sigma }^{i}_{L}(t-1))^{2}} \\(\hat{\sigma }^{i}_{L}(t))^{2} &= (\hat{\sigma }^{i}_{L}(t-1))^{2}+\frac{\frac{1}{SA-1}\sum _{s,a}\left(\widetilde{Q}^{i}_{L}(s,a)-Q(s,a)\right)^{2}-(\hat{\sigma }^{i}_{L}(t-1))^{2}}{t}\\(\hat{\sigma }^{i}_{LA}(t))^{2} &= (\hat{\sigma }^{i}_{LA}(t-1))^{2}+\frac{\frac{1}{SA(N-1)-1}\sum _{s,a,N}\left(\widetilde{Q}^{j,i}_{A}(s,a)-Q(s,a)\right)^{2}-(\hat{\sigma }^{i}_{LA}(t-1))^{2}}{t}\end{aligned}\) 
\((\hat{\sigma }^{i}_{L}(t))^{2}\)  represents the estimation of the variance of the learner-agent noise (\(\sigma ^{2}_{L}\) ) at step \(t\) , and \((\hat{\sigma }^{i}_{LA}(t))^{2}\)  represents the estimation of the agent-agent communication noise variance (\(\sigma ^{2}_{A}+\sigma ^{2}_{L}\) ) at step \(t\) , such that the formula for \(\hat{w}_{ii}(t)\)  is an adaptive version of the optimal weight. The estimation of \((\hat{\sigma }^{i}_{L}(t))^{2}\)  in each step is done by using a weighted average of the current estimate with the past one. The estimate given the current \(Q\) -function is an empiric average over all values in \(\widetilde{Q}^{i}_{L}\) . \((\hat{\sigma }^{i}_{LA}(t))^{2}\)  is calculated similarly, but uses all of the \(N-1\)  values in \(\lbrace \widetilde{Q}^{j,i}_{A}\rbrace \) .
<FIGURE>We use a version of Algorithm  that accumulates samples for each state-action instead of replacing them with new ones whenever the number of active samples increases. We use \(\epsilon _{a}=10^{-7}, \epsilon _{b}=0.1, k=9, k_{m}=3, \gamma =0.98,Q_{\max }=\frac{R_{max}}{1-\gamma }\)  and 30 Bellman iterations at most each time we perform them.
MATLAB Code is available at https://github.com/ravehor92/ALT2020-Raveh, and contains the following sub-codes, while the dependencies are shown in Figure REF .

main.m: The main code generating Figure REF .

Noisy_PAC.m: performs a parallel experiments of agents exploring in a grid world, for the same parameters.

VI.m: performs Value Iterations over a given Q function.

move:.m given a state and actions, progress to the next state and receive a reward.

Active_Samples.m: Number of active samples for a given state-action (samples that are currently being used for value iteration).

shadedErrorBar.m: Creates a shaded error bar.

