Large Scale Mapping of Indoor Magnetic Field by Local and Sparse Gaussian Processes

Authors: iad  abdul-raouf
vincent  gay-bellile
cyril  joly
steve  bourgeois
alexis  paljic
Venue: robot-learning_org_CoRL_2024_Conference
Published: 2024-06-05 21:34:45

Abstract: Magnetometer-based indoor navigation uses variations in the magnetic field to determine the robot's location. For that, a magnetic map of the environment has to be built beforehand from a collection of localized magnetic measurements. Existing solutions built on sparse Gaussian Process (GP) regression do not scale well to large environments, being either slow or resulting in discontinuous prediction. In this paper, we propose to model the magnetic field of large environments based on GP regression. We first modify a deterministic training conditional sparse GP by accounting for magnetic field physics to map small environments efficiently. We then scale the model on larger scenes by introducing a local expert aggregation framework. It splits the scene into subdomains, fits a local expert on each, and then aggregates expert predictions in a differentiable and probabilistic way. We evaluate our model on real and simulated data and show that we can smoothly map a three-story building in a few hundred milliseconds.

Main paper: 
 introduction: 
 Indoor positioning for robotics is an active and challenging field of research [1, 2] . While visionbased approaches [3] are the most common to address this issue, they struggle in environments with low or repetitive textures and when visual cues change over time [4] . Recently, a novel and promising approach has emerged that leverages the spatial anomalies of the indoor ambient magnetic field. [5] . These anomalies are sufficiently distinct to differentiate between locations [6] and stable enough to ensure that magnetic field maps remain valid over the years [4] .Indoor magnetic-field-based localization solutions, such as [4] built upon extended Kalman filtering and [1] built upon particle filtering, require a map of the field values at any position, and uncertainty estimation. Kalman filters additionally assume a differentiable observation model. Classic Gaussian Process (GP) regression [7] from magnetometer measurements can achieve a differentiable and probabilistic representation. However, it has a O(N 3 ) computational complexity, with N the number of data points, which becomes intractable for large datasets.To address this, Vallivaara et al. [8] have downsampled the training data and then learned map values using only training data in a small radius. Downsampling is suboptimal, and this local approach is discontinuous [9] . Instead, Solin et al. [10] and Menzen et al. [11] proposed sparse approximations of the GP. While these approaches scale well with increasing data density, they perform less effectively with increasing area size [12, 13] . Besides, to date, the fastest sparse magnetic map [11] is built upon structured kernel interpolation [14] , which may produce discontinuous predictions [9] . This work proposes a new approach to model the magnetic field that scales with increasing area size and data density. It also results in differentiable magnetic predictions that are efficiently computed, thus being compatible with any localization algorithm. We propose splitting the scene into subdomains and learning a Deterministic Training Conditional sparse GP [15] that has been modified to include prior knowledge from Maxwell equations. Then, neighboring experts' predictions are smoothly and efficiently aggregated into what we call a Local Bayesian Committee Machine, derived from the classic Bayesian Committee Machine [16] in a principled way.We thoroughly illustrate our method with several toy examples. We also compare it against the state of the art, on real and simulated data. Evaluations are performed regarding mean square error, mean standardized log loss, smoothness, and runtimes. In particular, we can learn a magnetic map in a three-story building from a trajectory of one kilometer in less than 100 milliseconds. 
