Modeling and interpolation of the ambient magnetic field by Gaussian processes

I. INTRODUCTION

• Magnetic maps can support indoor positioning without extra infrastructure or line-of-sight, but accurate mapping at magnetometer sampling rates makes GP computational cost a central challenge.The cited prior work downsampled data because GP complexity scales cubically with the number of observations.Source: §I · p. 1, §II · p. 2
• The method places a GP prior on a latent magnetic potential, uses a Hilbert-space representation, models Earth-field bias, and supports sequential online estimation and possible time-varying maps.The paper motivates temporal mapping with changes such as furniture being moved.Source: §I · p. 2
• The paper reports feasibility and accuracy in simulations and a small-scale experiment, and applicability in mapping experiments with a mobile robot and hand-held smartphone.For time-varying fields, it builds on a sequential GP approach solved with Kalman filtering.Source: §I · p. 2, §II · p. 3

II. RELATED WORK

• Earlier magnetic maps ranged from grid interpolation and spatial binning to GP models; Wahlström et al. incorporated Maxwell-based curl- and divergence-free priors.Some approaches used magnetic-field intensity, while others used full three-dimensional measurements.Source: §II · p. 2
• This paper builds on Wahlström et al. by using full three-dimensional magnetometer measurements and encoding physical knowledge through the magnetic-field potential.The introduction frames this as a GP prior on a latent, unobservable magnetic potential function.Source: §II · p. 2, §I · p. 2

III. THE AMBIENT MAGNETIC FIELD

• The paper assumes no free current in the indoor region of interest, making the magnetic field curl-free there and allowing it to be represented by a scalar potential.The text identifies metallic structures, rather than free currents in wires, as the major source of ambient-field variation in most indoor environments.Source: §III · p. 3
• Modeling the scalar potential instead of the magnetic field directly implicitly imposes physical constraints on the modeled field.A scalar potential is associated with the curl-free region; the paper excludes a vortex point where this condition fails.Source: §III · p. 3

A. Gaussian process regression

• GP regression treats functions as draws from a prior defined by a covariance function and uses noisy observations to compute the posterior process at test inputs.Because the prior and measurement model are Gaussian, the posterior is Gaussian.Source: A. Gaussian process regression · p. 4
• At test inputs, GP regression computes closed-form conditional means and variances; the marginal likelihood can also support inference of covariance-function and noise parameters.The covariance matrix uses pairwise covariances of training inputs, while k* contains covariances between a test input and training inputs.Source: A. Gaussian process regression · p. 4
• The covariance function encodes prior knowledge about the process; the squared exponential covariance is stationary and isotropic, with magnitude and length-scale parameters learnable from data.The paper names the squared exponential as also known as exponentiated quadratic, radial basis function, or Gaussian.Source: A. Gaussian process regression · p. 4

B. Interpolation of magnetic fields

• Interpolation predicts the latent, noise-free three-dimensional field at unmeasured locations from spatial measurements of the magnetic field corrupted by independent Gaussian noise.The measurements are vector-valued, so components may be modeled separately or linked by an assumed relationship.Source: B. Interpolation of magnetic fields · p. 4
• The magnetic-field mean lies around the local Earth field and can deviate indoors because of magnetic material; an additional covariance component can model the unknown mean.Source: B. Interpolation of magnetic fields · p. 4

C. Separate modeling of the magnetic field components

• The separate-component GP captures the general field variation along route D–E, with the simulated ground truth inside its 95% credibility interval, but encodes no magnetic-field physics.The model is described as flexible and conservative; the scalar-potential model is introduced as the alternative that incorporates physical knowledge.Source: Figure 4 · p. 5, C. Separate modeling of the magnetic… · p. 5
• Each separate GP models a smooth spatial component with a constant shift from zero mean, using a constant covariance term for the mean and a squared exponential term for small-scale variation.Source: C. Separate modeling of the magnetic… · p. 5
• Independent components have no encoded cross-component relationship, which can make hyperparameter estimation converge to local optima; practitioners may instead fix hyperparameters to reasonable values.Source: C. Separate modeling of the magnetic… · p. 5
• Jointly learning shared hyperparameters for three independent GPs requires inverting an n-sized matrix for marginal likelihood, rather than a 3n-sized matrix.Source: C. Separate modeling of the magnetic… · p. 5

• The model assumes the magnetic field is the gradient of a scalar potential, with magnetometer readings treated as its gradients corrupted by Gaussian noise.Source: D. Modeling the magnetic field as the… · p. 5
• The squared exponential covariance models small-scale anomalies, while a linear potential component represents the local Earth magnetic field.Source: D. Modeling the magnetic field as the… · pp. 5–6
• Because differentiation is linear and preserves Gaussianity, the model can use derivative measurements for GP regression and learn its four hyperparameters by maximizing marginal likelihood.Source: D. Modeling the magnetic field as the… · p. 6
• A GP prior on the potential yields a curl-free kernel whose samples obey the magnetic field's curl-free constraints, enabling physical structure in the model.The scalar-potential formulation is described as crucial to extending the model to efficient approximate GP inference.Source: D. Modeling the magnetic field as the… · p. 6
• The scalar-potential GP improves interpolation accuracy and produces a narrower 95% credibility interval than the independent-component GP in the route D–E example.Source: D. Modeling the magnetic field as the… · p. 6

