©2025 IEEE. Personal use of this material is permitted. Permission from IEEE must be obtained for all other uses, in any current or future media, including reprinting/republishing this material for advertising or promotional purposes, creating new collective works, for resale or redistribution to servers or lists, or reuse of any copyrighted component of this work in other works.

## This paper has been accepted for publication in the IEEE Robotics and Automation Letters . Please cite this paper as:

## @article{saavedra2025hef,

| title   | = {The Harmonic Exponential Filter for Nonparametric Estimation on Motion Groups},                          |
|---------|-------------------------------------------------------------------------------------------------------------|
| author  | = {Saavedra-Ruiz, Miguel and Parkison, Steven A. and Arora, Ria and Forbes, James Richard and Paull, Liam}, |
| year    | = 2025,                                                                                                     |
| journal | = {IEEE Robotics and Automation Letters},                                                                   |
| volume  | = {},                                                                                                       |
| number  | = {},                                                                                                       |
| pages   | = {1--8},                                                                                                   |
| doi     | = {10.1109/LRA.2025.3527346}                                                                                |

}

## The Harmonic Exponential Filter for Nonparametric Estimation on Motion Groups

Miguel Saavedra-Ruiz, Steven A. Parkison, Ria Arora, James Richard Forbes, and Liam Paull

Abstract -Bayesian estimation is a vital tool in robotics as it allows systems to update the robot state belief using incomplete information from noisy sensors. To render the state estimation problem tractable, many systems assume that the motion and measurement noise, as well as the state distribution, are unimodal and Gaussian. However, there are numerous scenarios and systems that do not comply with these assumptions. Existing nonparametric filters that are used to model multimodal distributions have drawbacks that limit their ability to represent a diverse set of distributions. This paper introduces a novel approach to nonparametric Bayesian filtering on motion groups, designed to handle multimodal distributions using harmonic exponential distributions. This approach leverages two key insights of harmonic exponential distributions: a) the product of two distributions can be expressed as the element-wise addition of their log-likelihood Fourier coefficients, and b) the convolution of two distributions can be efficiently computed as the tensor product of their Fourier coefficients. These observations enable the development of an efficient and asymptotically exact solution to the Bayes filter up to the band limit of a Fourier transform. We demonstrate our filter's performance compared with established nonparametric filtering methods across simulated and real-world localization tasks. 1

## I. INTRODUCTION

Recursive Bayesian filtering is a powerful approach for state estimation in robotics. It allows the fusion of different sensor streams while producing belief estimations over the state that can be used in downstream tasks such as planning, control, and mapping. To render the state estimation problem computationally tractable, certain assumptions are frequently invoked, including the presumption that the state distribution, as well as the motion and measurement noise, are unimodal and Gaussian distributed. However, there are many scenarios where a unimodal Gaussian belief fails to place a high probability (density) over the true state. Estimates involving uncertain data associations [1], propagating errors through rigid-body motions [2], and measurements that suffer from multipath errors [3] and symmetries [4] all result in distributions that are distorted or multimodal.

Gaussian filters have been extended to handle some nonGaussian distributions by modeling the error of the state in

Miguel Saavedra-Ruiz, Steven A. Parkison, Ria Arora, and Liam Paull are with the Department of Computer Science and Operations Research, Universit´ e de Montr´ eal, Montr´ eal, QC, Canada, and with Mila - Quebec AI Institute, Montr´ eal, QC, Canada.

James Richard Forbes is with the department of Mechanical Engineering, McGill University, Montr´ eal, QC, Canada.

This research was partially funded by an FRQNT B2X Scholarship [MS]. We acknowledge the support of the Natural Sciences and Engineering Research Council of Canada (NSERC), [funding reference numbers RGPIN2024-04956 and RGPIN-2023-03621], as well as CIFAR under the CCAI Chairs program.

1 Project page https://montrealrobotics.ca/hef/

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00907v3-figure-f3ddc3bc-7ab6-4a61-a708-14086ad7df22.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00907v3-figure-6ab56187-cd24-4ade-bc3f-567962a08949.jpg)

(b) Convolution

Fig. 1: The product (a) and convolution (b) operations between two distributions play an important role in recursive Bayesian filtering. The examples here are shown with harmonic exponential distributions on the circle S 1 .

Lie groups using parametric distributions [5], [6]. However, to propagate errors, these methods rely on approximations or parametric assumptions to maintain a tractable solution.

Nonparametric filtering offers an alternative to Gaussianbased filters, such as particle filters [7]-[9] and histogram filters [10]-[12]. However, these approaches are not free of limitations. Particle filters are susceptible to particle deprivation [1, Chapter 4], where the correct state is no longer represented within the particle set. Similarly, histogram filters require resource-intensive computations, including convolutions during motion propagation, and are only accurate up to the quantization error when computing probability bins.

This paper draws on seminal ideas in Fourier analysis for optimal filtering on rotation groups [13]-[16] and proposes a novel nonparametric, multimodal filtering approach for motion groups, leveraging harmonic exponential distributions [17]. We name this method the harmonic exponential filter (HEF) and it relies on three key benefits of harmonic exponential distributions. First, the prediction step exploits the Fourier convolution theorem, yielding an efficient tensor product instead of a convolution during motion propagation. Second, the update step of the filter becomes an element-wise addition. Finally, the Fourier coefficient representation of the harmonic exponential distribution represents more naturally smoothlyvarying density functions when compared to histogram binning. These steps are illustrated in Fig 1.

In summary, our main contributions are:

- · An asymptotically exact Bayesian filter that operates on a compact subset of motion groups, for band-limited prior, motion, and measurement functions.
- · A nonparametric filtering approach for noncommutative motion groups that captures multimodal distributions.

- · An efficient motion step implementation for motion groups that leverages the Fourier convolution theorem.
- · An experimental validation of our approach against classical filtering methods.

We show that our method outperforms the baselines in three non-Gaussian filtering tasks using one simulated and two realworld datasets.

## II. RELATED WORK

Various approaches have extended Bayesian filtering beyond the unimodal Gaussian assumptions of the Kalman filter (KF). For instance, the KF has been adapted to operate with a mixture of parameterized distributions instead of a unimodal Gaussian [18]. Gilitschenski et al. [19] adapted a mixture model to use direction statistics on the manifold of rotation groups. However, a key challenge with such approaches is the growth of mixture components with each iteration. Common approaches to overcome this are based on distance metrics to merge mixture components [20] or resampling [21].

The KF has also been extended to exploit symmetries of the system by characterizing the state of the robot using Lie groups [22]. If the state error is modeled in the Lie algebra of motion groups with a Gaussian belief and then mapped back to the group with the exponential map, a non-Gaussian posterior is obtained on the manifold [2], [5]. However, to propagate the error state over time and attain a tractable solution, these methods rely on second-order approximations. The invariant extended Kalman filter (IEKF) [6] also models the state error in the Lie algebra with a Gaussian density, which restricts it to unimodal distributions.