background on gaussian process regression: 
 This section introduces the necessary background on Gaussian process regression, sparse approximations, and local expert aggregations. Computational complexities are discussed with respect to data input x density and domain volume instead of the usual dataset size N and number of latent variable M . It brings a new light on the scalability of each model.Let us consider a vector field f (x) ∈ R d defined for any x ∈ Ω ⊂ R d . It must be estimated from a training dataset D = {(x 1 , y 1 ), . . . , (x N , y N )}, where each observation y i of f (x i ) (noted f i ) is corrupted by a Gaussian additive white noise. A classic approach is to model f as a zero mean Gaussian process indexed by x [7] f ∼ GP(0, κ), y = f (x) + , ∼ N (0, σ 2 noise I d ),where κ is the d × d covariance kernel such that κ(x, x ) = cov(f (x), f (x )). The choice of a specific covariance function encodes the a priori knowledge about the underlying process. It is often defined with respect to hyperparameters, which are learned jointly with σ noise from data or manually tuned to sensible values. For any test input x * , the prior and posterior of f (x * ) (noted f * ) are Gaussians and admit a closed-form expressionEQUATIONwhere Σ = (K f ,f + σ 2 noise I nd ) -1 . The vectors y and f are the concatenation (y 1 , . . . , y n ) and (f 1 , . . . , f n ) respectively. We also used the standard shorthand notation x |Ω| 2 ). It is intractable for dense data inputs or in a large environment.To tackle the intractability of exact regression, a large portion of the literature describes sparse approximations [9, 17] . We overview the deterministic training conditional (DTC) sparse approximation formalized by [15] . The starting point is a set of latent variables u = (u 1 , . . . , u M ) called inducing variables in the literature, where M N . They are values of the GP (like any f i and f * ) at the latent inputs z 1 , . . . , z M ∈ R d . The latent variables are involved in two successive approximations of the joint prior.EQUATIONis a conditional independence assumption between the test variable f * and training variables f . The training conditional is further approximatedEQUATIONwhere K u,u , K f ,u and K u,f are again shorthand notation for covariance matrices computed by block (e.g. cov(f i , u j ) = κ(x i , z j )). The approximation (4) assumes the training conditional covariance cov(f |u) is almost 0. Marginalization of u from the approximated p(f * , f , u) create an approximation of p(f * , f ). By injecting it into the classic GP framework, we get the DTC posteriorEQUATIONwhereS = (σ -2 noise K u,f K f ,u + K u,u ) -1 .Training has a computational complexity O(d 3 N M 2 ) and an additional O(d 3 M 2 ) operations are required at each test input x * .Prediction quality depends on the number of latent variables. Experiments show that using too few of them in large input spaces smooths out the model [12, 13] . Thus, defining M proportionally to the environment size: M = ρ z |Ω| is natural. The training and test complexities become respectivelyO(d 3 |Ω| 3 ρ x ρ 2 z ) and O(d 3 |Ω| 2 ρ 2 z ).The key takeaway is that sparse approximations scale well with increasing input density ρ x , but not so much when the environment size |Ω| increases. Therefore, it is tempting to model the vector field on a partition of the input space [10] .As discussed above, splitting the data and fitting a (sparse) Gaussian process on each subset is natural. It yields a discontinuous model, and prediction quality may be poor around the frontier, even if sub-datasets overlap [11] . To improve upon this naive idea, we can take inspiration from the literature on the aggregation of Gaussian experts Liu et al. [9] . We consider the classic Bayesian Committee Machine (BCM) introduced by Tresp [16] .Let D 1 , . . . , D J be a partition of the dataset D. We assume a Gaussian model, i.e., the prior p(f * ) and the posterior of each expert p(f * |D i ) are Gaussians. The BCM introduces a conditional independence approximation between the D i . Combined with Bayes theorem, it follows that the BCM posterior has the form of a product of expert regularized by a prior termEQUATIONThe BCM posterior is Gaussian with closed-form mean and covariance [16] . It is a smooth aggregation that automatically decreases the influence of weak experts (i.e., experts with high uncertainty). Compared to the naive approach, the quality is improved but slower since each prediction requires a prediction from all the experts. 
