Resource-Efficient Cooperative Online Scalar Field Mapping via Distributed Sparse Gaussian Process Regression

To speed up GPR, sparse variational GPR (SVGPR) is developed [11], [12], which employs a set of $M\ll N$ pseudo-points based on KL divergence to summarize the $N$ sample data thereby reducing computational costs to $\mathcal{O}(NM^{2})$. Meanwhile, many Gaussian process aggregation methods are proposed for parallel computation, like product-of-experts (PoE) [13] and robust Bayesian committee machine [14]. Then, SVGPR combined parallel computation [15] and stochastic optimization [16] methods are developed to further reduce computational costs. Moreover, Refs. [17] and [18] provide recursive frameworks for standard GPR and sparse variational GPR to handle online streaming data. Ref. [19] combine these sparse variational, mini-batches and streaming methods and apply them to Bayesian robot model learning. However, these methods both require a central node, in which SVGPR [11, 12, 17, 18, 19] requires global information of the observations for sparse approximations, and the parallel methods [13, 14, 15, 16] require a central node for model aggregations. Therefore, these methods cannot be easily extended to distributed multi-robot systems since all observations and models are required to communicate across networks, which increases inter-robot communications and causes repeated computations about the same observations on different robots.

For dynamic networks without a center, distributed Gaussian process aggregation methods based on average consensus are proposed in [20], [21] and [22]. Then, Ref. [23] applies an online distributed Gaussian process method to scalar field mapping. These methods avoid communicating observations on networks. However, these methods do not provide error bounds with centralized aggregation methods and do not consider the sparse computation, causing $\mathcal{O}(N^{3})$ computational costs. Ref. [24] provides a local online Gaussian process evaluation method to reduce computation and inter-robot communications by only computing and communicating partial observations. Ref. [25] employs the SVGPR for multi-robot systems, in which each robot obtains sparse subsets based on its local observations. Unfortunately, these approaches have limited global accuracy because the online evaluation only uses local information.

This letter proposes a distributed sparse online Gaussian process regression method for cooperative online scalar field mapping. The main contributions of this letter are as follows.

• •

  Propose a distributed recursive Gaussian process based on dynamic average consensus for streaming data and cooperative online mapping. The error bounds between distributed consensus results and the centralized results are guaranteed theoretically (Theorem 1).

• •

  Propose a novel distributed online evaluation method for Gaussian process. Compared with [24] and [25], each robot online evaluates observations based on distributed consensus results instead of only local information to obtain better global accuracy without additional communications.

• •

  Develop a resource-efficient cooperative online scalar field mapping method via distributed sparse Gaussian process regression. The performances of the proposed algorithms are evaluated by real online light field mapping experiments.

(Network Model) Consider a team of $p$ robots, labeled by $i\in V=\{1,...,p\}$. The communication networks at time $t$ can be represented by a directed graph $\mathcal{G}(t)=(V,E(t))$ with an edge set $E(t)\subset V\times V$. We consider that $(i,j)\in E(t)$ if and only if node $i$ communicates to node $j$ at time $t$. Define the set of the neighbors of robot $i$ at time $t$ as $\mathbb{N}_{i,t}=\{j\in V|(i,j)\in E(t)\}$. The matrix $A(t):=[a_{ij,t}]_{i,j=1}^{p,p}\in\mathbb{R}^{p\times p}$ represents the adjacency matrix of $\mathcal{G}(t)$.

(Observation Model) Each robot $i$ independently senses the fields $f:\mathcal{X}\rightarrow\mathcal{Y}$ and gets the measurements $y_{i,t}\in\mathcal{Y}$ with zero-mean Gaussian noise $e_{i,t}$ at time $t$, where $\mathcal{X}\subset\mathbb{R}^{n_{\mathbf{x}}}$ is the input feature space for $f$ and $\mathcal{Y}\subset\mathbb{R}$ is the corresponding output space. This letter considers 2D cases ($n_{\mathbf{x}}=2$) and it can also be used in 3D cases as [25]. The observation model is given by

where $\mathbf{x}_{i,t}\in\mathcal{X}$ is the position of robot $i$ at time $t$. Let the local datasets of robot $i$ at time $t$ be in the form $\mathbf{D}_{i,t}:=\{\mathbf{X}_{i,t},\mathbf{Y}_{i,t}\}$, where $\mathbf{X}_{i,t}=\{\mathbf{x}_{i,k}\}_{k=1}^{t},~{}\mathbf{Y}_{i,t}=\{y_{i,k}\}_{k=1}^{t}$ are the sets of observations from begin to current time $t$.

