Computer Science
Computer Science-2021.gz:143223
Computer Science
1827b13a3e7d78b0a5fabe771bfb9d8b40fa57f2
Lightweight Distributed Gaussian Process Regression for Online Machine Learning

In this paper, we study the problem where a group of agents aim to collaboratively learn a common static latent function through streaming data. We propose a lightweight distributed Gaussian process regression (GPR) algorithm that is cognizant of agents' limited capabilities in communication, computation and memory. Each agent independently runs agent-based GPR using local streaming data to predict test points of interest; then the agents collaboratively execute distributed GPR to obtain global predictions over a common sparse set of test points; finally, each agent fuses results from distributed GPR with agent-based GPR to refine its predictions. By quantifying the transient and steady-state performances in predictive variance and error, we show that limited inter-agent communication improves learning performances in the sense of Pareto. Monte Carlo simulation is conducted to evaluate the developed algorithm.


I. INTRODUCTION
Networks of agents can access large amount of streaming data online in many applications, e.g., emergency response [1], autonomous driving [2], and precision agriculture [3]. Machine learning has been increasingly adopted to extract reliable and actionable information from big data and enable agents to adapt and react in uncertain and dynamically changing environments. Nowadays machine learning algorithms have achieved remarkable performances in terms of consistency [4], generalization [5] and robustness [6]. However, limited resources challenge implementation of the algorithms on physical agents, e.g., mobile robots.
Gaussian process regression (GPR) [7] is an efficient nonparametric statistical learning model. GPR models a target function as a sample from a Gaussian process prior specified by a pair of mean and covariance functions only depending on inputs. With proper choice of prior covariance function, also referred as kernel, and mild assumptions on the target function, GPR is able to consistently approximate any continuous function [8]. With optimal sampling in the input space and covariance functions obeying Sacks-Ylvisaker conditions of order r, the generalization error of GPR diminishes at the rate of O(n −(2r+1)/(2r+2) s ) (Section V.2, [9]), where n s is the number of samples. GPR is able to quantify learning uncertainties since it predicts the target function in the form of a posterior probability distribution. GPR has demonstrated powerful capabilities in various applications, e.g., optimization [10], motion planning [11], and trajectory estimation [12].
GPR scales as O(n 3 s ) in computational complexity and O(n 2 s ) in memory (page 171, [7]), which prohibits applications Zhenyuan  with large datasets. There are multiple sparse approximation methods for large datasets. One major class of approximation methods, which is also referred as global approximation, tackles the computational complexity by achieving the sparsity of the Gram matrix. Methods include using a subset of data to approximate the whole training dataset, designing a sparse kernel, and sparsifying the Gram matrix. The best possible result can be achieved by global approximation algorithms is O(m 3 s ) in computational complexity and O(m 2 s ) in memory, where m s n s is the number of inducing points or the size of a subset of training data. More details about global approximation can be found in the recent survey paper [13]. In the community of Geostatistics, Nearest-neighbor GPR [14] is applied [15] [16], where the predictions are made only using the training data of the nearest input. It requires only O(n s ) in both memory and (worst-case) computation.
Centralized implementation of GPR is not suitable for networks of agents due to poor scalability in data size, high cost in communication and memory, and fragility to singlepoint failures. There have been studies on distributed GPR over server-client architecture, which is also referred to as divide-and-conquer approach or local approximations [13]. In the server-client architecture, a server acts as the centralized entity that partitions a dataset and assigns each subset of the data to computing units (clients). The clients perform training independently and send their learning results to the server for post-processing. These methods speed up the training process and are able to scale to arbitrarily large datasets. Communication budget constraint is considered in [17] by reducing the dimensionality of transmitted data to approximate the whole dataset. Sparse approximation of full GPR is used in [17] to further relieve the communication overhead. Notice that the server-client architecture requires each client being well-connected with the server, and is not robust to the failure of the server. Paper [18] decentralizes sparse approximations of full GPR for fixed datasets over complete communication graphs. A distributed algorithm is also proposed to deal with fixed and sparse graphs. However, this paper considers offline learning with static datasets on the agents and does not provide theoretic guarantee on the distributed algorithm.
Our work is related to multi-agent regression using kernel methods and basis functions. Papers [19] [20] [21] study offline learning, where all training data is provided before learning, using kernel methods. Papers [22][23] [24] [25] study online learning, where training data is collected successively by mobile robots, using basis functions. In particular, they approximate the unknown functions with a linear combination of a finite number of known basis functions. This reduces arXiv:2105.04738v4 [cs.LG] 2 Dec 2021 the problem into a parameter estimation problem. From the perspective of regression, the problem investigated in [22] is equivalent to selecting the centers of a finite number of basis functions defined by Voronoi partition. In contrast, this paper considers online learning of abstract agents using Gaussian processes, where the unknown function is modelled as a sample from a distribution of functions.
Contribution statement. We consider the problem where a group of agents aim to collaboratively learn a common static latent function through streaming data. We propose Lightweight Distributed Gaussian Process Regression (LiDGPR) algorithm for the agents to solve the problem. More specifically, each agent independently runs agent-based GPR using local streaming data to predict the test points of interest; then the agents collaboratively execute distributed GPR to obtain global predictions over a common sparse set of test points; finally, each agent fuses the results from distributed GPR with those from agent-based GPR to refine its predictions. Our analysis of the transient and steady-state performances in predictive variance and error reveals that through communication agents whose data samples have lower dispersion (or observation noise has lower variance) help improve the performance of the agents whose data samples have higher dispersion (or observation noise has higher variance). The improvements in learning performances are in the sense of Pareto, i.e., some agents' performances improve without sacrificing other agents' performances. In summary, our major contributions are two-fold: • We develop LiDGPR that is cognizant of agents' limited capabilities in communication, computation and memory. • We analyze the predictive mean and variance of LiDGPR and quantify the improvements of the agents' learning performances resulted from inter-agent communication.
Monte Carlo simulation is conducted to evaluate the developed algorithm. In addition to the preliminary version of this paper [26], this version includes a new set of theoretical results on the steady-state performance of predictive errors, identifies the factors that affect the improvement of learning performances, provides the proofs of the theoretical results, and discusses how the algorithm is cognizant of limited resources. Notations: We use lower-case letters, e.g., a, to denote scalars, bold letters, e.g., a, to denote vectors; we use uppercase letters, e.g., A, to denote matrices, calligraphic letters, e.g., A, to denote sets, and bold calligraphic letters, e.g., A, to denote spaces. For any vector a, we use a i to denote the i-th entry of a. For any matrix A, we denote a ij as the entry at i-th row j-th column. Denote I n ∈ R n×n the n-by-n-dimensional identity matrix, 1 n ∈ R n the n-dimensional column vector with all 1's, i.e., [1, · · · , 1] T , and 0 n analogously.
We use superscript (·) [i] to distinguish the local values of agent i, and (·) max ((·) min ) denote the maximum (minimum) of the local values, e.g., a max max i∈V a [i] . We denote superscript (·) T the transpose of a vector or matrix, bracket [·] E the column vector with elements satisfying event E. Denote E a [·] the expectation taken over the distribution of random variable a, and P {·} a distribution. We use O(·) to denote the conventional Big O notation, i.e., O(g(t)) represents the limiting behavior of some function f (t) if lim t→∞ f (t) g(t) = a for some constant a > 0.
We use ( ) to denote element-wise comparison between two vectors, i.e., for any a, b ∈ R n , a ( )b if and only if a i ( )b i for all i = 1, · · · , n. Operation |a| takes the absolute values element-wise on vector a, |A| returns the cardinality of set A, a ∞ max j |a j | for any vector a. Define the distance metric ρ(z, z ) z − z , the point to set distance as ρ(z, Z) inf z ∈Z ρ(z, z ). Define proj(z, Z) {z ∈ Z|ρ(z, z ) = ρ(z, Z)} the projection set of point z onto set Z. Denote the supremum of a function η as η Z sup z∈Z |η(z)|.
II. PROBLEM STATEMENT Network model. Consider a network of agents represented by a directed time-varying communication graph G(t) (V, E(t)), where V {1, · · · , n} represents the agent set, and E(t) ⊆ V × V denotes the edge set at time t. Notice that (i, j) ∈ E(t) if and only if agent i can receive messages from agent j at time t. Define the set of the neighbors of agent i We make the following standard assumptions [27] about the network topology: Assumption II.1. (Periodical Strong Connectivity). There exists positive integer b 1 such that, for all time instant t 0, the directed graph (V, E(t)∪E(t+1)∪· · ·∪E(t+(b−1))) is strongly connected.
This guarantees the information of each agent can reach any other agents in the network within finite time.
Assumption II.2. (Balanced Communication). It holds that 1 T n A(t) = 1 T n and A(t)1 n = 1 n , for all t 0. In the consensus literature, the first part of Assumption II.2 is called column stochasticity and is a standard sufficient condition to reach consensus. The second part is called row stochasticity and is needed to guarantee average consensus.
That is, each agent assigns nontrivial weights on information from itself and its neighbors.
Observation model. At each time instant t, each agent independently observes the outputs of a continuous common static latent function η : Z → Y with zero-mean Gaussian noise, where Z ⊆ R nz is the compact input space for η. The observation model is given by Y is the observation of agent i, and e [i] (t) is independent Gaussian noise. Note that we do not assume that input z [i] (t) follows any distribution, which is a standard assumption in statistical learning [4]. We let η(Z) return a column vector [η(z)] z∈Z , and similarly for other functions. Input: network of agents: V; test inputs: Z * ; common inputs: Z agg ; adjacency matrix: A(t); prior mean function: µ; kernel function: k; noise variance: (σ

4:
for t = 1, 2, · · · do 5: for i ∈ V do 6: Zagg|D(t) ) 2 ) 10: end for 11: end for 12: end procedure For notational simplicity, it is assumed that the output space Y ⊆ R because multi-dimensional observations can always be decomposed as aggregation of one-dimensional observations. Problem Statement. The objective of this paper is to design a distributed algorithm for the agents to learn the common static latent function η via streaming data The challenges of the problem stem from the fact that the training dataset is monotonically expanding due to incremental sampling while the agents have limited resources in communication, computation and memory.
The followings are examples of potential applications of this formulation. One example can be a group of mobile robots deployed in a vast open area to collaboratively monitor a static signal, such as temperature or wind field (see the case study in Section VI). Other examples includes the learning of the dynamics of a moving target using a network of static sensors [28]. In addition to robotic applications, this formulation also applies to profit predictions in marketing and wheat crop prediction [29].