method: 
 We present a probabilistic method to learn the magnetic field at a large scale without discontinuities in the prediction from magnetic vectors y i ∈ R 3 and their respective positions x i ∈ R 3 . For that, we split the domain into smaller subdomains. On each, we fit a sparse approximation that we call Gradient-DTC (G-DTC) as described in section 3.1. It embeds prior knowledge from Maxwell equations to reduce the computational complexity and improve the quality of the predictions. G-DTC scales well with increasing spatial densities of the magnetic observation. Then, in section 3.2, we develop a local expert aggregation technique, built upon BCM, that scales with respect to the environment size. Assuming local influence of magnetic data, each local expert in (6) makes a negligible contribution far from its training domain. In our local approximation, we modify the BCM posterior in a principled way so that the contribution of distant experts disappears seamlessly.Inspired by [10, 18] , we wish to inject prior knowledge of magnetic field physics in the DTC formulation to increase the model fidelity and reduce the computational cost. According to Maxwell's equations, the magnetic field is curl-free, provided no free current exists in the domain Ω [10] . Under such conditions, a scalar potential φ exists such that H = -∇φ. Instead of u representing magnetic field variables like f , the inter-domain Gaussian process method [19, 17] allows us to mix f and u of different natures. Using this approach, we will manage to do the bulk of the calculation with the scalar latent variables u 1 , . . . , u M , representing -φ(z 1 ), . . . , -φ(z M ). It will reduce the complexity related to the observation dimension d.By definition, if f and u are two Gaussian processes defined over Ω f and Ω u , and if for every x 1 , . . . , x n ∈ Ω f and z 1 , . . . , z m ∈ Ω u the vector (f (x 1 ), . . . , f (x n ), u(z 1 ), . . . , u(z m )) is Gaussian, then f and u are said to be jointly Gaussian, or to form an inter-domain Gaussian process. It is completely defined by their respective mean and kernel functionµ f (x), κ f (x, x ), µ u (z), κ u (z, z ),as well as the inter-domain kernel κ uf (z, x) = κ f u (x, z) . [10] suggest to define u as a Gaussian process with differentiable mean µ u = 0 and twice differentiable kernel:EQUATIONwhere l SE and σ SE are respectively the characteristic length-scale and amplitude of the process. Under such regularity conditions, f = ∇u is a well defined Gaussian process and f, u are jointly Gaussian [20] with explicit expressions for the mean and covariance kernelsEQUATIONEQUATIONEQUATIONWe directly substitute the overloaded kernel κ u , κ f , κ f u in (5) to compute the covariance matrices and get what we call G-DTC (G for gradient). Using the overloaded kernel reduces the size of K f * ,u , K f ,u , K u,f and K u,u . Indeed, each block computed from κ f u (x, x ) is a vector of size d × 1, and The remaining question is whether the DTC sparse process is a good approximation in this gradient setting. Assuming a differentiable scalar function is perfectly known, then its gradient must also be perfectly known. With enough latent variable u, the DTC hypothesis cov(f |u) 0 from (4) should be valid. However, Fig. 1 shows that the posterior is not properly learned when we test this hypothesis on a toy example. In this example, all latent inputs z i ∈ R 3 are set along the Ox axis. Therefore, the variations of u are modeled in this direction only. Only the first component of f = ∇u can be estimated confidently, and it is indeed possible to check that cov(f |u) = 0. This weakness is shared among state-of-the-art maps using the latent scalar potential [10, 11] . They map the field in 3D even though some of their datasets have positions on a 2D plane.κ u (x,To solve this issue, we set the latent inputs z i so that the scalar potential variations are modeled in all directions. The z i are vertices of a cubic grid of step δ = l SE /2. We keep only the vertices at a distance R or less from at least one training x i . The parameter R should be as small as possible to reduce the number of latent inputs but large enough such that even isolated x i have enough z i in their neighborhood. We set R = δ 2 + (δ/2) 2 + (δ/2) 2 1.23 • δ as justified in the appendix F.As explained previously, the number of latent variables in sparse GP must grow with the domain size. Therefore, we split the domain and fit a G-DTC for each subdomain. Then, we aggregate each local expert prediction with a local version of the BCM. To reduce BCM test runtime, we seek to consider only experts 'nearby' the prediction location x * [21, 22] , where the 'proximity' is defined by the Gaussian process kernel κ. The subdomain Ω i is 'far' from x * if κ(x, x * ) is small for all x ∈ Ω i [21] . Using the classic squared exponential (SE) kernel ( 8), we would like to ignore all the experts farther than a distance l max set to a few l SE .In the BCM framework, it is easy to see from ( 6) that discarding the prediction of the i-th expert is equivalent to replacing the expert posterior by the prior (it cancels out with the denominator). Thus, farther than l max from Ω i , we approximate the posterior of the i-th expert by the prior. Inside Ω i , we keep the posterior with no further approximations. Everywhere in between, we approximate the posterior of each expert by a geometric mean (GM) of the posterior and the priorEQUATIONwhere 0 ≤ β(r i ) ≤ 1, and r i is the distance separating x * from Ω i . When β(r i ) = 1, there is no approximation; when β(r i ) = 0 the posterior is replaced by the prior. We may drop the argument r i and note it β i for notional simplicity.Let A(x * ) = {i|β(r i ) = 0} be the set of index of active experts. By injecting the geometric approximation (11) in the BCM posterior (6) , we define the Local Bayesian Committee Machine (LBCM) posteriorEQUATION)Then, as developed in the appendix E, the predictive equations areEQUATIONwhereΛ = (1 -i∈A(x * ) β i )cov(f * ) -1 + i∈A(x * ) β i cov(f * |D i ) -1 . The notation E(f * |D i ) and cov(f * |D i )are for the mean and covariance of the exact posterior p(f * |D i ). When all the experts are active, this posterior is precisely what was called the Robust Bayesian Committee Machine (RBCM) in [23] . The difference between our LBCM and their RBCM is the choice of β i . Theirs are not functions of the distance r i . They set β i equal to the differential entropy between the prior and the posterior of expert i since their goal was to increase the robustness in the presence of weak experts.On the contrary, we design our β i to save computation time by imposing that its support is known and bounded around Ω i . We use the cubic hermit spline β(r) = 2(r/l max ) 3 -3(r/l max ) 2 + 1 where 0 ≤ r ≤ l max . It is defined such that it is differentiable, even at the transition points 0 and l max .The optimal l max depends on your computational budget and quality requirements. In appendix C we illustrate β and derive the exact expression of p GM . It is a Gaussian, which allows us to compute its Kullback-Leibler Divergence from the true local expert posterior and justify that l max = 2l SE is a reasonable choice.Our final model defines a G-DTC expert on each of the J subdomains of volume V , and it aggregates their predictions with LBCM. The training and test complexities are O(d|Ω|V 2 ρ x ρ 2 z ) and O(d|A(x * )|V 2 ρ 2 z ) respectively, where |A(x * )| is the number of active experts. Computing predictions of the magnetic field became independent from the input density and the scene volume.LBCM is based on BCM, whose performance depends on the domain partition. With smaller subdomains, each expert is faster. However, the conditional independence assumption in ( 6) is challenged at the frontier between domains. The predictions are better at the interior of a domain, and therefore, the domains must be large enough with respect to l SE . You should also use 'bulky' subdomains to minimize the contact surface. For instance, prefer cubes over flat cuboids. Furthermore, if possible, it is better to set boundaries such that there is little or no data near them. For instance, you could cut horizontally between two floors in a multi-story building. We create a regular partition of input space using boxes of shape [0, L 1 ) × • • • × [0, L d ). 