Unlike the various KF approaches, our method is nonparametric and can fit any distribution in a broader family of distributions while preserving system symmetries. A wellknown nonparametric filtering approach is the histogram filter, which approximates the posterior belief by discretizing the state space and assigning a probability value to each bin [10]. These methods can be adapted to the topology and action of the group [23] and have been widely studied to model distributions on rotation groups [11], [12]. Our method resembles histogram filters as we require discrete samples of the manifold to represent the posterior. However, our formulation operates in both rotation and motion groups and is more computationally efficient during the propagation step than histogram filters. What is more, our approach is exact up to the band limit of a discrete Fourier transform.

Another family of nonparametric and multimodal filters is the particle filter which represents the posterior belief using a discrete set of particles [7]. These methods can model general distributions, but are known to suffer from particle deprivation, particularly in high-dimensional state spaces [1, Chapter 4]. Research efforts have been devoted to addressing this problem by proposing better sampling schemes [8], [9], and by changing how particles evolve in the prediction step to better match the posterior distribution [24]. In contrast, our method does not use particles to model the posterior belief and does not suffer from particle deprivation.

Kernel density estimates are a prevalent method to represent nonparametric probability density functions. They can be adapted to handle belief propagation through probabilistic graphical models [25], and the product of distributions within a filtering framework [26]. Similarly, learning-based methods have been used to approximate sampling from a non-Gaussian posterior in SLAM by modeling conditional distributions in a Bayes tree with normalizing flows [27].

Harmonic analysis for optimal filtering has been widely studied on rotation groups. Willsky [13] derived state estimation equations for filtering on the circle using infinite Fourier coefficients and proposed finite-dimensional suboptimal implementations [14]. Filtering with 'Fourier densities' has since expanded to groups like SO N ( ) [16]. Similarly to our work, Lo et al. [15] highlighted the advantages of rotational Fourier densities for the update step on SO (3) . Fourier densities have also been applied to non-linear filtering, with Brunn et al. [28], [29] formulating Bayesian prediction and update steps on d-dimensional Euclidean spaces. Unlike these methods, our approach extends these ideas beyond compact rotation groups or Euclidean spaces and applies them to motion groups.

Pfaff et al. [30] proposed a method conceptually similar to ours, called Fourier filtering, which was applied to the unit circle and sphere [31], [32]. They also showed their method applied to the problem of multimodal density estimation on the hypertorus, and addressed arbitrary transition densities in the prediction step [33], [34].

Although these approaches resemble ours in that probability densities are represented by Fourier coefficients, our approach extends this concept in several meaningful ways. First, since we model density functions in the log-density space, the product of distributions becomes an element-wise addition instead of the more intricate tensor product needed for Fourier filtering. Second, by representing distributions with the Fourier coefficients derived from the log-density function, the corresponding probability density function is guaranteed to be non-negative [35]. Lastly, our approach can represent densities on compact subsets of motion groups SE N ( ) , particularly any matrix Lie group and not just on a circle or a sphere (the latter of which is not a Lie group).

## III. HARMONIC EXPONENTIAL FILTERING

Fourier analysis provides a way to parameterize wellbehaved functions in R N (i.e. one that is square integrable L 2 on V ⊆ R N and has a finite number of discontinuities in any finite interval in R N ) using an infinite sum of trigonometric basis functions. This approach has proven to be useful for various engineering applications for several reasons. Moreover, the convolution theorem shows that computing integrals in spatial domains can be replaced by an efficient tensor product in the frequency domain. A generalized form of Fourier analysis can be used on compact subsets of motion groups. The topic is presented in depth by Chirikjian and Kyatkin in [36, Chapters 8,10]. This section presents an overview of our approach. First by providing a brief introduction to Fourier analysis on motion groups, particularly SE (2) (Sec. III-A), then showing how it is used in the family of harmonic exponential distributions to represent a varied set of densities (Sec. III-B), and finally how these concepts can be combined to perform Bayesian filtering (Sec. III-C).

## A. Harmonic Analysis on Motion Groups

Motion groups are a specific category of Lie groups [22] that represent rigid-body motions within the Euclidean space R N . This classification includes other important subgroups, such as rotation groups SO N ( ) and rigid-body-transformation groups SE N ( ) . Particularly, these are Lie groups because they satisfy the group axioms (see [36, Chapter 7]) and are a differentiable manifold. A group consists of a set G and a binary operator ◦ that satisfies the closure property; ∀ g, h ∈ G , their composition g ◦ h also belongs to G .

Conducting Fourier analysis on functions supported on Lie groups, such as SE (2) , requires a collection of carefully chosen basis functions. These basis functions are in turn the matrix elements of irreducible unitary representations (IURs). A definition of IURs is provided below.

A matrix representation for SE (2) , denoted R λ ( g ) , is an infinite-dimensional matrix function that satisfies the homomorphism property R λ ( g ◦ h ) = R λ ( g R ) λ ( h ) , ∀ g, h ∈ G = SE (2) , where λ ∈ R + indexes the representation. Here, g ◦ h represents composition of group elements, while R λ ( g R ) λ ( h ) denotes standard matrix multiplication.

A unitary representation U λ mn ( g ) , with matrix elements m,n ∈ Z , is a representation that satisfies ⟨ U λ ( g φ , U ) 1 λ ( g φ ) 2 ⟩ = ⟨ φ , φ 1 2 ⟩ for all m,n ∈ Z . Here, the inner product is defined as ⟨ φ , φ 1 2 ⟩ = 1 2 π ∫ 2 π 0 φ 1 ( ψ ) H φ 2 ( ψ d ) ψ where ( ) · H denotes the Hermitian transpose and φ i is a function defined on S 1 . An irreducible unitary representation (IUR) is a unitary representation that cannot be block-diagonalized further into another representation. A representation can be block-diagonalized if it is equivalent to the direct sum of other representations [36]. An IUR is represented as T λ ( g ) ∈ R m n × .

IURs are the foundation for performing harmonic analysis on motion groups. If we have a function, f : G → R , that is square integrable, and an IUR of motion group G T , λ ( G ) , then we can perform a generalized Fourier transform to compute the frequency components of f as

$$\mathcal { F } [ f ] _ { m n } ^ { \lambda } = \int _ { G } f ( g ) T _ { m n } ^ { \lambda } ( g ) ^ { \mathrm H } d ( g ) = \eta _ { m n } ^ { \lambda }, \quad ( 1 )$$