V. EFFICIENT GP MODELING OF THE MAGNETIC FIELD

• Naive GP inference scales as O(n^3) and is described as impractical for large datasets, around n > 1000; the proposed approach uses differential structure and basis projection to reduce this burden.Source: §V · p. 6
• The framework distinguishes batch estimation, sequential updates of a static field, and spatio-temporal estimation when the field changes over time.The online extension is intended to update estimates as data arrives without repeating batch computations after each sample.Source: §V · p. 6

A. Reduced-rank GP modeling

• Truncating the covariance eigenfunction expansion at m ≪ n gives GP regression complexity O(nm^2), hyperparameter-learning complexity O(m^3), and memory use O(nm).The eigenfunctions are independent of covariance hyperparameters and can often be solved in closed form for many domains.Source: A. Reduced-rank GP modeling · p. 6
• For cuboid domains with Dirichlet boundary conditions, the basis functions and eigenvalues have analytic expressions and need only be evaluated once because the basis is hyperparameter-independent.Source: A. Reduced-rank GP modeling · pp. 6–7
• The reduced-rank method uses Laplace-operator eigenfunctions, which fit the scalar-potential formulation because observations are gradients of the latent potential.Source: A. Reduced-rank GP modeling · p. 7

• The reduced-rank prediction equations use basis functions evaluated at training and prediction locations, with the spectral-density matrix and vectorized observations.Source: B. Batch estimation · p. 7
• The batch workflow evaluates the gradient basis at observed locations, forms the regression quantities, and solves the GP regression problem for predictions.Hyperparameters are optimized using the marginal likelihood; its remaining covariance-function dependence is in the diagonal spectral-density matrix.Source: B. Batch estimation · pp. 7–8
• The batch algorithm takes spatial locations and magnetic readings, prediction locations, domain boundaries, and approximation degree m, and returns marginal means and variances of the predicted field.The scalar potential can also be returned using the potential basis at prediction locations.Source: B. Batch estimation · p. 8
• Prediction scales as O(nm^2), marginal-likelihood evaluation as O(m^3), and memory as O(nm), compared with O(n^3) for the naive full GP.Source: B. Batch estimation · p. 8

C. Sequential estimation

• The sequential formulation is a linear Gaussian estimation problem and matches the batch solution within numerical precision.The formulation is widely known as the Kalman filter.Source: C. Sequential estimation · p. 8
• The algorithm initializes the state from the GP prior, updates it for each observation, and can return a current prediction after each iteration.The algorithm treats hyperparameters as known and evaluates predictions at a test location.Source: C. Sequential estimation · p. 8
• Hyperparameter optimization is not included; the authors suggest optimizing with an initial batch and re-optimizing later if necessary.The marginal likelihood can be evaluated through the recursion.Source: C. Sequential estimation · p. 8
• Sequential computation scales as O(m^3) per update and O(nm^3) overall, with O(m^2) memory if intermediate results are not stored.Source: C. Sequential estimation · p. 8

D. Spatio-temporal modeling

• The sequential model can track dynamic field changes with virtually no additional computational burden, and the spatio-temporal model has the same computational costs as sequential modeling.Source: D. Spatio-temporal modeling · pp. 8–9
• The spatio-temporal GP places a prior on scalar potential ϕ(x, t), extending the model to location and time.Source: D. Spatio-temporal modeling · p. 8
• A temporal covariance function specifies assumptions about temporal behavior, with ℓtime controlling the temporal length-scale.Source: D. Spatio-temporal modeling · p. 8
• The temporal model assumes a continuous but not necessarily differentiable process, identified as an Ornstein–Uhlenbeck process.Source: D. Spatio-temporal modeling · pp. 8–9
• Temporal effects are assigned only to the anomaly component; the bias is tracked as static.Source: D. Spatio-temporal modeling · p. 9

VI. EXPERIMENTS

• The four experiments cover simulated data, a small-scale proof of concept, handheld smartphone mapping in a building, and online robot mapping with real-time tracking of field changes.Source: §VI · p. 9
• The simulation study permits comparison with known ground truth, while the later experiments demonstrate mapping in practical settings.Source: §VI · p. 9
• The methods were implemented in MATLAB on an Apple MacBook Pro with a 3.1 GHz Intel Core i7 and 16 GB RAM.Source: §VI · p. 9

A. Simulated experiment

