Modeling and interpolation of the ambient magnetic field by Gaussian processes
Arno Solin, Manon Kok, Niklas Wahlström, Thomas B. Schön, Simo Särkkä
AACCEPTED TO IEEE TRANSACTIONS ON ROBOTICS 1
Modeling and Interpolation of the AmbientMagnetic Field by Gaussian Processes
Arno Solin, Manon Kok,
Member, IEEE , Niklas Wahlstr¨om, Thomas B. Sch¨on,
Senior Member, IEEE, and Simo S¨arkk¨a,
Senior Member, IEEE
Abstract —Anomalies in the ambient magnetic field can beused as features in indoor positioning and navigation. By usingMaxwell’s equations, we derive and present a Bayesian non-parametric probabilistic modeling approach for interpolation andextrapolation of the magnetic field. We model the magnetic fieldcomponents jointly by imposing a Gaussian process (GP) prioron the latent scalar potential of the magnetic field. By rewritingthe GP model in terms of a Hilbert space representation,we circumvent the computational pitfalls associated with GPmodeling and provide a computationally efficient and physicallyjustified modeling tool for the ambient magnetic field. The modelallows for sequential updating of the estimate and time-dependentchanges in the magnetic field. The model is shown to work well inpractice in different applications: we demonstrate mapping of themagnetic field both with an inexpensive Raspberry Pi poweredrobot and on foot using a standard smartphone.
Index Terms —Gaussian process, magnetic field, Maxwell’sequations, mapping, online representation
I. I
NTRODUCTION
Magnetic material causes anomalies in the ambient mag-netic field. In indoor environments, large amounts of suchmagnetic material are present in the structure of buildingsand in furniture. Our focus is on building maps of the indoormagnetic field these structures are inducing. These maps areconstructed by interpolating three-dimensional magnetic fieldmeasurements obtained using magnetometers. An illustrationof a map obtained using our proposed method is available inFigure 1.Magnetic maps of indoor environments can be used inindoor positioning and navigation applications (see, e.g. , [1]).In these applications, sensors providing position informationthat is accurate on a short time-scale—but drifts on longertime-horizons—are typically combined with absolute positionmeasurements. For example, data from wheel encoders andinertial sensors—which give accurate position informationonly on a short time-scale—can be combined with ultra-wideband, Wi-Fi, or optical measurement equipment such ascameras (see, e.g. , [2, 3]). The downside of these sourcesof absolute position is that they typically rely on additional
Arno Solin is with the Department of Computer Science, Aalto University,02150 Espoo, Finland, and with IndoorAtlas Ltd., Helsinki, Finland, e-mail:arno.solin@aalto.fi.Manon Kok is with the Department of Engineering, University of Cam-bridge, CB2 1PZ Cambridge, United Kingdom, e-mail: [email protected] Wahlstr¨om and Thomas B. Sch¨on are with the Department ofInformation Technology, Uppsala University, SE-751 05 Uppsala, Sweden,e-mail: { niklas.wahlstrom, thomas.schon } @it.uu.se.Simo S¨arkk¨a is with the Department of Electrical Engineering and Automa-tion, Aalto University, 02150 Espoo, Finland, e-mail: simo.sarkka@aalto.fi. M a gn it ud e ( µ T ) (a) Interpolated magnetic field strength x -component y -component z -component − −
20 0 20 (b) Vector field components
Fig. 1: Interpolated magnetic field of the lobby of a buildingat the Aalto University campus. The marginal variance (un-certainty) is visualized by the degree of transparency.infrastructure or require certain conditions to be fulfilledsuch as line-of-sight. The advantage of using the magneticfield for positioning is that it can be measured by a smalldevice, without additional infrastructure and without line-of-sight requirements. Furthermore, magnetometers are nowadayspresent in (almost) any inertial measurement unit (IMU) orsmartphone. A requirement for localization using the ambientmagnetic field as a source of position information is thataccurate maps of the magnetic field can be constructed withinreasonable computational complexity which is the focus of thiswork. This can also be regarded as a step towards simultaneouslocalization and mapping (SLAM) using magnetic fields inwhich localization is done while building the map (see, e.g. ,[4, 5]). a r X i v : . [ c s . R O ] M a r CCEPTED TO IEEE TRANSACTIONS ON ROBOTICS 2
We interpolate the magnetic field using a Bayesian non-parametric approach where prior knowledge about the proper-ties of magnetic fields is incorporated in a Gaussian process(GP) prior. GPs (see, e.g. , [6, 7]) are powerful tools forBayesian non-parametric inference and learning, and they pro-vide a framework for fusing first-principles prior knowledgewith quantities of noisy data. This has made them populartools in signal processing, machine learning, robotics andcontrol [8–10].The contributions of this paper are three-fold. First, wemodel the ambient magnetic field using a Gaussian processprior in which we incorporate physical knowledge about themagnetic field. This extends the work by Wahlstr¨om et al. [11]by presenting an approach where the GP prior is a latent(unobservable) magnetic potential function. Second, we use acomputationally efficient GP implementation that allows us touse the large amounts of data provided by the magnetometer.To circumvent the well-known computational challenges withGPs (see, e.g. , [7]), we rewrite the model in terms of a Hilbertspace representation introduced by Solin and S¨arkk¨a [12]. Weextend the approach to allow for modeling of the bias causedby the Earth magnetic field. Third, we use this method incombination with the sequential approach introduced in S¨arkk¨aet al. [13]. This allows for online updating of the magneticfield estimate. It also opens up the possibility to focus onthe spatio-temporal problem in which the magnetic field canchange over time, for instance due to furniture being movedaround. An extensive evaluation of the proposed method isdone using both simulated and empirical data. The simulationstudy and a small-scale experiment illustrate the feasibilityand accuracy of the approach and allow for comparison withother methods. Experiments with a mobile robot and witha hand-held smartphone show the applicability to real-worldscenarios.This paper is structured as follows. The next section coversa survey of existing work, which also provides additionalmotivation for the approach. Section III provides a briefbackground of the properties of magnetic fields relevant to thiswork. The Gaussian process regression model is constructed inSection IV, which is then extended to explicit algorithms forbatch and sequential estimation in the next section. Section VIcovers the experiments. The experimental results and someadditional comments regarding the methodology are discussedat the end of the paper.II. R
ELATED WORK
Spatial properties of the magnetic field have been of interestin a large variety of research domains. For instance, themagnetic field has been extensively studied in geology (see, e.g. , [14]) but also in magnetospheric physics, geophysics,and astrophysics. In all of these domains, interpolation of themagnetic field is of interest (see, e.g. , [15–19] for examplesof magnetic field interpolation in the respective areas).In recent years interest has emerged in using the magneticfield as a source of position information for indoor positioning[20]. Feasibility studies have been conducted, focusing bothon the time-varying nature of the magnetic field and on the amount of spatial variation in the magnetic field. Liet al. [21] report experiments showing that the magnetic fieldin a building typically shows large spatial variations and smalltime variations. This is also supported by the experimentalstudy reported by Angermann et al. [22] in which significantanomalies of the ambient magnetic field are reported. Theseexperiments give confidence that the magnetic field providessufficient information for localization purposes. However, Liet al. [21] also report significant temporal changes in themagnetic field in the vicinity of mobile magnetic structures,in their case an elevator.A number of approaches have been reported on buildinga map of the ambient magnetic field for indoor localizationpurposes. Le Grand and Thrun [23] propose a method to builda map of the magnetic field by collecting magnetometer datain a grid and linearly interpolating between these points. Thismap is subsequently used for localization with a particle filtercombining magnetometer and accelerometer measurementsfrom a smartphone. Robertson et al. [24] present a SLAMapproach for pedestrian localization using a foot-mountedIMU. They use the magnetic field intensity which they modelusing spatial binning. Frassl et al. [25] discuss the possibilityof using more components of the magnetic field (for instancethe full three-dimensional measurement vector) in the SLAMapproach instead. Vallivaara et al. [26, 27] present a SLAMapproach for robot localization. They model the ambient mag-netic field using a squared exponential GP prior for each of themagnetic field components. Wahlstr¨om et al. [11] incorporateadditional physical knowledge by making use of Maxwell’sequations resulting in the use of curl- and divergence-free GPpriors instead.As can be concluded, there exists a wide range of existingliterature when it comes to modeling the ambient magneticfield. The amount of information that is used differs betweenthe approaches. For instance, some approaches use full three-dimensional magnetic field vectors while others only use aone-dimensional magnetic field intensity. Furthermore, theamount of physical information that is included differs. In thispaper, we build on the approach by Wahlstr¨om et al. [11] anduse the full three-dimensional magnetometer measurements.We include physical knowledge in terms of the magnetic fieldpotential.As discussed above, GPs have frequently been used in mod-eling and interpolation of magnetic fields. GP regression hasalso successfully been applied to a wide range of applications(see, e.g. , [28–33]). Furthermore, it has previously been usedfor SLAM (see, e.g. , [34–38]). One of the challenges in usingGPs is the computational complexity (see, e.g. , [7]), whichscales cubically with the number of training data points. Con-sidering the high sampling rate of the magnetometer and thefact that each observation contains three values, a large numberof measurements is typically available for mapping. Becauseof these computational challenges, the data in Wahlstr¨omet al. [11] was downsampled.Attempts to speed up GP inference have spawned a widerange of methods which aimed at bringing GP regression todata-intensive application fields. These methods (see [39] fora review) typically build upon reducing the rank of the Gram
CCEPTED TO IEEE TRANSACTIONS ON ROBOTICS 3 + =
Earth magnetic field Ferromagnetic object Distorted field
Fig. 2: A ferromagnetic object deflects the Earth’s magneticfield and introduces distortions in the field.(covariance) matrix and using the matrix inversion lemma tospeed up matrix inversion. For stationary covariance functions,the spectral Monte Carlo approximation by L´azaro-Gredillaet al. [40], the Fourier features by Hensman et al. [41], orthe Laplace operator eigenbasis based method introduced bySolin and S¨arkk¨a [12] can be employed. For uniformly spacedobservations, fast Fourier transforms can provide computa-tional benefits [42, 43]. As will be shown later on in thispaper, the Laplace operator approach by Solin and S¨arkk¨a [12]falls natural to modeling of the magnetic field in terms of amagnetic field potential.All approaches on mapping of magnetic fields discussedabove assume that the magnetic field is constant over time.However, as shown by Li et al. [21] significant temporalchanges in the magnetic field occur in the vicinity of mobilemagnetic structures. For GP models evolving in time, spatio-temporal GP models (see, e.g. , [9]) can be solved efficientlyusing Kalman filtering methods [13, 44–47]. In this paper,we will take the approach of S¨arkk¨a et al. [13] to compose aspatio-temporal GP prior for the model and solve the inferenceproblem by a sequential Kalman filtering setup. This allowsfor online estimation of the magnetic field estimate and canbe used to allow for time variations in the magnetic field.Combining models from physics with GPs has also beenstudied under the name
Latent force models by ´Alvarez et al.[48, 49]. The connection of these models with spatio-temporalKalman filtering was studied, for example, in the work byHartikainen and S¨arkk¨a (see also [13, 50, 51]). However,the Kalman filtering approach itself dates back to Curtainand Pritchard [52] in the same way that GPs date back toO’Hagan [6].III. T
HE AMBIENT MAGNETIC FIELD
On a macroscopic scale, magnetic fields are vector fields ,meaning that at any given location, they have a direction and astrength (magnitude). These properties are familiar from every-day life: the force created by permanent magnets attracting andrepelling ferromagnetic materials is used in various utensils,and the compass aligning itself with the direction of the Earth’smagnetic field has proved invaluable for mankind during thepast centuries. The Earth’s magnetic field sets a backgroundfor the ambient magnetic field, but deviations caused by thebedrock and anomalies induced by man-built structures deflectthe Earth’s magnetic field. This makes the magnetic field varyfrom point to point, see Figure 2.We describe the magnetic field with a function H ( x ) , where H : R → R . For each point in space x , there will be anassociated magnetic field H ( x ) . Such a vector field can be Ω Fig. 3: Illustration of a vector field with non-zero curl. Thevortex point makes it non-curl-free as the vector field curlsaround it. However, the subset Ω excludes the vortex pointand the vector field is curl-free in this region. To this regiona scalar potential ϕ can be associated, here illustrated withshading.visualized with field lines, where points in space are associatedwith arrows. The principles under which the magnetic fieldis affected by structures of buildings are well known andgoverned by the very basic laws of physics (see, e.g. , [53–55]).In this work, we make use of the fact that the magnetic field H is curl-free ∇ × H = (1)provided that there is no free current (current in wires forexample) in the region of interest (see [56] for more details).This assumption is valid in most indoor environments wherethe major source for variations in the ambient field is causedby metallic structures rather than free currents in wires.One property of curl-free vector fields is that the lineintegral along a path P only depends on its starting point A and end point B , and not on the route taken (cid:90) P H ( x ) · d x = ϕ ( A ) − ϕ ( B ) , (2)where ϕ : R → R . This can be rewritten by interpreting ϕ as a scalar potential H = −∇ ϕ. (3)Figure 3 illustrates the curl-free property and the scalar po-tential. Domain Ω is curl-free and has an associated scalarpotential, while the entire domain is not curl-free due to thevortex point. Note that in a non-curl-free vector field no suchscalar potential exists since a line integral around the swirl isnon-zero. For the magnetic field H , the swirl corresponds to awire of free current pointing perpendicular to the plane, whichwe assume is not included in the region of interest.The relation (3) is the key equation that we will exploit inour probabilistic model of the ambient magnetic field. We willchoose to model the scalar potential ϕ instead of the magneticfield H directly. This implicitly imposes the constraints on themagnetic field that the physics is providing. This model willbe explained in the next section. CCEPTED TO IEEE TRANSACTIONS ON ROBOTICS 4
IV. M
ODELING THE MAGNETIC FIELD USING G AUSSIANPROCESS PRIORS
In this section, we introduce our approach to modeling andinterpolation of the ambient magnetic field. We use a Bayesiannon-parametric model in which we use knowledge about thephysical properties of the magnetic field as prior information.We tackle the problem of interpolating the magnetic field usingGaussian process (GP) regression. In Section IV-A we firstgive a brief background on GPs. After this, we introduce theproblem of modeling the magnetic field in Section IV-B. Acommonly used GP model of the magnetic field is introducedin Section IV-C. In Section IV-D, we subsequently introduceour proposed method for modeling the magnetic field in whichwe encode the physical properties that were presented in theprevious section.
A. Gaussian process regression
In GP regression [7] the model functions f ( x ) are assumedto be realizations from a Gaussian random process priorwith a given covariance function κ ( x , x (cid:48) ) . Learning amountsto computing the posterior process at some test inputs x ∗ given a set of noisy measurements y , y , . . . , y n observed at x , x , . . . , x n , respectively. This model is often written in theform f ( x ) ∼ GP (0 , κ ( x , x (cid:48) )) ,y i = f ( x i ) + ε i , (4)where the observations y i are corrupted by Gaussian noise ε i ∼ N(0 , σ noise ) , for i = 1 , , . . . , n . Because both theprior and the measurement model are Gaussian, the posteriorprocess will also be Gaussian. Hence, the learning problemamounts to computing the conditional means and covariancesof the process evaluated at the test inputs.Prediction of yet unseen process outputs at an input location x ∗ amounts to the following in GP regression: p ( f ( x ∗ ) | D ) =N( f ( x ∗ ) | E [ f ( x ∗ )] , V [ f ( x ∗ )]) . The conditional mean andvariance can be computed in closed-form as (see [7]) E [ f ( x ∗ )] = k T ∗ ( K + σ noise I n ) − y , V [ f ( x ∗ )] = κ ( x ∗ , x ∗ ) − k T ∗ ( K + σ noise I n ) − k ∗ , (5)where K i,j = κ ( x i , x j ) , k ∗ is an n -dimensional vector withthe i th entry being κ ( x ∗ , x i ) , and y is a vector of the n observations. Furthermore, due to Gaussianity, the marginallikelihood (evidence) of the covariance function and noiseparameters can also easily be computed, allowing for Bayesianinference of the parameters as well [7].The choice of a specific covariance function encodes the a priori knowledge about the underlying process. One of themost commonly used covariance functions, which will alsofrequently be used in the next sections, is the stationary andisotropic squared exponential (also known as exponentiatedquadratic, radial basis function, or Gaussian). Following thestandard notation from Rasmussen and Williams [7] it isparametrized as κ SE ( x , x (cid:48) ) = σ SE exp (cid:18) − (cid:107) x − x (cid:48) (cid:107) (cid:96) SE (cid:19) , (6) where the hyperparameters σ SE and (cid:96) SE represent the mag-nitude scale and the characteristic length-scale, respectively.These can be learned from data, for instance by maximizingthe marginal likelihood. B. Interpolation of magnetic fields
In this work we tackle the problem of interpolating themagnetic field to spatial locations from where we do not haveany measurements. In other words, we will tackle the problemof predicting the latent (unobservable noise-free) magneticfield f ( x ∗ ) (such that f : R → R ) at some arbitrary location x ∗ given a set of measurements D = { ( x i , y i ) } ni =1 of themagnetic field. Here, the measurements y correspond to theH-field corrupted by i.i.d. Gaussian noise.Two important things need to be noted with regard to theinterpolation of the magnetic field. First, note that the mea-surements of the magnetic field are vector-valued (contrary tothe scalar observations in (4)). This raises the question of howto deal with the different magnetic field components. They caneither be treated separately as will be done in Section IV-C, ora relation between the different components can be assumed,as is the case in the method we propose in Section IV-D.Secondly, note that the function describing the magnetic fieldis not zero-mean, contrary to the GP model in (4). Instead, itsmean lies around a local Earth magnetic field. This dependson the location on the Earth but can also deviate from theEarth’s magnetic field in indoor environments due to magneticmaterial in the structure of the building. The unknown meancan be modeled as an additional part of the covariance function κ ( x , x (cid:48) ) [7].An illustration of GP regression for magnetic fields isprovided in Figure 4, where noisy readings of a magnetic fieldhave been collected along route A–D (comprising D ), and themagnetic field along route D–E (comprising the predictionlocations x ∗ ) needs to be inferred from the measurements.Each component varies around a local magnetic field dueto magnetic material in the vicinity of the sensor. Differentinterpolation techniques can be used based on different priorknowledge that can be incorporated in the GP. Two differentinterpolation results are shown: one based on independentmodeling of each vector field components (with shared hy-perparameters) and another based on associating the GP priorwith the scalar potential of the (curl-free) vector field. Theseare based on the models we will introduce in the coming twosections. C. Separate modeling of the magnetic field components
The most straightforward approach to GP modeling ofvector-valued quantities is to model each of the field com-ponents as an independent GP. This approach has been widelyapplied in existing literature (see, e.g. , [26, 27, 57–59]). Foreach of the three magnetic field components d ∈ { , , } , thismodel can be written as f d ( x ) ∼ GP (0 , κ const. ( x , x (cid:48) ) + κ SE ( x , x (cid:48) )) ,y d,i = f d ( x i ) + ε i,d , (7)where the observations y d,i , i = 1 , , . . . , n , are corrupted byindependent Gaussian noise with variance σ noise . The non-zero CCEPTED TO IEEE TRANSACTIONS ON ROBOTICS 5
AB CDE (a) The route
D E
IndependentGPs
D E
ScalarpotentialGP
A B C H x H y H z (b) The data and the GP prediction for D–E Fig. 4: A simulated example of the interpolation problem. (a) Training data has been collected along the route A–D, but themagnetic field between D–E is unknown. (b) The noisy observations of the magnetic field between A–D, and GP predictionswith 95% credibility intervals. Both the independent GP modeling approach (with shared hyperparameters) and the scalarpotential based curl-free GP approach are visualized. The simulated ground truth is shown by the solid lines.mean of the magnetic field is handled by a constant covariancefunction (see [7]) κ const. ( x , x (cid:48) ) = σ const. , (8)where σ const. is a magnitude scale hyperparameter. The small-scale variation in the field is modeled by a squared exponentialcovariance function (6). Hence, the model (7) encodes theknowledge that the realizations are expected to be smoothfunctions in space with a constant shift from zero mean.The model has four hyperparameters: two magnitude scaleparameters ( σ const. and σ SE ), a length-scale parameter ( (cid:96) SE ),and a noise scale parameter ( σ noise ). Assuming that the com-ponents are completely separate, each component has fourhyperparameters to learn. The resulting model is flexible,as it does not encode any relation between the vector fieldcomponents. In practice, this might lead to problems in hy-perparameter estimation, with parameter estimates convergingto local optima and magnetic field components behavingvery differently with respect to each other. Therefore, thehyperparameters are often fixed to reasonable values—insteadof learned from data (see, e.g. , [57]).A more sensible approach for separate, but not completelyindependent, modeling of the magnetic field measurements,models them as realizations of three independent GP priorswith joint learning of the shared hyperparameters (see also[57]). Note that for this model, the covariance in the GPposterior is independent of the outputs y and only depends onthe input locations x (which are shared for all components in y ). Hence, calculating the marginal likelihood only requiresinverting a matrix of size n (not n ). For this model, theexpression for evaluating the log marginal likelihood functionfor hyperparameter optimization can be written as L ( θ ) = − log p ( y | θ , D ) = 32 log | K θ + σ noise I n | + 12 tr (cid:2) y ( K θ + σ noise I n ) − y T (cid:3) + 3 n π ) , (9) where y ∈ R × n and K θ ∈ R n × n .Figure 4b shows the results of predicting the magnetic fieldbehavior along the route D–E for the GP prior modelingthe three magnetic field components separately but with jointlearning of the shared hyperparameters. The colored patchesshow the 95% credibility intervals for the prediction with themean estimate visualized by the white line. The simulatedground-truth (solid colored line) falls within the shown inter-val, and the model captures the general shape of the magneticfield variation along the path. The strengths of this model arethat it is flexible and that the assumptions are conservative.The weaknesses on the other hand are evident: The modeldoes not incorporate physical knowledge of the magnetic fieldcharacteristics. In the next section we will instead explorethis knowledge by modeling the magnetic field as derivativemeasurements of a scalar potential. D. Modeling the magnetic field as the gradient of a scalarpotential
Following our choices in Section III, we assume that themagnetic field H can be written as the gradient of a scalarpotential ϕ ( x ) according to (3). Here, ϕ : R → R and x ∈ R is the spatial coordinate. We assume ϕ ( x ) to be a realizationof a GP prior and the magnetic field measurements y i ∈ R to be its gradients corrupted by Gaussian noise. This leads tothe following model ϕ ( x ) ∼ GP (0 , κ lin. ( x , x (cid:48) ) + κ SE ( x , x (cid:48) )) , y i = −∇ ϕ ( x ) (cid:12)(cid:12) x = x i + ε i , (10)where ε i ∼ N( , σ noise I ) , for each observation i =1 , , . . . , n . The squared exponential covariance function (6)in (10) allows us to model the magnetic field anomaliesinduced by small-scale variations and building structures. Thelocal Earth’s magnetic field contributes linearly to the scalarpotential as κ lin. ( x , x (cid:48) ) = σ lin. x T x , (11) CCEPTED TO IEEE TRANSACTIONS ON ROBOTICS 6 where σ lin. is the magnitude scale hyperparameter. To simplifythe notation in the next sections, we introduce the notation f ( x ) for the gradient field evaluated at x .Derivative measurements can straightforwardly be used inGaussian process regression because nabla ( ∇ ) is a linearoperator and Gaussianity is preserved under linear opera-tions [7]. Hence, the model (10) allows us to learn a map ofthe magnetic field and make predictions at unseen locationsusing the methods described in Section IV-A. We learn the fourhyperparameters of the model—the magnitude scale parame-ters ( σ lin. and σ SE ), the length-scale parameter ( (cid:96) SE ), and thenoise scale parameter ( σ noise )—by maximizing the marginallikelihood.The model (10)—where we place a GP prior on thescalar potential ϕ ( x ) and the magnetometer measurements arederivative measurements—can in fact equivalently be writtenas a GP prior on the magnetic field H ( x ) as H ( x ) ∼ GP ( , σ const. I + K curl ( x , x (cid:48) )) , (12)and magnometer measurements being of the standard form (4).Here, K curl ( x , x (cid:48) ) is the so-called curl-free kernel (see, e.g. ,[60–62]). This kernel ensures that any sample drawn from thisGP prior obeys the curl-free constraints (1). The equivalenceof (10) and (12) is shown in Wahlstr¨om [56]. In Jidlinget al. [63] a more general approach to force GP priors toobey linear constraints is presented. For our work, however,the model formulation through the scalar potential is crucial, asit will enable us to easily extend the model to an approximateform for efficient GP inference in the next section.The second set of predictions for the route D–E in Figure 4bshows the interpolation outcome from the model (10). Incomparison to the independent GP model, the scalar potentialbased GP prior provides additional information to the modelby tying the vector field components to each other. Thisimproves the estimates in terms of accuracy and makes the95% credibility interval more narrow.V. E FFICIENT GP MODELING OF THE MAGNETIC FIELD
Gaussian processes are convenient tools for assigning flexi-ble priors to data—as we saw in the previous section. However,the main problem with the model in the previous section isits high computational cost. The approach scales as O ( n ) (recall that each observation is three-dimensional, meaning n becomes large very quickly). This computational complexityrenders the approach more or less useless in practice, whenthe number of observations becomes large (say n > ).This is a fundamental restriction associated with the naiveformulation of GP models involving the inversion of thecovariance matrix. Using special structure of the problemand/or approximative methods, this high computational costcan often be circumvented. In this section we present anapproach which both uses the special differential operatorstructure and projects the model on a set of basis functionscharacteristic to the covariance function. We first present themethod for spatial batch estimation, and then extend it to atemporal dimension as well. Existing GP methods for mapping and interpolation of themagnetic field have been considering only batch estimation,where the data is first acquired and then processed as a batch.In this section, we aim to extend this to an online method,enabling the GP regression estimate of the magnetic field tobe updated when new data is acquired. We denote such a dataset as D n = { ( x i , y i ) | i = 1 , , . . . , n } , and thus D i denotesall the data that has been observed up to time instance t i .Considering time as part of the data stream enables us tothink of three distinctive setups for estimation of the magneticfield: • Batch estimation of the magnetic field, where the datais first acquired and then the field is estimated at once. • Sequential updating of the field estimate, where weassume all the measurements to be of the same staticmagnetic field. • Spatio-temporal estimation of the time-dependent mag-netic field, where we assume the field to change overtime.In the next sections we will present how these scenarios can becombined with the scalar potential based GP scheme withoutrequiring to repeat the batch computations after each sample.
A. Reduced-rank GP modeling
A recent paper by Solin and S¨arkk¨a [12] presents anapproach that is based on a series expansion of stationarycovariance functions. The approximation is based on thefollowing truncated series: κ ( x , x (cid:48) ) ≈ m (cid:88) j =1 S ( λ j ) φ j ( x ) φ j ( x (cid:48) ) , (13)where S ( · ) is the spectral density of the covariance func-tion κ ( · , · ) , φ j ( x ) is the j th eigenfunction of the negativeLaplace operator and λ j is the corresponding eigenvalue.The efficiency of this approach is based on two properties:(i) the eigenfunctions are independent of the hyperparametersof the covariance function, and (ii) for many domains theeigenfunctions and eigenvalues can be solved beforehand inclosed-form. Truncating this expansion at degree m (cid:28) n allows the GP regression problem to be solved with a O ( nm ) and the hyperparameters to be learned with a O ( m ) timecomplexity. The memory requirements scale as O ( nm ) .Our interest lies in modeling the magnetic field in compactsubsets of R , allowing us to restrict our interest to domains Ω comprising three-dimensional cuboids (rectangular boxes)such that x ∈ [ − L , L ] × [ − L , L ] × [ − L , L ] ⊂ R (recallthat a stationary covariance function is translation invariant).In this domain, we can solve the eigendecomposition of theLaplace operator subject to Dirichlet boundary conditions (cid:40) −∇ φ j ( x ) = λ j φ j ( x ) , x ∈ Ω ,φ j ( x ) = 0 , x ∈ ∂ Ω . (14)The choice of the domain and boundary conditions is arbi-trary, but for regression problems with a stationary covariancefunction the model reverts back to the prior outside theregion of observed data, so the Dirichlet boundary condition CCEPTED TO IEEE TRANSACTIONS ON ROBOTICS 7 does not restrict the modeling if Ω is chosen suitably. Thisparticular choice of domain and boundary conditions yieldsthe following analytic expression for the basis functions andthe corresponding eigenvalues: φ j ( x ) = (cid:89) d =1 √ L d sin (cid:18) πn j,d ( x d + L d )2 L d (cid:19) , (15) λ j = (cid:88) d =1 (cid:18) πn j,d L d (cid:19) , (16)where the matrix n ∈ R m × consists of an in-dex set of permutations of integers { , , . . . , m } ( i.e. , { (1 , , , (1 , , , . . . , (1 , , , . . . , (2 , , , . . . } ) . The basisfunctions are independent of the hyperparameters, and thusonly need to be evaluated once.The Laplace operator eigenbasis approximation methodcan be combined with the independent GP approach, theindependent GPs with shared hyperparameters, and the scalarpotential GP approach. Our presentation will be specific tothe scalar potential model presented in Section IV-D, but asimilar setup can be constructed for the other methods as well.The approximation method is even better suited for the scalarpotential model because the approximation is based on theeigendecomposition of the Laplace operator. This eigenbasisfalls naturally to the problem formulation in which the latentpotential field is observed through gradients.The covariance in our model (10) consists of a squaredexponential and a linear covariance function. The former isstationary and the approach from [12] can straightforwardly beapplied. The latter is not stationary, but there is also no need toapproximate it. Consequently, we consider the approximation κ ( x , x (cid:48) ) = κ lin. ( x , x (cid:48) ) + κ SE ( x , x (cid:48) ) ≈ σ lin. x T x + m (cid:88) j =1 S SE ( λ j ) φ j ( x ) φ j ( x (cid:48) ) . (17)Note that adding the linear covariance function, the approxi-mation will no longer revert to zero on the domain boundary,but will instead revert to the scalar potential of the local Earthmagnetic field.The computational benefits come from the approximateeigendecomposition of the Gram (covariance) matrix, K i,j = κ ( x i , x j ) (see [12] for derivations and discussion). It can nowbe written out in terms of the basis functions and spectraldensities: K ≈ ΦΛΦ T . The basis functions, which span thesolution, are collected in the matrix Φ ∈ R n × (3+ m ) , with thefollowing rows Φ i = (cid:0) x T i , φ ( x i ) , φ ( x i ) , . . . , φ m ( x i ) (cid:1) , (18)for i = 1 , , . . . , n . Accordingly, we define the correspondingmeasurement model matrix projecting the derivative obser-vations onto the basis functions. Analogously, we define thematrix ∇ Φ ∈ R n × (3+ m ) as the following block-row matrix: ∇ Φ i = (cid:0) ∇ x T i , ∇ φ ( x i ) , ∇ φ ( x i ) , . . . , ∇ φ m ( x i ) (cid:1) , (19)for i = 1 , , . . . , n . Similarly we define Φ ∗ and ∇ Φ ∗ asvectors evaluated at the prediction input location x ∗ defined analogously to Equations (18) and (19), respectively. Thematrix Λ is defined by Λ = diag( σ lin. , σ lin. , σ lin. , S SE ( λ ) , S SE ( λ ) , . . . , S SE ( λ m )) . (20)For three-dimensional inputs, the spectral density function ofthe squared exponential covariance function (6) is given by S SE ( ω ) = σ SE (2 π(cid:96) SE ) / exp (cid:18) − ω (cid:96) SE (cid:19) , (21)where the hyperparameters σ SE and (cid:96) SE characterize the spec-trum. B. Batch estimation
We first tackle the batch estimation problem which providesthe approximative solution to the GP regression problem inEquation (5) for the scalar potential GP.Following the derivations of Solin and S¨arkk¨a [12], predic-tions for interpolation and extrapolation of the magnetic fieldat yet unseen input locations x ∗ are given by: E [ f ( x ∗ )] ≈ ∇ Φ ∗ ([ ∇ Φ ] T ∇ Φ + σ noise Λ − ) − [ ∇ Φ ] T vec( y ) , V [ f ( x ∗ )] ≈ σ noise ∇ Φ ∗ ([ ∇ Φ ] T ∇ Φ + σ noise Λ − ) − [ ∇ Φ ∗ ] T , (22)where vec( · ) is the vectorization operator which converts amatrix to a column vector by stacking its columns on top ofeach other, such that the × n matrix is converted into avector of size n . The basis functions ∇ Φ and ∇ Φ ∗ need tobe evaluated by Equation (19), and Λ by (20).For this model, the expression for evaluating the logmarginal likelihood function for hyperparameter optimizationcan be written as L ( θ ) = 12 log | K θ + σ noise I n | +12 vec( y ) T ( K θ + σ noise I n ) − vec( y ) + 3 n π ) , (23)where the quantities can be approximated by log | K θ + σ noise I n | ≈ (3 n − m ) log σ noise + m (cid:88) j =1 [ Λ θ ] j,j + log | σ noise Λ − θ + [ ∇ Φ ] T ∇ Φ | , (24) vec( y ) T ( K θ + σ noise I n ) − vec( y ) ≈ σ noise (cid:2) vec( y ) T vec( y ) − vec( y ) T ∇ Φ ( σ noise Λ − θ + [ ∇ Φ ] T ∇ Φ ) − [ ∇ Φ ] T vec( y ) (cid:3) , (25)where the only remaining dependency on the covariancefunction hyperparameters is in the diagonal matrix Λ definedthrough the spectral density in Equation (20). In a softwareimplementation, Cholesky decompositions can be employedfor numerical stability in the calculation of determinants andmatrix inverses. For optimizing the hyperparameters, gradientbased optimizers can be employed (see [12] for details onderiving the partial derivatives).Algorithm 1 describes the step-by-step workflow for apply-ing these equations in practice. The inputs for the method are CCEPTED TO IEEE TRANSACTIONS ON ROBOTICS 8
Algorithm 1:
Algorithm for batch estimation of the scalarpotential GP magnetic field with the reduced-rank ap-proach.
Input: D = { ( x i , y i ) } ni =1 , x ∗ , Ω , m . Output: E [ f ( x ∗ )] , V [ f ( x ∗ )] . Use Eq. (19) to evaluate the basis functions ∇ Φ from x i s and Ω . Use Eq. (23) to optimize hyperparameters θ = { σ lin. , σ SE , (cid:96) SE , σ noise } . Use Eq. (19) to evaluate the basis functions ∇ Φ ∗ from x ∗ s and Ω . Solve the GP regression problem by Eq. (22).the data D (spatial points and the magnetic field readings),the test points x ∗ to predict at, the domain boundaries Ω ,and the approximation degree parameter m (controlling theaccuracy of the Hilbert space approximation, see Eq. (13)).The algorithm returns the marginal mean and variance of thepredicted magnetic field at x ∗ . The scalar potential could bereturned instead by using Equation (18) in step three.The computational complexity of this method scales as O ( nm ) for prediction and O ( m ) for evaluating the marginallikelihood during optimization, which makes it computation-ally appealing in comparison to the computational complexityof the naive full GP solution which is O ( n ) . The memoryrequirements scale as O ( nm ) . C. Sequential estimation
Many applications require online (sequential) estimates ofthe magnetic field. The following formulation provides the same (within numerical precision) solution as the batch esti-mation solution (22) in the previous section. The inferencescheme in the previous section is in practice the solutionof a linear Gaussian estimation problem. Sequential solutionsfor this type of problems have been extensively studied, andthis mathematical formulation is widely known as the
Kalmanfilter . The connection between Kalman filtering and Gaussianprocess regression has recently been studied, for example, in[13, 46, 47].Reformulation of a batch problem to a sequential algorithmis discussed (with examples) in the book by S¨arkk¨a [64].Following this formulation (and notation to large extent) wemay write the following recursion: Initialize µ = and Σ = Λ θ (from the GP prior). For each new observation i = 1 , , . . . , n update the estimate according to S i = ∇ Φ i Σ i − [ ∇ Φ i ] T + σ noise I , K i = Σ i − [ ∇ Φ i ] T S − i , µ i = µ i − + K i ( y i − ∇ Φ i µ i − ) , Σ i = Σ i − − K i S i K T i . (26)This means that for a test input location x ∗ we get predictionsfor the mean and the variance of the magnetic field which aregiven by E [ f ( x ∗ ) | D i ] ≈ ∇ Φ ∗ µ i , V [ f ( x ∗ ) | D i ] ≈ ∇ Φ ∗ Σ i [ ∇ Φ ∗ ] T , (27) Algorithm 2:
Algorithm for sequential modeling of thescalar potential GP magnetic field estimate. Alternative (a)corresponds to the sequential model, and (b) to the spatio-temporal modeling approach.
Input: D n , x ∗ , Ω , m , θ . Output: E [ f ( x ∗ )] , V [ f ( x ∗ )] . Initialize µ = and Σ = Λ θ from Eq. (20). for i = 1 , , . . . , n do Evaluate ∇ Φ i by Eq. (19) from x i . (a) Perform an update by Eq. (26). (b)
Perform an update by Eqs. (31–32). Evaluate the current prediction at x ∗ by Eq. (27). end for and conditional on the data observed up to observation i .Writing out the conditioning on D was stripped in the earliersections for brevity.Here we do not consider optimization of the hyperparam-eters θ . The marginal likelihood can be evaluated throughthe recursion, but in an online setting we suggest optimizingthe hyperparameters with some initial batch early in the datacollection and then re-optimizing them later on if necessary.Algorithm 2 presents the scheme of how to apply theequations in practice. The inputs are virtually the same asfor the batch algorithm, but now the hyperparameters θ areconsidered known. The current estimate can be returned aftereach iteration loop. The computational complexity of thesequential modeling approach scales as O ( m ) per updateand thus has a total computational complexity of O ( nm ) ,but a memory scaling (if intermediate results are not stored)of O ( m ) . D. Spatio-temporal modeling
The sequential model allows for extending the modelingto also track dynamic changes in the magnetic field withoutvirtually any additional computational burden. Let the data D n = { ( t i , x i , y i ) } ni =1 now also comprise a temporal vari-able t which indicates the time when each observation wasacquired.We present the following spatio-temporal model for trackingchanges in the ambient magnetic field. The spatio-temporalGP prior assigned to the scalar potential ϕ ( x , t ) , dependingon both location and time, is defined as follows: ϕ ( x , t ) ∼ GP (0 , κ lin. ( x , x (cid:48) ) + κ SE ( x , x (cid:48) ) κ exp ( t, t (cid:48) )) , (28)where the additional covariance function κ exp ( t, t (cid:48) ) defines theprior assumptions of the temporal behavior. This covariancefunction is defined through κ exp ( t, t (cid:48) ) = exp (cid:18) − | t − t (cid:48) | (cid:96) time (cid:19) , (29)where (cid:96) time is a hyperparameter controlling the length-scaleof the temporal effects. In the temporal domain, this modelis also known as the Ornstein–Uhlenbeck process (see, e.g. ,[7]). The assumption encoded into it is that the phenomenon is CCEPTED TO IEEE TRANSACTIONS ON ROBOTICS 9 continuous but not necessarily differentiable. Therefore it pro-vides a very flexible means of modeling the changing ambientmagnetic field. Also note that in Equation (28) the temporaleffects are only associated with the anomaly component, thebias being tracked as a static component.Following the derivations in Hartikainen and S¨arkk¨a [44],we can write down the dynamic state space model associatedwith the time evolution of the spatio-temporal GP prior model(28) A i = blkdiag( I , I m exp( − ∆ t i /(cid:96) time )) , Q i = blkdiag( , I m [1 − exp( − t i /(cid:96) time )]) , (30)where ∆ t i = t i +1 − t i is the time difference between twoconsecutive samples and denotes a × zero matrix.For the time update (Kalman prediction step) we may thuswrite ˜ µ i = A i − µ i − , ˜ Σ i = A i − Σ i − A T i − + Q i − , (31)and the modified measurement update (Kalman update step) S i = ∇ Φ i ˜ Σ i [ ∇ Φ i ] T + σ noise I , K i = ˜ Σ i [ ∇ Φ i ] T S − i , µ i = ˜ µ i + K i ( y i − ∇ Φ i ˜ µ i ) , Σ i = ˜ Σ i − K i S i K T i . (32)Algorithm 2 features the workflow of applying the method(option (b)). In practice the only additional input is includingthe temporal length-scale in θ . The computational costs arethe same as for the sequential model.VI. E XPERIMENTS
The experiments in this paper are split into four parts. Wefirst demonstrate the feasibility of the scalar potential approachwith simulated data, where comparison to known ground truthis possible. After this we present a small-scale proof-of-concept demonstration of the approach. The third experimentis concerned with mapping the magnetic field in a buildingusing a handheld smartphone. The final experiment uses aninexpensive mobile robot for online mapping and real-timetracking of the changing magnetic field. The methods wereimplemented in Matworks M
ATLAB R (cid:13) on an Apple MacBookPro (3.1 GHz Intel Core i7, 16 GB RAM). A. Simulated experiment
As a first part of our experimental validation of the model,we present a simulation study. This will be used to illustrateour method and to quantify its performance through MonteCarlo simulations. In the simulations, we assume that themagnetic field measurements can indeed be modeled usinga scalar potential as argued in Section III. Hence, we simulatethe magnetometer data from the GP that models the magne-tometer measurements as gradients of a scalar potential fieldas discussed in Section IV-D. Because we are interested insimulating a large number of data points, simulating this datausing a full GP approach is computationally very expensive.We therefore simulate the data using the computationally efficient approach described in Section V-B. We use a largenumber of basis functions ( m = 4096 ). This has been shownto be a good approximation of the true model in [12]. We willalso show that this is a good approximation in Section VI-Bbased on experimental data.The magnetometer is assumed to move in a three-dimensional volume with x, y, z ∈ [ − . , . such that x = ( x, y, z ) . The training data used for training the GP israndomly uniformly distributed over this volume. A validationdata set is used to assess the predictive power of the trainedGP. This validation data is a three-dimensional meshgrid overthe same volume and consists of n val. = 9261 positions x val. with a true magnetic field f true ( x val. ) . The magnetic field ispredicted at these points using the trained GP, leading to f train. ( x val. ) , after which the quality of the GP solution can beassessed in terms of the root mean square error (RMSE). Weset the domain Ω in the GP model to L = L = L = 0 . and simulate using the following hyperparameters σ const. =0 . , σ SE = 1 , (cid:96) SE = 0 . , and σ noise = 0 . (see Sec. V-Band IV-D for more details on the notation).We compare our proposed method with two other ap-proaches. The first is the approach by Vallivaara et al. [26, 27],where the magnetic field is modeled using independent GPsfor all three components. Each GP consists of a constant and asquared exponential kernel and has its own hyperparameters.The second approach, considered in Kemppainen et al. [57],models the magnetic field similarly but with shared hyperpa-rameters. For details, see also Section IV-C.In a first set of Monte Carlo simulations, we analyze theperformance of the three different approaches depending onthe number of (randomly distributed) training data points, thatis in terms of the sparseness of the magnetometer data andthe amount of interpolation that is needed for prediction. Forall three approaches we use a large number of basis functions( m = 4096 ). They are hence expected to approach the perfor-mance of the full GP solution. To exclude problems with localminima in the hyperparameter optimization—which will bethe topic of a second set of simulations—the hyperparameteroptimization is started in the values used for simulating thedata. The results from 30 Monte Carlo simulations are shownin Figure 5a. Naturally, the more data is used for training theGP, the smaller the RMSE becomes. The GP that models themagnetic field measurements as gradients of a scalar potentialfield, outperforms the other two approaches, independent ofthe amount of training data used. This can be understoodfrom the fact that this approach incorporates most physicalknowledge.In a second set of Monte Carlo simulations, the sensitivityto the initialization of the hyperparameter optimization isanalyzed for the three different methods. Only one simulateddata set is used but the hyperparameter optimization is startedin 30 randomly selected sets of hyperparameters θ a for avarying length of the training data. The hyperparametersare assumed to lie around the estimates that are obtainedusing the same optimization strategy as above for 8000 datapoints. Hence, these sets of hyperparameters θ atrue are knownto results in small RMSE values as depicted in Figure 5a.The superscript ‘a’ on θ atrue is used to explicitly denote that CCEPTED TO IEEE TRANSACTIONS ON ROBOTICS 10 − Number of data points, n R M S E Scalar potential GPSeparate GP modelsShared GP model (a) Randomized over data sets − Number of data points, n R M S E Scalar potential GPSeparate GP modelsShared GP model (b) Randomized over initializations
Fig. 5: The average RMSE (with standard deviations) from 30 Monte Carlo simulations as a function of the number of simulateddata points used for training the GP. The constantly lower error for the Hilbert space scalar potential GP in comparison to theseparate GP models and shared hyperparameter GP model is explained by the additional prior physical knowledge encodedinto the model. In (b), the shared hyperparameter GP shows a slight advantage over the fully independent models.
Optical marker Xsens IMU Magnets
Fig. 6: Setup of the experiment: On left the sensor board thatwas used for collecting the data. The right side figure showsthe magnetic environment used in the experiment shown frombelow comprising small magnets to ensure sufficient excitationof the magnetic field.these sets of hyperparameters actually differ between the threedifferent approaches. For the approach where the magneticfield components are modeled using independent GPs, eachcomponent results in a set of hyperparameters θ atrue . Forsimplicity, in this approach, θ atrue is chosen to be the meanof these three sets of hyperparameters.For each of the three approaches, the initial parameters θ are then assumed to deviate from θ atrue by at most as θ = θ atrue [1 + 0 . − , . (33)The Monte Carlo simulation results are depicted in Figure 5b.As can be seen, the approach which models the three magneticfield components using separate GPs suffers most from localminima. Our proposed model using a scalar potential stilloutperforms the other two. B. Empirical proof-of-concept data Ω in the GP model is set as L = L = 0 . and L = 0 . . Theactual two-dimensional movement is performed in a rectangleof cm × cm and is hence well within the domain Ω . Thelearned scalar potential and its marginal variance (uncertainty)are shown in Figure 7b. Because the training data covers CCEPTED TO IEEE TRANSACTIONS ON ROBOTICS 11 . . (a) Training data (left) and validation data (right) − . − . . . · − · − . · − (b) Learned scalar potential (left) and its uncertainty (right) − − (c) The predicted magnetic field: H x , H y , and H z . . (d) Norm of the magnetic field Fig. 7: Illustration of the magnetic field data and the results from the GP approach from Algorithm 1 for the experimentdiscussed in Section VI-B. In all figures, the y -axis is − . , . . . , . m. The x -axis is − . , . . . , . m. The units in the fieldsurface plots are arbitrary due to normalization.almost the whole displayed area, the learned scalar potentialhas very low uncertainty. It is also possible to compute thepredicted magnetic field which is shown in Figure 7c. Forcompleteness, the intensity of these predicted magnetic fieldmeasurements are shown in Figure 7d. Although this quantityis only indirectly related to the outcome of the GP approach,it is frequently used in the remaining sections because of itseasy and intuitive visualization.By predicting the measurements at the locations of thevalidation data set, it is also possible to compute the RMSEon the validation data. In Figure 8 we visualize the RMSE as afunction of the number of basis functions used in Algorithm 1.We also compare to the RMSE from a full GP approach,that is using the same GP prior but without the Hilbert spaceapproximation scheme to speed up the inference. To allow forcomparison with a full GP approach—which suffers from ahigh computational complexity for large data sets—the datahas been downsampled to 5 Hz. As can be seen, already foraround m = 1000 basis functions, the quality of the estimatesfrom Algorithm 1 approaches that of the full GP approach. C. Mapping the magnetic field in a building
The third experiment was concerned with estimating a mapof the magnetic environment inside a building by only usinga smartphone for the data collection. The mapped venue islocated on the Aalto University campus, and a floor plan sketchis shown in Figure 9. For practical reasons, we limited ourinterest to the lobby which is approximately 600 m in size.For the measurements, we used an Apple iPhone 4 and itsbuilt-in 9-dof IMU (3-axis AKM AK8975 magnetometer). Allsensors were sampled at 50 Hz, and the data was streamedonline to a laptop computer for processing and storing. Thephone was held at waist-height and pointed towards theheading direction. m R M S E Full GPHilbert space approximation GP
CCEPTED TO IEEE TRANSACTIONS ON ROBOTICS 12
Domain boundaryTraining path N o r m o f e rr o r( µ T ) Fig. 9: A training (black) and validation (colored) free-walking path that was used in the experiment. Trajectorieswere collected by a mobile phone, and the magnetometer datawas corrected for gravitation direction and heading using theinertial sensors in the device. The domain boundaries for thereduced-rank method are shown by the dashed line.302 meters long, respectively. The number of magnetometerdata samples acquired along the paths were 9868, 10500, and12335. Prior to each acquisition, the phone magnetometer wascalibrated by a standard spherical calibration approach. Thecombined size of the training data set was n = 32703 . Forvalidation, we collected a walking path passing through thevenue (length 54 m, n = 2340 ).We considered a batch interpolation problem of creating amagnetic map of the lobby. The map was assumed static overtime, and we applied Algorithm 1 to the training data with m = 1024 basis functions. The optimized hyperparameterswere σ lin. ≈
575 ( µ T ) , σ SE /(cid:96) SE ≈
373 ( µ T ) , (cid:96) SE ≈ . m,and σ noise ≈ .
53 ( µ T ) . The coefficients of the inferred linear(bias) model corresponding to the linear covariance functionwas ( − . , . , − . .Figure 1 shows the interpolated magnetic field magnitude( (cid:107) f (cid:107) ) and the vector field components. The component-wisemaximum of the associated marginal standard deviation fieldsis visualized by the degree of transparency (with a standarddeviation of µ T being fully transparent). The overall shapeof the estimate agreed even when the model was trainedseparately with each of the training paths. To most part, thestrong fluctuations in the magnetic field are located near wallsor other structures in the building. The strong magnetic fieldin the open area in the lower right part of the floorplan wasidentified to most likely be due to a large supporting structureon the lower floor-level. We also used the model for predictingthe measurements along the validation path (see Fig. 9). Thecomponent-wise RMSEs were (2 . µ T , . µ T , . µ T ) andmean absolute errors (1 . µ T , . µ T , . µ T ) . Figure 9shows the norm of the error along the validation path. Themeasurement noise level of the magnetometer is in the mag-nitude of µ T and the uncertainty in the PDR estimatecontributes to the remaining variance.
Optical marker SmartphoneTrivisio IMUInvensense IMUDiddyBorg robot board
Fig. 10: The robot was built on a DiddyBorg robotics boardand controlled by a Raspberry Pi single-board computer. Thethree sensors (Invensense, Trivisio, and smartphone) providingmagnetometer readings were mounted on the top. The ref-erence locations were provided by a Vicon optical trackingsystem.
D. Online mapping
CCEPTED TO IEEE TRANSACTIONS ON ROBOTICS 13 µ T (a) Final interpolated magnitude map x -component y -component z -component − − (b) Vector field(c) Snapshot 1 (d) Snapshot 2 (e) Snapshot 3 (f) Snapshot 4 (g) Snapshot 5 Fig. 11: The GP interpolation task the robot was faced with. The final interpolation outcome of the magnitude field is shownin (a), and the different vector field components are shown in (b). Snapshots along the temporally updating field estimate areshown in (c–g) together with the path travelled since the previous update. The marginal variance (uncertainty) is visualized bythe degree of transparency.provided by the robot, but the interest in this experiment wasrather to focus on the interpolation of the magnetic field, notthe path estimation.The task was to map the magnetic field inside a markedregion roughly 6 m × m = 1024 .The length-scale, magnitude, and noise variance hyperpa-rameters were learned from the first two data sets ( n = 17980 vector valued observations) by maximizing with respect tomarginal likelihood (see Alg. 1). The obtained values were (cid:96) SE ≈ . m, σ SE /(cid:96) SE ≈
287 ( µ T ) and σ noise ≈ .
27 ( µ T ) .The linear model magnitude scale parameter was fixed to σ lin. = 500 ( µ T ) . The noise model is not only capturing thesensor measurement noise, but the entire mismatch betweenthe data and the model. This explains the rather large noisevariance. We also checked, that the hyperparameter estimatesremained stable when optimized using the rest of the data.Figure 11 shows the results for the static magnetic fieldexperiment. The estimate was updated continuously five timesa second, and we show five snapshots of the evolution ofthe magnetic field estimate in Figures 11c–11g (vector fieldmagnitude shown in figures). These snapshots also show thepath travelled since the previous snapshot. The alpha channelacts as a proxy for uncertainty; the marginal variance ofthe estimate is giving the degree of transparency. The finalmagnitude estimate—after iterating through all the n = 43029 observations—is shown in Figure 11a together with the vectorfield components in Figure 11b.The frontal part of the mapped region shows strong mag-netic activity, whereas the parts further back do not showas strong fields. Inspection of the venue suggested metallicpipelines or structures in the floor to blame (or thank) forthese features. In this particular case most parts of the effectis seen in the x -component. We repeated the reconstructionwith data collected from the Trivisio and smartphone sensors,and the results and conclusions remained unchanged. As a sup-plementary file to this paper, there is a video demonstratingthe online operation which has been sped-up × . CCEPTED TO IEEE TRANSACTIONS ON ROBOTICS 14
The last part of the experiment was dedicated to dynamical(time-dependent) modeling of the magnetic field. We usedall the data from the first part of the experiment to traina sequential model and used that as the starting point forchanging the field ( t = 0 ). During acquisition of data while therobot was driving around, we brought in two metallic toolboxshelves: first a larger toolbox shelf on wheels (Fig. 12c) andthen a smaller box (Fig. 12d). We acquired altogether some s of data ( n = 15513 ) of the changed environment.For encoding the assumptions of a changing magneticfield, we used Algorithm 2. The additional hyperparametercontrolling the temporal scale was fixed to (cid:96) time = 1 hour, thusencoding an assumption of slow local changes. This choice isnot restrictive, because the data is very informative about theabrupt changes.Figure 12a shows the evolution of the magnetic field compo-nents for one fixed location (indicated by a cross marker in thefigures). The two toolboxes induce clear changes in the localanomaly field, but the effects are restricted to the immediatevicinity of the boxes. Thus the spatio-temporal model onlygains information about the changed field, when the robotpasses by the location of interest. This effect is clearly visiblearound t = 20 s, and later on around t = 70 s and t = 130 s.After this the estimate stabilizes and only drifts around forthe remaining time. Even though the changes in the fieldcomponents appear clear in Figure 12a, they are only around µ T and thus only account for a variation of about 2% inthe scale of the entire field visualized in Figure 11b. Inducingmore noticeable changes in the magnetic field would requiremoving around larger structures (say an elevator). Yet, evenchanges this small can be tracked by the modeling approach.VII. D
ISCUSSION
In recent years, interest has emerged in mapping of themagnetic field by Gaussian processes for robot and pedestrianlocalization. This paper has aimed at presenting a new efficientmethod for mapping, but also at providing a study of bestpractices in using GPs in this context. Thus, we went throughthree different approaches for formulating GP priors for themagnetic field (independent, shared hyperparameters, and curl-free/scalar potential). These three models differ in the amountof prior knowledge encoded in the model, and the moreinformation available in the prior, the better the interpolationand extrapolation capabilities in the model—as long as thedata agrees with the assumptions. This was also demonstratedin Figure 4 and in Section VI-A.The methods presented in this paper are related to the use ofGaussian random field priors in inverse problems [65, 66]. Thisconnection has been explored from various points of views ( cf. [67]). However, the machine learning [7] way of interpretingthe GP priors indeed brings something new on the table—we are explicitly modeling the uncertainty in the field byusing a stochastic model, which has interesting philosophicalimplications. The formulation of the prior through a GPcovariance function provides both an intuitive and theoreticallyjustified way of encoding the information.The scalar potential approach is not the only way to buildthe model. It would also be possible to include disturbance − − − − − − − − x -component y -component z -componentTime (seconds) F i e l d c o m pon e n t s ( µ T ) (a) Magnetic field estimate at a fixed location(b) t = 0 s (c) t = 6 s (d) t = 27 s Fig. 12: (a) Evolution of the magnetic field at one spatial loca-tion over a time-course of 200 seconds. (b–d) The blue metallictoolboxes are brought in at the beginning of the experiment.The abrupt changes in the field estimate corresponds to timeinstances when the robot has passed the toolboxes and gainedinformation about the changed environment.in the model, which would model the effect of free currentsin the area or its boundaries. Furthermore, more complicatedassumptions of the temporal time-changing behavior of themagnetic field could be included in the temporal covariancefunction. For example, various degrees of smoothness orperiodicity could be included in the framework.This paper has been considering the ‘M’ (mapping) partin SLAM. The ‘L’ (localization) part based on these mapsis presented in Solin et al. [1]. The online mapping schemepresented in Algorithm 2 opens up for possibilities for simul-taneously building the map and localization within the map.As seen in the experiments presented both in this work and inSolin et al. [1], this appears feasible and provides an interestingdirection for further research.VIII. C
ONCLUSION
Small variations in the magnetic field can be used asinputs in various positioning and tracking applications. Inthis paper, we introduced an effective and practically feasibleapproach for mapping these anomalies. We encoded priorknowledge from Maxwell’s equations for magnetostatics into aBayesian non-parametric probabilistic model for interpolationand extrapolation of the magnetic field.The magnetic vector field components were modeled jointlyby a Gaussian process model, where the prior was associateddirectly with a latent scalar potential function. This ensuresthe field to be curl-free—a justified assumption in free spaces.This assumption couples the vector field components and
CCEPTED TO IEEE TRANSACTIONS ON ROBOTICS 15 additionally encodes the assumption of a baseline field withsmooth small-scale variations. We also presented connectionsto existing formulations for vector-valued Gaussian processmodels.In addition to constructing the model, we also presenteda novel and computationally efficient inference scheme forinterpolation and extrapolation using it. We built upon aLaplace operator eigenbasis approach, which falls natural tothe formulation of the model. The inference scheme ensuresa linear computational complexity with respect to the numberof observations of the magnetic field. We also extended themethod to an online approach with sequential updating of theestimate and time-dependent changes in the magnetic field.We presented four experiments demonstrating the feasibilityand practicality of the methods. A simulated experimentshowed the benefit of including additional knowledge fromphysics into the model, and a simple proof-of-concept exam-ple demonstrated the strength of the approximation schemein solving the model. Two real-world use cases were alsoconsidered: we mapped the magnetic field in a building onfoot using a smartphone, and demonstrated online mappingusing a wheeled robot.A
CKNOWLEDGEMENTS
This work was supported by grants from the Academyof Finland (266940, 273475, 308640), by CADICS, a Lin-naeus Center, and by the project Probabilistic modelling ofdynamical systems (contract number: 621-2013-5524), fundedboth by the Swedish Research Council (VR), by the SwedishFoundation for Strategic Research under the project Coop-erative Localization and by the EPSRC grant
Autonomousbehaviour and learning in an uncertain world (Grant number:EP/J012300/1).We are grateful for the help and equipment provided by theUAS Technologies Lab, Artificial Intelligence and IntegratedComputer Systems Division (AIICS) at the Department ofComputer and Information Science (IDA), Link¨oping Uni-versity, Sweden. We also acknowledge the computationalresources provided by the Aalto Science-IT project. Finally, wewould like to thank IndoorAtlas Ltd. for providing expertiseand for lending us equipment for the measurements.R
EFERENCES [1] A. Solin, S. S¨arkk¨a, J. Kannala, and E. Rahtu, “Terrain navi-gation in the magnetic landscape: Particle filtering for indoorpositioning,” in
Proceedings of the European Navigation Con-ference , 2016.[2] O. J. Woodman, “Pedestrian localisation for indoor environ-ments,” Ph.D. dissertation, University of Cambridge, UnitedKingdom, 2010.[3] J. D. Hol, “Sensor fusion and calibration of inertial sensors,vision, ultra-wideband and GPS,” Ph.D. dissertation, Link¨opingUniversity, Sweden, 2011.[4] J. J. Leonard and H. F. Durrant-Whyte, “Simultaneous mapbuilding and localization for an autonomous mobile robot,” in
Proceedings of the International Workshop of Intelligent Robotsand Systems (IROS) , vol. 3. IEEE, 1991, pp. 1442–1447.[5] H. Durrant-Whyte and T. Bailey, “Simultaneous localizationand mapping: Part I,”
IEEE Robotics & Automation Magazine ,vol. 13, no. 2, pp. 99–110, 2006. [6] A. O’Hagan, “Curve fitting and optimal design for prediction(with discussion),”
Journal of the Royal Statistical Society.Series B (Methodological) , vol. 40, no. 1, pp. 1–42, 1978.[7] C. E. Rasmussen and C. K. I. Williams,
Gaussian Processes forMachine Learning . MIT Press, 2006.[8] N. A. Cressie,
Statistics for Spatial Data . New York: Wiley-Interscience, 1993.[9] N. Cressie and C. K. Wikle,
Statistics for Spatio-Temporal Data .Hoboken, NJ: John Wiley & Sons, 2011.[10] M. P. Deisenroth, D. Fox, and C. E. Rasmussen, “Gaussianprocesses for data-efficient learning in robotics and control,”
IEEE Transactions on Pattern Analysis and Machine Intelli-gence , vol. 37, no. 2, pp. 408–423, 2015.[11] N. Wahlstr¨om, M. Kok, T. B. Sch¨on, and F. Gustafsson, “Mod-eling magnetic fields using Gaussian processes,” in
Proceedingsof the 38 th International Conference on Acoustics, Speech andSignal Processing (ICASSP) , 2013, pp. 3522–3526.[12] A. Solin and S. S¨arkk¨a, “Hilbert space methods for reduced-rankGaussian process regression,”
ArXiv preprint arXiv:1401.5508 ,2014.[13] S. S¨arkk¨a, A. Solin, and J. Hartikainen, “Spatiotemporal learn-ing via infinite-dimensional Bayesian filtering and smoothing,”
IEEE Signal Processing Magazine , vol. 30, no. 4, pp. 51–61,2013.[14] M. N. Nabighian, V. J. S. Grauch, R. O. Hansen, T. R. LaFehr,Y. Li, J. W. Peirce, J. D. Phillips, and M. E. Ruder, “Thehistorical development of the magnetic method in exploration,”
Geophysics , vol. 70, no. 6, pp. 33ND–61ND, 2005.[15] A. Guillen, P. Calcagno, G. Courrioux, A. Joly, and P. Ledru,“Geological modelling from field data and geological knowl-edge: Part II. Modelling validation using gravity and magneticdata inversion,”
Physics of the Earth and Planetary Interiors ,vol. 171, no. 1, pp. 158–169, 2008.[16] P. Calcagno, J.-P. Chil`es, G. Courrioux, and A. Guillen, “Ge-ological modelling from field data and geological knowledge:Part I. Modelling method coupling 3D potential-field interpola-tion and geological rules,”
Physics of the Earth and PlanetaryInteriors , vol. 171, no. 1, pp. 147–157, 2008.[17] F. Mackay, R. Marchand, and K. Kabin, “Divergence-freemagnetic field interpolation and charged particle trajectoryintegration,”
Journal of Geophysical Research: Space Physics ,vol. 111, no. A6, p. A06208, 2006.[18] B. K. Bhattacharyya, “Bicubic spline interpolation as a methodfor treatment of potential field data,”
Geophysics , vol. 34, no. 3,pp. 402–423, 1969.[19] V. Springel, “Smoothed particle hydrodynamics in astro-physics,”
Annual Review of Astronomy and Astrophysics ,vol. 48, pp. 391–430, 2010.[20] J. Haverinen and A. Kemppainen, “Global indoor self-localization based on the ambient magnetic field,”
Robotics andAutonomous Systems , vol. 57, no. 10, pp. 1028–1035, 2009.[21] B. Li, T. Gallagher, A. G. Dempster, and C. Rizos, “Howfeasible is the use of magnetic field alone for indoor position-ing?” in
Proceedings of the International Conference on IndoorPositioning and Indoor Navigation (IPIN) , 2012, pp. 1–9.[22] M. Angermann, M. Frassl, M. Doniec, B. J. Julian, andP. Robertson, “Characterization of the indoor magnetic field forapplications in localization and mapping,” in
Proceedings ofthe International Conference on Indoor Positioning and IndoorNavigation (IPIN) , 2012, pp. 1–9.[23] E. Le Grand and S. Thrun, “3-axis magnetic field mapping andfusion for indoor localization,” in
IEEE Conference on Mul-tisensor Fusion and Integration for Intelligent Systems (MFI) ,2012, pp. 358–364.[24] P. Robertson, M. Frassl, M. Angermann, M. Doniec, B. J. Julian,M. Garcia Puyol, M. Khider, M. Lichtenstern, and L. Bruno,“Simultaneous localization and mapping for pedestrians usingdistortions of the local magnetic field intensity in large indoorenvironments,” in
Proceedings of the International Conference
CCEPTED TO IEEE TRANSACTIONS ON ROBOTICS 16 on Indoor Positioning and Indoor Navigation (IPIN) . IEEE,2013, pp. 1–10.[25] M. Frassl, M. Angermann, M. Lichtenstern, P. Robertson, B. J.Julian, and M. Doniec, “Magnetic maps of indoor environmentsfor precise localization of legged and non-legged locomotion,”in
Proceedings of the IEEE/RSJ International Conference onIntelligent Robots and Systems (IROS) , 2013, pp. 913–920.[26] I. Vallivaara, J. Haverinen, A. Kemppainen, and J. R¨oning, “Si-multaneous localization and mapping using ambient magneticfield,” in
Proceedings of the IEEE Conference on MultisensorFusion and Integration for Intelligent Systems (MFI) , 2010, pp.14–19.[27] ——, “Magnetic field-based SLAM method for solving thelocalization problem in mobile robot floor-cleaning task,” in
Proceedings of the 15 th International Conference on AdvancedRobotics (ICAR) , 2011, pp. 198–203.[28] S. T. O’Callaghan and F. T. Ramos, “Gaussian process occu-pancy maps,”
The International Journal of Robotics Research ,vol. 31, no. 1, pp. 42–62, 2012.[29] M. Smith, I. Posner, and P. Newman, “Adaptive compression for3D laser data,”
The International Journal of Robotics Research ,vol. 30, no. 7, pp. 914–935, 2011.[30] S. Kim and J. Kim, “Hierarchical Gaussian processes for robustand accurate map building,” in
Proceedings of AustralasianConference on Robotics and Automation , 2015, pp. 117–124.[31] F. Ramos and L. Ott, “Hilbert maps: Scalable continuousoccupancy mapping with stochastic gradient descent,”
The In-ternational Journal of Robotics Research , vol. 35, no. 14, pp.1717–1730, 2016.[32] R. Senanayake, L. Ott, S. O’Callaghan, and F. T. Ramos,“Spatio-temporal Hilbert maps for continuous occupancy rep-resentation in dynamic environments,” in
Advances in NeuralInformation Processing Systems 29 . Curran Associates, Inc.,2016, pp. 3925–3933.[33] T. Vidal-Calleja, D. Su, F. De Bruijn, and J. V. Miro, “Learningspatial correlations for Bayesian fusion in pipe thickness map-ping,” in
Proceedings of the IEEE International Conference onRobotics and Automation (ICRA) , 2014, pp. 683–690.[34] C. H. Tong, P. Furgale, and T. D. Barfoot, “Gaussian processGauss–Newton for non-parametric simultaneous localizationand mapping,”
The International Journal of Robotics Research ,vol. 32, no. 5, pp. 507–525, 2013.[35] B. Ferris, D. H¨ahnel, and D. Fox, “Gaussian processes for signalstrength-based location estimation,” in
Proceedings of Robotics:Science and Systems , 2006.[36] S. Barkby, S. B. Williams, O. Pizarro, and M. V. Jakuba,“Bathymetric particle filter SLAM using trajectory maps,”
TheInternational Journal of Robotics Research , vol. 31, no. 12, pp.1409–1430, 2012.[37] T. D. Barfoot, C. H. Tong, and S. S¨arkk¨a, “Batch continuous-time trajectory estimation as exactly sparse Gaussian processregression,” in
Proceedings of Robotics: Science and Systems ,2014.[38] S. Anderson, T. D. Barfoot, C. H. Tong, and S. S¨arkk¨a,“Batch nonlinear continuous-time trajectory estimation as ex-actly sparse Gaussian process regression,”
Autonomous Robots ,vol. 39, no. 3, pp. 221–238, 2015.[39] J. Qui˜nonero-Candela and C. E. Rasmussen, “A unifying viewof sparse approximate Gaussian process regression,”
Journal ofMachine Learning Research , vol. 6, pp. 1939–1959, 2005.[40] M. L´azaro-Gredilla, J. Qui˜nonero-Candela, C. E. Rasmussen,and A. R. Figueiras-Vidal, “Sparse spectrum Gaussian processregression,”
Journal of Machine Learning Research , vol. 11, pp.1865–1881, 2010.[41] J. Hensman, N. Durrande, and A. Solin, “Variational Fourierfeatures for Gaussian processes,”
Journal of Machine LearningResearch , accepted.[42] C. J. Paciorek, “Bayesian smoothing with Gaussian processesusing Fourier basis functions in the spectralGP package,”
Jour- nal of Statistical Software , vol. 19, no. 2, pp. 1–38, 2007.[43] J. Fritz, I. Neuweiler, and W. Nowak, “Application of FFT-based algorithms for large-scale universal kriging problems,”
Mathematical Geosciences , vol. 41, no. 5, pp. 509–533, 2009.[44] J. Hartikainen and S. S¨arkk¨a, “Kalman filtering and smoothingsolutions to temporal Gaussian process regression models,”in
Proceedings of the International Workshop on MachineLearning for Signal Processing (MLSP) . IEEE, 2010, pp. 379–384.[45] S. Reece and S. Roberts, “An introduction to Gaussian processesfor the Kalman filter expert,” in
Proceedings of the 13th Inter-national Conference on Information Fusion (FUSION) , 2010,pp. 1–9.[46] M. Osborne, “Bayesian Gaussian processes for sequential pre-diction, optimisation and quadrature,” Ph.D. dissertation, Uni-versity of Oxford, United Kingdom, 2010.[47] M. F. Huber, “Recursive Gaussian process: On-line regressionand learning,”
Pattern Recognition Letters , vol. 45, pp. 85–91,2014.[48] M. A. ´Alvarez, D. Luengo, and N. D. Lawrence, “Linear latentforce models using Gaussian processes,”
IEEE Transactions onPattern Analysis and Machine Intelligence , vol. 35, no. 11, pp.2693–2705, 2013.[49] M. A. ´Alvarez and N. D. Lawrence, “Latent force models,” in
Proceedings of the 12 th International Conference on ArtificialIntelligence and Statistics , ser. JMLR W&CP, vol. 5, 2009, pp.9–16.[50] J. Hartikainen and S. S¨arkk¨a, “Sequential inference for latentforce models,” in
Proceedings of the International Conferenceon Uncertainty in Artificial Intelligence , 2011, pp. 311–318.[51] S. S¨arkk¨a and J. Hartikainen, “Infinite-dimensional Kalmanfiltering approach to spatio-temporal Gaussian process regres-sion,” in
Proceedings of the 15 th International Conference onArtificial Intelligence and Statistics , ser. JMLR W&CP, vol. 22,2012, pp. 993–1001.[52] R. F. Curtain and A. J. Pritchard,
Infinite Dimensional LinearSystems Theory . Springer–Verlag, 1978.[53] J. D. Jackson,
Classical Electrodynamics , 3rd ed. New York:Wiley, 1999.[54] J. Vanderlinde,
Classical Electromagnetic Theory . Dordrecht:Springer, 2004.[55] D. J. Griffiths and R. College,
Introduction to Electrodynamics .Upper Saddle River, NJ: Prentice hall, 1999.[56] N. Wahlstr¨om, “Modeling of magnetic fields and extendedobjects for localization applications,” Ph.D. dissertation,Link¨oping University, Sweden, 2015.[57] A. Kemppainen, J. Haverinen, I. Vallivaara, and J. R¨oning,“Near-optimal exploration in Gaussian process SLAM: Scalableoptimality factor and model quality rating,” in
Proceedings ofthe 5 th European Conference on Mobile Robots (ECMR) , 2011,pp. 283–289.[58] J. Jung, T. Oh, and H. Myung, “Magnetic field constraintsand sequence-based matching for indoor pose graph SLAM,”
Robotics and Autonomous Systems , vol. 70, pp. 92–105, 2015.[59] A. Viseras Ruiz and C. Olariu, “A general algorithm forexploration with Gaussian processes in complex, unknown envi-ronments,” in
Proceeding of the IEEE International Conferenceon Robotics and Automation (ICRA) , 2015, pp. 3388–3393.[60] E. J. Fuselier, Jr, “Refined error estimates for matrix-valued ra-dial basis functions,” Ph.D. dissertation, Texas A&M University,2007.[61] L. Baldassarre, L. Rosasco, A. Barla, and A. Verri, “Vectorfield learning via spectral filtering,” in
Machine Learning andKnowledge Discovery in Databases , ser. Lecture Notes inComputer Science. Springer, 2010, vol. 6321, pp. 56–71.[62] M. A. ´Alvarez, L. Rosasco, and N. D. Lawrence, “Kernels forvector-valued functions: A review,”
Foundations and Trends inMachine Learning , vol. 4, no. 3, pp. 195–266, 2012.[63] C. Jidling, N. Wahlstr¨om, A. Wills, and T. B. Sch¨on, “Linearly
CCEPTED TO IEEE TRANSACTIONS ON ROBOTICS 17 constrained Gaussian processes,” in
Advances in Neural Infor-mation Processing Systems (NIPS) 30 , 2017.[64] S. S¨arkk¨a,
Bayesian Filtering and Smoothing . CambridgeUniversity Press, 2013.[65] A. Tarantola,
Inverse Problem Theory and Methods for ModelParameter Estimation . SIAM, 2004.[66] J. Kaipio and E. Somersalo,
Statistical and ComputationalInverse Problems . Springer, 2005.[67] S. S¨arkk¨a, “Linear operators and stochastic partial differentialequations in Gaussian process regression,” in
Artificial NeuralNetworks and Machine Learning – ICANN 2011 , ser. LectureNotes in Computer Science, vol. 6792. Springer, 2011, pp.151–158.
Arno Solin received his M.Sc. (Tech.) degree (withdistinction) in Engineering Physics and Mathematicsand D.Sc. (Tech.) degree (with distinction) in com-putational science from Aalto University, Finland, in2012 and 2016, respectively.Dr. Solin is an Academy of Finland PostdoctoralResearcher with Aalto University, and Technical Ad-visor of IndoorAtlas Ltd. His research interests arein probabilistic inference for temporal, spatial, andspatio-temporal models with applications in sensorfusion for tracking and navigation, brain imaging,and machine learning problems.
Manon Kok received M.Sc. degrees in AppliedPhysics and in Philosophy of Science, Technologyand Society (2009 and 2007, respectively), both fromthe University of Twente. From 2009 until 2011she worked as a Research Engineer at Xsens Tech-nologies. She received the PhD degree in AutomaticControl from Link¨oping University in 2017.Dr. Kok is currently a Postdoc in the MachineLearning Group of the Computational and BiologicalLearning Lab at the University of Cambridge. Herresearch interests are in the fields of probabilisticinference for sensor fusion, signal processing and machine learning.
Niklas Wahlstr¨om received the M.Sc. degree inApplied Physics and Electrical Engineering in 2010and Ph.D. degree in automatic control in 2015, bothfrom Link¨oping University, Sweden.Dr. Wahlstr¨om is since 2016 a Researcher atthe Department of Information Technology, UppsalaUniversity. His research interests include machinelearning, deep learning, sensor fusion, localizationand mapping, especially using magnetic sensors.
Thomas B. Sch¨on received the M.Sc. degree inApplied Physics and Electrical Engineering in 2001,the B.Sc. degree in Business Administration 2001and the Ph.D. degree in Automatic Control in 2006,all from Link¨oping University.Dr. Sch¨on is Professor of the Chair of AutomaticControl in the Department of Information Technol-ogy at Uppsala University. His main research interestis nonlinear inference problems, especially withinthe context of dynamical systems, solved usingprobabilistic methods, more specifically sequentialMonte Carlo, particle MCMC and graphical models. He is active within thefields of machine learning, signal processing and automatic control. He is aSenior Member of the IEEE.