where d g ( ) is the bi-invariant volume element of G . Matrix elements of the Fourier coefficients η λ ∈ R m n × , and IUR T λ ( g ) ∈ R m n × are indexed by their row and column indices m and n . Exploiting the fact that T g ( ) is a unitary representation, the generalized inverse Fourier transform is

$$\cdots & \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \\ & \mathcal { F } ^ { - 1 } [ \eta ] ( g ) = \sum _ { n, m \in \mathbb { Z } } \int _ { 0 } ^ { \infty } \eta _ { m n } ^ { \lambda } T _ { n m } ^ { \lambda } ( g ) \lambda d _ { \lambda } = f ( g ), \quad ( 2 ) \quad \begin{array} { c } \eta \text{ aus} \\ \eta \end{array} \cdots \\ \cdots & \quad \cdot \quad \cdot \quad \cdot \quad \cdot \end{array} \cdot$$

with F F [ -1 [ η ] ] = η and F -1 [ F [ f ]] = f . Chirikjian and Kyatkin provide a proof of the Fourier inversion formula for SE (2) in [36, Section 10.4.2] and offer a detailed explanation of the fast Fourier transform (FFT) over 2D motion groups in [36, Section 11.2].

In this work we use the IUR of SE (2) presented in [36], [37]. For g x, y, θ ( ) ∈ SE (2) with ( x, y ) ∈ R 2 and θ ∈ [0 , 2 π ] parameterizing an element of SE (2) , a variable ψ ∈ [0 , 2 π ) ,

Fig. 2: The relationship between probability density p g ( ) , logdensity ln ( p g ( )) , and their spectral coefficients η and M . We model the log-density ln ( p g ( )) using η . To recover the p g ( ) , we compute the inverse FFT of η and exponentiate the result. The spectral coefficients M are obtained by taking the FFT of p g ( ) and are used to compute convolutions.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00907v3-figure-d534179e-6ed1-4550-a58a-a8ad94d1644f.jpg)

bi-invariant volume dµ ≜ dxdy dθ dψ and complex unit i , the Fourier coefficients η λ mn are computed as

$$\stackrel { - } { \mathfrak { e }, \ g } \stackrel { - } { \mathfrak { e } } h \quad \eta _ { m n } ^ { \lambda } = \int _ { x, y, \theta, \psi } f ( x, y, \theta ) e ^ { i \psi n } e ^ { i ( x \lambda \cos \psi + y \lambda \sin \psi ) } e ^ { - i m ( \psi - \theta ) } d \mu.$$

To make the integral in (1) tractable, we assume that f is supported on compact subsets of the group and use a discrete subset of function samples. For SE (2) , the orientation component SO (2) is compact, and for translation, we consider a predefined compact subset of R 2 , as described in [36]. Additionally, we assume f is bandlimited with λ ∈ (0 , λ max ] and truncate the infinite sums in (2) to -k ≤ m,n ≤ k . We will omit indices and denote the generalized Fourier transform for compact subsets of motion groups as η = F [ f ] .

## B. Harmonic Exponential Distributions

Probability distributions, which must abide by the probability axioms, are non-negative functions, p g ( ) ≥ 0 , ∀ g ∈ G , which integrate to 1, ∫ G p g dg ( ) = 1 . Cohen and Welling [17] presented the harmonic exponential distribution as a class of probability distributions supported on groups and homogeneous spaces. Its density function is defined as

$$p _ { \eta } ( g ) \stackrel { \triangle } { = } p ( g \, | \, \eta ) = \frac { 1 } { Z _ { \eta } } \exp ( \eta \cdot T ( g ) ), \quad \quad ( 3 )$$

where η are the (spectral) parameters, T g ( ) is the IUR of group G Z , η is the normalizing constant, and ( ) · is a shorthand notation for η T · ( g ) ≜ F -1 [ η ]( g ) . Both notations will be used interchangeably to remain consistent with [17]. The parameters η and Z η can be computed via harmonic analysis from ϕ g ( | η ) ≜ exp( η · T g ( )) as

$$\cdot \cdot \cdot \cdot \cdot \cdot \\ Z _ { \eta } = \int _ { G } \phi ( g \, | \, \eta ) T _ { 0 0 } ^ { 0 } ( g ) \, d ( g ) = \mathcal { F } [ \phi ] _ { 0 0 } ^ { 0 }. \quad ( 4 \mathbf b )$$

$$\int _ { \Delta } \int _ { \Delta } \int _ { \Delta } \int _ { \Delta } \int _ { \Delta } \int _ { \Delta } \int _ { \Delta } \int _ { \Delta } \int _ { \Delta } \int _ { \Delta } \int _ { \Delta } \int _ { \Delta } \int _ { \Delta } \int _ { \Delta } \int _ { \Delta } \int _ { \Delta } \int _ { \Delta } \int _ { \Delta } \int _ { \Delta } \int _ { \Delta } \int _ { \Delta } \int _ { \Delta } \int _ { \Delta } \int _ { \Delta } \int _ { \Delta } \int _ { \Delta } \int _ { \Delta } \int _ { \Delta } \int _ { \Delta } \int _ {\Delta } \int _ { \Delta } \int _ { \Delta } \int _ { \Delta } \int _ { \Delta } \int _ { \Delta } \int _ { \Delta } \int _ { \Delta } \int _ { \Delta } \int _ { \Delta } \int _ { \Delta } \int _ { \Delta } \int _ { \Delta } \int _ { \Delta } \int _ { \Delta } \int _ { \Delta } \int _ { \Delta } \int _ { \Delta } \int _ { \Delta } \int _ { \Delta } \int _ { \Delta } \int _ { \Delta } \int _ { \Delta } \int _ { \Delta } \int _ { \Delta } \int _ { \Delta } \int _ { \Delta } \int _ { \Delta } \int _ { \Delta } \left ( 4 \mathbf a \right )$$

with F [ ϕ ] 0 00 being the zeroth Fourier coefficient and T 0 00 ( g ) = 1 . Representing the probability density with Fourier coefficients in log-density space, ensures that the probability density function is non-negative for any set of parameters η , and the normalizing constant Z η ensures that it integrates to 1. From this point forward, we use the shorthand notation p a ( g ) ≜

p η a ( g ) . The harmonic exponential distribution provides an elegant solution for the product of two distributions, achieved through the element-wise addition of their Fourier coefficients:

$$p _ { a } ( g ) p _ { b } ( g ) & = \frac { 1 } { Z _ { \eta } ^ { a } Z _ { \eta } ^ { b } } \exp ( ( \eta _ { a } + \eta _ { b } ) \cdot T ( g ) ) \quad \text{$density$} \\ & = p ( g | \eta _ { a } + \eta _ { b } ). \quad \text{$p_{z_{t}|x_{t}$} } \text{functio} \\ \text{Convolving two distributions is inherently more compli-coeffs}$$