experiments: 
 In this section, we compare our approach against the state-of-the-art [10] on both real and simulation data. We first illustrate how G-DTC can take advantage of the spatial distribution of the dataset on a manifold, such as when a robot trajectory is constrained in narrow spaces or limited to 2D navigation on the floor. Then, we show that LBCM combined with G-DTC map efficiently and smoothly the magnetic field of large areas.Datasets. Here we introduce three small simulated datasets called Simu1D, Simu2D and Simu3D illustrated Fig. 4 and a Real data in Corridor were collected by the magnetometer and the IMU in an SBG-Ellipse-N rigidly fixed to 4 FLIR Blackfly S cameras. The IMU and the four cameras are used in a graph SLAM [3] to obtain accurate localization of the magnetic observations. The dataset comprises one training and one test trajectory of approximately one kilometer each. They contain 15600 and 16600 observations, respectively, in 1400 m 3 of corridors and staircases in a three-story building. Our model uses a 0 mean prior, so the training observations are preprocessed to subtract their empirical mean, which is then added back to the predictions.Metrics. We evaluate all models with the Mean Squared Error (MSE) and the Mean Standardized Log Loss (MSLL). Both are standard metrics described in [7] . The MSE evaluates how close the posterior mean is from the ground truth. The MSLL also takes into account the quality of the uncertainty. The lower both are, the better the model is. All runtimes are measured on the same laptop from a C++ implementation. We used one core of an 11th Gen Intel® CoreTM i7-11800H with 16 GB of RAM. Prediction times refer to processing the entire test dataset, not just one x * .G-DTC And LBCM Setup. G-DTC depends on the hyper-parameters l SE , σ SE , and σ noise . They have the same physical meaning as in a full GP using the covariance kernel (10) . Therefore, on Simu1D, Simu2D, and Simu3D we use the available ground truth. Despite our best efforts, unaccounted-for localization errors remain on Corridor. We reflected it in the sensor noise by setting σ noise = 4 µT. Then, we randomly selected 1000 training observations in a small area and optimized a GP log marginal likelihood to tune l SE and σ SE . We get l SE = 1.35 m and σ SE = 6.9 µT. LBCM depends on the domain partition and l max . For reasons that will become clear after the experiment 4.2, we set the side length of the boxes to be equal to 3l SE , except when the data covers a multi-story building. In this case, we exploit the natural floor partition and decrease the box height to the height of a story ( 3 meter). Finally, we take l max = 2l SE as justified section C.Competing Method Setup. Solin et al. [10] introduced a sparse approximation of local experts and a strategy to partially smooth the predictions at the subdomains' boundaries. We abbreviate both contributions G-Hilbert and Overlap respectively. Overlap requires to define a space partition, and we use the same as for LBCM. Furthermore, each local expert extends its domain by a range l train to include training data from neighbor experts, such that the "training boxes" overlap. We set l train = 0.5l SE as in [11] . G-Hilbert relies on basis functions analogous to our latent variables. Three of them are responsible for modeling the unknown mean of the GP prior. In this work, all datasets are preprocessed to be compatible with our 0 mean GP prior assumption. Thus, we remove these three additional basis functions for a fair comparison. Because of mathematical details, G-Hilbert must be defined on a larger domain than the training boxes. We extend each training box by 2l SE as in [11] . Finally, G-DTC and G-Hilbert share the same values for l SE , σ SE , and σ noise .In this first experiment, we illustrate G-DTC efficiency when training inputs are located on or near a manifold of lower dimension. For that, we evaluate G-DTC and G-Hilbert on Simu1D, Simu2D, Simu3D and summarize the results in table 1. On Simu1D, G-DTC has similar prediction quality and is one order of magnitude faster than G-Hilbert. Increasing the manifold's dimension in Simu2D and then in Simu3D decreases our speed advantage. In this box of size 3l SE × 3l SE × 3l SE , we could have at most 8 × 8 × 8 = 512 latent inputs. Our worst case is the same as the normal case for G-Hilbert. Suppose we use the same number of latent variables; by looking at the posterior equations, we expect G-DTC to have the same fit time and twice larger prediction time than G-Hilbert. This experiment confirms it on Simu3D. For robotic indoor navigation, 1D and 2D input distributions are the most common, which is where our method shines. Notice that we consider the intrinsic dimension and that in this regard, even a trajectory in a staircase is 1D.We proved the pertinence of G-DTC on a small domain. Now, we evaluate the scalability and regression quality of LBCM combined with G-DTC. We compare LBCM against Overlap and against the naive aggregation method that uses the space partition as is, without any smoothing strategy. We also compare our full model G-DTC & LBCM against [10] for completeness.G-DTC & Naive G-DTC & Overlap G-DTC & LBCM G-Hilbert & Overlap [10] MSE (µTFig. 2 shows the map built from our model, and table 2 exposes the quantitative results. Surprisingly, the MSE and MSLL are similar for LBCM, Overlap and the naive approach. These metrics do not convey the discontinuity issues visible in Fig. 3 that arise at the box boundaries. We see that overlapping training data is insufficient and that LBCM is the only smooth aggregation method.It matters because such discontinuities happen every few meters at each domain transition. Using LBCM comes at the cost of slightly slower predictions because it aggregates a few neighbor experts.On the other hand, Overlap is slower to train since overlapping boxes are bigger and contain more inputs and latent inputs. Furthermore, overlapping the domains creates 186 experts instead of 140, with additional experts on the training trajectory's side, top, or bottom. Compared to ours (G-DTC & LBCM), the approach of [10] is two orders of magnitude slower to fit and one order of magnitude slower to predict, mainly because G-Hilbert does not exploit the intrinsic smaller dimension of Corridor. Finally, our method combines efficiency with prediction quality and differentiability. 
discussion: 
 Limitations: Our method demonstrated its scalability on large datasets. However, several reasons limit our approach's deployment:For now, the method fits the map to the data offline. Adding measures incrementally to DTC (called Projected Process in [24, 25] ) is possible, and we could extend it to G-DTC.Furthermore, we assumed that the localization x of each measurement y is perfect, which is false in practice and deteriorates predictions. We solved it by increasing the observation noise, which is an easy fix for model discrepancies [1] , but it would be more rigorous to explicitly model input noise [26, 27] .Moreover, we use a low number of latent points to improve runtimes, but it can lead to underconfident G-DTC predictions on dense datasets, especially when the latent point layers are thin, such as in Simu1D and Simu2D.Another limitation is in the number of new parameters introduced. There is l max and the subdomain size and shape for LBCM, and G-DTC requires a step size δ and a range R for its partial grid. We described some heuristics to tune them throughout the article, but they remain empirical.Finally, we describe how to scale the fit and predict operations on large datasets. Still, we do not provide a way to scale the tuning of our model hyper-parameter σ noise , σ SE , and l SE . Instead, we learned them from the full GP on a subset of data. It is sensible because they share the same physical meaning for both models, but it is suboptimal.We introduced an efficient approach that creates a smooth map of the magnetic field values and associated uncertainties. We first introduced G-DTC to model the magnetic field, which scales well with respect to increasing data density. It includes physics knowledge from Maxwell equations and can exploit the intrinsic lower dimensions of inputs that naturally arise when collecting measurements along a trajectory. It is faster than other smooth and sparse GP approximations of the magnetic field. We then introduced LBCM to scale with the scene size. Unlike the other approaches, it aggregates local experts in a differentiable way. Efficiency and differentiability are key properties required to develop graph SLAM localization algorithms. For that, online update strategies and robustness to training localization error may be the focus of future research. The three simulated datasets Simu1D, Simu2D, Simu3D displayed in Fig. 4 in fact, just one dataset of 3100 data points split into four parts. Each simulated dataset comprises 1000 data points for training, and all share the same 100 test data points. The 3100 observations were sampled jointly from a Gaussian process prior using the kernel (10) where l SE = 1 and σ SE = 1. Gaussian white noise was also added to the training observation of each dataset, with σ noise = 0.1.This appendix provides useful identities for manipulating Gaussian distributions and computing the Kullback-Leibler divergence.Let p(x) = N (x|µ 0 , Σ 0 ) be the density of a multivariate Gaussian vector of dimension D, and α > 0 a positive scalar. After re-normalization, the power of a Gaussian remains GaussianEQUATIONFurthermore, if q(x) = N (x|µ 1 , Σ 1 ) is another density of a multivariate Gaussian vector, then after re-normalization, the following product is Gaussian as wellEQUATIONwhereµ 2 = Σ 2 (Σ -1 0 µ 0 + Σ -1 1 µ 1 ), Σ 2 = (Σ -1 0 + Σ -1 1 ) -1 .Their quotient also admits a Gaussian expressionEQUATIONwhereµ 3 = Σ 3 (Σ -1 0 µ 0 -Σ -1 1 µ 1 ), Σ 3 = (Σ -1 0 -Σ -1 1 ) -1 .Notice that Σ 3 can be non-positive, making this quotient unsuitable for a random variable density. However, if Σ 3 is definite positive, this expression can be normalized into a proper Gaussian density N (µ 3 , Σ 3 ).Working with Gaussian densities is convenient since, among many other reasons, their Kullback-Leibler divergence admits a closed-form expressionEQUATIONwhere the notation |A| and tr(A) mean respectively the determinant and trace of a matrix A. Section 3.2 we approximated the BCM posterior to introduce LBCM, which happens to have the same posterior as RBCM [23] , up to the choice of β. To the best of our knowledge, it is the first time RBCM is viewed as an approximation of each posterior by a geometric mean between the posterior and prior. We use this theoretical reinterpretation to derive a new metric on the approximation quality with respect to β.Assuming the BCM is a good approximation for the full GP posterior, a reasonable and simple way to evaluate the quality of LBCM is to evaluate the approximation separately for each expert. We consider the Kullback-Leibler Divergence (KLD) of the geometric approximation compared to the true posterior:D KL (p(f * |D i ) || p GM (f * |D i )).To compute it, we need an explicit expression for p GM (f * |D i ). After normalization, identities ( 14) and ( 15) tell us that the power of a Gaussian density and the product of two Gaussian densities are also Gaussian. Combining these two results, we getEQUATIONwhereΛ = β i cov(f * |D i ) -1 + (1 -β i )cov(f * ) -1 .Then, the KLD between two Gaussian distributions admits a closed-form expression.We can now experiment with the β function to choose one for which the KLD remains small everywhere. For differentiability, we restrict our-self to the cubic hermit spline β(r) = 2(r/l max ) 3 -3(r/l max ) 2 + 1, where r is the distance separating x * from the expert domain. β is defined up to a hyper-parameter l max that we will tune experimentally in the remainder of this appendix.For clarity of the exposition, we start by analyzing a 1D toy example illustrated in 6] , the KLD between the posterior and prior is almost 0, which means the posterior is well approximated by the prior. Therefore, we set l max = 2l SE , and we verify that the KLD between the posterior and the geometric approximation remains small everywhere.Then, we conduct the same analysis to tune l max on real 3D magnetic data. , where 4.05m = 3l SE and 3m is the height of one story in the building. Then, we fit a GP on the training data using the kernel defined in (10) , which embeds prior knowledge from magnetic field physics. The second subfigure of Fig. 6 shows that the KLD between the posterior and the prior is almost zero at a distance of 2l SE or more from the training data. Therefore, we use l max = 2l SE in the cubic hermit spline formulae given above, which defines the β function in the geometric approximation. We then display the KLD between the posterior and the geometric approximation to show that it remains small everywhere, as intended. Here, we introduce some well-known results about definite positive symmetric matrices that will be useful to manipulate covariances and derive the LBCM posterior in appendix E. We use the standard notation S 0 and S 1 S 2 to indicate that the symmetric matrices S and S 1 -S 2 are definite positive.Theorem 1 Let S 1 , S 2 be two symmetric positive matrices with at least one definite positive, and α 1 , α 2 two strictly positive scalars, then α 1 S 1 + α 2 S 2 0.Theorem 2 Let S be a symmetric definite positive matrix, then S -1 0.Theorem 3 Let S 1 , S 2 be two symmetric definite positive matrices, then S 1 S 2 if and only if S -1 2 S -1 1 .Theorem 4 Let S be a symmetric definite positive matrix and V a (possibly rectangular) matrix with full column rank, then A SA 0.In this section, we detail the derivation of the LBMC posterior. It is based on an approximation of the BCM that we can also quickly derive here for completeness.In general, if we have J experts trained on a partition D 1 , . . . , D J of the dataset D, then from Bayes theorem:EQUATIONThe BCM introduces a conditional independence approximation between the D i , which we inject in (19) , and then we use Bayes a second time on each expert posteriorEQUATIONThen, we recall that we inject our geometric approximation in (20) to get the LBCM posteriorEQUATION)where A(x * ) = {i|β i = 0}. This expression comprises a quotient, products, and powers of Gaussian densities. We straightforwardly combine all the Gaussian identities from appendix B. With a zero mean Gaussian prior, the LBCM posterior has the following exponential formEQUATIONwhereEQUATIONFor p LBCM to be a proper Gaussian density N (µ, Λ -1 ), all there is left is to check that Λ -1 is symmetric definite positive. Symmetry is stable by sum, inverse, and multiplication by a scalar, so Λ -1 is symmetric. Showing positive definiteness is slightly more cumbersome. According to theorem 2, we can directly consider Λ. We by rewriting itEQUATIONwhich is a sum of matrices weighted by positive coefficients that would be positive definite if each matrix is positive definite as well (theorem 1). Using the inversion results (theorem 2), we have cov(f * ) -1 0. And according to 3, the matrix cov(f* |D i ) -1 -cov(f * ) -1 is positive definite if and only if cov(f * ) cov(f * |D i ).Intuitively, the prior covariance is larger than the posterior one. Thus, cov(f * ) cov(f * |D i ) should be valid for any sensible expert model, but to finish the formal proof, we need to limit ourself to specific examples. For instance, using full GP experts with kernel κ SE ,EQUATIONwhere Σ = (K f ,f + σ 2 noise I nd ) -1 . The matrix K f ,f * has full column rank and using theorem 1 and 2 we have Σ 0. From theorem 4, it follows that cov(f * ) cov(f * |D i ) holds for each GP expert and LBCM is well defined. If, instead, we use the DTC or the G-DTC sparse approximation defined from this GP covEQUATIONwhereS = (σ -2 noise K u,f K f ,u + K u,u ) -1 . We can see that S -1 = σ -2 noise K u,f K f ,u + K u,u K u,u , thus from theorem 3, we have K -1u,u -S 0. It follows from theorem 4 that cov(f * ) cov(f * |D i ) holds for each experts. Therefore, LBCM is again well-defined.G-DTC is based on extracting a partial grid of latent inputs near the training dataset. Here, we give efficient computation techniques for the partial grid, describe some subtleties about latent inputs near the subdomains border, and motivate the chosen value for the parameter R.Given inputs x 1 , . . . , x n and a radius R, we retrieve all the vertices z 1 , . . . z M from a cubic grid of step δ, that are at a distance R or less from at least one x i . To do so efficiently, for each x i , we generate all the vertices z in a box of center x i and side length 2R (the generation process is described below). Then, we loop through all the vertices of this box to remove those at a distance greater than R from x i . Each remaining z is inserted in O(log(M )) into a set data structure to avoid duplicates. Thus, the partial grid is created in O((R/δ) 3 N log(M )) operations, where (R/δ) 3 is proportional to the number of vertices in a ball of radius R.All that is left is to describe the vertices generation process in a box of center x = (α 1 , . . . , α d ) and side length 2R. Let p 0 = (γ 1 , . . . , γ d ) be the origin of the cubic grid of step δ, then it is possible to describe any vertex z by its grid macro coordinates v ∈ Z d : As a side note, when we split the space Ω in Ω 1 , . . . , Ω J to create the data partition, we allow the partial grid of expert i to go outside Ω i . In other words, each expert's partial grid can overlap at the subdomain boundaries.EQUATIONNow we can focus on the ideal value of R in the case of 3D inputs (d = 3). In particular, we want R large enough such that any input x has neighbor vertices spread in 3D, and we say that R = δ 2 + (δ/2) 2 + (δ/2) 2 is the smallest of such values that work for any input x. For any given input x, LBCM is based on efficiently retrieving to which expert domain x belongs and which are the neighbor experts. We must also compute the distance r from x to any neighbor expert domains. A regular partition of the space into boxes is convenient to compute these quantities in constant time. EQUATION)where a denotes the integral part of a.To get all the neighbor boxes that intersect the ball of center x and radius l, we start by finding all boxes that intersect the box of center x and shape [0, l) × • • • × [0, l). It creates a small list of candidate boxes. Then loop through all of them, compute the distance r according to (29), and keep only the ones such that r ≤ l.All that is left is to compute the candidates list.  
 
 Appendix:  