Let the unknown fields $f:\mathcal{X}\rightarrow\mathcal{Y}$ be the target regression functions. A robot can use Gaussian process regression [26] to predict the latent function value $f(\mathbf{x}^{*})$ at the test data points $\mathbf{x}^{*}\in\mathcal{X}$ using the training datasets $\mathbf{D}=\{\mathbf{X},\mathbf{Y}\}$.

In GPR, the Gram function $K(\mathbf{X},\mathbf{X}):=[\kappa(\mathbf{x}_{i},\mathbf{x}_{j})]_{i,j=1}^{N,N}$ is constructed from the kernel function $\kappa:\mathcal{X}\times\mathcal{X}\rightarrow\mathbb{R}$. A commonly used kernel function is the squared exponential (SE) kernel defined as

with hyper-parameters $\{\sigma_{f},\bm{\Sigma}_{\eta}\}$. Then, the posterior distribution is derived as $\rho(\mathbf{x}^{*}):=\mathcal{P}(f(\mathbf{x}^{*})|\mathbf{X},\mathbf{Y},\mathbf{x}^{*})\sim\mathcal{N}\left(\bm{\mu}(\mathbf{x}^{*}),\bm{\Sigma}(\mathbf{x}^{*})\right)$, where

In this letter, we assume the hyper-parameters can be typically selected by maximizing the marginal log-likelihood offline [26] and focus on GP predictions (2) for field mapping, like Refs. [20], [21], [23], [24] and [25].

In the online mapping mission, the data arrives sequentially. At any time $t$, the new data point $\{\mathbf{x}_{t},y_{t}\}$ is added to datasets $\mathbf{D}_{t}=[\mathbf{D}_{t-1};\mathbf{x}_{t},y_{t}]$ in streaming set. In this case, the Gaussian process prediction (2) will suffer from slow updating. To speed up the online mapping, we employ the recursive Gaussian process update [17],

where $\bm{\alpha}_{t}\in\mathbb{R}^{t}$ and $\mathbf{C}_{t}\in\mathbb{R}^{t\times t}$ are recursive updated variables. Suppose a new data point $\{\mathbf{x}_{t+1},y_{t+1}\}$ has been collected, the recursive updates are designed as,

with the scalars $q_{t+1}$, $r_{t+1}$ and $\gamma_{t+1}$ given by

Upon comparison with (3), $1/\gamma_{t+1}$ can be interpreted as the recursive predictive variance of the new data point $\mathbf{x}_{t+1}$ given the observations at $\mathbf{D}_{t}$ without noise. This property and the variable $\mathbf{Q}$ will be used in Section IV for sparse approximation.

At time $t$, each robot $i$ first samples a new training pair $\{\mathbf{x}_{i,t},y_{i,t}\}$ in the fields. Then, robots recursive update local Gaussian process posterior mapping $\rho_{i,t}^{[L]}(\mathbf{x}^{*})$ by (3) and (4), where $\mathbf{x}^{*}$ are the test map grid positions. With some abuse of notations, denote $\rho_{i,t}^{[L]}$ and $\rho_{i,t}^{[D]}$ as the local and distributed Gaussian posterior mapping of robot $i$ at time $t$, respectively.

Next, we employ discrete-time dynamic average consensus methods [27] for distributed map aggregation. The first-order iterative aggregation process can be written as

where $\bm{\xi}_{i,t}$ denotes the local consensus state of robot $i$ at time $t$ and $\Delta\mathbf{r}_{i,t}=\mathbf{r}_{i,t}-\mathbf{r}_{i,t-1}$ is the reference input, which determines the converged consensus results. Many aggregation schemes can be used to design this reference input. We employ the PoE method in this letter. The distributed Gaussian process update can be designed as

where $\bm{\xi}_{i,t+1}^{[k]}$ denotes the $k$-th element of the local consensus state and $\hat{\bm{\Sigma}}_{i,t}^{[L]}:=\bm{\Sigma}_{i,t}^{[L]}+\sigma_{n}^{2}$ is the local Gaussian posterior variance with the correction biases $\sigma_{n}^{2}$ for avoiding singularity. Due to the theoretical foundation of consensus algorithms, the distributed aggregation mean $\bm{\mu}_{i,t}^{[D]}$ will converge to the centralized PoE results

which will be discussed latter. In particular, robots are only required to communicate the local consensus state $\bm{\xi}_{i,t}$ with constant size, instead of the online datasets $\mathbf{D}_{i,t}$ with growing size.

For the derivation of distributed aggregation error bounds, we introduce the following assumptions.

The Assumption 1 describes the periodical strong connectivity between robots and Assumption 2 defines the parameters selection of adjacency matrix, which are common for average consensus methods [27].