Convolving two distributions is inherently more complicated. However, harmonic analysis provides a connection between convolutions and spectral coefficients,

$$( p _ { a } * p _ { b } ) ( g ) & = \int _ { G } p _ { a } ( h ) p _ { b } ( h ^ { - 1 } \circ g ) d ( h ) & ( 6 ) & \frac { \sqrt { \text{Algorit} } } { \text{Im} } \\ & = \mathcal { F } ^ { - 1 } [ \mathcal { F } [ p _ { b } ] \mathcal { F } [ p _ { a } ] ] ( g ) = \mathcal { F } ^ { - 1 } [ M _ { b } M _ { a } ] ( g ) & ( p _ { u _ { t } }, \text{me} \\ \text{with } \mathcal { F } [ n _ { n _ { t } } ] & = \int _ { n _ { t } } ( a \varsigma T ( a \varsigma d ( a ) } \, = \, M _ { n _ { t } } \text{ in} \sigma \text{ the Fourier} \quad \text{// } \mathbb { P }$$

with F [ p a ] = ∫ G p a ( g T ) ( g d g ) ( ) = M a being the Fourier coefficients of the density p a , instead of a log-density function as in (4a). Due to the nonconmutative nature of SE (2) , encoded through T g ( ) , the order of the Fourier coefficients is reversed in (6). Note that M a and M b are composed via matrix products: M λ = M M λ b λ a [36, Section 10.4]. The relation between these spaces is shown in Fig. 2.

## C. Bayesian Filtering

Following the notation of Thrun et al. [1, Chapter 2], equations (5) and (6) allow us to incorporate the harmonic exponential distribution into a Bayesian filter:

$$\begin{array} { c } \cdot & & \cdot \\ b e l ( x _ { t } ) \overset { \triangle } { = } p ( x _ { t } | z _ { 1 \colon t }, u _ { 1 \colon t } ) & = \frac { \cdot } { p ( z _ { t } | x _ { t } ) p ( x _ { t } | z _ { 1 \colon t - 1 }, u _ { 1 \colon t } ) } & & \cdot \\ p ( z _ { t } | z _ { 1 \colon t - 1 }, u _ { 1 \colon t } ) & & \cdot \\ & = \gamma p ( z _ { t } | x _ { t } ) \overline { b e l } ( x _ { t } ), & & \cdot \\ & & \quad \text{The He} \\ \text{where the predicted belief over the state $x_{t}$ after incorporating} & & spaces t \end{array}$$

where the predicted belief over the state x t after incorporating the motion model but before incorporating the measurement z t is denoted as bel ( x t ) ≜ p x ( t | z 1: t -1 , u 1: t ) , and the posterior belief after incorporating z t as bel ( x t ) .

The prediction step involves defining the motion model, which represents the likelihood of the current pose x t ∈ SE N ( ) given the previous pose x t -1 ∈ SE N ( ) and control input u t , i.e., p x ( t | u , x t t -1 ) . Defining the motion model as x t = x t -1 ◦ u t , where u t is a random variable with a known density, allows us to compute the relative motion as x -1 t -1 ◦ x t = u t . Therefore, the relative motion likelihood is p u t ( x -1 t -1 ◦ x t ) , giving us the prediction step

$$p _ { u _ { t } } ( x _ { t - 1 } \ominus x _ { t } ), \, \text{giving us up current program} \, \stackrel { \dots } { \text{intracta} } } \, \overline { b e l } ( x _ { t } ) & = \int _ { S E ( N ) } p ( x _ { t } | u _ { t }, x _ { t - 1 } ) b e l ( x _ { t - 1 } ) \, d x _ { t - 1 } \quad \text{ in Sec.} \\ & = \int _ { S E ( N ) } b e l ( x _ { t - 1 } ) p _ { u _ { t } } ( x _ { t - 1 } ^ { - 1 } \circ x _ { t } ) \, d x _ { t - 1 } \quad ( 8 ) \, \text{ certain} \\ & = \mathcal { F } ^ { - 1 } [ \mathcal { F } [ p _ { u _ { t } } ] \mathcal { F } [ b e l _ { t - 1 } ] ] ( x _ { t } ). \\ \text{The update step of a Bayes}, \, \text{ filter is the product of} \quad \text{When}$$

The update step of a Bayes' filter is the product of two probability distributions, p z ( t | x t ) bel ( x t ) . If we perform Fourier analysis in the domain x t of the measurement loglikelihood function we can get a set of Fourier coefficients η z t | x t ≜ F [ ln ( p z ( t | x t )] , and the measurement likelihood function becomes p z ( t | x t ) = 1 /Z η exp( η z t | x t · T x ( t )) . With that and using (5), the update step simplifies to

$$b e l ( x _ { t } ) = \gamma \exp \left ( ( \eta _ { z _ { t } | x _ { t } } + \eta _ { b \overline { e l } _ { t } } ) \cdot T ( x _ { t } ) \right ), \quad ( 9 ) \text{ \ \ } \text{$\ \ } \text{$\ \ } \text{$\ \ } \text{$\ \ } \text{$\ \ } \text{$\ \ } \text{$\ \ } \text{$\ \ } \text{$\ \ } \text{$\ \ } \text{$\ \ } \text{$\ \ } \text{$\ \ } \text{$\ \ } \text{$\ \ } \text{$\ \ } \text{$\ \ } \text{$\ \ } \text{$\ \ } \text{$\ \ } \text{$\ \ } \text{$\ \ } \text{$\ \ } \text{$\ \ } \text{$\ \ } \text{$\ \ } \text{$\ \ } \text{$\ \ } \text{$\ \ } \end{cases}$$

where the normalizing constant γ is computed from the zeroth Fourier coefficient as shown in (4b).

We use bel t : SE N ( ) → R ≥ 0 to represent the probability density function for the posterior at time t . The predicted belief bel t , motion model p u t , and measurement likelihood p z t | x t are similarly denoted. Note that these definitions are functions, not point estimates, allowing us to compute Fourier coefficients for our density functions as η z t | x t ≜ F [ln( p z t | x t )] and η bel t ≜ F [ln( bel t )] . The harmonic exponential Bayes filter is presented in Algorithm 1.

## Algorithm 1 The Harmonic Exponential Bayes Filter

Input : Previous posterior belief bel t -1 , motion likelihood p u t , measurement likelihood p z t | x t and IUR T .

$$<_C++_> // Prediction Step$$

$$\text{our} \quad \text{$ - $. $ - $ - $. $ - $. $ - $. $ - $. $ - $. $ - $. $ - $. $ - $. $ - $. $ - $. $ - $. $ - $. $. $ - $. $ - $. $. $ - $. $. $ - $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $ $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $ - $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $ $. $. $.. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $.. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $.. $. $. $ $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $ $. $. $.$. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $.$. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $ $. $. $. $ $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $.. $. $. $.. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $ $. $. $. $ \. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $.
$. \text{$ \text{$ - $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $........................................................................................................................................................................................................
\text{$ - $. $ - $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $ $. $.. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $.. $. $ $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $ $. $.$. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $ $. $. $ $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $.. $. $.. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $. $.$$

## D. Analysis of the Harmonic Exponential Filter

We now point out several advantages inherent to the HEF. One is efficiency compared to a histogram-based approach. The HEF requires a transition between two different Fourier spaces to compute a convolution, which inherently adds complexity. Nevertheless, it can still be more efficient than a standard convolution implementation in 3D. If we consider a discretized space of n s × n s × n s in 3D, a convolution has complexity O n ( 6 s ) , while a SE (2) Fourier transform using FFT methods is of O n ( 3 s log n 3 s ) . The previous calculation assumes m,n ≈ n L ≈ n s and λ ≈ n λ ≈ n s . Under these assumptions, the tensor product needed for (6) is O n n ( λ 3 L ) and the complexity of (2) is also O n ( 3 s log n 3 s ) .

Table I shows a runtime comparison between a histogrambased approach and the HEF over varying sample grid sizes. As shown, the histogram approach quickly becomes intractable. A workaround we use for the benchmark algorithm in Sec. IV is to convolve with a smaller kernel, eliminating areas with near zero probability, but that is only an option for certain densities, such as low covariance Gaussians.

Another advantage of representing density functions with Fourier coefficients over probability bins in histogram filters is the fidelity with which they can represent signals. When computing the Fourier coefficients for the harmonic exponential distribution, FFT operations behave as a lowpass filter on the signals, leading to information loss in the distribution. However, binning a probability distribution also causes a loss of information. To compare the two, we computed the Kullback-Leibler (KL) divergence D KL ( P || Q ) = ∫ X p x ( ) log ( p x ( ) q x ( ) ) dx , between the two approaches and a mixture of two von Mises distributions on the Lie group S 1 .

|              | Convolution Runtime (s)   | Convolution Runtime (s)   | Convolution Runtime (s)   |
|--------------|---------------------------|---------------------------|---------------------------|
|              | Grid size (10 , 10 , 8)   | Grid size (20 , 20 , 8)   | Grid size (40 , 40 , 8)   |
| Hist (PP)    | 0.2284                    | 15.7969                   | 997.2984                  |
| Hist (Scipy) | 0.0375                    | 2.2434                    | 146.4707                  |
| HEF (Ours)   | 0.0014                    | 0.0083                    | 0.0684                    |

TABLE I: Runtime comparison of a histogram-based approach and the HEF for convolving two functions. Hist (PP) is a direct convolution written in pure python, Hist (Scipy) uses scipy.signal.convolve() , and HEF (ours) implements (6) in Python using the NumPy FFT library.

This analysis is shown in Fig. 3, where the KL divergence is used as a measure of distance between two distributions.

Fig. 3: Analysis of the harmonic exponential distribution and histogram distributions in modeling a mixture of two von Mises distributions on S 1 . (a) depicts how each method reproduces the underlying mixture. (b) shows how D KL evolves for each method as the number of parameters used increases, and shows the effect of the concentration parameter κ on the fidelity of each method, as measured by D KL .

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00907v3-figure-1e40548a-0cfa-415e-ab35-89d9127fe964.jpg)

The von Mises distribution is a probability distribution on the circle S 1 whose density function is given by

$$p ( g | \mu, \kappa ) = \frac { \exp \left ( \kappa \cos ( x - \mu ) \right ) } { 2 \pi I _ { 0 } ( \kappa ) },$$

where I 0 ( ) · is the Bessel function and κ plays a role similar to the inverse of the covariance in a Gaussian distribution. We can also define a mixture of two von Mises distributions, that is, p g ( ) = 1 2 p g µ , κ ( | 1 1 ) + 1 2 p g µ , κ ( | 2 2 ) . Our analysis in Fig. 3 shows that harmonic exponential distributions have a smaller D KL across a range of situations. When both approaches use few parameters, the harmonic exponential distribution does a better job of representing the underlying density. Increasing the number of parameters improves the fidelity of both approaches; however, the harmonic exponential distribution maintains its advantage in terms of D KL .

The parameter κ illustrates the 'low-pass filter' effect of the FFT: as κ increases, the distribution gains high-frequency

Fig. 4: Ability of the HEF and other filtering approaches to represent the banana distribution. (a) shows the prior and motion model, which are modeled with a rectangular and Gaussian distribution, respectively. Both EKF and IEKF use a Gaussian prior. (b) presents the posterior distribution of the four approaches after five prediction steps. The HEF, PF, and IEKF are able to capture the banana distribution.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00907v3-figure-5390db03-376a-4582-9b93-d3e3e1ca3fc3.jpg)

components, leading to a rise in D KL , though, as Fig. 3 shows, the histogram approach performs worse. If the underlying density were discontinuous, like a sum of Haar wavelets, the histogram approach might outperform the HEF in terms of D KL . However, since many robotics problems involve smoothly varying functions, the HEF is better suited to handle continuous densities, as shown in this paper.

Another benefit of the HEF is its expressivity in capturing diverse distributions, such as the banana shape distribution presented by Thurn et al. [1, Chapter 5]. This distribution appears when the uncertainty in the heading of an agent increases due to the accumulation of error during multiple prediction steps. The standard EKF fails to capture this distribution, as it models uncertainty with error ellipses, which is inconsistent with the density of the banana distribution. Long et al. [2] showed that a Gaussian density on the Lie algebra can effectively represent the banana distribution, as shown in Fig. 4. However, a Gaussian density in the Lie algebra can only model unimodal distributions.

We compare the expressivity of our method against four other well-established filtering approaches by taking five prediction steps in a horizontal straight path to the right. Interestingly, the HEF, particle filter (PF) and IEKF are able to capture the banana distribution. Note that all filters start from a rectangular prior, as shown in Fig. 4a, except for the EKF and IEKF due to its Gaussian assumptions. Although histogram filters (HistF) are highly expressive, they fail to capture the banana distribution, as 3D convolutions fail to model the effect of rotations in SE (2) .

## IV. EVALUATION

To evaluate our method, we use one simulated and two real-world datasets where the robot state is represented in

SE (2) . Our assessment focuses on two main metrics. The first metric is the absolute trajectory error (ATE), calculated using two different estimators of the robot position: the mode and mean of the posterior belief. The ATE, measured in meters, is computed along the entire trajectory. However, as this metric does not fully capture the non-Gaussian nature of multimodal estimators, we also calculate the negative log-posterior (NLP) of the ground truth position under the posterior belief. We use the equation NLP = -1 T ∑ t log bel ( x GT t ) where T is the total number of time steps and bel ( x GT t ) is the posterior belief of the true state position, x GT t , computed at time t . The result of the HEF is the posterior belief over the state space which can be used in a state-belief planner [38]. NLP helps to determine how useful our approach would be when used with these systems.

To compute the FFTs, we define a 50 × 50 × 32 grid in SE (2) , with 50 × 50 Cartesian samples over ( x, y ) ∈ [ -0 5 . , 0 5] . 2 and 32 angular samples over θ ∈ [0 , 2 π ) . The experiment maps are scaled to fit this grid. All HEF experiments use λ max = 35 , where the inverse transform (2) is approximated using 111 samples, m = 669 , and n = 32 . We approximate the mean belief relying on the same samples used to calculate the FFTs, that is, µ ≈ ∑ n bel ( x n ) x n . The mode of the HEF is also approximated using the FFT samples: x ∗ = argmax x ∈ x ,x 1 2 ,...,x N F -1 [ η bel ]( x ) . In these assessments, our proposed HEF presented in Sec. III is compared against the following filtering approaches: an EKF implementation that assumes unimodal Gaussian state and error, a histogram filter (HistF), and a particle filter (PF). For the sake of comparison, HEF and HistF use 80 000 , parameters, that is, a 50 × 50 × 32 grid and PF uses 80 000 , particles. To determine the NLP of the PF we approximate the density function using a histogram whose resolution is the same one used for the HEF and HistF. The mode of the PF is estimated by taking the heaviest weighted particle before resampling. For all our experiments we use as a motion model the one of a differential drive robot corrupted with Gaussian noise. To estimate the motion and measurement noise of each filter, we sweep over different noise values and took those that minimize the NLP.

## A. Simulation

The first evaluation is performed in a simulation environment with range-only measurements. This environment includes landmarks arranged in a line. The robot travels in a counterclockwise circle around the landmarks, which alternate in providing a range measurement at each time step, while assuming known correspondence. We chose to simulate range measurements due to their highly non-Gaussian measurement likelihood, which produces a ring around each landmark. All filters start from an equally-weighted bi-modal prior, as shown in Fig. 5. The EKF uses the weighted mean and covariance of this mixture as its prior. We present the results of this experiment in Table II, where metrics are computed over 10 different random seeds. Overall, our method outperforms the baselines, obtaining the lowest ATE and NLP in this multimodal localization experiment.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00907v3-figure-74523e28-bb44-4324-82e5-82607d99e023.jpg)

Fig. 5: Sample results from the simulated environment. (a) shows the posterior belief for the baselines and the HEF. (b) depicts the steps of the HEF: the predicted belief (top) is bel t from (8), the range-only measurement likelihood (middle) is p z t | x t , and the posterior belief (bottom) is bel t from (9).

|            | ATE w/ mode (m)   | ATE w/ mean (m)   | NLP (-)       |
|------------|-------------------|-------------------|---------------|
| EKF        | 0.101 ± 0.03      | 0.101 ± 0.03      | 271.2 ± 23.06 |
| HistF      | 0.151 ± 0.02      | 0.132 ± 0.02      | 17.04 ± 0.57  |
| PF         | 0.131 ± 0.03      | 0.102 ± 0.02      | 5.521 ± 1.15  |
| HEF (ours) | 0.089 ± 0.03      | 0.082 ± 0.02      | 2.242 ± 1.11  |

TABLE II: Quantitative results for the range-only simulator. Reported metrics include the absolute trajectory error (ATE) computed using the mode and the mean of the posterior belief, and the negative log-posterior (NLP). Lower values are better. We report the mean and standard deviation of these metrics over 10 different random seeds.

## B. UWB Range Dataset

For our first real-world evaluation, we used the Clearpath Husky unmanned ground vehicle (UGV) shown in Fig. 6a. The UGV is equipped with an ultra-wideband (UWB) tag that communicates with stationary UWB beacons in the surrounding room. Using a ranging protocol between the UGVmounted UWB tag and the stationary UWB beacons, range measurements between each are computed. The experiment consists of two runs, each with five beacons that were mapped prior to the experiment. For each run, we removed beacons one by one ignoring their range measurements, producing five trials per run. We assume known correspondence.

The motion of the agent was determined by integrating the IMU and wheel encoders of the Husky platform. The motion noise was estimated to have a standard deviation of 3.08 cm and 0.57 degrees per step. Finally, the ground truth was provided by a NaturalPoint OptiTrack motion capture system. We repeat each run with a different number of beacons to demonstrate the robustness of our proposed HEF to handle unimodal and multimodal distributions. The layout of the

(a) UWB.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00907v3-figure-295df206-ef77-40fe-a01e-619a6ef2d109.jpg)

Fig. 6: Setup for the UWB and doors datasets. (a) shows the layout for the two UWB runs with a Clearpath Husky and five ultra-wideband (UWB) beacons. This image also provides a top-down view of the beacon configuration, showing how the measurement likelihood shifts from Gaussian to highly nonGaussian as beacons are removed. (b) presents a top-down view of the laboratory used in the doors dataset.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00907v3-figure-415694c7-8f59-4236-9443-5a25ee70790c.jpg)

beacons used in the experiments is shown in Fig. 6a.

|            | 5 Beacons       | 5 Beacons       | 5 Beacons   | 1 Beacon        | 1 Beacon        | 1 Beacon   |
|------------|-----------------|-----------------|-------------|-----------------|-----------------|------------|
|            | ATE w/ mode (m) | ATE w/ mean (m) | NLP (-)     | ATE w/ mode (m) | ATE w/ mean (m) | NLP (-)    |
| EKF        | 60161           | 60161           | 861e8       | 0.903           | 0.903           | 11.030     |
| HistF      | 0.221           | 0.210           | -2.653      | 1.157           | 1.202           | -0.364     |
| PF         | 0.304           | 0.210           | 3.311       | 3.731           | 2.639           | 9.253      |
| HEF (ours) | 0.199           | 0.194           | -2.984      | 1.616           | 1.076           | -0.533     |

TABLE III: Quantitative results for the UWB dataset, run I. Reported metrics include the absolute trajectory error (ATE) computed using the mode and the mean of the posterior belief, and the negative log-posterior (NLP). Lower values are better. Note that the EKF is highly uncertain about the ground truth position of the agent during the experiments. 2

The results of this evaluation are presented in Table III. In all runs and with different beacon configurations, the HEF is comparable to or better than the baseline filters even in the presence of a highly non-Gaussian measurement likelihood. Although the EKF obtained the best ATE among all filters in the first run with one beacon, its NLP is the largest among all filters, indicating that throughout the experiment the filter is highly uncertain about the ground truth position. As all the results of the second run closely resembled those of the first run, they were not included in the table.

## C. Doors Dataset

In our final real-world experiment, we draw inspiration from the semantic localization problem discussed in Thurn et al. [1, Chapter 1]. The experiment involves localizing an agent navigating anticlockwise in a corridor using doors detected from RGB images as landmarks. All doors share similar semantic features, such as color or handle. Therefore, we do

TABLE IV: Quantitative results for the doors dataset. Reported metrics include the absolute trajectory error (ATE) computed using the mode and the mean of the posterior belief, and the negative log-posterior (NLP). We also report the standard deviation of the translation error.

|            |   ATE mode | w/ (m)   |      |   w/ (m) (-) | ATE mean   |      |     NLP |
|------------|------------|----------|------|--------------|------------|------|---------|
| EKF        |      2.224 | ±        | 1.91 |        2.224 | ±          | 1.91 | 521.479 |
| HistF      |      9.739 | ±        | 7.78 |        8.806 | ±          | 7.12 |  16.719 |
| PF         |      5.66  | ±        | 4.78 |        3.541 | ±          | 2.82 |   4.605 |
| HEF (ours) |      1.525 | ±        | 1.34 |        1.559 | ±          | 1.36 |   0.986 |

not assume known correspondence for which of the 21 doors generated a given measurement in this problem.

We premapped the environment using gmapping [39] and assumed a bearing-only measurement model that enforces a multimodal distribution. Bearing measurements are computed using the centroid of the doors. Bearing measurements are taken only for samples in the free space of the occupancy grid map. As in Sec. IV-B, we integrate the IMU and wheel encoders to estimate the motion of a Jackal platform. The motion noise was estimated to have a standard deviation of 3.51 cm and 2.50 degrees per step. The ground truth position was estimated using the laser-based SLAM implementation of gmapping. The estimated entropy of the posterior distribution over the ground truth's pose using gmapping was 2 232 . nats. The layout of the experiment, as well as an image of a door within the dataset is shown in Fig. 6b.

As there is no known correspondence, we employ a greedy data association strategy, computing the measurement likelihood using only the most likely door for each measurement. We present the results of this experiment in Table IV. Our HEF surpasses by a large margin all baselines in this challenging task, showcasing its ability to handle multimodal distributions and the importance of symmetries for robot localization. A snapshot of the filters is shown in Fig. 7.

During these experiments, our proposed approach took an average of 1.14s to run for each timestep. Our approach was implemented in python and ran on an i7-12700H CPU with 16GB of RAM.

## V. CONCLUSION

In this paper, we have presented a novel approach to nonparametric Bayesian filtering using the harmonic exponential distribution. We showed how our method was capable of capturing multimodal distributions while preserving the structure of non-commutative motion groups by relying on harmonic analysis [36]. Harmonic analysis on groups and homogeneous spaces has attracted new interest due to the growing field of geometric deep learning. As work in this field shows, FFT operations and related convolution computations are fully differentiable [40]. This leads to the possibility of using the harmonic exponential filter as a fully differentiable filter for end-to-end training of measurement and motion models. Another promising avenue for exploration is to use harmonic exponential distributions in a SLAM framework or

Fig. 7: Sample results from the doors dataset. (a) shows the posterior belief for the baselines and the HEF. (b) features the steps of the HEF: the predicted belief (top) is bel t from (8), the bearing-only measurement likelihood (middle) is p z t | x t , and the posterior belief (bottom) is bel t from (9).

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00907v3-figure-f8b3507c-8820-43f8-b9dc-0b3b863e2282.jpg)