III. PRELIMINARIES
In this section, we provide necessary background on GPR. Let η : Z → R be the target function, where Z ⊆ R nz . Given input z(t) ∈ Z at time t, the corresponding output is: y(t) = η(z(t)) + e(t), e(t) ∼ N (0, σ 2 e ), where e(t) is the Gaussian measurement noise. Let training data be in the form D (Z, y), where Z {z(1), · · · , z(n s )} is the set of input data and y [y(1), · · · , y(n s )] T is the column vector aggregating the outputs. GPR aims to estimate the function over a set of test data points Z * ⊂ Z using D by modelling η as a sample from a Gaussian process prior. Define kernel function k : R nz × R nz → R that is symmetric and positive semi-definite; i.e., k(z, z )f (z)f (z )dν(z)dν(z ) 0 for all f ∈ L 2 (Z, ν), where ν denotes a measure (page 80, [7]). By modeling η(·) as a sample from the Gaussian process prior specified by mean function µ(·) and kernel k(·, ·), the training outputs y and the test outputs η(Z * ) are jointly distributed as: where K(Z, Z * ) returns a matrix such that the entry at the i th row and the j th column is k(z(i), z * (j)), z * (j) ∈ Z * , and analogously for K(Z, Z) and K(Z * , Z * ).

IV. LIGHTWEIGHT DISTRIBUTED GPR
In this section, we propose the Lightweight distributed GPR (LiDGPR) algorithm which allows the agents to collaboratively learn the static latent function subject to limited resources. As shown in Figure 1, LiDGPR is composed of three parts: (i) agent-based GPR (Algorithm 2), where the agents make their own predictions of η over a given set of points of interest Z * ⊆ Z using local streaming data {z [i] (1), · · · , z [i] (t)} aggregates local input data and y [i] (t) [y [i] (1), · · · , y [i] (t)] T aggregates the outputs; (ii) distributed GPR (Algorithm 3), where the agents integrate their predictions with those of their neighbors on a pre-defined common set Z agg ⊂ Z * and for z * ∈ Z * do 3: choose z 8: end procedure estimate the predictions on this set given the global training datatset D(t) ∪ i∈V D [i] (t); (iii) fused GPR (Algorithm 4), where the agents refine the predictions on Z * by fusing the results from distributed GPR with those from agent-based GPR. The formal statement of LiDGPR is presented in Algorithm 1. For each iteration t, each agent i collects data online and updates local dataset , and then sequentially executes agent-based GPR, distributed GPR, and fused GPR.

