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: CoRL 2024 (2024)
Keywords: Gaussian process regression
magnetic field maps
indoor localization
TL;DR: We introduce a new large-scale magnetic field map model based on a combination of a sparse approximation of Gaussian process regression and a local expert aggregation method.

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.

Large Scale Mapping of Indoor Magnetic Field by Local and Sparse Gaussian Processes
Iad ABDUI-RAOUFL Vincent GAY-BELLILE' JOLY?
Steve BOURGEOIS
Alexis PALJIC?
Université Paris-Saclay CEA, List, F-91120, Palaiseau, France
"Centre de Robotique, MINES PerisTech, Université PSL, 75006 Paris, France iad. abdulraouf@@ea..f.
Abstract: Magnetometer-baseed indoor navigation uses variations in the magnetic field to determine the robot's location. For that, magnetic map of the environ- ment has to be built beforehand from a collection of localized magnetic measure- ments. Existing solutions built on sparse Gaussian Process (GP) regression do not scale well to large environments, being either slow or resulting in discontin- uous prediction. In this paper, we propose 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 en- vironments 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 each, and then aggregates expert predictions in differentiable and probabilistic way. We evaluate our model on real and simulated data and show that we can smoothly map three-story building in a few hundred milliseconds.
Keywords: Gaussian process regression, magnetic field maps, indoor localization
Introduction
Indoor positioning for robotics is active and challenging field of research [1. 2]. While vision- based approaches [3] are the most common to address this issue, they struggle environments with low 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 sufficiently distinct to differentiate between locations [6] and stable enough to ensure that magnetic field maps remain valid over the years [4].
localization solutions, such [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 Gaus- sian Process (GP) regression [7] from magnetometer measurements can achieve a differentiable and probabilistic representation. However, it has ((N3) computational complexity, with N the num- ber of data points, which becomes intractable for large datasets.
Indoor magnetic-field-baae
To address this, Vallivaara et al. [8] have downsampled the training data and then learned map values using only training data in 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 effec- tively 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. also results in differentiable magnetic predictions that are efficiently com- puted, 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 mod- ified to include prior knowledge from Maxwell equations. Then, neighboring experts" predictions
8th Conference on Robot Learning (CoRI. 2024), Munich, Germany.

---

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 three-story building from trajectory of one kilometer less than 100 milliseconds.
2 Background Gaussian Process Regression
This section introduces the necessary background on Gaussian process regression, sparse approxi- mations, and local expert aggregations. Computational complexities discussed with respect to data imput 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.
2.1 Full Gaussian Process Regression
Let us consider a vector field f(x) Rd defined for any x E C Bdd It must be estimated from a training dataset = {x] vi)....y (æ yN)), where each observation y, of f(x;) (noted fi) is corrupted by Gaussian additive white noise. A classic approach is to model f as zero mean Gaussian process indexed by x [7]
(1) where is the d x d covariance kernel such that k(2z," z²) = cov(f(x), (x')). The choice of 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 nolse from data or manually tuned sensible values. For any test input va, the prior and posterior of f(x.) (noted f-) are Gaussians and admit closed-form expression
J ~ GP(O, R), y f()) +c, e e~ NO, o Ta),
(2) where (Krr, + ond)ss--ndd- The vectors y and f the concatenation (gi ..... 9n and (fI fI). respectively. We also used the standard shorthand notation Ke,s cov(a, b). The covariances Ky.t and K5..t are matrices of size Nd x Nd, and d x Nd respectively. They are defined by block via the covariance kernel (e.g. cov(fs. ff) x(x,, 2,)).
p(f+) = -V Kr.t.). p(fx(D) = = -. KJ... Ky,tKK f.,f)
Training the model is the pre-computation of all quantities in (2) that depend on training data only. Training computational complexity is dominated by the matrix inversion in Odd No) operations. Then a prediction requires (( N²) operations for each input x... We find helpful to rewrite the complexities as function of the volume [ and the number of points per unit volume Pz= = N/[2]- Training and test complexities respectively become (dd p (2 and O (dd ()2)). is intractable for dense data inputs or large environment.
2.2 DTC Approximation
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 approx- imation formalized by [15]. The starting point is a set of latent variables u = (u. called inducing variables in the literature, where M < N. They are values of the GP (like any fi
u
and f.) at the latent inputs 21,..., 2M E Rdd. The latent variables are involved in two successive approximations of the joint prior.
p(J,, u) = p(s,, p(u) = p(Alu) p(flu)
(3)
is a conditional independence assumption between the test variable f. and training variables f. The training conditional further approximated
uf = NKK,, Kr.s - Kr.u Ku.) ~ NK
DIC
0),
(4)
where Ku,.. KF., and Kuf are again shorthand notation for covariance matrices computed by block (e.g. cov (fe, uy) = ((x,,2,)) The approximation (4) assumes the training conditional co- variance cov(f)u) is almost O. 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 posterior (f = NN
noise K K Kn, "_,, k..,f.)) and
(5)
where S = (o noise K Ku.eKr,.u Ku, uu Training has computational complexity an additional Od) operations are required at each test input T...
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 Pa [Q] natural. The training and test complexities become respectively The key takeaway is that sparse approximations scale well with increasing input density Pas but not so much when the environment size |S2] increases. Therefore, is tempting to model the vector field on partition of the input space [10].
Odd |2] P=-p") and O(d"
2.3 Bayesian Committee Machine
As discussed above, splitting the data and fitting (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-dataseta 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].
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:) are Gaussians. The BCM introduces conditional inde- pendence approximation between the D. Combined with Bayes theorem, follows that the BCM posterior has the form of a product of expert regularized by prior term
Let , D,
P(F D:) p(f-)--
f( BCM II..).))
PBCM(f+|D) &
(6)
The BCM posterior Gaussian with frosed-form mean and covariance [16]. is smooth aggrega-
tion 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 prediction from all the experts.
3 Method
We present probabilistic method to learn the magnetic field large scale without discontinuities in the prediction from magnetic vectors y; E R and their respective positions 2; E R For that, 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 principled way so that the contribution of distant experts disappears seamlessly.
3.1 Gradient-DTCC
Inspired by [10, 18], we wish to inject prior knowledge of magnetic field physics in the DTC for- mulation 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]. Un- der such conditions, a scalar potential & exists such that H = Vá. 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 ui...., UMs representing -0(21),.... -((MM) It will reduce the complexity related to the observation dimension
By definition, if ʃ and u are two Gaussian processes defined over Ωy and Q,, and if for every Ti, In E S2 and Z1,..., Zm E Sr the vector (f(x)),... [(x)), ((=1),.... u(z,)) Gaus- sian, then f and are said to be jointly Gaussian, or to form inter-domain Gaussian process. It is
3

---

Ganges
Garges
g-ditc mM
my
(a) input space
(b) G-DTC posterior
(c) partial grid
(d) G-DTC posterior
Figure I: Latent inputs degeneracy problem. (a) Training inputs (black) are placed on the Ox axis. Latent inputs (orange) are spread on this line. (b) With degenerate ID placement in 3D space, the posterior of the field's Oy and Oz components are equal to their prior (fy and fz curves are superimposed). Only the Or component has been learned. (c-d) With con-degenerated latent inputs, the three components of the posterior are learned properly.
completely defined by their respective mean and kernel function (4 (f,,x,) ('),,) (()),,,,) (z, as well as the inter-domain kernel kuy(2, x) = (f](x,=) [10] suggest to define u as a Gaussian process with differentiable mean puu 0 and twice differentiable kernel:
Ke (2,, x²) (sE x²) Se exp 00.5 (1 2/1l//s))))
where ISE and OSE are respectively the characteristic length-scale and amplitude of the process.
Under such regularity conditions, f Vu is well defined Gaussian process and f, u are jointly Gaussian [20] with explicit expressions for the mean and covariance kernels
(r(x)) (x)
0,
(8)
Kyu(1, x²) V Min (I,
- (SE x²) (x)(xrr'))))))) - x')/ KSE(I, x x')(TT (x - xx)(x= (x)=((')' (2g)/)))/23
(9)
K5(x, Vz V of (x.
(10)
We directly substitute the overloaded kernel Ku, Ky, Kfu 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 j.,. Ke,u: Ku,f and Kn Indeed, each block computed from y (I, x²) is a vector of size Ka( x²) is scalar. Compare it against previous block matrices, all computed from x f a') of size d x d. By mixing vectors of size d with scalars, the new training and test complexities are the same as in vanilla DTC except that they are proportional to d 3 instead of d 27.
1, and
The remaining question is whether the DTC sparse process is a good approximation in this gradient setting. Assuming 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) 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, axis. Therefore, the variations of u are modeled in this direction only. Only the first component of Vu can be estimated confidently, and is indeed possible to check that cov (ov(flu) (f o O. This weakness is shared among wate--of-the-an 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.
from (4)
latent inputs 2; E RR.3 are set along the Or
To solve this issue, we set the latent inputs z; so that the scalar potential variations are modeled in all directions. The 2; are vertices of cubic grid of step = Ise/2. We keep only the vertices at distance R or less from at least one training Z;. The parameter R should be as small as possible reduce the number of latent inputs but large enough such that even isolated a; have enough 2; in their neighborhood. We set = 882 + (8/2)² + (8/2)2 1.23 as justified in the appendix F.
3.2 Local Expert Aggregation
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 local version of the BCM. To reduce BCM test runtime, we seek to

---

consider only experts 'nearby' the prediction location 2,, 22], where the "proximity" is defined by the Gaussian process kernel k. The subdomain 92, is 'far' from v, if x(x, x.)) ) is small for all E Ω4 [21]. Using the classic squared exponential (SE) kernel (8), we would like to ignore all the experts farther than a distance mmx set to a few ISE.
3.2.1 A Local Bayesian Committee Machine
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 baaz from S2, we approximate the posterior of the i-th expert by the prior. Inside Ω2, keep the posterior with no further approximations. Everywhere in between, we approximate the posterior of each expert by geometric mean (GM) of the posterior and the prior
p(f+[D:) PGm (P)) & p(f.)D.) (
(11)
where 0 S. (rr)) S 1, and Ti is the distance separating 2, from S4,. When ((r,) = 1, there is no approximation; when (rr;) 0 the posterior is replaced by the prior. We may drop the argument and note it 33 for notional simplicity.
Let A(z.) = {ir)(r;) :) * 0} be the set of index of active experts. By injecting the geometric approx- imation (11) the BCM posterior (6), we define the Local Bayesian Committee Machine (LBCM) posterior
PLBCM x Il' -(f-)- p(f.)= ((f.)) 8.
ILEACE P(
(12)
Then, as developed in the appendix E, the predictive equations are
PLBCM (J.1) =N N NN
-N(N- A-1
Σ
Be cov (..
A-1
(13)
(EA(Z.)
EiEn(a.)) 8. cov The notation E(f- DD)) and cov(f.[2;) are for the mean and covariance of the exact posterior p(f) When all the experts 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 Br. Theirs are not functions of the distance rú. They set 3; equal to the differential entropy between the prior and the posterior of expert since their goal was increase the robustness in the presence of weak experts. On the contrary, we design our 3: to save computation time by imposing that its support is known and bounded around Ω2 We use the cubic hermit spline 3(r) 2(r 3 33r +11 where 0 < r < 1maxx. It is defined such that it is differentiable, even at the transition points O and 4max- The optimal Imax depends on your computational budget and quality requirements. In appendix we illustrate 3 and derive the exact expression of PGM- It is a Gaussian, which allows us to compute its Kullback-Le Leibler Divergence from the true local expert posterior and justify that Imax 21ge is a reasonable choice.
where A = (1 CiE
Our final model defines a G-DTC expert on each of the subdomains of volume V, and it ag- gregates their predictions with LBCM. The training and test complexities are V p²) and (d) A (z²2))²²>)) respectively, where [A(x.)] is the number of active experts. Computing predic- tions of the magnetic field became independent from the input density and the scene volume.
3.2.2 Domain Partition
LBCM is based on BCM, whose performance depends the domain partition. With smaller subdo- mains, each expert is faster. However, the conditional independence assumption in (6) is challenged the frontier between domains. The predictions better at the interior of a domain, and therefore, the domains must be large enough with respect to ISE. 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, L1) X -.. x [0, La Furthermore, we advise that min(L;) > Umax so that the number of active experts [A(x..) < 2d' with equality when is near the corner of sub- domain. In three dimensions AA(x.) <= 8. Using cuboids makes straightforward to determine the active experts in constant time and to compute the distances r. Details are in the appendix G.
5

---

40
T0
20
Figure 2: Test trajectory of Corridor and the norm of the map learned from the training data
4 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 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.
4.1 Experimental Protocol
Datasets. Here introduce three small simulated datasets called SimulD, Simu2D and Simu3D illustrated Fig. 4 and a large-scale real one called Corridor shown Fig. 2. All datasets are available in our repository httts://githube
The simulated datasets are in a 3 x 3 x m box, typical of local expert's domain. Each dataset comprises 1000 noisy training data points and 100 test data points. Simu/ is representative of dataset acquired with motion constrained by the environment, such as in a narrow corridor. All positions are within a tube of 0.5 meter in diameter. Simu2D represents approximately planar motion in open spaces. All positions are within box of height 0.5 meter. Simu3D represents aerial motion, with positions spread in 3D. Observations are generated jointly from zero mean Gaussian process as described in more detail in appendix A.
Real data in Corridor were collected by the magnetometer and the IMU in an SBG-EllipseNN 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 observa- tions, respectively, in 1400 m2 of corridors and staircases in a three-story building. Our model uses a 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 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 C++ implementation. We used one core of 11th Gen Intel CoreTM 17-118000 with 16 GB of RAM. Prediction times refer to processing the entire test dataset, not just one T..
G-DTC And LBCM Setup. G-DTC depends on the hyper-parameters 4gF, OSE. and Nenise- They have the same physical meaning as in a full GP using the covariance kernel (10). Therefore, on Simu/ Simu2D, and Simu3/) use the available ground truth. Despite our best efforts, unaccounted-forr localization errors remain on Corridor. We reflected in the sensor noise by setting small area and optimized a GP log marginal likelihood to tune ISE and osg. We get ISE = 1.35 m and aSg 6.9
Tnoiee 4 4'T. Then, we randomly selected 1000 training observations
LBCM depends on the domain partition and 4maxe- For reasons that will become clear after the experiment 4.2, we set the side length of the boxes to be equal to 3/Sg, except when the data covers multi-story building. In this case, we exploit the natural floor partition and decrease the box height to the height of a story (- meter). Finally, we take mmax 21. as justified section
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 space partition,

---

SimulD
Simu2D
Simu3D
G-DTC
G-Hilbert
G-DTC
G-Hilbert
G-DTC
G-Hillert
MSE
7.7x10-5
7.8x10-8
1.9x10-4
1.9x104
4.9x10
5.1x10-4
MSLL
-11.8
-12.3
-11.6
-11.7
-11.1
-10.9
Fit (ms)
8.2
131
29
132
111
135
Predict (ms)
2.2
28
8.7
28
43
28
M
91
512
211
512
477
512
Table Comparison of G-DTC and G-Hilbert using datasets defined near ID,2D, and 3D manifolds.
-2
20
20
20
10
10
10
-4
magnetic field (uT
0
0
y (mmtte)
-6
-10
-10
-10
f
-20
-0
-0
F
-8
-30
-30
30
expert training test
440
-40
40
-10
50
50
550
-2
(meter)
distance (m)
distance (m)
distance (m)
(a) inputs boxes
(b) Naive
(c) Overlap
(d) LBCM
Figure 3: Regression discontinuities at the borders. The dotted lines are the test observations used as ground truth. When the test trajectory crosses domain borders, simple aggregation methods are discontinuous, whereas LBCM is differentiable by definition.
and we use the same as for LBCM. Furthermore, each local expert extends its domain by range train to include training data from neighbor experts, such that the "training boxes" overlap. We set herain 0.5/gg as in [11]. G-Hilbert relies basis functions analogous 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 mean GP prior assumption. Thus, we remove these three additional basis functions for fair comparison. Because of mathematical details, G- Hilbert must be defined on larger domain than the training boxes. We extend each training box by 2|gr in [11]. Finally, G-DTC and G-Hilbert share the same values for (SF, USE, and Tnoise-
4.2 Evaluation of G-DTC Efficiency
In this first experiment, illustrate G-DTC efficiency when training inputs are located a manifold of lower dimension. For that, we evaluate G-DTC and G-Hilbert on SimuID, Simu2D, Simu3l) and summarize the results in table 1. On Simu/ G-DTC has similar prediction quality and is one order of magnitude faster than G-Hilbert. Increasing the manifold's dimension in Simu21) and then in Simu3I) decreases have at most 8 x 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 staircase is ID.
speed advantage. In this box of size 3lS 3/gE 31SE, could
4.3 Smooth Local Experts Aggregation In Large-Scale Regression
We proved the pertinence of G-DTC on a small domain. Now, we evaluate the scalability and re- gression quality of LBCM combined with G-DTC. We compare LBCM against Overlap and against

---

G-DTC & Naive
G-DTC & Overlap
G-DTC & LBCM
G-Hilbert & Overlap [10]
MSE (uT)²
1.45
1.14
1.17
1.17
MSLL
-5.38
-5.35
-5.37
-4.70
Fit (ms)
65
135
64
4061
Predict (ms)
108
170
578
4738
Table Aggregation methods evaluation large building with 16600 test predictions.
the naive aggregation method that uses the space partition as is, without any smoothing strategy. We also compare full model G-DTC & LBCM against [10] for completeness.
Fig. shows the map built from our model, and table 2 exposes the quantitative results. Surprisingly, the MSE and MSLL 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. matters because such discontinuities happen every few meters at each domain transition. Using LBCM comes the cost of slightly slower predictions because it aggregates 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 orders of magnitude slower 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 differeentiaiiityy
5 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 r 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 would be more rigorous to explicitly model input noise [26, 27).
Moreover, we low number of latent points to improve runtimes, but it can lead under- confident G-DTC predictions on dense datasets, especially when the latent point layers are thin, such in DD and Simu2D.
Another limitation is in the number of new parameters introduced. There is mmax and the subdomain size and shape for LBCM, and G-DTC requires step size 8 and 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 large datasets. Still, we do not provide a way to scale the tuning 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.
our model hyper-parameter Tnoises OSE, and İSR. Instead, we
We introduced an efficient approach that creates 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 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 differentiable way. Efficiency and differentiability are key properties required to develop graph SL.MM localization algorithms. For that, online update strategies and robustness to training localization error may be the focus of future research.
Conclusion:
8

---

Acknowledgments
We thank Olivier GOMEZ for its technical support in dealing with the implementation dependencies and cross- cross-platform building. We also thank the reviewers who gave valuable comments improve the publication's s clarity. This work was supported by the CEA List.
References
[1] M. Kok and A. Solin. Scalable Magnetic Field SLAM in 3D Using Gaussian Process Maps. 21 International Conference on Information Fusion (FUSION), 2018.
[2] Z. Wu, W. Wang. J. Zhang. Y. Wang, G. Peng, and D. Wang. MGLT: Magnetic-Lead Global Localization and Tracking in Degenerated Repetitive Environments. JEEE/ASME Transactions on Mechatronics, 2024.
[3] C. Campos, R. Elvira, J. J. G. Rodríguez, J. M. M. Montiel, and J. D. Tardós. ORB-SLAM3. An Accurate Open-Source Library for Visual, Visual-Inertial, and Multimap SLAM. JEEE Transactions Robotics, 37(6), 2021.
[4] Coulin, R. Guillemard, V. Gay-Bellile, C. Joly, and A. d. L. Fortelle. Tightly-Coupled Magneto-Visuall Inertalln Fusion for Long Term Localization in Indoor Environment. IEEE Robotics and Automation Letters, 7(2), 2022.
[5] G. Ouyang and K. Abed-Meraim. A Survey of Magnetic-Ficee Saaettic-Fiel--Basee Beeld-Based Indoor Localization. Electronics, 11(6), 2022.
[6] Abdul Raouf, V. Gay-Bellile Bourgeois, C. Joly, and Paljic. MagHT: A Magnetic Hough Transform for Fast Indoor Place Recognition. In 2023 IEEE/RSJ International Confer- ence on Intelligent Robots and Systems (IROS). IEEE, 2023.
[7] C. E. Rasmussen and C. K. I. Williams. Gaussian processes for machine learning. Adaptive computation and machine learning. MIT Press, 2006.
[8] Vallivaara, Haverinen, A. Kemppainen, and J. Roning. Simultancous localization and mapping using ambient magnetic field. In 2010 IEEE Conference on Multisensor Fusion and Integration. IEFE, 2010.
[9] H. Liu, Y.-S. Ong, X. Shen, and J. Cai. When Gaussian Process Meets Big Data: A Review of Scalable GPs. IEEE Transactions on Neural Networks and Learning Systems, 31(11), 2020.
[10] A. Solin, M. Kok, N. Wahlström, T. Schön, and Särkkä. Modeling and Interpolation of the Ambient Magnetic Field by Gaussian Processes. IEEE Transactions on Robotics, 34(4), 2018.
[11] C. Menzen, M. Fetter, and M. Kok. Large-scale magnetic field maps using structured ker- nel interpolation for Gaussian process regression. In 2023 26th International Conference on Information Fusion (FUSION), 2023.
[12] A. Datta, S. Banerjee, A. O. Finley, and A. E. Gelfand. Hierarchical Nearest-Neighbor Gaus- sian Process Models for Large Geostatistical Datasets. Journal of the American Statistical Association, 11 111(514), 2016.
[13] L. Wu, G. Pleiss, and J. P. Cunningham. Variational nearest neighbor Gaussian process. In Proceedings of the 39th International Conference on Machine Learning. PMLR, 2022.
[14] A. Wilson and H. Nickisch. Kernel Interpolation for Scalable Structured Gaussian Processes (KISS-GP). In Proceedings of the 32nd International Conference on Machine Learning. PMLR, 2015.
[15] R. Herbrich, C. E. Rasmussen, and J. Quiñonero-Candela. A Unifying View of Sparse Ap- proximate Gaussian Process Regression. Journal of Machine Learning Research, 2005.
[16] V. Tresp. A Bayesian Committee Machine. Neural Computation, 12(1 1), 2000.

---

[17] F. Leibfried, V. Dutordoir, S. T. John, and N. Durrande. A Tutorial on Sparse Gaussian Pro- cesses and Variational Inference, 2022. arXiv:2012.1396..
[18] N. Wahlstrom, M. Kok, T. B. Schon, and F. Gustafsson. Modeling magnetic fields using Gaussian processes. In 2013 IEEE International Conference on Acoustics, Speech and Signal Processing. IEEE, 2013.
[19] M. Lázaro-Gredilla and A. Figueiras- -idda.. Inter-domain Gaussian Processes for Sparse In- ference using Inducing Features. In Advances Neural Information Processing Systems, vol- ume 22. Curran Associates, Inc., 2009.
[20] R. J. Adler. The Geometry of Random Fields. SIAM, 1981.
[21] D. Nguyen-tuone J. Peters, and M. Seeger. Local Gaussian Process Regression for Real Time Online Model Learning. In Advances Neural Information Processing Systems, volume 21. Curran Associates, Inc., 2008.
[22] J. Lee, J. Feng, M. Humt, M. G. Müller, and R. Triebel. Trust Your Robots! Predictive Uncertainty Estimation of Neural Networks with Sparse Gaussian Processes. In Proceedings of the 5th Conference Robot Learning. PMLR, 2022.
[23] M. Deisenroth and J. W. Ng. Distributed Gaussian Processes. In Proceedings of the 32nd International Conference on Machine Learning. PMLR, 2015.
[24] M. Seeger, editor. Bayesian Gaussian Process Models: PAC-Bayesian Generalisation Error Bounds and Sparse Approximations. University of Edinburgh, 2003.
[25] L. Csató and Opper. Sparse On-Line Gaussian Processes. Neural Computation, 14(3), 2002.
[26] A. Mchutchon and C. Rasmussen. Gaussian Process Training with Input Noise. In Advances in Neural Information Processing Systems, volume 24. Curran Associates, Inc., 2011.
[27] A. C. Damianou, M. K. Titsias, and N. D. Lawrence. Variational Inference for Latent Variables and Uncertain Inputs in Gaussian Processes. Journal of Machine Learning Research, 2016.
10

---

X train 50
150
100
Z inducing
0.75
075
d:0
0.00
0.00
0:75
0.75
1.50
.50
1.30
1.50
250
0.75
0.75
-1.50 -0.7 0.00
00.75
50
-1.50
-0775
0.75
0.75
136
0.00
0.75
13650
0.00
0.75
15650
(a) SimulD)
(b) Simu21)
(c) Simu3D
Figure 4: Simulated dataset and G-DTC latent inputs. G-DTC uses a partial grid that creates fewer latent inputs (also called inducing inputs) SimulD and Simu2D because the training data cover only a fraction of the complete domain.
A Simulated Datasets
The three simulated datasets SimulD, Simu2D), Simu3D displayed in Fig. 4 are, in fact, just 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 ISE = and OSE = 1. Gaussian white noise also added the training observation of each dataset, with Ooise 0.1.
B Gaussian Identities
This appendix provides useful identities for manipulating Gaussian distributions and computing the Kullback-Leibler divergence.
Let p(x) (z Eo) be the density of a multivariate Gaussian vector of dimension D, and a positive scalar. After re-normalization, the power of Gaussian remains Gaussian
1 const p(x) =N(<ppp -x))-
(14)
Furthermore, if (x))= = (r) Σi) is another density of multivariate Gaussian vector, then after re-normalization, the following product is Gaussian as well
const p(=)g(z) (zlpp,, Ez),
(15)
where
#22 Σ2 H0 + Eg (Io + +
Their quotient also admits Gaussian expression
p(x) ((x)
x exp
- -- - H3)
1 (a)
"gg (x ps)
(16)
where
#43 Es 40 - Σ3= -E
Notice that Eg can be non- non-positive, making this quotient unsuitable for a random variable density. However, if Eg is definite positive, this expression be normalized into a proper Gaussian density ((H, Eg).
Working with Gaussian densities is convenient since, among many other reasons, their Kullback- Leibler divergence admits a closed-form expression
Σ En
DKL. (p)lg) 2
log
+ (uo - -11 14) " (No - pi) tr(
(17)
where the notation |A| and tr(A) mean respectively the determinant and trace of matrix A.