Decision: Summary of Strengths:
* Magnetic mapping is a relevant topic.
* The concepts presented are well formulated (G-DTC and LBCM) can be applied to other problems.
* The contribution lies in the multiple expert modelling and physics-based learning techniques.
* Evaluation on both simulated data and a real dataset collected by the authors.
* Demonstrated computational advantage and mapping accuracy (including uncertainty) compared to state-of-the-art methods.

Summary of Weaknesses:
* Limitations are not discussed adequately -> this has been addressed in the rebuttal and needs to be added to the paper
* The paper needs a proper related work section, where it differentiates clearly from other state-of-the-art approaches. Noting that similar concepts have been presented already in the literature.  -> A Related Work section should be added (some space can be saved in Section 2).
* Readability. -> to be addressed in the paper
* Reproducibility of the results. -> to be addressed as per discussion in the rebuttal

Reviews: review: title: title: Exposition seems to need a major revision
paper_summary: paper_summary: This paper focuses on magnetic maps of indoor environments, which can be used in indoor positioning and navigation applications. Here, magnetic field is measured by magnetometers for example. Main contributions of this paper are twofolds. First, the paper modifies a Determinsitic Training Conditional spasre GP to include prior knowledge from Maxwell equations. Second, the paper modifies aggregation methods of local Gaussian Processes (e.g., an ensemble of GPs that divides input space) for improved efficiency and continuous predictions. To create large scale map of indoor magnetic field, the paper combines these two techniques. Simulation and real dataset are used to evaluate the proposed methodologies.
main_review: main_review: Strenghts:
- The topic of magnetic mapping is interesting and has potential to complement visual SLAM.
- There are no severe grammar or spelling mistakes in the paper.
- Ideas of leveraging local GPs along with inducing point approximation seems sensible and worthwhile to try.