A. Agent-based GPR
To reduce computational complexity, we implement Nearest-neighbor GPR as agent-based GPR. Instead of feeding the whole training dataset to full GPR in (2), agent-based GPR only feeds the nearest input z    can be the average of the observations. The predictive mean and variance for each z * are given in Line 4 and 5 of agent-based GPR. Agent-based GPR returnsμ Z * |D [i]

B. Distributed GPR
Note that each agent only maintains D [i] (t), a portion of the global training dataset D(t). Besides collecting more data, information exchanges between the agents could enhance the learning performance upon agent-based GPR. However, limited communication budget prevents the agents from sharing D [i] (t), whose size monotonically increases. Hence, we develop distributed GPR where the agents communicate with the predictive means and variances over a common set Z agg .
In order to deal with large dataset using GPR, local approximation methods such as Product of Expert (PoE) [30] and Bayesian Committee Machine [31] are proposed to factorize the training process. We consider the following PoE aggregation model for predicting each z * ∈ Z agg : for z * ∈ Z agg do {Dynamic average consensus} 3: z * (t) 12: end for 13: Zagg|D(t) ) 2 15: end procedure which is consistent with full GPR (Proposition 2, [32]).
The two summations in (3) and (4) involve the global training dataset. To decentralize the computation, we consider the computation of the two summations as a dynamic average consensus problem and use FODAC in [33] for each z * ∈ Z agg . State λ [i] (t) is used as one of the criteria for applying fusion between agent-based GPR and distributed GPR in fused GPR. The dynamics of FODAC is shown in Line 4, Line 6 and Line 8 of distributed GPR respectively for each consensus state, where ∆r(t) r(t) − r(t−1) denotes the temporal change of the signal r. Specifically, θ  θ,z * (t) defined in Line 3 among the agents. In particular, agent i computes a convex combination of θ , and then adds the combination into the temporal change of r z * (t) are similar. The updated states are sent to each agent in N [i] (t) as in Line 13. Notice that consensus is not necessarily reached at each time t. We will show that consensus is reached in an asymptotic way in Section V-B.