a state-belief planner [38]. As Fourie et al. [41] showed, all that is needed to perform belief propagation in a Bayes tree is the product operation between two distributions.

## REFERENCES

- [1] S. Thrun, W. Burgard, and D. Fox, Probabilistic Robotics (Intelligent Robotics and Autonomous Agents) . The MIT Press, 2005.
- [2] A. Long, K. Wolfe, M. Mashner, and G. Chirikjian, 'The banana distribution is Gaussian: A localization study with exponential coordinates,' in Robotics: Science and Systems (RSS) , Sydney, Australia, July 2012.
- [3] D. Dardari, A. Conti, U. Ferner, A. Giorgetti, and M. Z. Win, 'Ranging with ultrawide bandwidth signals in multipath environments,' Proceedings of the IEEE , vol. 97, no. 2, pp. 404-426, 2009.
- [4] Y. Xiang et al. , 'Posecnn: A convolutional neural network for 6d object pose estimation in cluttered scenes,' 2018.
- [5] Y. Wang and G. S. Chirikjian, 'Nonparametric second-order theory of error propagation on motion groups,' The International Journal of Robotics Research , vol. 27, pp. 1258-1273, 2008.
- [6] A. Barrau and S. Bonnabel, 'The invariant extended Kalman filter as a stable observer,' IEEE Trans. on Autom. Control , vol. 62, no. 4, pp. 1797-1812, 2017.
- [7] A. Doucet et al. , 'Sequential Monte Carlo methods in practice,' in Statistics for Engineering and Information Science , 2001.
- [8] T. Li et al. , 'Adapting sample size in particle filters through KLDresampling,' Electronics Letters , vol. 49, pp. 740-742, 2013.
- [9] V. Elvira, L. Martino, D. Luengo, and M. F. Bugallo, 'Improving population Monte Carlo: Alternative weighting and resampling schemes,' Signal Processing , vol. 131, pp. 77-91, 2017.
- [10] D. Scott, Multivariate Density Estimation: Theory, Practice, and Visualization , 03 2015.
- [11] Z. Zhou et al. , 'A discrete quaternion particle filter based on deterministic sampling for IMU attitude estimation,' IEEE Sensors J. , vol. 21, no. 20, pp. 23 266-23 277, 2021.
- [12] K. Li et al. , 'Grid-based quaternion filter for SO(3) estimation,' in European Control Conference (ECC) , 2020, pp. 1738-1744.
- [13] A. Willsky, 'Fourier series and estimation on the circle with applications to synchronous communication-I: Analysis,' IEEE Trans. Inform. Theory , vol. 20, no. 5, pp. 577-583, 1974.
- [14] --, 'Fourier series and estimation on the circle with applications to synchronous communication-II: Implementation,' IEEE Trans. Inform. Theory , vol. 20, no. 5, pp. 584-590, 1974.
- [15] J. T. H. Lo et al. , 'Exponential Fourier densities on SO(3) and optimal estimation and detection for rotational processes,' SIAM Journal on Applied Mathematics , vol. 36, no. 1, pp. 73-82, 1979.
- [16] S. I. Marcus et al. , 'The use of harmonic analysis in suboptimal estimator design,' IEEE Trans. on Autom. Control , vol. 24, no. 4, pp. 530-532, 1979.
- [17] T. Cohen and M. Welling, 'Harmonic exponential families on manifolds,' in 32nd International Conference on Machine Learning , vol. 37. PMLR, 07-09 Jul 2015, pp. 1757-1765.
- [18] H. Sorenson and D. Alspach, 'Recursive Bayesian estimation using Gaussian sums,' Automatica , vol. 7, no. 4, pp. 465-479, 1971.
- [19] I. Gilitschenski, R. Sahoo, W. Schwarting, A. Amini, S. Karaman, and D. Rus, 'Deep orientation uncertainty learning based on a Bingham loss,' in International conference on learning representations , 2019.
- [20] J. R. Schoenberg et al. , 'Localization with multi-modal vision measurements in limited GPS environments using Gaussian sum filters,' in Proc. IEEE Int. Conf. Robot. and Automation , 2009, pp. 1423-1428.
- [21] R. Chen and J. S. Liu, 'Mixture Kalman filters,' Journal of the Royal Statistical Society , vol. 62, no. 3, pp. 493-508, 2000.
- [22] J. Sola, J. Deray, and D. Atchuthan, 'A micro Lie theory for state estimation in robotics,' arXiv preprint arXiv:1812.01537 , 2018.
- [23] G. Kurz et al. , 'Application of discrete recursive Bayesian estimation on intervals and the unit circle to filtering on SE(2),' IEEE Trans. Ind. Inform. , vol. 14, pp. 1197-1206, 2018.
- [24] F. A. Maken, F. Ramos, and L. Ott, 'Stein particle filter for nonlinear, non-Gaussian state estimation,' IEEE Robot. and Automation Lett. , vol. 7, no. 2, pp. 5421-5428, 2022.
- [25] E. B. Sudderth et al. , 'Nonparametric belief propagation,' Communications of the ACM , vol. 53, no. 10, pp. 95-103, 2010.
- [26] L. Song et al. , 'Kernel embeddings of conditional distributions: A unified kernel framework for nonparametric inference in graphical models,' IEEE Signal Processing Mag. , vol. 30, pp. 98-111, 2013.
- [27] Q. Huang et al. , 'Incremental non-Gaussian inference for SLAM using normalizing flows,' IEEE Trans. on Robot. , vol. 39, no. 2, pp. 14581475, 2022.
- [28] D. Brunn, F. Sawo, and U. D. Hanebeck, 'Efficient nonlinear Bayesian estimation based on Fourier densities,' in IEEE Int. Conf. on Multisensor Fusion and Integration for Intell. Syst. , 2006, pp. 317-322.
- [29] --, 'Nonlinear multidimensional Bayesian estimation with Fourier densities,' in 45th IEEE Conference on Decision and Control , 2006, pp. 1303-1308.
- [30] F. Pfaff, G. Kurz, and U. D. Hanebeck, 'Multimodal circular filtering using Fourier series,' in 18th International Conference on Information Fusion . IEEE, 2015, pp. 711-718.
- [31] --, 'Filtering on the unit sphere using spherical harmonics,' in IEEE Int. Conf. on Multisensor Fusion and Integration for Intell. Syst. , 2017, pp. 124-130.
- [32] F. Pfaff, K. Li, and U. D. Hanebeck, 'Fourier filters, grid filters, and the Fourier-interpreted grid filter,' in International Conference on Information Fusion , 2019, pp. 1-8.
- [33] F. Pfaff et al. , 'Multivariate angular filtering using Fourier series,' IEEE Trans. Inform. Theory , vol. 11, pp. 206-226, 12 2016.
- [34] F. Pfaff, G. Kurz, and U. D. Hanebeck, 'Nonlinear prediction for circular filtering using Fourier series,' in 19th International Conference on Information Fusion , 2016, pp. 1821-1828.
- [35] J. J. Fern´ andez-Dur´ an, 'Circular distributions based on nonnegative trigonometric sums,' Biometrics , vol. 60, no. 2, pp. 499-503, 06 2004.
- [36] G. S. Chirikjian and A. B. Kyatkin, Engineering Applications of Noncommutative Harmonic Analysis: With Emphasis on Rotation and Motion Groups , 1st ed. CRC Press, 2001.
- [37] G. Chirikjian and Y. Wang, 'Engineering applications of the motiongroup Fourier transform,' Modern Signal Processing MSRI Publications Volume , vol. 46, 01 2003.
- [38] R. Platt Jr et al. , 'Belief space planning assuming maximum likelihood observations.' in Robotics: Science and Systems (RSS) , vol. 2, 2010.
- [39] G. Grisetti, C. Stachniss, and W. Burgard, 'Improved techniques for grid mapping with Rao-Blackwellized particle filters,' IEEE Trans. on Robot. , vol. 23, no. 1, pp. 34-46, 2007.
- [40] T. S. Cohen, M. Geiger, J. K¨ ohler, and M. Welling, 'Spherical cnns,' in 6th International conference on learning representations , 2018.
- [41] D. Fourie et al. , 'A nonparametric belief solution to the Bayes tree,' in Proc. IEEE/RSJ Int. Conf. Intell. Robots and Syst. , 2016.
Miguel Saavedra-Ruiz Steven A. Parkison Ria Arora James Richard Forbes Liam Paull
The Harmonic Exponential Filter for Nonparametric Estimation on Motion Groups
Bayesian estimation is a vital tool in robotics as it allows systems to
update the robot state belief using incomplete information from noisy sensors.
To render the state estimation problem tractable, many systems assume that the
motion and measurement noise, as well as the state distribution, are unimodal
and Gaussian. However, there are numerous scenarios and systems that do not
comply with these assumptions. Existing nonparametric filters that are used to
model multimodal distributions have drawbacks that limit their ability to
represent a diverse set of distributions. This paper introduces a novel
approach to nonparametric Bayesian filtering on motion groups, designed to
handle multimodal distributions using harmonic exponential distributions. This
approach leverages two key insights of harmonic exponential distributions: a)
the product of two distributions can be expressed as the element-wise addition
of their log-likelihood Fourier coefficients, and b) the convolution of two
distributions can be efficiently computed as the tensor product of their
Fourier coefficients. These observations enable the development of an efficient
and asymptotically exact solution to the Bayes filter up to the band limit of a
Fourier transform. We demonstrate our filter's performance compared with
established nonparametric filtering methods across simulated and real-world
localization tasks.