Weakness:

- (major) The paper does not distinguish existing ideas properly. The paper needs related work section to locate the contributions within the state-of-the-art.

For example, embedding maxwell equation into GP prior is an existing idea [1]. Leveraging only few local neighboring experts for more efficient local GP aggregation is also a known idea [2]. There are several ways to handle discontinuity of local GPs, e.g., [3-4] and the paper misses to differentiate itself from existing works clearly. I think a thourough related work section, especially with a focus on new techniques, can help the paper.

- (major) The paper writing can improve significantly.

It was difficult to parse section 3.1 and 3.2. I could only understanding by reading [5], where several figures and intuition are prsented clearly. Looking back, section 3.1 alone will be difficult to comprehend by many readers. Secondly, the connection from the mapping to the use of GPs is poorly written. There should be one subsection devoted instead of a sentence in section 3. Lastly, a direct derivation of equation 12 and 13 should be presented in the paper or supplementary.

- (minor) Table 1 and table 2 seem to need standard deviations rather than only reporting the mean values.

[1] Modeling magnetic fields using gaussian processes
[2] Trust Your Robots! Predictive Uncertainty Estimation of Neural Networks with Sparse Gaussian Processes
[3] Efficient Computation of Gaussian Process Regression for Large Spatial Data Sets by Patching Local Gaussian Processes
[4] Patchwork kriging for large-scale gaussian process regression
[5] Modeling and interpolation of the ambient magnetic field by gaussian processes
questions: questions: One of the key advantage of local GPs is also in the data uncertainty term. Instead of having one constant noise, each local GP can represent more complex noise terms, which intuitively, will make more sense in mapping applications (each magnetometer might depending on the environment produce different noise characteristics). Why not exploit such characteristics?
limitations: limitations: 1
review_summary: review_summary: My recommendation is a weak reject due to the quality of the presentation. I can guess what is meant, and I could look up related literature to guess what the contributions are, but having to guess make it difficult to judge the main contribution of the paper.
title: title: Interesting research problem, good solution that improves on state of the art, unclear limitations
paper_summary: paper_summary: This paper presents a method for large scale mapping of indoor magnetic fields. The method includes splitting the domain into subdomains, fitting sparse approximations on them, and a local expert aggregation technique. The approach is validated on real and simulated data.
main_review: main_review: The paper addresses an important problem with an original approach, where the infamous intractability of Gaussian Processes is intelligently addressed. The exposition is sufficiently clear given the page limit and the required mathematical background. I recommend making the source code available online, as the method might contain non-obvious implementation details, for reproducibility. The experimental campaign is convincing, although it would have been nice to see an experiment comparing the current approach vs a traditional mapping approach to corroborate the claims made in the introduction. Without this experiment, I am not sure this paper will convince people to try this method instead of just tweaking their favorite traditional mapping library.