C. Fused GPR
Fused GPR aims to refine predictions of η(Z * ) by integrating agent-based GPR with distributed GPR. The goal is to obtain an estimate of the predictive distribution P {η(z * )|D(t)} for each z * ∈ Z * . Note that distributed GPR obtains new

13:
end for 14: Z * |D(t) ) 2 15: end procedure estimates of η(Z agg ) by combining results from each agent through convex combination. It can return results with more uncertain predictions, and these predictions should be ignored. The set of inputs predicted by distributed GPR with lower uncertainty is defined as Z agg (t) is the set of inputs z agg ∈ Z agg , where the two variance estimates from distributed GPR, the estimates of (σ Notice that for all z agg ∈ Z [i] However, P {η(z agg )|D(t)} and P {η(z * )|η(z agg ), D(t)} are unknown but can be estimated using the results from distributed GPR and agent-based GPR respectively. In particular, the results from distributed GPR are used to estimate P {η(z agg )|D(t)} since the estimate has lower variance (uncertainty). The results from agent-based GPR are used to estimate P {η(z * )|η(z agg ), D(t)} because distributed GPR is limited to Z agg and z * may not be in Z agg . The product of P {η(z * )|η(z agg ), D(t)} and P {η(z agg )|D(t)}, which yields P {η(z * ), η(z agg )|D(t)}, then contains information from the local agent and those from the other agents in the network. The overall process can be interpreted as a fusion of global information with local information, where improvement is ex-pected because more information is provided. After integrating over η(z agg ), we obtain the estimate of P {η(z * )|D(t)}. The detailed procedure is broken down into the following steps.

D. Choice of the kernel
In this paper, we assume the following properties of the kernel k used in LiDGPR algorithm.
3) (Monotonicity). It holds that κ(r) is monotonically decreasing as r increases and κ(0) = σ 2 f . Remark IV.2. In GPR, kernel can be interpreted as the prior correlation between function evaluations. For a continuous function, it is reasonable to assume bounded correlation and the correlation is negatively related to the distance between two inputs. One example that satisfies Assumption IV.1 above is the class of squared exponential kernels having the form [7]). To obtain the theoretic guarantees in Section IV-E, σ 2 f is chosen for initialization as follows. Let σ 2 When , µ χ and σ 2 χ converge to positive constants, + has growth rate O(σ 2 f ), which gives the left hand side of (5) diminishing at O( 1 (5) is satisfied when σ 2 f is sufficiently large. A distributed way to choose a single σ 2 f is as follows. By using the Floodset algorithm (page 103, [34]), each agent i sends (σ e ) 2 to its neighbors. By Assumption II.1, within n(b − 1) iterations, each agent obtains a copy of (σ Then all the values in (5) can be calculated. Therefore, each agent can identify σ [i] f to satisfy (5). Finally, the agents employ maximum consensus [35] to compute σ f = max i∈V {σ

E. Performance guarantee
In this section, we present the performance of predictive mean and variance returned by LiDGPR. The main results are summarized in Theorem IV.3 and Theorem IV.8, and their proofs are presented in Section V-B and Section V-C.
Part of the performance is quantified in terms of the dispersion of local data defined as We can interpret dispersion as a measurement of how dense the sampled data is distributed within a compact space. For notational simplicity, we introduce shorthand ρ Z z ρ(z, Z). Theorem IV.3 shows that LiDGPR makes predictions with lower uncertainty than agent-based GPR.
Theorem IV.3. (Uncertainty reduction). Part I: Suppose Assumption IV.1 holds. For all z * ∈ Z and i ∈ V, the predictive variance by agent-based GPR is bounded as Part II: Suppose Assumptions II.1, II.2, II.3 and IV.1 hold. For all z * ∈ Z and i ∈ V, there exists a non-negative sequence γ σ,z * (t) such that the predictive variance by LiDGPR is In particular, if Z [i] σ,z * (t) = 0; otherwise : We provide the steady-state results assuming that the dispersion is diminishing. Lemma 6 in [36] shows that dispersion does go to zero under uniform sampling.
To ensure the improvement on prediction accuracy, we need to assume that the prior covariance function of η is correctly specified. Note that any non-zero mean Gaussian process can be decomposed into a deterministic process plus a zero-mean stochastic process such that GPR can be performed over the zero-mean stochastic process (page 27, [7]). Therefore, without loss of generality, we assume η follows a zero-mean Gaussian process for notational simplicity.
That is, the target function η is completely specified by a zero-mean Gaussian process with kernel k. This assumption is common in the analysis of GPR (Theorem 1, [37]).
Furthermore, we need to assume that the state transition matrix induced by A(t) is constant.
One example that satisfies this assumption is each entry of A(t) being constant 1 n , which is a complete graph. Furthermore, we assume η is Lipschitz continuous.
Assumption IV.7. There exists some positive constant η ∈ R such that sup z,z ∈Z |η(z) − η(z )| η ρ(z, z ). Theorem IV.8 below compares the predictive errors of agent-based GPR with those of LiDGPR. , > σ max e , the predictive error resulted from agent-based GPR is bounded as The two theorems indicate that LiDGPR leverages interagent communication to improve transient and steady-state learning performance; meanwhile, no agent suffers from degraded learning performance. This improvement of learning performance is achieved by the fact that the agents whose data samples have higher dispersion (or observation noise has higher variance) benefit from those with data samples having lower dispersion (or observation noise having lower variance) via communication. Next we elaborate on the fact.
Theorem IV.3 indicates that agent i benefits in variance prediction when Note that ρ is closely related to dispersion d [j] (t) and recall the monotonicity property of κ. Hence, LiDGPR enables agents whose data samples have higher dispersion (data sparsely sampled) and observation noise has higher variance to benefit from those with data samples having lower dispersion (data densely sampled) and observation noise having lower variance.
Steady-state improvement. If lim inf σ,z * (t) > 0, it indicates that agent i obtains improvement in steady-state learning performance in predictive variance. From Corollary IV.4, we can see that lim inf z * |D(t) −η(z * )) 2 − (μ z * |D [i] (t) − η(z * )) 2 ] are positively related to κ(ρ agg * (t) z * ) 2 and κ(ρ Zagg z * ) respectively. By monotonicity of κ in Assumption IV.1, these terms indicate that the benefit brought by communication decays as z * is moving away from z agg * (t) and Z agg respectively. That is, a denser set Z agg could induce larger improvements.

F. Discussion
Relevance: The two theorems indicate that both prediction uncertainties and prediction errors reduce as local dispersion d [i] (t) reduces. This provides insights on data sampling such that the agents should sample in a way that agg * (t) z * ) 2 and κ(ρ Zagg z * ) in Theorem IV.3 and Theorem IV.8 show that the improvement of learning performances obtained from communication decreases as the test point z * is moving away from Z agg . This can guide the design process of Z agg such that if the test points in Z * are known a priori, Z agg should be allocated such that sup z * ∈Z * min zagg∈Zagg ρ(z * , z agg ) is minimized; otherwise Z agg should be designed such that sup z * ∈Z min zagg∈Zagg ρ(z * , z agg ) is minimized.
Complexities related to Z * and Z agg . Nearest-neighbor GPR vs. full GPR. Part I of Theorem IV.3 and Theorem IV.8 characterize the steady-state errors of agentbased GPR. Paper [8] shows that σ 2 e ) 2 and µ z * |D [i] (t) → η(z * ) almost surely as t → ∞ for full GPR. Part e ) 2 , hence the variance for noisy prediction (page 19 [7]) equals e ) 2 , and Theorem IV.8 indicates The discrepancy can be caused by the fact that full GPR in [8] makes prediction using all the data in the dataset while Nearest-neighbor GPR only uses the data of the nearest input. Full GPR has computational complexity O(t 3 ) while Nearest-neighbor GPR has the same computational complexity as nearest neighbor search, which is O(t) for the worst case [38]. This is the trade-off between learning accuracy and computational complexity. Note that both full GPR and Nearest-neighbor GPR have the same steady-state errors under noise-free condition, i.e., (σ

V. PROOFS
In this section, we present the derivation of Line 11-12 in fused GPR and the proofs of Theorem IV.3 and Theorem IV.8. Table I shows the symbols that are used in multiple important results and the relation among them.

Symbol
Meaning/definition Equivalencě Predictive mean from agent-based GPR of agent iμ Predictive mean from fused GPR of agent ǐ Predictive variance from agent-based GPR of agent iσ2 Predictive variance from fused GPR of agent i θ,z * (t) Reference signal for consensus state θ in distributed GPR r Reference signal for consensus state ξ in distributed GPR λ,z * (t) Reference signal for consensus state λ in distributed GPŘ Real-valued component of r Recall thatP [i] {η(z * )|η(z agg ), D(t))} is given by applying identities of joint Gaussian distribution (page 200, [7]) tõ This givesP [i] {η(z * )|η(z agg ), D(t))} as a Gaussian distribution with mean and variancẽ Notice that the mean and variance ofP [i] {η(z agg )|D(t)} is given in distributed GPR. Then we have the product z * |zagg,D(t) ) 2 ).
Replacing z agg with z First of all, we introduce some properties of functions f 1 : f +x . These will be used in later analysis. Lemma V.1. It holds that