• The simulation assumes magnetometer measurements are gradients of a scalar potential; a full-GP simulation was considered computationally expensive, so an efficient approach was used.The study used m = 4096 basis functions, described as a good approximation of the true model.Source: A. Simulated experiment · p. 9
• Prediction quality was assessed by RMSE against the true field on a 9,261-position validation mesh spanning the simulated volume.Source: A. Simulated experiment · p. 9
• The comparison included independent GPs with separate hyperparameters and independent GPs with shared hyperparameters.Source: A. Simulated experiment · p. 9
• The first Monte Carlo study varies the number of randomly distributed training points to assess data sparsity and interpolation demand.Source: A. Simulated experiment · p. 9
• A second Monte Carlo study varied hyperparameter-optimization initialization across 30 randomly selected initial sets and training-data lengths.Source: A. Simulated experiment · p. 9
• The scalar-potential GP has lower RMSE than separate-component and shared-hyperparameter GPs across training-data amounts, with more training data reducing RMSE.The authors attribute its advantage to incorporating more physical knowledge.Source: A. Simulated experiment · p. 9, Figure 5 · p. 10
• The scalar-potential model still outperforms the alternatives in the initialization study, while separate-component GPs suffer most from local minima.Source: A. Simulated experiment · p. 10

B. Empirical proof-of-concept data

• The proof-of-concept used magnetometer measurements from a sensor board slid over tables arranged above irregularly placed magnets, with separate training and validation datasets.Each dataset comprised approximately three minutes sampled at 100 Hz; optical measurements supplied positions and orientations for transforming readings into the lab frame.Source: B. Empirical proof-of-concept data · p. 10
• Training data covering nearly the displayed area yielded a low-uncertainty scalar potential, from which the approach predicted the magnetic field and its intensity.The intensity is only indirectly related to the GP output, but was used for intuitive visualization.Source: B. Empirical proof-of-concept data · pp. 10–11
• At around m = 1000 basis functions, Algorithm 1’s validation estimates approached those of the full GP, compared after downsampling data to 5 Hz.Figure 8 plots validation RMSE against the number of basis functions.Source: B. Empirical proof-of-concept data · p. 11, Figure 8 · p. 11

C. Mapping the magnetic field in a building

• For predictions of measurements along the validation path, component-wise RMSEs were (2.35 µT, 3.05 µT, 2.71 µT), and mean absolute errors were (1.72 µT, 2.42 µT, 2.03 µT); the validation walk was 54 m long with 2,340 samples.Source: C. Mapping the magnetic field in a… · p. 12
• The estimated map’s overall shape agreed when trained separately on each training path; validation-path component-wise RMSEs were (2.35, 3.05, 2.71) µT, with mean absolute errors (1.72, 2.42, 2.03) µT.The lobby was approximately 600 m²; training data totaled n = 32,703, and validation used a 54 m path (n = 2,340).Source: C. Mapping the magnetic field in a… · pp. 11–12

D. Online mapping

• Sequential reduced-rank estimation used m = 1024 basis functions; calculations were performed offline, though the algorithm was described as fast enough for real-time online estimation.The first experiment assumed a stationary magnetic environment.Source: D. Online mapping · p. 13
• The static-field estimate updated five times per second across 43,029 observations, with snapshots showing the evolving estimate and uncertainty.Transparency represented marginal variance; the final magnitude and vector components were shown separately.Source: D. Online mapping · p. 13, Figure 11 · p. 13
• The static map showed stronger magnetic activity at the front of the region than farther back, and reconstructions with Trivisio and smartphone sensors left the results and conclusions unchanged.The authors suspected metallic pipelines or structures in the floor as the source of the strong activity.Source: D. Online mapping · p. 13
• For dynamic modeling, the robot collected about 300 seconds of data after two metallic toolbox shelves were introduced, using a temporal scale of ℓtime = 1 hour.The authors described this as encoding an assumption of slow local changes, while noting that the data was informative about abrupt changes.Source: D. Online mapping · p. 14
• Toolbox-induced changes were about 2 µT, roughly 2% of the field scale, and were tracked when the robot passed the affected location.The estimate changed around t = 20 s, 70 s, and 130 s, then stabilized; effects were restricted to the boxes’ immediate vicinity.Source: D. Online mapping · p. 14

VII. DISCUSSION

• More physical information in the GP prior improved interpolation and extrapolation when the data agreed with its assumptions, as illustrated by the compared prior formulations.The discussion considers independent components, shared hyperparameters, and a curl-free scalar-potential model.Source: §VII · p. 14
• The mapping method models uncertainty through a stochastic field model, while the GP covariance function provides a theoretically justified way to encode prior information.Source: §VII · p. 14

VIII. CONCLUSION

• The vector-field components were modeled jointly through a GP prior on a latent scalar potential, enforcing a curl-free field—a justified assumption in free spaces.The prior also encodes a baseline field with smooth small-scale variations.Source: §VIII · pp. 14–15
• The Laplace-operator eigenbasis inference scheme has linear computational complexity in the number of magnetic-field observations and was extended to sequential, time-dependent estimation.The discussion notes that disturbance effects and richer temporal behavior could also be incorporated into the model.Source: §VII · p. 14, §VIII · p. 15
• Four experiments demonstrated feasibility and practicality, including a physics-informed simulation, an approximation proof of concept, smartphone mapping on foot, and online mapping with a wheeled robot.Source: §VIII · p. 15