The average consensus methods generally suppose the derivative of the consensus state is bounded [20], [21] and [27]. The derivative bounds are related to the consensus errors and can be used for parameter selections of the communication network in practice, like communication periods and adjacency matrix. However, the derivative bound is hard to get by sample if the field distribution is unknown in practice. This letter employs the character of the recursive Gaussian process to loosen this requirement to the bounded observations (Assumptions 3) for easier parameter selections on online mapping. Moreover, since the Gaussian process predictions are determined by the observations and the prior kernel function, it is bounded as the observations.

The error bounds given in Theorem 1 do not require the specific GP kernels and hyper-parameters and thus can be used in other distributed GP systems for communication parameter selections.

For a robot $i$, the recursive variable scalar $q_{i,t+1}$ and $r_{i,t+1}$ (5) for can be interpreted as the variance-weighted prediction at point $\mathbf{x}_{i,t+1}$ using datasets $\mathbf{D}_{i,t}$. Ref. [24] use the sum of the changes in posterior mean at all current sample points $\{\mathbf{x}_{i,1},...,\mathbf{x}_{i,t}\}$ to denote the score for the observation $\{\mathbf{x}_{i,t+1},~{}y_{i,t+1}\}$, which has been derived as

where $[t+1]$ denotes the $(t+1)$-th item of the vectors and $\mathbf{diag}(\cdot)$ denotes the diagonal vectors of the matrix. However, for a multi-robot system, each robot $i$ has different datasets $\mathbf{D}_{i,t}$ and recursive Gaussian process variables $\bm{\alpha}_{i,t},~{}\mathbf{C}_{i,t}$. The score (12) can only evaluate the points based on local information, which causes repeated information and limited global performances. A 1-D toy example is shown in Fig. 2. Two robots sample the 1-D toy function and compress observations. Compression based on only local information tends to retain repeated information in common areas because of the same local metrics.

In Section III, we have obtained the distributed aggregation results based on consensus methods. And then we employ the the distributed aggregation results $\rho_{i,t}^{[D]}$ to improve the compression performances.

Bhattacharyya-Riemannian distance is employed to evaluate the sample points using the distributed aggregation results. For two Gaussian distributions $\nu_{1}(\mathbf{x})=\mathcal{N}(\bm{\mu}_{1},\bm{\Sigma}_{1})$ and $\nu_{2}(\mathbf{x})=\mathcal{N}(\bm{\mu}_{2},\bm{\Sigma}_{2})$ over feature space $\mathcal{X}$, Bhattacharyya-Riemannian distance is defined as

where $\bar{\bm{\Sigma}}=\frac{1}{2}(\bm{\Sigma}_{1}+\bm{\Sigma}_{2})$, and $\{\lambda_{1},...,\lambda_{n_{\lambda}}\}$ is the eigenvalues of $\bm{\Sigma}_{2}^{-1}\bm{\Sigma}_{1}$. This metric (13) modifies the Hellinger distance by employing Riemannian item $d_{\mathcal{R}}$, which is efficient for the Gaussian variance representation [28].

Denote $\mathbf{S}_{i,k}:=\{\mathbf{x}_{i,k},~{}y_{i,k}\}\subset\mathbf{D}_{i,t}$ as the evaluated observation. We design a distributed sparse metrics as

where $k_{\phi}\in(0,1)$ is a constant weight for the local score and distributed score and $\Omega(\cdot)$ is the normalization function for all evaluated observations. This letter uses the simple min-max normalization. $\epsilon_{i,k}$ denotes the local score (12) for robot $i$ about the evaluated observation $\mathbf{S}_{i,k}$. $\widetilde{\rho}_{i,t,k}^{[L]}$ is the local recursive Gaussian process distribution removing the evaluated observation, which will be introduced latter.

In Section III, we have proved the distributed mean will converge to the centralized PoE results. Therefore, the distributed sparse metrics (14) can find observations of sufficient global utility.

Then, we employ the distributed sparse metrics (14) to select the information-efficient subsets with $M$ points and design the distributed sparse Gaussian process.

Reorder the elements of ${\bm{\alpha}}_{i,t},~{}\mathbf{C}_{i,t},~{}\mathbf{Q}_{i,t}$ so that they are consistent with $\mathbf{S}_{i,k}$ being the last element in the current datasets $\mathbf{D}_{i,t}$, and define

where $\mathbf{C}^{[i,j]}$ denotes the elements of the matrix $\mathbf{C}$ at $i$-th row and $j$-th column and $\tilde{k}=[1,\dots,k-1,k+1,\dots,M+1]$ denotes an index vector from $1$ to $M+1$ except $k$. The recursive variables removing points $\mathbf{S}_{i,k}$ are given as [17]