Proof:
It is obvious that f 1 is convex. Then Jensen's inequality (page 77, [40] It is obvious that f 2 is concave. By Jensen's inequality and concavity, we have Applying Jensen's inequality utilizing the monotonicity and the convexity of inverse function f 3 (x) = 1 x for x > 0, we can also obtain 1

1) Variance analysis of agent-based GPR:
In this section, we present the proof of Theorem IV.3 Part I.
Proof of Theorem IV.3 Part I: Pick any z * ∈ Z. By monotonicity of κ in Assumption IV.1, Line 5 in agent-based GPR gives the predictive varianceσ 2 Combining this with the decomposition property of κ in Assumption IV.1 giveš The definition of local dispersion . Combining this with the monotonicity of κ in Assumption IV.1 gives κ(d [i] (t)) κ(ρ Applying the boundedness of κ in Assumption IV.1 to (6), we converges to its lower bound.
2) Variance analysis of distributed GPR: First, we define the following notations. We let operators ∆, sup, max and min be applied element-wise across the vectors: First of all, we introduce several properties in Lemma V.3.
Proof: We prove the claims one-by-one: Proof of Claim V.3.1: Recall the boundedness property in Assumption IV.1 requires that κ(·) > 0. Therefore it follows from Part I of Theorem IV.3 that Combining this with the definition of r ξ,z * (t) on Line 5 in distributed GPR, gives Combining with initial condition r
Since we have initial condition r 1 |Zagg| , the claim holds for t = 1. Suppose it holds for t = m. Then for t = m + 1, according to distributed GPR Line 6, we have where the inequality follows from the row stochasticity in Assumption II.2. This proves the upper bound of the induction. By Claim V.3.2, ∆r z * (m), following from (7) we have ξ,z * (0) + ∆r The proof for the induction is completed. Second, we find the upper bound of ξ,z * (τ ) 0 for all τ 0 and all j ∈ V. Hence Given initial condition r ξ,z * (−1) = 0 and recall the definition of operator ∆ where ∆r(t) r(t) − r(t−1), it follows that ξ,z * (t − 1). Therefore  ξ,z * (t − 1) n( which, by (8), is also an upper bound for ξ,z * (τ ). Proof of Claim V.3.4: By Claim V.3.2, we have ∆r min,ξ (t) 0 |Zagg| , and hence The definition of r ξ,z * (t) on Line 5 in distributed GPR gives

Applying the upper bound in Theorem IV.3 Part I renders
Combining the definition of r ξ,z * (t) and Theorem IV.3 Part I renders ξ,z * (s) ( Therefore, we have Based on the monotonicity of κ in Assumption IV.1 and since lim > 0 for all t 0, we can apply manipulation Then we can further obtain e ) 4 . Plugging the above inequality back into the definition of ∆r max (t) and apply the monotonicity of κ in Assumption IV.1, it follows that The proof of the lemma is completed.
We define the subsequence {t j } as follows: t −1 1, where ceil(x) min{x ∈ Z | x x}, and for all j 1, The proposition below characterizes the convergence of predictive variance from distributed GPR. for all i ∈ V and t 1, the convergence of (σ z * |D(t) ) 2 in distributed GPR to (σ (agg) z * |D(t) ) 2 is characterized by: For t < t 0 : For t t 0 : where l(t) is the largest integer such that t l(t) j=−1 t j . Proof: Line 9 in distributed GPR indicates that ξ z * |D(t) ) −2 , and combining Line 5 in distributed GPR with (4) gives 1 n n j=1 r ξ,z * (t) . (9) The upper bound of ξ ξ,z * (t) is found by the following two claims, whose proofs are at the end.  Hence we havē Proof of Claim V.4.2: Outline: Write t = τ 0 + τ for some τ 0 , τ > 0. We first derive a general form of ||δ m (τ 0 + τ )|| ∞ . Then we prove the cases when t < t 0 and t t 0 respectively. First, we give the general form. Applying inequality (B.1) in [33] in vector form, we have Second, we prove the case of t < t 0 . Let τ 0 = 0, τ = t in (10). We first find a uniform upper bound for ω(0, t). Recall that Claim V.3.4 shows that δ r ξ (t) ∆ r max (t) and Claim V.3.5 shows that∆r max (t) is element-wise non-increasing. It follows that δ r ξ (t) ∆ r max (1) for all t 1 and hence Since 0 < ζ 1, we can further write that ||ω(0, t)|| ∞ 2(nb − 1) ∆ r max (1) ∞ 1 ζ . Then (10) becomes Note that the time-dependent term on the right hand side of (11) is exponentially decreasing. Suppose t 0 is the smallest integer such that where t 0 can be obtained as defined. Hence we have which proves Claim V.4.2 for t < t 0 , and for t t 0 , we have Finally, we prove the case of t t 0 . Write τ 0 = t 0 . Then (10) becomes Applying analogous algebra as of ω(0, t) gives Using this and (12) as the upper bound for ||δ m (t 0 )|| ∞ , we can rewrite (13) as Similarly, let t 1 be the smallest integer such that Using similar manipulation as t 0 renders t 1 as defined and By similar logic, we have t j as defined for all j 1 and ||δ m (t)|| ∞ 4(nb − 1) ∆ r max (t l(t)−1 ) ∞ 1 ζ . Corollary V.5 shows that the predictive variance from distributed GPR converges to that from aggregated method.
Proof of Theorem IV.3 Part II: By Line 4 in fused GPR, it is obvious that if Z σ,z * (t) = 0. Now we consider the case when Z σ,z * (t) and its uniform lower bound; verfication of the selection of g(z * , t); derivation of the growth factor of γ σ,z * (t). First, we show the expression of γ σ,z * (t) and derive its uniform lower bound. According to Line 12 in fused GPR, we have (σ Line 2 of fused GPR rules that (σ agg * (t)|D(t) . Obviously, γ σ,z * (t) 0. Second, we verify the selection of g(z * , t). We verify that the selection of g(z * , t) is valid by showing (σ z * |D(t) ) 2 > 0. We analyze each factor of γ 1. Line 2 of fused GPR rules that (σ . By Claim V.3.3 and Line 9 in distributed GPR, σ,z * ,2 < 1. By definition, σ σ,z * (t) by analyzing the growth factor of γ σ,z * ,3 (t) and σ [i] min,z * (t) 2 respectively.
We first consider γ The upper bound given in Proposition V.4 can be written as Then Proposition V.4 gives agg * (t)|D(t) Hence we have lower bound .
The upper bound and lower bound of (σ (agg) z * |D(t) ) 2 is given below, whose proof can be found at the end of the proof.
Claim V.7.1. For each z * ∈ Z * , the aggregated variance returned from (4) can be characterized as Denote φ [i] σ 2 f + (σ ) 2 in the inequality above gives agg * (t) agg * (t) The boundedness of κ in Assumption IV. 1 gives Applying this upper bound to the denominator of the lower bound above gives Now we characterize the rest of factors of γ agg * (t)) 2 = κ(ρ f . Combining this with the lower bound of σ [i] min,z * (t) 4 above gives Combining the lower bounds of all the factors gives The definition in (16) and Claim V.3.6 renders h(t) → 0. This renders the Big O notion.
Proof of Claim V.7.1: Using equality (6), we can characterize (σ (agg) z * |D(t) ) 2 in (4) as Taking the inverse and applying Lemma V.1 by substituting The lower bound provided in Part I of Theorem IV.3 and equation (6) give where the last inequality follows from the fact that C. Proof of Theorem IV.8 In this section, we present the theoretical results that leads to Theorem IV.8. We first present the error between the predictive mean of agent-based GPR and the ground truth, which is the result of Theorem IV.8 Part I, in Section V-C1. Secondly, we characterize the predictive mean returned from distributed GPR in Proposition V.9 in Section V-C2. Lastly, we finish the proof for Part II of Theorem IV.8 in Section V-C3.
1) Mean analysis of agent-based GPR: In this section, we provide the proof of Part I of Theorem IV.8.
Proof of Part I of Theorem IV.8: By Assumption IV.5 and decomposition and monotonicity properties in Assumption IV.1, Line 4 of agent-based GPR becomesμ z * |D [i] . It implies thať By boundedness of κ in Assumption IV.

Combining this with triangular inequality gives
Now we analyze the upper bound of each term on the right hand side of (17). Recall that z Utilizing the Lipschitz continuity of η in Assumption IV.7 gives The observation model (1) gives y . Therefore by Chebyshev inequality (page 151, [42]), for all > 0, we have P {|y Note that |η(z * )| η Z . Applying these two inequalities to (17) gives . The proof is completed by using inequality and the monotonicity property of κ in Assumption IV.1.
depends on latent function η andě e ) 2 ) independent over agent i ∈ V and input z * ∈ Z. It is obvious that lim 2) Mean analysis of distributed GPR: Before presenting the results, we derive the solution to the consensus state θ [i] z * (t), i ∈ V, in terms of input signal ∆r [i] θ,z * (t). We also show the decompositions of r [i] θ,z * (t) and θ [i] z * (t), which separate the two terms into real-valued parts and stochastic parts.
Term 2. By definitions,ř Proof: Since function f 2 in Lemma V.1 is strictly increasing, we have Combining these two statements with (37) gives Taking the inverse of (37) gives where the last inequality follows from Lemma V.1. This gives where the equality follows from (4). This is equivalent to where the equality follows from Corollary V.5. Given (38) and the inequality above, since z agg ∈ Z agg is arbitrary, we have lim We now proceed to finish the proof of Part II of Theorem IV.8.
Proof of Theorem IV.8 Part II: By Line 3-4 in fused GPR, it is obvious that if lim agg (t) = ∅. Lemma V.14 corresponds to the sufficient condition for Proposition V.13, which is the sufficient condition for lim