---

beta
train posterior
Du (post|jorio)) (post]joM)
GM
0
postertor
2
-10
prior
1
-20
2
-30
-2
16
11
31
Figure 5: Illustration of the geometric approximation on 1D toy dataset. On all figures, dotted vertical lines delimit the transitions from / = to 3 = 0. The geometric approximation KLD is low everywhere, with a small peak toward the end of the transition zone. Small peaks translate into small posterior mean errors, remaining within the la-error curves of the true posterior. It also translates into higher uncertainties.
C Choosing The Beta Function
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 B. To the best of our knowledge, it is the first time RBCM viewed as an approximation of each posterior by geometric mean between the posterior and prior. We use this theoretical reinterpretation to derive a new metric on the approximation quality with respect to 3.
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-Leible Divergence (KLD) of the geometric approximation compared to the true posterior: DKL. (p(f- (D.) (cM(f.|P:)) To compute it, we need an explicit expression for (f Di)) 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 get
PGM (f. -A (1:00~~88,0o0//00 (f. (F + (f)))oo(()- A-!) (18)
where A (f (1- (f))oo((f))) Then, the KLD between two Gaussian distribu- tions admits a closed-fone form expression.
We can now experiment with the 3 function to choose one for which the KLD remains small ev- erywhere. For differentiability, we restrict our-self to the cubic hermit spline ((r) 2(r/Imax)
3(r//max)2 ) + 1, where is the distance separating from the expert domain. 3 is defined up to hyper-parameter 4max that we will tune experimentally in the remainder of this appendix.
For clarity of the exposition, we start by analyzing a ID toy example illustrated in Fig 5. The dataset comprises 400 observations on the segment 52, = [0, 6). They are sampled from Gaussian process prior defined by KSE and additive white noise. crpper-parametry are set to [sg 1. OSE = 13.4 and Ouoise 1.4. The KLD between the prior and the posterior is high in the training domain because the posterior incorporates lot of information from training data. However, at a distance 2|sF. from [0, 6), the KLD between the posterior and prior is almost 0, which means the posterior well approximated by the prior. Therefore, we set /max 21sE, and we verify that the KLD between the posterior and the geometric approximation remains small everywhere.
Then, we conduct the same analysis to tune mmx on real 3D magnetic data. Fig 6 displays the dataset extracted from Corridor. We select all the observations made inside cube of size 4.05 4.05 m3, where 4.05m = 31sg 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. shows that the KLD between the posterior and the prior is almost zero at a distance of 2/SE or more from the training data. Therefore, we use mmax = 2/gE in the cubic hermit spline formulae given above, which defines the 3 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.
12

---

beta
21
0.
45
0.6
0.4
-8
02 -8
00
(a) training data and beta
(b) KLD post. prior
(c) KLD post. vs geometric
Figure 6: Illustration of the geometric approximation on a real dataset. All figures are top views of the 3D space. The red and green squares delimit the transitions from 3 to 3 O. The figure on the leftt displays the training dataset of local expert, well as 3. The figure in the middle shows KLD between the posterior and the prior (saturated to for visualization purposes, but goes up to 43), which is almost zero at distance of 21sg more from the training data. The figure the right shows KLD between the posterior and the geometric approximation. is small everywhere and peaks at 3 toward the end of the transition zone.
D Definite Positive Symmetric Matrices
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 S2 to indicate that the symmetric matrices S and S S2 are definite positive.
Theorem Lei S1, S2 be two symmetric positive matrices with at least one definite positive, and 01, 02 two strictly positive scalars, then 01 + 0:22 0.
Theorem Let S be symmetric definite positive matrix, then S- >
Theorem Let S1, S2 be two symmetric definite positive matrices, then S1 > Sy if and only if ST > -ST
Theorem Let S be a symmetric definite positive matrix and V a (possibly rectangular) matrix with full column rank, then A SA > O.
E Derivation of LBCM Posterior
In this section, detail the derivation of the LBMC posterior. It is based approximation of the BCM that we can also quickly derive here for completeness.
In general, if we have J experts trained on partition Dı,..., DJ of the dataset D. then from Bayes theorem:
pf
D)
Hayes
p(P)S.)p(())
Py1S-)p(f.).
(19)
The BCM introduces a conditional independence approximation between the D;, which we inject in (19), and then we use Bayes second time each expert posterior
p(f. D.).
P(A.
BCM p(D.,...,PrlAA)) DrI)))p(s.) & p(J.) II ((D.A-)
Bayes & p(S.) II
(20)
x(-))
p(f.)=-1
Then, we recall that we inject our geometric approximation in (20) to get the LBCM posterior
(..C))(.lPD) I1' ((.
IliE (LEEA(..)P(..(...) P(S.
(21)
P(1.)-11
p -1+2, -1 B.
13

---

where A(x.) [i]3, #0}. This expression comprises 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 form
PLBCM(Sa|D) x exp (-2(r - 1)) =) A(x (1)
(22)
where
== A-1
Σ
Be cocovo (E(ID))
tEA(x.)
(23)
A (1- 8.)cov(fa)-
E (EA(Z.)
cov(ff
)EA(x.)
For PLBCM to be proper Gaussian density Nu,, A[*), all there is left is to check that A is symmetric definite positive. Symmetry stable by sum, inverse, and multiplication by a scalar, so A is symmetric. Showing positive definiteness is slightly more cumbersome. According to theorem 2, we can directly consider A. We start by rewriting it
A = = +
Σ (EA(Z.)
Be (cov(fx - cov(f.)=
(24)
which 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 cfv(f.) O. And according 3, the matrix cov(f. -ccofff- is positive definite if and only cov(f.) > cov (f.
Intuitively, the prior covariance is larger than the posterior one. Thus, cov (f-) > cov(J- should be valid for any sensible expert model, but to finish the formal proof, we need limit ourself to specific examples. For instance, using full GP experts with kernel KōŞE.
cov - cov = K EKr.j.
(25)
where Σ = (Kr,r + naa) 2 we have > O. From theorem 4, it follows that cov(f.) cov(f.|D) 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
1
The matrix Kr,s. has full column rank and using theorem I and
cov(f-) - cov = KT,,, (K S)Ku.,.
(26)
where S = (a noise K KurrKr,u + Ku,u) -1, We can see that S-1 o moise Ku,kKr, Ku Ku thus from theorem 3, we have K-1 S > O. It follows from theorem 4 that cov(f.) cov(f. holds for each experts. Therefore, LBCM is again well-defined.
F Partial Grid
G-DTC is based on extracting partial grid of latent inputs the training dataset. Here, 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 Ii,..., Tn and a radius R, we retrieve all the vertices .I.... ZM from a cubic grid of step 5, that are at distance R or less from at least one To do efficiently, for each we generate all the vertices z in a box of center I; and side length 2R (the generation process is described below). Then, we loop through all the vertices of this box to remove those distance greater than R from 2, Each remaining 2 is inserted in O (log(AM)))) (M) into a data structure to avoid duplicates. Thus, the partial grid is created in (((R)8) 8 'N log log(M)) operations, where (R/8)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 = (ai, side length 2R. Let Po = (Yi...,, .a) be the origin the cubic grid step 8, then it is possible to describe any vertex & by its grid macro coordinates v E zd
and
= --x-tp
(27)
14

---

0.0
L
(mmeeer)
2 00.0
0.0
085
0.2 (nveter)
0.2
(mmtter
0.1
0.9
(a) R
(b) R
(c) R = 3/2 x
(d) Worse case
Figure 7: Vertices selection from one cell of the cubic grid. The selected vertices (orange) are extracted from a cubic grid (gray) at distance I or less than an input x (blue). We represent one grid cell only for clarity. (a) R is too small, and only
vertex is selected. (b) R is still too small because the four selected vertices are coplanar (beware, the 4 gray vertices on the back are outside the (chere)(c) R is larger, and the six selected vertices not coplanar anymore.
(d) In the worst case, the input located the middle of a face, and R should be greater than the red segment length /8a + (6/2)² + (6/2)² = 33/2 x 6.
The lowest and highest macro coordinate values all the vertices inside the box are respectively
R 7
a R Va
Umin
(28)
+ R n
O.d' R Va'
Umax
where [a] and [a] denotes the ceiling and integer part of scalar respectively. Then [Umin, Umax,1 x ..... x [[min,d', Umaxxd" is the list of grid macro coordinates of the vertices inside the box.
As a side note, when we split the space Q Q,....,22, , to create the data partition, we allow the partial grid of expert to go outside Qz In other words, each expert's partial grid can overlap the subdomain boundaries.
Now we can focus on the ideal value of R in the case of 3D inputs (d' 3). In particular, we want large enough such that any input x has neighbor vertices spread in 3D, and we say that R / 82 + (8/2)² + (8/2)² the smallest of such values that work for any input v. Fig. illustrates it.
G
Box Partition
For any given input x, LBCM is based efficiently retrieving to which expert domain 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.
Let (71, ... Na) be the center of a box of shape [0, L2) x [0, Lvr.. Then its distance from any = (a) (.)...., axr is
Ta
max [a: - 7v
(29)
where the notation |a] stands for the absolute value of a scalar a.
Fach box can be identified by its macro coordinates' h E zdd which represents how many boxes should be skipped in each direction from origin box b (0...., 0), defined arbitrarily. If is the center of the origin box, then the macro coordinates of the box containing are
(1
a1 71
ad'' 7d )
b(x) =
(30)
L
La
15

---

where |a] denotes the integral part of a.
To get all the neighbor boxes that intersect the ball of center and radius I, start by finding all boxes that intersect the box of center I and shape [0, x ... x [0, It creates a small list of candidate boxes. Then loop through all of them, compute the distance according (29), and keep only the ones such that r S L.
All that is left is to compute the candidates list. For that, compute b(x X) and b(x + X) where À = (1, ....... 1). Then [b(x- - A)1, b(x + A)1] x .... [b(x A)ar, b(x + X)ar] is the list of candidate macro coordinates, where the double brackets [[,, b] stands for the set of integers between and b included.
16


Reviews:
[{"review_id": "BKZ9DMVwwS", "reviewer": "Reviewer_XHZx", "summary": "", "questions": "", "limitations": "", "rating": 2, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": 2, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "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.", "summary_of_recommendation": "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.", "technical_quality": 2, "clarity_of_presentation": 1, "potential_impact": 2, "robotics_focus": 3, "quality_of_limitations": 1, "questions_for_rebuttal": "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?", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "yrvjwMwJaO", "reviewer": "Reviewer_KGGy", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": 4, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "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.", "summary_of_recommendation": "I did not go through the Appendix, but the paper seems solid, hence my recommendation.", "technical_quality": 4, "clarity_of_presentation": 4, "potential_impact": 3, "robotics_focus": 3, "quality_of_limitations": 1, "questions_for_rebuttal": "1. What are the limitations of the work?", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "84qlpLxq3w", "reviewer": "Reviewer_MA22", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": 3, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "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.", "summary_of_recommendation": "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..", "technical_quality": 4, "clarity_of_presentation": 5, "potential_impact": 4, "robotics_focus": 3, "quality_of_limitations": 1, "questions_for_rebuttal": "- 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?\n- 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?\n- 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?\n\nMinor comments:\n- Broken reference to appendix in line 198 of manuscript\n- \"two one order of magnitude slower\" in line 279\n- should refer to table 1 in line 256\n- mathcal Ox missing in first line of Figure 1 caption", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: 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
Decision: Accept