I could not find the limitations section, which should have been included (see the official call for papers https://www.corl.org/contributions/call-for-papers). For example, one major limitation is that the current method might not work well online (which is briefly mentioned in the conclusions).
 

Minor comments:

1. Missing reference to the appendix in Section 3.2.2
2. Fix the citation at Line 95
3. Wrong citation in the last column of Table 2?
questions: questions: 1. What are the limitations of the work?
limitations: limitations: 1
review_summary: review_summary: I did not go through the Appendix, but the paper seems solid, hence my recommendation.
title: title: A valuable contribution to robotics research for large scale mapping, applicable beyond magnetic field mapping.
paper_summary: paper_summary: This paper presents a methodology for large scale mapping of indoor magnetic fields. The application to robotics is the use-case of magnetometer-based navigation and localisation, which require a magnetic map of the environment a priori. The methodology is based on a physics-informed sparse Gaussian process regression framework over multiple localised domains, with a novel expert aggregation technique to recover large scale modelling from the local regression. Aggregation is performed such that the result is differentiable and smooth. The entire framework demonstrates computational advantage compared to state-of-the-art alternatives, with competitive mapping accuracy and uncertainty estimation.
main_review: main_review: This is a well written paper with clear structure and presentation of ideas and results.

The paper is significant as it contributes to multiple expert modelling and physics-based learning techniques. The concepts in this paper of Gradient-DTC (G-DTC) and Local Bayesian Committee Machines (LBCM) can be applied outside of the application presented (magnetic field mapping) and to other mapping problems. e.g. G-DTC may prove useful for mapping any curl-free field via sparse Gaussian processes, LBCM may be used in any probabilistic/Bayesian multiple-expert learning framework. Discontinuities at borders of submaps is a common hurdle in large-scale mapping problems in robotics, and LBCM is naturally differentiable in these regions, lending itself well to the problem. 

Strengths:
- strong theoretical contributions that are applicable to a broad range of mapping problems
- well formulated paper 
- evaluation on both simulated data and on a real dataset collected by the authors
- demonstrated computational advantage compared to state-of-the-art methods
- competitive mapping accuracy and uncertainty estimation with state-of-the-art methods

Weaknesses:
- no explicit limitations section and only minor discussion of limitations in text
- a stronger robotics application would be an extension to online mapping as a physical robot collects measurements with a magnetometer.
questions: questions: - The corridor dataset the authors have collected could be another valuable contribution to research -- have you considered making it and the details of its collection available?
- The choice of the $l_{max}$ for the beta function as discussed in Appendix B does not seem rigorous. Fig 5. does demonstrate that $2l_{SE}$ is suitable, but is this true of all datasets and partitions? How transferrable is this analysis to higher dimensional problems and other domains?
- Beyond $l_{max}$, hyperparameter choice and sensitivity in general could be discussed more. i.e. how the domain is partitioned, are boxes always recommended? Does this change with domain shape?

Minor comments:
- Broken reference to appendix in line 198 of manuscript
- "two one order of magnitude slower" in line 279
- should refer to table 1 in line 256
- mathcal Ox missing in first line of Figure 1 caption
limitations: limitations: 1
review_summary: review_summary: I believe this is a good theoretical contribution to addressing mapping problems via mixed-expert learning frameworks, and to incorporating physics-based knowledge into sparse Gaussian process regression. Uncertainty in my recommendation comes from a lack of explicit discussion regarding the limitations of the methodology. Some may also believe this paper has a weaker focus on learning for a physical robotics problem -- it is evaluated on real data but there is no true "hardware" involved in this problem, though it is certainly a problem relevant to robotics..