VI. SIMULATION
In this section, we conduct Monte Carlo simulation to evaluate the developed algorithm. For the algorithms introduced below, we use (NN) to denote the version of the algorithm related to Nearest-neighbor GPR and (full) to denote the version related to full GPR. We compare LiDGPR (NN), i.e., Algorithm 1, with five benchmarks: (i) agent-based GPR (NN), i.e., Nearest-neighbor GPR (Algorithm 2); (ii) agent-based GPR (full), i.e., Algorithm 2 is replaced by (2) and hence   e ) 2 ), t 1. The set Z * of test points are uniformly separated over Z, and |Z * | = 1600. We use 25% of the test points for the set Z agg , i.e., |Z agg | = 400. The points in Z agg are uniformly separated. The kernel is k(z, z ) = σ 2 f exp(−2 z − z 2 ), where σ 2 f is chosen by following the procedure under Remark IV.2 such that (5) is satisfied. The prior mean is µ(z) = 0 for all z ∈ Z.
The performances of the robots are similar, and we present (a) Predictive variance (b) Predictive error Figure 3: Predictive variance and error of cNN-GPR the figures for robot 1 due to space limitation. Let the predictive error at z * ∈ Z * be the distance between predictive mean and the ground truth of η at z * , where the distance adopts 2-norm. For example, the predictive error at z * of agent-based GPR (NN) is (μ z * |D [i] (t) −η(z * )) 2 . When robots' trajectories and η are those in Figure 2, Figure 3 shows the predictive variance and predictive error over Z * of cNN-GPR.
We can see that the predictive variances and errors are smaller near the trajectories of the robots. Figure 4 shows the predictive variances and predictive errors of agent-based GPR (NN) and LiDGPR (NN) over Z * of robot 1. We can see that by only communicating a portion of the testing sets, LiDGPR (NN) improves the learning performances over agent-based GPR (NN) with reduced predictive variances and errors. The red dots in Figures 4c and 4d are the points of Z agg , and the "holes" indicate that the improvements take place around the trajectories (training data) of the other robots, which corresponds to the term κ(ρ Notice that the differences in predictive variances and errors between cNN-GPR and centralized full GPR are small, while the diminishing rates are comparable. This shows that cNN-GPR has small performance loss compared to the benefit in reducing the complexities in computation and memory. Comparing the curves of LiDGPR (NN) with agent-based GPR (NN) and agent-based GPR (full), we can see that LiDGPR (NN) not only compensates the information loss of using agent-based GPR (NN) to approximate agent-based GPR (full), but also gains extra information from the other robots. The wall clock time for prediction using LiDGPR (NN) versus t, the number of local data points, after linear leastsquare fitting, has slope 2.13e-6 second per test point per data point and the bias is 1.19e-3 second per test point. Recall that Section IV-F indicates that agent-based GPR (NN) has complexity linear in the number of data points. Over the simulation, the average is 1000 test-point predictions per second, or 1 kHz, with standard deviation 153 predictions per second.

VII. CONCLUSION
We propose the algorithm LiDGPR which allows a group of agents to collaboratively learn a common static latent function through streaming data. The algorithm is cognizant of agents' limited resources in communication, computation and memory. We analyze the transient and steady-state behaviors of the algorithm and quantify the improvement brought by interagent communication. Simulations are conducted to confirm the theoretical findings. Possible future works include analysis with different implementations of agent-based GPR and the consideration of dynamic latent functions.