where $\mathbf{\Gamma}_{i,k}=\mathbf{q}_{i,k}\mathbf{c}^{T}_{i,k}+\mathbf{c}_{i,k}\mathbf{q}^{T}_{i,k}$, so as to compute $\widetilde{\rho}_{i,t,k}^{[L]}$ using the reduced $\hat{\bm{\alpha}}_{i,k},~{}\hat{\mathbf{C}}_{i,k},~{}\hat{\mathbf{Q}}_{i,k}$ by (3).

The entire cooperative online field mapping method is shown in Algorithm 1 and the pipline is presented in Fig. 3. At the time $t$, each robot samples at the field and updates local Gaussian process predictions $\rho_{i,t}^{[L]}$ by (3). Next, robots share local consensus state $\bm{\xi}_{i,t}$ (6) with neighbors and update distributed Gaussian process predictions $\rho_{i,t}^{[D]}$ by (7). Moreover, robots use the distributed sparse metrics (14) to greedily select the information-efficient subsets (Algorithm 2). The recursive variables $\bm{\alpha}_{i,t},~{}\mathbf{C}_{i,t},~{}\mathbf{Q}_{i,t}$ keep lower $M$ dimensions by (17) for fast Gaussian process predictions.

A comparison of time, memory, and communication complexity is presented in Table I. Recall that the robot number is $p$, and each robot samples $N$ training points and selects $M$ sparse points. The numbers of test points $\mathbf{x}^{*}$ and communication edges are $\mathcal{T}$ and $|E|$.

Computation and memory: for each robot, the recursive sparse update of local mapping yields $\mathcal{O}\left(N(M^{2}+M^{3})\right)$ time and $\mathcal{O}(M^{2})$ memory complexity, which is significantly less than distributed recursive GP with $\mathcal{O}(N^{3})$ time and $\mathcal{O}(N^{2})$ memory complexity. Compared with local sparse metrics (12), the proposed distributed sparse metrics (14) have better global accuracy with additional $\mathcal{O}\left(NM^{3}\right)$ computation. Since $M\ll N$, the additional computation is acceptable for online mapping.

Communication: compared with communicating original observations $\mathcal{O}(p^{2}N)$ or compressed observations $\mathcal{O}(p^{2}M)$ [24], [25], the proposed methods only communicate consensus state $\bm{\xi}$ and have $\mathcal{O}(|E|\mathcal{T})$ communication complexity, where $p-1\leq|E|\leq p(p-1)/2$. Therefore, the proposed methods have lower communication costs for large-scale multi-robot systems where $p$ is large. Moreover, the proposed methods can avoid repeated computations based on the same observations for GP mapping in different robots.

Select a $50\times 50$ mesh grid in the workspace as the test points to construct the light field mapping. The size of sparse training subsets is set as $M=13\times 5$. The averaged distributed predictions of all 5 robots for all test points in the experiment and the centralized PoE results are shown in Fig. 4. The communication network is dynamic and Robots $4$ and $5$ are out of communication ranges of Robots $1,~{}2$ and $3$ at first time. When sample number $n>166$, Robots $4$ and $5$ begin sharing information in the communication network and the distributed aggregation results of all robots converge to the centralized PoE results.

Fig. LABEL:fig:field represents the progress of this experiment and the corresponding distributed GP results of Robot $3$ during time $t=35$ ($n=1056$) seconds. Fig. LABEL:fig:field(a) shows the experiment snapshots and the robot number. Fig. LABEL:fig:field(b) shows the results of the distributed online GP variance predictions for Robot $3$ with the exploration trajectories of all robots. The black points denote the sparse subsets for GP predictions. Fig. LABEL:fig:field(c) shows that the GP mean predictions for Robot $3$ gradually converges to the real light field in Fig. LABEL:fig:field(a). Compared with the local GP results of Robot $3$ shown in Fig. 6, the distributed aggregation methods efficiently construct light field maps for unvisited areas.

To analyze the map accuracy with different GP prediction methods, $10\times 10$ equal interval points in the workspace are sampled by hand as the test points. Fig. 7(a) shows the mean RMSE of the distributed predictions on the test points for all robots. Compared with the local data compression [24], the proposed distributed compression improves the prediction accuracy. Fig. 7(b) shows the mean online prediction computational time for all robots on $20$ repeated experiments. The proposed sparse GP has an efficient constant prediction time with the increasing sample $N$. In addition, the data compression is run all the time here to illustrate the constant computation time. In practice, it is not necessary when the sample number $N$ is small.