Max Entropy Moment Kalman Filter for Polynomial Systems with Arbitrary Noise
Authors: Sangli Teng
Harry Zhang
David Jin
Ashkan Jasour
Ram Vasudevan
Maani Ghaffari
Luca Carlone
Venue: NeurIPS 2025 poster ()
Keywords: Bayes Filter
Kalman Filter
Moment Problem
Convex Optimization
Robotics
Primary area: probabilistic_methods
TL;DR: A deterministic asymptotically optimal Bayes Filter for nonlinear non-Gaussian system.

Abstract: Designing optimal Bayes filters for nonlinear non-Gaussian systems is a challenging task. The main difficulties are: 1) representing complex beliefs, 2) handling non-Gaussian noise, and 3) marginalizing past states. To address these challenges, we focus on polynomial systems and propose the Max Entropy Moment Kalman Filter (MEM-KF). To address 1), we represent arbitrary beliefs by a Moment-Constrained Max-Entropy Distribution (MED). The MED can asymptotically approximate almost any distribution given an increasing number of moment constraints. To address 2), we model the noise in the process and observation model as MED. To address 3), we propagate the moments through the process model and recover the distribution as MED, thus avoiding symbolic integration, which is generally intractable. All the steps in MEM-KF, including the extraction of a point estimate, can be solved via convex optimization. We showcase the MEM-KF in challenging robotics tasks, such as localization with unknown data association.

Max Entropy Moment Kalman Filter for Polynomial Systems with Arbitrary Noise
Sangli Teng" sang 1itgiiii @ @ ch.ichhoe edu
Harry Zhang harryzemit edu
David Jin jindavidemit.en edu
Ashkan Jasour jasourO caltech.
Ram Vasudevan" canv@uml ch. edu
Maani Ghaffari" ma maanig jinigjjen 18 umi ch edu
Luca Carlone 1carlone@mit.ede
Abstract
Designing optimal Bayes filters for nonlinear non-Gaussian systems is challeng- ing task. The main difficulties are: 1) representing complex beliefs, 2) handling non-Gaussian noise, and 3) marginalizing past states. To address these challenges, we focus on polynomial systems and propose the Max Entropy Moment Kalman Filter (MEM-KF). To address 1), represent arbitrary beliefs by Moment- Constrained Max-Entrops Distribution (MED). The MED can asymptotically approximate almost any distribution given increasing number of moment straints. To address 2), we model the noise in the process and observation model MED. To address 3), propagate the moments through the process model and recover the distribution as MED, thus avoiding symbolic integration, which is generally intractable. All the steps MEM-KF, including the extraction of point estimate, can be solved via convex optimization. We showcase the MEM-KF KF in challenging robotics tasks, such as localization with unknown data association.
Introduction
The Kalman Filter (KF) is an optimal estimator for linear Gaussian systems in the sense that provably computes minimum-variane estimate of the system's state. From the perspective of the Bayes Filter, the KF recursively computes Gaussian distribution that captures the mean and covariance of the state. However, the optimality of the KF breaks down when dealing with nonlinear systems or non-Gaute Gaussian noise. A plethora of extensions of the KF has been developed to deal with nonlinearity and non- non-Gaussianity. For example, linearization-base Kalman Filter (EKF) [1]. linearize the process and measurement models and propagate the covariance using standard KF. Other techniques, such the Unscented Kalman Filter (UKF) [2], rely deterministic sampling to better capture the nonlinearities of the system but then still fit a Gaussian distribution to the belief. Unfortunately, these extensions do not enjoy the desirable theoretical properties of the KF: they are not guaranteed to produce an optimal estimate and -in the nonlinear or non-Gaussian case- their mean and covariance might be far from describing the actual belief distribution. Consider for example system with non-Gaussian noise sampled from a discrete
methods, such as the Extended
Uniiersity Michigan, Ann Arbor, MI 48109, USA Massachusetts Institute of Technology, Cambridge, MA, 02139, USA
$ Equaal contributions
Team 347T-Robotic Aerial Mobility, Jet Propulsion Lab, Pasadena, CA, 91109, USA
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

Generalized Moment Kalman Filter (Polynomial Non-Gaussian System)
Robot localization with unknown data association
Kalman Filter (Linear Gaussian System)
k++1
Az +
wk
3k+1 f (# Ubs wc) h (yt, 2x) Uk
yk
Ca 24
Figure 1: We generalize the classical Kalman Filter (KF) by proposing the Max Entropy Moment Kalman Filter (MEM-KF), which operates on a polynomial system with arbitrary noise. In the Bayes Filter framework, the KF predicts and updates a Gaussian distribution that can be characterized by its mean and covariance. The MEM-KF predicts and updates max entropy distribution that can be uniquely characterized by its higher-order moments. The proposed algorithm evaluated a robot localization problem with unknown data association. this case, the robot does not know which landmark it observes, which can be modeled by multi-modal distribution.
distribution over {-1, 1;:
clearly, approximating this noise as a zero-mean Gaussian fails to capture the bimodal nature of the distribution. For this reason, existing nonlinear filters, such as the EKF and the UKF, can hardly capture the noise distribution precisely with Gaussian assumptions.
Thus, a natural question is, how can we provide provably optimal state estimation for nonlinear systems corrupted by arbitrary noise? When the noise is non-Gaussian, and the system is nonlinear, the estimation problem becomes less structured for analysis, and traditional methods only focus on the first two moments (i.e., mean and covariance) describe the noise and the belief, which are too crude to capture the true underlying distributions. We argue that by accounting for higher-order moments, we can recover the state distribution more precisely. In this work, we focus on polynomial systems, a broad class of potentially nonlinear systems described by polynomial equations (Fig. 1). These systems can model dynamics on matrix Lie groups [3] and broad set of measurement models [4]. To handle arbitrary noise, we let the user specify an arbitrary number of moments to better describe the process and measurement noise.
Contribution. We propose the Max Entropy Moment Kalman Filter (MEM-KF) with three key contributions. Our first contribution is to formulate the state estimation problem as Kalman-type filter using the Moment-Constraing The MED can asymptotically approximate any moment determinate distribution. Given the moments, the coefficients of the MED can be recovered by convex optimization.
Max-Entropy Distribution (MED) as the belief representation.
Our second contribution is formulate the prediction and update step of the MEM-KF and solve them via convex optimization. In the prediction step, we propagate the moments through the process model and recover the belief MED. In the update step, we formulate the sensor model as MED and show how to obtain the posterior distribution.
Our third contribution is to extract the optimal point estimate via convex optimization. Due to the polynomial nature of the MED, we apply semidefinite moment relaxation techniques to extract the optimal point estimate given the belief representation. The global optimality of the solution (resp. non-trivial suboptimality bound) can be certified by using a rank condition (resp. the relaxation gap) of the moment relaxation. We showcase the resulting MEM-KF in challenging roboties tasks, such as localization with unknown data association.
Compared to the traditional marginalization step in the Bayes Filter that requires symbolic integrations and is generally intractable, the proposed method only requires numerical integrations to evaluate the gradients for the convex optimization. Although the numerical integration is still relatively expensive, MEM-KF provides a general solution for filtering polynomial systems with arbitrary noise.
Related Work
Bayes Filter. The Bayes filter recursively estimates the posterior probability of a system's 's state by predicting the prior from the process model and updating it with new measurements using Bayes' rule [5]. For a nonlinear process model, even a simple distribution can be skewed and transformed to a multi-model distribution that is hard to approximate. As a solution, the prior can be approximated by parametric functions, such as a Gaussian mixture [6. 7]. Despite that, the belief can be asymptotically approximated with larger Gaussian sums, the prediction still assumes linear process models [6] or relies on linearization [7]. Other approximation of belief includes point mass [8, 9] and piece-wise constant [10, 11] that are nonsmooth. Contrary to the deterministic method that approximates the

---

prior by parametric functions, sampling- campling-based methods such as partiele filters [12] have been applied to approximate the predicted distribution. Though the sampling-based methods are more versatile for nonlinear non- non-Gaussian systems, these methods generally lack formal guarantees. We recommend the readers to [5] for more comprehensive review.
Kalman Filtering. When the process and measurement models are linear, and the noises are Gaussian, the Kalman Filter [13] becomes provably optimal and returns the minimum-varianee estimate among all estimators. To deal with nonlinear systems, the Extended Kalman Filter (EKF) linearizes the process and measurement models around the estimated states [14] and propagates the covariance through the linearized system. However, the performance of EKF deteriorates as the linearization points deviate from the true states, as the process model becomes more nonlinear and the noise increases [1, 15]. To tackle this problem, symmetry presseripg lincarization has been applied in homogeneous space [16-18] that can avoid the inconsistency caused by imperfect state estimations. As an alternative to linearizations, the Unscented Kalman Filter ( (UKF) [2] applies a deterministic sampling approach to handle nonlinearities. However, for non- non-Gaussian noise, UKF lacks systematic way to handle arbitrary noise. We note that both EKF and UKF adopt the Kalman gain and only use the mean and covariance of the noise distribution. Shimizu et al. [19] and Jasour et al. [20] are closer to our work, where moment propagation is derived for filtering problems by assuming knowledge of the characteristic function of the Gaussian noise. In our methods, we only utilize moment information instead of the entire distribution, e.g., the characteristic function, which is inaccessible in typical applications. Though [21] utilized higher-order moments to derive a minimal variance filter, the prediction applies approximations that is not exact.
Semidefinite Relaxation Semidefinite relaxations have been critical to designing certifiable algo- rithms for robot state estimation. These relaxations transform polynomial optimization problems into semidefinite programs (SDPs) that are convex, using Lasserre's hierarchy of moment relaxations [22]. Subsequently, SDP relaxation-basee certifiable algorithms have been extended to various geometric perception problems, including pose graph optimization [23, 24], rotation averaging [25, 26), trian- gulation [27, 28), 3D registration [29-32], absolute pose estimation [33], and category- perception [34, 35]. Compared to belief, e.g., [6-11] that are hard to extract or certify the globally optimal point estimation, the optimal point estimation of MED can be solved by SDPs.
object
3 Preliminaries
We review key concepts in polynomial systems. Let R[x] be the ring of polynomials with real coefficients, where x := [ .2. Fn. Nº I, 0; S S (i.e., is a vector of integers that sums up to no more than r). A monomial of degree up to r is denoted as x := γी xg aa a E Nr. For a polynomial p(æ) == E, x? its degree deg p is defined as the largest 1//// with + O.
Given an integer r, we define the set N == {a e
Let s(r, n) = | ni + )) We then define the basis function (x) : R=" + Ro using all
the entries of the canonical basis (ar) 12, Liy a. 21322 2 x x
Given a probability distribution p(x:) over support K,, we can compute the moment matrix of prob- ability p(x) as M, (..) = JK. (x)x))((zz dx) a E N where we have S. := [xc x:" ((xx()) da to denote the moment of æ². The inverse of this process is the moment problem that tries to identify the distribution given the moment sequence. A special case is the moment-determinntt distribution that describes wide range of well-behaved probabilities that are neither long-tailed nor have higher-order moments growing too fast:
Definition 1 (Moment-determing distribution). A distribution p(x) is said to be moment- determinate if it can be uniquely determined by all the moments with order from 0 to oo.5
We introduce the following equality-constraiin Polynomial Optimization Problem (POP) that will be applied to extract the optimal point estimation given an MED:
WWit tinite order there can be infinite number of distributions with the specified moments.
3

---

Definition 2 (Equality-onnstriined constrained POP). An gquality-constranne polynomial optimization problem (POP) an optimization problem with polynomial cost and constraints g(x))(())))))))))))))))) E R[x]
= inf c(x)
s.t. g,(x;) 0,
Vj E (1....,m)}
(POP)
The moment relaxation that solve (POP) via SDPs is presented in Appendix A. Some introductory examples of the moment matrix and polynomial systems are provided in Appendix and C. In the perception problems, such relaxed SDPs usually retrieve globally optimal solutions to the original (POP), meaning that the relaxation is often exact [28, 36, 37].
Problem Formulation
We consider the following polynomial dynamical system:
Æk++1 f uk; wx) h 2x) = 22k
E K, A Ak,
where ag is the deterministic state to estimate, yı: are the measurements, and ug is the control input, all defined discrete time We assume that the state 2& restricted to the domain K,, e.g., the set of 2D poses. Both the process model f and the observation model h are vector-valued real polynomials. We take the standard assumption that the process noise wy and measurement noise U/ are identically and independently distributed across time steps. For polynomial systems, we make the following assumption:
Assumption 1 (K)) The domain K. is either the Euclidean space K. = R (i.e., the state is uncon- strained), or can be described by polynomial equality constraints = {x g(x) = E R")
Assumption 1 relatively mild and captures broad set of roboties problems where the variables belong to semi-algebraic sets (e.g., rotations, poses); see, for instance [38, 39]. Given the system in eq. (I), formally define the filtering problem that recursively estimates the belief p(æ)):
Problem 1 (Recursive State Estimation). Assume the state time k is solely dependent on the previous state Tk. action Ug. and measurement yk, we have the process and observation model as:
p( 1kk, Ug, '.- ,20) p( [xx ux).
p(yx Lkk yk 30) p(yx (xx). (2)
Our goal is to recursively estimate the belief p(æx) == p(æx uk. uo yk. ...,, yo) for the system in (1) dependent on all the information from process and observation models until time k.
In the Bayes Filter, we apply the prediction and update steps fuse the dynamics and measurements:
p (*x+1) & [ p(a mxx us (p(x))d
(Prediction)
p(ax)) x p(yx (x (#+).
(Update)
As the system (1) nonlinear and corrupted by arbitrary noise, solving Problem via Bayes filter is challenging in representing the belief p(x) and the noise. When non-trivial noise belief is involved, marginalizing the distribution in the prediction steps is also generally intractable. In the following sections, we leverage the polynomial structure of the system to address these challenges.
5
Moment- Connet--onntrraiee Max-Entropy Distribution
In this section, we introduce the MED and show how to recover the distribution given the moment constraints using convex optimization.
5.1 Asymptotic Property of MED
Though a moment-determinate distribution can be uniquely determined by its moments, for computational purposes, we expect to describe distributions with finite number of moments. Given a finite number of moments, the underlying distribution is usually not unique Therefore, to guarantee uniqueness, define the notion of Moment- Conteennnon Max- Ent-ooy Distribution.

---

which is the distribution with maximum entropy among all the distributions that match the given moments. We denote the set of distributions that match set of given moments of order up to r as P. = {p(x) Sxx p(ar)dæ= = T.. Va E ) and its limit as Px == lim,- Pr.
Thus, the MED can be obtained by maximizing the entropy functional over P,:
Se p(x) In ((x)dd s.t. p(z) Pr. p())
(x) arg max
(Pr)
We denote the solution its limiting problem (Poo) with feasible set Poo as p (2). For moment- determinate distribution in Definition 1, Poo (x) can be uniquely determined by Poo The asymptotic property of (x) can be summarized by:
Theorem (Asymptotic approximation of Poo is defined by the moments of a moment determinate distribution. the solution P (æ) is unique and pł (x) converges to P (x) norm," i.e., lim (1);(x) P (x) O.
Determinate Distribution [40]). Suppose that
Proof. The complete proof shown [40].
We consider MED as the belief representation for our filtering approach. The resulting approach avoids the integral in the (Prediction) step -which generally intractable- by directly operating over the moments, and then recovers the corresponding distribution via (P,). The approach guarantees convergence to the true belief the limit, i.e. - oc, thanks to Theorem I.
Remark 1.
Remark 2. For p(x) with compact support. p(æ) is guaranteed be moment determinate, i.e.. P (x) be unique; see Proposition 12.17 of [41]. For an engineering problem where states are not possible to be infinite, additional assumptions on the compactness of support are reasonable.
5.2 Recover MED via Convex Optimization
Now we briefly review the structure of the solution of (P,) and how we can retrieve it via convex optimization. The full derivation by [42] in a normed linear space is summarized Appendix D-E.
By the first-order optimality condition, we can show that the solution to (P,) takes the following form with the coefficients Aς be determined:
p(x)= exp Ea Aax=")
Via the Legendre-Fenchel dual of (P+), the optimal coefficients can be obtained by minimizing the thermal dynamics potential A(A) of p(x) [42]:
min
A A = Se exp (- Ea ) + Eo Z0
(4)
We observe the gradients,
0A(A) 8Ag
exp
exp (- E. Ja 3z N:
(5)
equal zero when p(x) satisfies all the moment constraints. Thus, we minimize A(X) search p(æ) that meets the first-order optimality constraints. We can further verify the convexity of ((A) by verifying that the following Hessian positive definite:
H()) = aD(A) aAB2y
[ exp (-
Σ Ja T E CANERNNA NI
We find that H()) is exactly the moment matrix generated by p(æ), which is positive semidefinite as p(x) is nonnegative. We note that numerical integration is required to evaluate the Hessian and gradients. For an n-dimensional system with moments up to order r are matched, computing the exact gradients and Hessians requires evaluating |N]| and |Ng, terms of moments, respectively.
For nonempty defined by g(r) in Assumption 1, additional linear constraints is needed (4) to consider the ideal generated by g(æ), i.e., Z[g(=)) = {E; detailed derivation is summarized in Appendix F.
R(x]). The
6 [400 considers the nonorm [|p(=)]]1 So [p[x]]dx.
5

---

Main MEM-KF Algorithm
In this section, we introduce the Max Entropy Moment Kalman Filter (MEM-KF) algorithm. In the prediction steps, instead of directly marginalizing the distribution, we predict the moments via the dynamics first, and then recover the MED matching those moments. In the update steps, we formulate the observation model as MED and update the parameters via Maximum A Posterior estimation.
6.1 Extended Dynamical System
To generate higher-order moments of the states to better approximate the belief, we need an extended version of system (1) that exposes these moments in the state-space model. To achieve this goal, we apply (-)- to both sides of (I) to generate monomials of degree up to for 2R, Wh and ux:
o, (#k+1) (f (#k, Uk wx))
{ Φ (x(+1) A(wg, un ()))((()))) C( = (1(4))
(6)
Φ, (h(yk, xx)) (Uk)
Intuitively, (6) combines the entries of the original process and measurement models into higher-order polynomials by taking powers and products of the entries of ʃ and h. Since each polynomial of línear combination of the canonical basis Φb, (x:), ,e can rewrite both the extended process or the observation models as linear functions of Φb, (æx). Note that as the Φs (2 on the RHS of the dynamics may have different degrees as LHS due the nonlinearity, thus we have the order possibly # s, # We also note that the first term of (-,(-) (-) is always so we have the dynamics as homogeneous linear equalities. An example of the extended dynamical system is presented in Appendix G.
degree up can be written
6.2 Initialization
We denote the belief parameters after the prediction step with superscript (-) For recursive estimator, we set the initial state as a MED (o) = exp (- Ea x0, with coefficients Ag- One typical choice is the Gaussian distribution, max-entropy distribution with constraints on the mean and covariance.
6.3 Prediction Steps
It is possible to formulate the prediction step marginalizing the joint density between zµ and #4+1
p (= 1) xx fK p(#+++ Øx, un )p(x Au dakx.
(7)
However, the integration on the RHS is generally intractable for arbitrary distribution. Even p(0 Dx. Ug) and p( (<x))>)))))))))))))))))>>))))( MED, p (xx++1 is not guaranteed to be MED, that is, not claim MED is a valid conjugate prior in general.
To mitigate this issue, we instead use the extended dynamics to predict the moments for p (œx++1 and then recover the distribution. Consider the extended dynamics and take the mean value on both sides:
Φ (0x+1) = A Wk: ux (x))..())) + (x++1) = B[A(ww (wy uy (x)
(8)
Under the assumption that the noise distribution is independent of the state distribution, the expectation of the cross term of and w can be separated, and the above equation can be simplified to:
Ta,k+1 A( ug 3
with the moment sequence specified by a E NP; Ng. Given the moments sequence 1o can apply (4) recover the distribution represented by A i.c., p( (
Remark According to Theorem 1, (œx+1) converges the true distribution norm as r op.
which, in this limiting case, marginalizes the distribution (7) without the symbolic integration.
6.4 Update Steps
In the update steps, apply Maximum A Posterior estimation to update the belief. Consider the
extended measurement model (()) (x) Φ,,(v) Given the moments of Φ, (v), i.c., Ca we model the noise distribution MED parameterized by μl:
p(ulu)) exp (- Ea 2)
(10)

---

We note that the noise relates to the observation model via = h(y, Thus, we have the likelihood function by substituting the with hyy
ply æ. v) exp (- Ea Mahhy x)
exp (- Ea vax"),
(11)
where v is the coefficients corresponding the canonical basis (x).. With batch of m measure-
ments, and the prior distribution p( A ) the posterior distribution æ can be obtained by:
p(=|2) x p(= I[ ply x, p1) exp (- E., >,x"),
(12)
with À being the summation of the coefficients À = = + Σ. Vi- To predict the belief of the next step, normalize the likelihood function via numerical integration to compute the moments.
6.5 Extraction of Optimal Point Estimation
Algorithm I Max Entropy Moment Kalman Filter
Given the belief at each step of MEM-KF, we extract the optimal point estimation with the highest probability density by:
Require: Initial distribution parameterized by X dy- namics model #x+1 f(xx, wx), observation model h(y, ax) =
maX p(a X) A minlex -In p(x)))
// Obtain the observation medel μa (4)
(13) Given Assumption and (3), (13) is (POP) with the polynomial cost function being the neg- ative log-likelihood. Thus, we instead solve the semidefinite moment relaxation of (POP) (see details in Appendix A) to extract the optimal point estimate that has the highest density value. Under the condition that the resulting SDP has rank-1 solution, the optimal point estimate can be extracted and certified by the rank condition. We summarize MEM-KF in Algorithm 1.
0
time k= I.....,,NN do do
// Compute the moments at h k-1
- J p( da
// Propagation of moments
Fa, + A(w.. 4
// Reconstruct the distribution
A- Ta.k
// Update with measurements
(1)
ye. us A < A +
// Normalize the distribution At + Ak,no - In Se p(æ|Ax)ds
Remark (KF subset of MEM-KF). The KF is a special case of MEM-KF with linear pro- cess and measurement models, and the moments up to the second order are matched.
// Optimal point estimation
ai. X; (SDP) minzex Inp(a/Ax)
end for
7 Numerical Experiments
7.1 Asymptotic Property of MED
We first showcase the asymptotic property of MED in approximating complicated distributions. We apply a line search method to implement (4) with the initial step size determined by exact first-order gradients (5) and the Hessian approximated by the BFGS algorithm. Although Newton's 's method is introduced in [42], we find obtaining the exact Hessian matrix extremely expensive for high-order systems. We adopt SPOTLESS [43] to process the polynomials and an adaptive Gaussian quadrature method [44] for numerical integrations. We apply the MOSEK [45] to solve the (SDP) to extract the point estimations. More results on nonempty K. like (E(() can be found the Appendix. H.2.
In the first case, we consider distribution on IR generated by the following sampling strategy:
sin w +1 v cos w + 62
(14)
where w ~ (10 x) is a uniform distribution between 0 and and é1, 62 ~ (0, 0.01) is Gaussian noise. We draw N = 1e5 samples from the distribution to compute mean values, an unbiased estimation of the underlying moments, i.e., Oa N Σ Us E N..
Then recover the distribution 2 as p(=|A). The samples and MED with moment up order are shown in Figure 2. We observe that as more moments are matched, the distribution asymptotically approaches the sample distribution. The result is expected according Theorem 1. When r the MED is the Gaussian distribution. Due to the nonlinear transformation, the samples are concentrated on the curve with higher curvatures, which starts emerging p(æ) A) when 2 4. We can further verify that the result of 12 is sufficiently good by comparing it with the heat map of the samples.
7

---

P-222 r-d F-0 T-8
P to 12
Sangles
Hoat Mag Sample
DO
The MED with moment matched up order The heat map and samples are shown on the RHS of the plot. MED asymptotiallyy approaches the sample distribution as the order increases.
Figure 2:
Table 1: Covariance of rank-1 point estimation with N accumulated samples using noise model (14) and (15) with moments up to order r matched. The normalized log-likelihood of the belief, i.c., In p(æ), ), is shown in (-) for small N. We note that increasing consistently results in higher density.
50
(BLUE)
1.58
1.54)
0,806 (-0.844)
0.539
(109(-00.000 0.439)
0.326
0.163
0.0784
0.0319
0.0301
(-0.840)
1.58 (0.034)
(0.665)
0.0303
(14)
4.48
4.32 (0.8922
0.161)
1,80 (0.885)
1.83 (2.03)
0.696 (1.65)
0.661 (2.81)
0.131
0.150
7.09ce3
0.0210
9.30c
3.54c-
2.72c
8.04
4.27e
1.34e
4.04 (1.26)
1.63 (3.04)
0.628 (3.96)
0.0586
1 1.00c-3
3.08c-4
1.07
4. 4.60e-5
(1.72)
1.72 (4.25)
(5.92)
3.78c
(15)
2 2(BLUE)
2.07 (-1.88)
4.47 (0.073)
1,01 11.18)
2.05 (1.42)
0.706 (-U.779)
(.14(2.22)
0.420
0.272
0.0128
0.208
4.39c
0.0999
1.90c
0.04U7
9.38e
0.0199
7.2 Update Steps
Now we verify the update step of MEM-KF. We consider a basic linear observation model y = R², where the moments of known in advance. Then, we consider the non-Gaussian noise generated by the following nonlinear sampling strategy:
24 1 11
41 92 Bernonlli (1.0) (0.)))),.....2(((
(15)
242 - 1 + 62
We compare the proposed method with the Best Linear Unbiased Estimator (BLUE)[46), which is the minimum variance estimator using only mean and covariance. We compare the asymptotic property of the update steps with accumulated measurements. We accumulate N measurements to obtain the
negative likeiiiieodd,, i.e., In p( yı y2. ... yN) = Scalar + Ex= E. (yy x)
The normalized belief distribution is illustrated in Figure 3. After receiving the first measurement, notice that the belief computed by r = retains the multi-modality caused by the noise distribution and exhibits four different peaks. After more measurements the true state with dramatically smaller covariance. On the other hand, the BLUE (r 2) provides a unimodal (Gaussian) approximation of the belief and even with more data, still exhibits relatively large covariance around the true state. The interested reader can also find results with noise (14) are shown in Figure 6 in the Appendix. H.1 which exhibit similar patterns.
incorporated, the belief converges to
We further sample 1000 trials in each noise with different N to compute the empirical variance of the point estimation extracted by (SDP). The results are presented in Table 1. We observe that with sufficiently many measurements accumulated (N 2 5), increasing consistently lowers the covariance and outperforms BLUE dramatically. With small N, the multi-modality of the belief introduces additional bias when computing the point estimation that makes the covariance higher than BLUE. However, note that the point estimation is not full belief representation and can represent the multi-modality. The log-likelihood of the normalized density for the cases with N < 3 computed to show that increasing can provide point estimations with higher likelihood.
7.3 Localization with Unknown Data Association on SE(2)
Now we consider localization problem where robot with pose:
Z - RR E SE((2), including rotation R -s and position =
Py
moves in the plane while collecting aange-andd-berrin measurements. 7 The problem is modeled by the following dynamical system:
Z W.U.W E SE(2) yv = R. [L,, + - Px) R
(16)
"Note that the unknown rotation has to satisfy the equality constraint g(x) = 6² + 5² - - 0. due to the structure of the Special Orthogonal group SO(2).
8

---

N
--3 --4 N
N
Figure 3: Applying the update step of MEM-KF to the linear system corrupted by noise (15). With larger r, the posterior distribution converges faster to the truth states. The case for noise (14) is presented in the Appendix.
IaEKF
UKP
UXPM
MEM-XY
RXP
1sRKF
UKK
IXPM
MFM-XF
Figure 4: Localization with unknown data association. In the first row, the landmark association follows the distribution in (17). In the second row, the measurement corner. The performance baselines degrades due to incorrect data associations, while MEM-KF can consistently track the ground truth.
from the landmarks the lower
where the pose Z, input U, and noise W all belong SE(2). U models the wheel odometry of differential-drive wheel robot, and W can be generated by the associated noise distribution in the Lie algebra [47]. Due the SO(2) group constraints on R, the system is nonlinear, We denote with L, E R² the (known) position of landmark i I and with Lix the position of the landmark observed time while is the measurement noise.
In our case study, we consider the robot receiving observations from four landmarks with unknown data associations. Due the unknown data associations, for given yk. do not know which landmark is actually measured, so the index ik is unknown. The work [48] decides the data association between observations and landmarks by taking the ik variables. For filtering problem, [48] is equivalent to choosing the most likely association and leads to an intractable estimation problem. To mitigate this issue, assume each measurement P1 _... = PN = N. Thus, such "sensor model" with the max entropy association strategy can be formulated as
of the landmarks with equal probability
Rey + = = == Li +
(17)
where = E I and u E (O, The RHS can modeled as single noise term Uz that can be approximated by MED following the same procedure for (15). We consider 4 in our case, which sufficiently good to represent the four peaks caused by the landmarks. For r = representing the four landmarks Gaussian distribution with extremely large covariance leads to little innovation in the update steps.
We present the result of the proposed method and compare it with the baselines, including the EKF, UKF, and their variants on SE(2) manifolds, i.e., the InEK [16] and UKFM [49]. We consider two scenarios, such that 1) the underlying association distribution is identical to the max entropy association, i.e., the measurement comes from four landmarks with the equal probability, and 2) the robot receives more measurements from the landmark on the lower-left corner. The results of the two

---

N-14 X-Jrdd
N
NEMKP
X
2
Figure 5: Comparison of particle filter with particles with MEM- KF. In the first row, the landmark association follows the same distribution in (17). In the second row, the measurement comes more from the landmarks at the lower left corner. We that the particle filter can better track the trajectories with more particles. We find that hoth the particle filter with sufficient number of particles and MEM-KF estimate the trajectories with similar quality. The statistics are presented in Figure 8.
cases are shown in Figure 4. We can see that MEM-KF consistently outperforms the baselines in both scenarios. The baselines 2) have clear bias in the trajectory, due the fact that they greedily choose the most likely landmark the originator of a measurement.
The MEM-KF is also compared with the particle filter in Figure 5. When sufficiently many particles are added, both methods can successfully estimate the trajectories with similar estimation error. However, we note that MEM-KE KF deterministic and provides a smooth belief representation that is analytical and convenient for the extraction of the globally optimal point estimation via convex optimization. Such property not possible for the partiele filter that is based on random sampling. More statistics about the comparison are presented in the Appendix. H.3.
8 Conclusions
In this work, we presented the Max Entropy Moment Kalman Filter (MEM-KF) for optimal state esti- mation of polynomial systems corrupted by arbitrary noise. By the Moment-Constraang Constrained Max- Suxt-oy Distribution (MED), we can asymptotically approximate almost any belief and noise distributions given an increasing number of moments the underlying distributions. We then leverage the system dynamics to predict future moments overcome the difficulties in the marginalization steps in the prediction steps. The update step is formulated as a Maximum A Posteriori estimation problem with measurements modeled MED. All the steps in MEM-KF, as well the extraction of the optimal point estimation, can be exactly solved by convex optimization. We showcase the MEM-KF in challenging roboties tasks, such as localization with unknown data association.
Limitations and Future Work The proposed method relies on numerical integration, which is computationally intensive for real-time deployment. The clussian-quadratte applied to evaluate the gradients (5) in (4) still scales exponentially w.r.t the dimension of the systems. Given this reason, the average time consumption for the localization case in the update is 6.87 seconds and 34.51 seconds in recovering MED from moments. Future work will consider GPU-based numerical integration that scales better for high-dimensiooal systems and real-time deployment. Possible implementation can be [50], which includes a GPU version of the deterministic integration methods in [51]. To further improve the scalability of MEM-KF, we can leverage the sparsity of polynomials by considering only a partial amount of moments [52].
Societal Impact This paper provides optimal filtering framework that can asymptotically approxi- mate any moment-determine This paper pushes the boundary of reliable and certifiable state estimations, especially for robotics.
belief for non- aun-Gussiian polynomial systems with smooth functions.
A Acknowledgment
This work was partially funded by the ONR RAPID program and Carlone's NSF CAREER. Maani Ghaffari and Ram Vasudevan are sponsored by AFOSR grant FA9550-23-1-04000 (MURI).
10

---

References
[1] Yongkyu Song and Jessy W Grizzle. The extended Kalman filter as local asymptotic observer for nonlinear discrete-time systems. In Proceedings of the American Control Conference, pages 3365-3369. IEEE, 1992. 1,3
[2] Eric A Wan and Rudolph Van Der Merwe. The unscented Kalman filter for nonlinear estimation. In Proceedings of the IEEE Adaptive Systems for Signal Processing, Communications, and Control Symposium, pages 153-158. leee, 2000. 1,3
[3] Sangli Teng. Optimization-basee Robot Control and State Estimation on Matrix Lie Groups. PhD thesis, 2025. 2
[4] Heng Yang and Luca Carlone. Certifiably optimal outlier-robust geometric perception: Semidef- inite relaxations and scalable global optimization. IEEE Transactions Pattern Analysis and Machine Intelligence, 2022. 2
[5] Zhe Chen et al. Bayesian filtering: From kalman filters to partiele filters, and beyond. Statistics, 182(1):1-66 69, 2003. 2, 3
[6] Harold W Sorenson and Daniel L Alspach. Recursive bayesian estimation using gaussian sums. Automatica, 7(4):465-479, 1971. 2,
[7] Daniel Alspach and Harold Sorenson. Nonlinear bayesian estimation using gaussian sum approximations. JEEE transactions on automatic control, 17(4):439-448, 2003.
[8] Richard S Bucy. Bayes theorem and digital realizations for non-linear filters. Journal of the Astronautical Sciences, 17:80, 1969.
[9] Richard S Bucy and Kenneth D) Senne. Digital synthesis of non-linear filters. Automatica, 7(3): 287-298, 1971.
[10] Genshiro Kitagawa. Non-gaussian state-space modeling of nonstationary time series. Journal of the American statistical association, 82(400): 1987.
[11] Stuart C Kramer and Harold W Sorenson. Recursive bayesian estimation using piece-wise constant approximations. Automatica, 24(6))789-801 1988. 2,
[12] Petar M Djuric, Jayesh H Kotecha, Jianqui Zhang. Yufei Huang, Tadesse Ghirmai, Mónica Bugallo, and Joaquin Miguez. Particle filtering. IEEE signal processing magazine, 20(5):19-38, 2003. 3
[13] Rudolph F. Kalman. A New Approach to Linear Filtering and Prediction Problems. Journal of Basic Engineering, 82(1))35545, 03 1960.
[14] Sebastian Thrun. Probabilistic robotics. Communications of the ACM, 45(3):522255 2002.
[15] Guoquan Huang, Anastasios I Mourikis, and Stergios I Rooumeliotis. Observability-basee rules for designing consistent EKF SLAMM estimators. International Journal of Robotics Research, 29(5):002-528, 2010. 3
[16] Axel Barrau and Silvère Bonnabel. The invariant extended Kalman filter stable observer. IEEE Transactions on Automatic Control, 62(4): 2016. 3,
[17] Axel Barrau and Silvere Bonnabel. Invariant Kalman filtering. Annual Review of Control, Robotics, and Autonomous Systems, 1:237-257, 2018.
[18] Pieter van Goor, Tarek Hamel, and Robert Mahony. Equivariant filter (eqf). IEEE Transactions on Automatic Control, 2022. 3
[19] Yutaka Shimizu, Ashkan Jasour, Maani Ghaffari, and Shinpei Kato. Moment- based Kalman filter: Nonlinear Kalman filtering with exact moment propagation. arXiv preprint 20X19.230..09130, 2023.

---

[20] Ashkan Jasour, Allen Wang, and Brian C Williams. Moment Moment-based exact uncertainty propaga- tion through nonlinear stochastic autonomous systems. arXiv preprint arXIv:2101.12990, 2021.
3
[21] Sangli Teng, Harry Zhang, David Jin, Ashkan Jasour, Maani Ghaffari, and Luca Carlone. Gmkf: Generalized moment kalman filter for polynomial systems with arbitrary noise. arXiv preprint arXiv:2403.04712 2024. 3
[22] Jean B Lasserre. Global optimization with polynomials and the problem of moments. SIAM Journal on optimization, 1133)796-817, 2001. 3, 16
[23] Luca Carlone, Giuseppe C Calafiore, Carlo Tommolillo, and Frank Dellaert. Planar pose graph optimization: Duality, optimal solutions, and verification. JEEE Transactions on Robotics, 32 (3):545-565, 2016. 3
[24] David M Rosen, Luca Carlone, Afonso S Bandeira, and John J Leonard. SE-Sync: A certifiably correct algorithm for synchronization over the special euclidean group. International Journal of Robotics Research, 38(2-3):95-125, 20109.
[25] Johan Fredriksson and Carl Olsson. Simultaneous multiple rotation averaging using Lagrangian duality. In Asian Conference on Computer Vision, pages 245-253 Springer, 2012. 3
[26] Anders Eriksson, Carl Olsson, Fredrik Kahl, and Tat-Jun Chin. Rotation averaging with the chordal distance: Global minimizers and strong duality. IEEE Transactions on Pattern Analysis and Machine Intelligence, 43(1):256-268, 2019.
[27] Chris Aholt, Sameer Agarwal, and Rekha Thomas. qeqp approach to triangulation. In Proceedings of the European Conference Computer Vision, pages 654-667. Springer, 2012.
3
[28] Diego Cifuentes. A convex relaxation to compute the nearest structured rank deficient matrix. SIAM Journal on Matrix Analysis and Applications, 42(2):708-772, 2021, 3,
[29] Jesus Briales and Javier Gonzalez-limenez. Convex global 3D registration with Lagrangian duality. Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 4960-4969, 2017.
[30] Kunal N Chaudhury. Yuehaw Khoo, and Amit Singer. Global registration of multiple point clouds using semidefinite programming. SIAM Journal on Optimization, 25(1):46--501, 2015.
[31] Haggai Maron, Nadav Dym, Itay Kezurer, Shahar Kovalsky, and Yaron Lipman. Point reg- istration via efficient convex relaxation. ACM Transactions on Graphics (TOG), 35(4):1-12, 2016.
[32] José Pedro Iglesias, Carl Olsson, and Fredrik Kahl. Global optimality for point set registration using semidefinite programming. In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 8287-8295, 2020. 3
[33] Sérgio Agostinho, João Gomes, and Alessio Del Bue. Cvxpnpl: A unified convex solution to the absolute pose estimation problem from point and line correspondences. Journal of Mathematical Imaging and Vision, 65(3)492--512 2023.
[34] Jingnan Shi, Heng Yang, and Luca Carlone. Optimal pose and shape estimation for category- level 3D object perception. In Proceedings of the Robotics: Science and Systems Conference, 2021.
[35] Heng Yang and Luca Carlone. In perfect shape: Certifiably optimal 3D shape reconstruction from 2D landmarks. In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 621-6300 2020. 3
[36] Diego Cifuentes, Sameer Agarwal, Pablo A Parrilo, and Rekha R Thomas. On the local stability of semidefinite relaxations. Mathematical Programming, pages 1-35, 2020. 4
[37] Diego Cifuentes, Corey Harris, and Bernd Sturmfels. The geometry of sdp-exactness in quadratic optimization. Mathematical programming, 182(1-2):399442 2020. 4
12

---

[38] Luca Carlone. Estimation contracts for outlier-robust geometric perception. Foundations and
Trends (FnT) in Robotics, 2023.
[39] Sangli Teng, Ashkan Jasour, Ram Vasudevan, and Maani Ghaffari Jadidi. Convex Geometric Motion Planning on Lie Groups via Moment Relaxation. In Proceedings of the Robotics: Science and Systems Conference, Daegu, Republic of Korea, July 2023. doi: 10.15607/RSS. 2023.XIX.058. 4
[40] Jonathan M Borwein and Adrian Stephen Lewis. Convergence of best entropy estimates. SIAM Journal on Optimization, 1(2):1912055, 1991. 5, 14
[41] Konrad Schmüdgen et al. The moment problem, volume 9. Springer, 2017. 5
[42] Lawrence R Mead and Nikos Papanicolaou. Maximum entropy in the problem of moments.
Journal of Mathematical Physics, 25(8):2404-24111 1984. 5, 7, 18
[43] Mark M. Tobenkin, Frank Permenter, and Alexandre Megretski. SPOTLESS: Conic and Polynomial Programming Toolbox. Massachusetts Institute of Technology, 2013. https:
/ //github.cc con// wttaa//spotttobooo/oooooo wwwpollbxxxpptteess.... spotless.
[44] Robert W Johnson. Algorithm 988: Amgkq: an efficient implementation of adaptive multivariate gauss-kronrod quadrature for simultaneous integrands in octave/matlab. ACM Transactions on Mathematical Software (TOMS), 44(3):1-19, 2018. 7
[45] MOSEK ApS. MOSEK Modeling Cookbook. 2022. URL https:////csss... mosek. com/ MOSEKModel htte:/ooee.nnncooooooooo--eeeeee
paf.
[46] Jeffrey Humpherys, Preston Redd, and Jeremy West. A fresh look at the Kalman filter. SIAM review, 54(4):801-88) -823, 2012. 8
[47] Andrew W Long, Kevin C Wolfe, Michael Mashner, Gregory S Chirikjian, al. The banana distribution is gaussian: A localization study with exponential coordinates. Robotics: Science and Systems VIII, 265:1, 2013, 9, 22
[48] Yihao Zhang, Odin A Severinsen, John Leonard, Luca Carlone, and Kasra Khosoussi. Data- association-free landmark-based slam. In 2023 IEEE International Conference on Robotics and Automation (ICRA), pages 8349-8355 IEFE, 2023. 9
[49] Martin Brossard, Silvere Bonnabel, and Jean- San-Philippe Condomines. Unscented Kalman filtering on Lie groups In Proceedings of the IEEE/RSJ International Conference on Intelligent Robots and Systems, pages 2485-2491. IEEE, 2017.
[50] Ioannis Sakiotis. Portable Numerical Library for the Calculation of Tulti-Dimensional Integrals. PhD thesis, Old Dominion University, 2024. 10
[51] Thomas Hahn. Cuba-a library for multidimensional numerical integration. Computer Physics Communications, 168(2):78 -5,7, 2005. 10
[52] Victor Magron and Jie Wang. TSSOS: a julia library to exploit sparsity for large-scale polyno- mial optimization. arXiv preprint arXiv:2103.00915 2021. 10
[53] Jean Bernard Lasserre. introduction polynomial and semi-algebraic optimization, vol- ume 52. Cambridge University Press, 2015. 16
13

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
We developed the MEM-KF that provides asymptotically optimal state estima- tion for nonlinear non-Gan Gaussian systems. the limiting case, the belief maintained in the filter converges to the true distributions.
Justification:
2. Limitations
Question: Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: We have included this part in the conclusions.
3. Theory Assumptions and Proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?
Answer: [Yes]
Justification: We assume the system is polynomials. We then referred the work [40] to show that the induced max entropy distribution can asymptotically approximate any belief that is moment-dettrmmintee.
4. Experimental Result Reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Answer: [Yes]
The algorithm can be reproduced providing that the 1) polynomial basis is computed correctly, 2) proper numerical integration method, e.g., Gaussian quadrature, is available, The algorithm block reveals all the information to reproduce the result.
Justification:
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [Yes]
Justification: We have provided the core algorithm for the main algorithm in the appendix.
6. Experimental Setting/eetails
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: The experiment classical problem in robotics that involves multimodal distributions to show the main asymptotic property of the MEM-KE.
7. Experiment Statistical Significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [Yes]
Justification: We have provided the evaluation of the experiment on substantial number of trials and provided the box plot.
8. Experiments Compute Resources
14

---

For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed reproduce the experiments?
Question:
Answer: [Yes]
Justification:
The experiments can be implemented on a laptop.
9. Code Of Ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics https: /ttps:////uuissspp httpr//s./.cpplllc/chccceooeoolo
Answer: [Yes]
Justification: The research respect the NeurIPS Code of Ethies.
10. Broader Impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [Yes]
Justification: We have included this part in the conclusions. This paper is foundational research on state estimation that does not have a direct path to negative societal impact.
Safeguards
Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Answer: [NA]
Justification: This paper is for nonlinear state estimation that does not have this concern.
12. Licenses for existing assets
Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [Yes]
Justification:
We have properly cited the references and the owners of the submitted materials.
13. New Assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [Yes]
Justification: All the data and code are well documented.
14. Crowdsourcing and Research with Human Subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification: Not used.
15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: Not used.
15

---

A Polynomial Optimization and Semidefinite Relaxations
In general, POPs are hard non-convex problems. In order obtain convex relaxation, first rewrite the POP as function of the moment matrix. The key idea here is that we can write polynomial c(x)) of degree up to 2r as linear function of the moment matrix (x),,, which contains all the monomials of degree up to 2r:
<(x) (r(CM, (C, M,),(
(18)
where we observe that (x). by construction, is a positive semidefinite rank-1 I matrix. We then relax the problem by dropping the rank-one constraint the moment matrix M, while enforcing:
X = 0, X1,1
(Semidefinite relaxation)
Since the moment matrix M, (z) contains multiple repeated entries, we further add moment con- straints that enforee these repeated entries be identical:
( = 0.
(Moment constraints)
Finally, enforce the equality constraints, using localizing constraints in the form:
(G,, X) = 0.
(Localizing constraints)
These steps lead to the standard semidefinite moment relaxation of (POP), see [22, 53].
Definition (Moment Relaxation of (POP) [22, 53)).
ه ;= in,(C, X) X
s.t.
X 0, 11,1 =1
(SDP)
B X) 0, (G, = Vi,j-
If rank-1 solution is obtained from (SDP), its solution X is valid moment matrix. In such case, extract X (.) from the first column of X.
B Auxiliary Matrices
This appendix formally defines the localizing matrices (By}y [53] that locate the monomials in moment matrix and its orthogonal complement { (B,,. 1y.0. These matrices are used the definition of the moment constraints for the moment relaxations for (POP).
We first introduce the matrix B, which retrieves the monomial x from the moment matrix M, (æ),
ie., (B,, Mz)) x7:
By = C. y=0+
1
Eeace C, = [[(o, 8) a +
(19)
where Ca [ is the canonical basis corresponding to the n-tuples No and C, is normalizing factor. We note that the moment matrix can be written as:
M, EC,B,'" 7
(20)
Then introduce B which is used to enforce that repeated entries (corresponding to the same monomial x²) ) in the moment matrix are identical. We first define the symmetric matrix, for any given a. No:
Ea.9 e Egee
(21)
Then we construct 7,a for some given 7 E NJ, and for a E No (with a < 7) as:
B y,ax E.,3 Ea.S,
(22)
where 3 is such that + 3 and B < a, and where ã, B E No are chosen such that a + 3 = and o is the smallest vector (in the lexicographic sense) that satisfies ā + B 7. An example of matrices B,, and B for simple problem is given in Appendix C.
16

---

C Examples of Moment Matrix and Constraints
For =2, 2, we have the moment matrix:
22
I
2j42
2 22
r
M2(2)=
22
T122
x2
31
31 22
31
22
2122
22
2122
For = [2, [2,0], which corresponds to zi.
we have:
0
1
0
0
1
0
0
B,
C. =3.
(23)
CT
1
0
0
0
0
0
0
Now let us compute B for [1,0] $ y- We note that / y
[1, 0] and choose
[0, 0] and
[2, 0]. The corresponding vectors in the canonical basis are:
Ca = 0
1
0
0
0
0] es
10
0 0 0],
(24)
ea =
10 2001 0],
from which we obtain:
B 7,a
Ea,
Ea,9 = eaey T
eze (exe] + e
0
0
0
+1
0
0
0
0
0
0
0
(25)
+1
0
0
0
0
0
0
By adding the constraint:
(B 7, X) 0,
X
(26)
enforce that the entries corresponding to 2 are identical.
D Derivation of MED
Consider the optimization (P,), given the multiplier À E R! we have the Lagrangian:
L(p(z), A) = / p(x) In dxx) + Σ da / x0p(d)d Z.).
(27)
For fixed p(æ), the dual problem to (P,) becomes:
max min p(m)
((p(x),
(28)
Via taking the variation of p(x) on L(p(x:), A), we have the first-order optimality condition:
V C(p(() A) p(æ) p(x)
+ Inp(x)) + Σ >xx="
(29)
where V the Fréchet derivative in normed linear space defined all possible distributions, For simplicity, we merge the constant 1 produced by the entropy term and let the (29) equal to zero. Then we have the stationary point (3):
p(x) = exp (- Σ )a
17

---

E Derivation of Thermal Dynamics Potential
Now we introduce the thermal dynamics potential term that can minimize to obtain the coefficients of (3). We mainly refer to the result in [42] to introduce the thermal dynamics potential and the equivalent potential-like term for the optimization. Consider the integrand:
(+())
exp (-
Σ dax
(30)
K
NENN
the normalized version of (3) can be expressed as:
exp (- Ea )
p" (x)
(31)
(0()
By substituting the p" (a:) to the inner minimization of minp(a) ((p(x), A),
we have the Legendre-Fente Fenchel dual as (a:), A)
1 p¹ (x)) (x) In lnpppp(p (x))da Σ dal /fc " T.)
E>al/
{ p'(x)) In p(x) dæ + Eval Aaal / x0p"(x0) Ta.)
((()
[ p¹ (1x))pp())da + Σda / x² (x0)dx o( p³ la Σ Aata
(32)
K
K
/ (n
Σ ta
In
Σ -I
We note that the coefficients Xag corresponding the normalization constraints can be extracted from (A) and we have:
T(A) In / (-
xx exp (~do)dæ + (Ao +
Σ
Auta)
aENN (ovo
aENp/(oo>
(33)
= In jk exp (-
Σ
d de +
Σ
AaTa
wEN://oo)dds
ooN:/(o0)
Given the fact that p(x) is normalized the stationary point:
S
x (-
Σ
X, ( exp (- (-Ao)dæ 1,
(34)
oENy/ [oo)
we have the no-nomlllzze thermal dynamies potential:
A(A) = //f exp
Σ
" - 1)+
Σ
(35)
aEN
aENT
The fact that A(A) = т(A) can be verified by extract the term Dag and Ag = In J exp (- EaErp/ Jao xx )dx
F Recover Distribution on Quotient Ring
Now, we consider nonempty IC that introduces an ideal structure to the systems, which is common in robotics. Consider the ideal generated by the polynomial equality constraints g(x):
I I[g(x)] =
Σ
gi (1))
e R[x
(36)
18

---

The integration of any term in [[g(x)) will be zero, i.,
[ p(æ) Σ (t)) = O, Va,
(37)
which can be checked by the fact that (x) = 0, Va E K. Thus, we need to incorporate these constraints when recovering À to avoid degenerate search direction or reducing the systems to the quotient ring [xxx//Igg(x))
Now, we proceed construct the ideal from the canonical basis (x) by linear algebra techniques. Consider the elements of [[<]/[[g(<>) with finite degree r; each element can be expressed as linear combination of the terms in the set:
B {œ" (x) deg x" (x S Vi, Va).
(38)
Let the vector of all elements in B as b(x), we can express b(x) linear combination of (2),, ie, b(x) (x).
(39)
Via QR decomposition of B i.e
R 0
B = 1Q Q
(40)
we can extract the null space of B, i.e, the span of column space of Q1, to formulate the basis of the quotient ring. The R is an upper diagonal matrix. Then, we conclude that
(xx) Q (x)
(41)
is a basis for the quotient ring. In this case, the max entropy distribution can be expressed in terms of (æ)):
p(z) exp (- A (x)),
(42)
and then we can minimize the thermal dynamics potential on the quotient ring:
A(x,) = f exp (-11,(())) + Σ Nny.
(43)
By the definition of X, and Q we have the equivalence relation of
Q1 E colsp(Q):
(44)
Thus, we can alternatively perform constrained optimization with moments on R[æ] by incorporating equality constraints to avoid the degenerate search directions:
min X
A))
(45)
s.t. T==
As the constraints are linear, the programming remains convex. When À is fixed, the gradients and Hessian can be computed by numerical integration and used to determine the optimal search direction for the optimal coefficients. As the convex optimization is not agnostic to the initial guess, we can initialize À by an MED with finite integrals over K., for example, the Gaussian distribution.
G Example of Extended Polynomial System and Its Affine Form
We provide an example of the extended polynomial system and transform to the affine form as in (6). Consider linear measurement model with state IR?, measurements y, and measurement noise w:
V1 Ij W1,
(46)
1/2 T2
W2-
The extended system 2 is:
Φ x) = (w)
(47)
19

---

where
1
1 91 31
11 31
10
Ør (y- x)
1/2 J2
1/2 22
(w)(
W2
(y1
21 (8/2 22))
9192 + F1
3 12 y231
W1
w2
(11 11 )
2y1 21
(3/ T2)
2 22
22
w02
(48)
We then have the extended system
A(y) v2
(49)
where
0
001
D1
0
1
0
0
0
102
(x))
322
A(y)
1
0
(50)
W1 W12
31
O
0
w
2y/1
02
0
0
2y2
0
H More Numerical Results
H.1
Update Steps
We compare the update steps in the linear system (??) using the noise model (14) with different orders of moments constraints in Figure 6. We find that with higher moments, the density is more concentrated and converges faster to the ground truth.
H.2 Prediction Steps
We only apply the prediction steps to the system on the matrix Lie group:
Z Z W& SE(2),
(51)
and compare the recovered distributions with Monte Carlo simulations. We find that with larger r. better recover the distribution. As marginalizing the distribution p(æs. A requires symbolic integration that is generally intractable, we apply the Langevin dynamics on SFF(2) to sample and visualize the position part of MED. The MED is sampled by Langevin dynamics and visualized in y plane in Figure 7.
H.3 Comparison with Particle Filter
We compare MEM-KF with particle filters with N 1e1, 1c2, 1e3, 1e4 and 1e5 particles. For both cases in Section 7.3, we sample the max entropy association strategy to avoid overconfidence in the observation models. We consider the threshold for resampling as 0.5 among all the trials. The estimated trajectories of all the cases are shown in Figure 5. The box plot statistics for the position error is shown in Figure 8.
find that as more particles are considered, the particle filters converge to better solutions. In terms of numerical results, MEM-KF has compatible solutions with slightly larger median error than the particle filters with 2 le3. However, we note that MEM-KF a deterministic method that maintains smooth analytical belief representation that convenient for the extraction of global point estimation. Such property is not possible for particle filters that are based on random sampling.
20

---

N
N= 2
N 3
N=
N
.V
Figure 6: From the to bottom, using noise model with moments with order up to r = 2, 4, 6, 10 and 12 matched. As more measurements accumulated, the belief converges the state. With larger r. the belief converges much faster the Dirac measure.
21

---

r
2
7 7=
4
3.5
3.5
Monte Carlo Simulation Max Entropy Distribution
2.5
2.5
0.5
0.5
-0.5
-0.5
2
2
r
r
3.5
3.5
2.5
2.5
0.5
-0.5
-0.5
2
2
Figure 7: The propagation of the helief SF((2 with only prediction steps. We visualize the position part in system dynamics, and red dots
the - y plane. The blue dots represent Monte Carlo simulations represent the MED sampled by Langevin dynamics using the log likelihood of the recovered distribution. We find that with r the banana shape [47] is not fully captured by the MED. With larger r. the MED better captures the underlying distribution, especially the tail of the banana shape.
22

---

10-1
10-3
10-3
10
10
lnl
Nate2
Nade3
Notes
Notes
AEIM-RFF
Notel
Notel
Nodled
Nodes
ARM-KF
Figure 8: Location error for localization tasks with unknown data association. (a)The underlying real association follows the same distribution in (17). (b)The measurement comes more from the landmarks at the lower left corner We compare the 2-norm between the estimated position and the ground truth for all 100 trajectories. We compared the particle filters with N particles with the proposed MEM-KF KF with polynomial order up to r 4. We find that the particle filter with a sufficient number of particles and proposed JEM-KF can estimate trajectories with similar quality.
23


Reviews:
[{"review_id": "O5cJMBOUpi", "reviewer": "Reviewer_5jz7", "summary": "This paper considers the problem of filtering, i.e., state estimation, in non-linear, non-Gaussian systems. This is a well-studied, classical problem with contributions dating back decades and literature too vast to summarize in any meaningful way. Most contributions to this problem are, in some sense, inspired by attempts to approach the ideal of Kalman filtering in linear Gaussian models. This is done either through computationally intensive approaches such as particle filtering or Markov chain Monte Carlo or through careful design of approximate, yet analytically fast, update steps of the filtering distribution while keeping model expressivity. This paper is a contribution to the latter. The authors introduce a novel update scheme utilizing MED in order to perform filtering in a non-linear, non-Gaussian state space model. The paper assumes a belief (i.e. state) distribution with some fixed number of moments having maximum entropy. The key idea is that by thinking of the belief distribution in terms of its moments only, the prediction step can be viewed in terms of predicting the moments directly using an auxiliary dynamical system. Then the corresponding belief MED is computed by solving an optimization problem given the predicted moments. The methodology is tested on several examples, ranging from simple to more complex.", "questions": "Is there some relatively fast way to approximate the gradients that are used to find the predicted belief distribution? It would be good to see how well this method does in larger-scale problems when approximations are used.\n\nCan you extend this approach to handle a variable number of moments to save computation time? Behave like a linear Gaussian system when necessary, and use more moments when necessary?\n\nCan you consider dynamics in high-dimensional state spaces where we have a higher-moment model for some of the latent variables?", "limitations": "Yes.", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 3, "strengths_and_weaknesses": "Strengths: the paper is well-motivated and addresses a still-relevant problem. I like the use of simple and transparent examples to illustrate the methodology, making it easy to replicate and test it in practice. The comparison with particle filtering as a sanity check is appreciated. The limitations (i.e., reliance on numerical integration) are clearly stated. The approach appears to be theoretically grounded. Reliance only on convex optimization despite the potential for arbitrarily complex distributions is an advantage. The approach seems well-suited for problems with complex dynamics in low-dimensional state spaces.\n\nWeaknesses: reliance on numerical integration and optimization steps for every update. No good understanding of the computation time is provided for this method, even with comparison to e.g., particle filtering. It seems that the scalability issue may be harder to overcome and parallelization is not developed in any detail other than mentioning the use of GPUs, but it not clear how this can be done in practice. Perhaps the authors can mention some potential approximate schemes that can be used?", "was_revised": true, "final_justification": "The authors have made some nice elaborations on my questions, and I have adjusted my score, though I am still not convinced about the practical constraints in overcoming the computational burdens in this method that would allow wide deployment.", "rebuttal": "We thank the reviewer for the great suggestions and appreciations of the novelty of this work. \n\n**Time**\n\nWe note that all the UKF, EKF, UKFM, and InEKF can be run in real-time in the cases considered in this work, both prediciton and update can be completed less than $1e-4s$. For the PF with $1e5$ particles, the prediction and update steps with resampling and CPU parallelization require less than $0.2s$. \n\nFor the MEM-KF in the localization example, the mean (std.dev) time (secs) is $6.0165e^{-4}$ $(2.0275e^{-4})$ in the prediction, $6.8753$ ($8.7022$) in the update, and $34.5102 $ ($26.2063$) in recovering the MED from moments. The Laptop is equipped with Intel i7-11850H CPU. The convex optimization is run until the KKT violation is below $1e-4$. \n\nWe will add the table to the revised version. We admit that the cost for integration is still high for real-time deployment at this time. The complexity and possible GPU implementation is shown below: \n\n**The complexity:**\n\nWe show some theoretical analysis of the scalability of the proposed method. The complexity is mainly determined by the dimension of the system $n$ and the order of moments $r$. The number of the monomials for such system is $$s(r, n):= |\\mathbb{N}_r^n| = \\begin{pmatrix}\n    n+r \\\\ \n    n\n\\end{pmatrix} = \\frac{(n+r)!}{n!r!}.$$ \nFor fixed $r$, we have $s(r, n) \\sim \\mathcal{O}(n^r)$ by eliminating the $n!$ in the denominator: $$ s(r, n) = \\frac{(n+1)(n+2)\\cdots(n+r)}{r!}, \\frac{(n+r)^r}{r!} \\ge s(r, n) \\ge \\frac{n^r}{r!}.$$\nFor fixed $n$, we have $s(r, n) \\sim \\mathcal{O}(r^n) $ by eliminating the $r!$ in the denominator: $$ s(r, n) = \\frac{(r+1)(r+2)\\cdots(r+n)}{n!},  \\frac{(r+n)^n}{n!} \\ge s(r, n) \\ge \\frac{(r+1)^n}{n!}.$$\nThus, the terms to evaluate scales polynomially w.r.t the order and exponentially w.r.t the dimension. \n\nFor the numerical integration, the number of quadrature scales exponentially w.r.t. $n$. Say that the number of regions is $c^n$.  \n\nCompute the moments: Requires evaluating $s(r, n)$ terms of moments by numerical integration. These require $s(r, n) \\times c^n $ operations. \n\nPropagation of moments: We need to compute $x_{\\alpha, k+1}=Ax_{\\alpha, k}$,  which requires $\\mathcal{O}(s(r, n)^2)$ multiplications and $\\mathcal{O}(s(r, n)^2)$ additions. To see this, we find each entry of $\\bar{x}_{\\alpha, k+1}$ needs $s(r, n)$ of multiplication and add them together.\n\nReconstruction of distribution: We need to run the convex optimization, that each gradient descent to compute the moments of $\\phi(x)$, which involves numerical integrations of $s(r, n)$ terms of the $ n$-dimensional system.  \n\nUpdate with measurements: This step requires $s(r, n)$ operation of addition to update the coefficients given measurements and prior belief.\n\nNormalization: This step requires integrating the unnormalized belief once. \n\nIn summary, the overall complexity scale exponentially w.r.t the dimension $n$ and polynomially w.r.t the order $r$. We will include these discussions in the revised version. \n\n**The GPU implementation:**\n\nWe agree that the integration is the main bottleneck of the proposed method, given the complexity of deterministic numerical integration. However, we note that the proposed algorithm is suitable for parallelization.\n\nFor numerical integration that uses adaptive quadrature, the idea is to partition the domain into grids of different sizes based on the tolerance and shape of the functions to integrate. The region with more peaks requires more grids, namely $N_g$. In each region, the point to evaluate to integration requires $N_p$ points. Thus, the operation will require computing the $ N_g * N_p * s(r, n)$ number of operations to evaluate the integrations. Such operations only require matrix multiplications that are suitable for parallelization. Usually, the $N_g$ scales exponentially w.r.t the dimension. \n  \nIn our original implementation, we used a recursive adaptive quadrature algorithm ``Robert W Johnson. Algorithm 988: Amgkq: an efficient implementation of adaptive multivariate gauss-kronrod quadrature for simultaneous integrands in octave/matlab. ACM Transactions on Mathematical Software (TOMS), 44(3):1–19, 2018.'' that supports the CPU parallelization to compute integration of vector value functions on grids.\n\nFor future work, a possible GPU implementation is presented in the following dissertation and open-sourced : ``A Portable Numerical Library for the Calculation of Multi Dimensional Integrals, Ioannis Sakiotis''. \n\n**GPU memory expectation**\n\nWe consider the sparse grid shown in ``Cuba library for multidimensional numerical integration''. Note that the thesis with GPU implementation considers denser tensor-product-based grids. We here consider a less conservative one. \n\nSay that the term to integrate is $s(r, n)$. The number of partitioned regions approximately satisfies $100 * 1.3^n$. We consider the table in Section. 5 of the above paper. Each point needs $n$-number to represent its value. With the degree-9 rule, the points in each region is $C=[153, 273, 453, 717, 1105, 1689, 2605, 4117, 6745]$ from dimension 4 to 12. Thus, we can approximate the memory in GB by $C_{n+r}^r * (100*1.3^n) * C(n) * n / (1024^3) * 4$. \n\nThus, by approximation, we find that for $n = 6, r = 10$, the memory consumption is 39GB for the input using Float32. The dimension of the moment is approximately $8000$. The $r=10$ can represent highly nonlinear distributions, which can represent highly complicated distributions. \n\nFor $n=8, r=4$, the memory consumption is $13.3$ GB. While for $n=8, r=6$, the memory becomes $80.7 GB$. \n\nThus, with full moment space, we may expect dimensions up to **n=8** with **r=4** to be deployed on an Nvidia Orin for real-time computation. Then the actual time depends on the FLOPS of the GPU. \n\n\n**Q1:**\n\nYes, we can apply Monte Carlo integration to approximate the gradients and perform the optimization. Based on our experiments, such stochastic integration is faster than deterministic integrations when setting a smaller budget for function evaluation. Then the convex optimization becomes a stochastic gradient descent. However, such stochastic optimization is less stable, as not all coefficients of the MED ensures the convergence of integrations. The noisy gradients make the optimization easier to diverge. In the localization case, we set the budget to $1e7$ to ensure a relatively stable convergence of the convex optimization. \n\nThough MCMC is one way to approximate the gradients, it violates our initial idea to develop a completely deterministic and certifiable algorithm to solve the nonlinear filtering problems. If we have such particles, it would be more reasonable to use them as in PF.\n\nOther than MCMC, there are possible sparse representations that may help to reduce the burden, which can be referred to in the following paper with sparsity pattern in $A$ matrix in the dynamics taken into consideration. We consider this as an interesting future direction. \n\nVictor Magron and Jie Wang, Sparse Polynomial Optimization: Theory and Practice, World Scientific Press, 2023\n\n**Q2:**\n\nWe thank the reviewer for this great suggestion. \n\nYes, an approximation strategy is to set a threshold of covariance to decide whether to switch to Gaussian. The logic is that MED with $r>2$ always has lower entropy than the Gaussian distribution with the same covariance. Thus, when the covariance is smaller than a threshold, the distribution can be replaced by such a Gaussian distribution without the risk of overconfidence.\n\nIn the localization case, after the update steps, we first integrate the moments with $r\\le2$ to check the covariance. If the covariance $trace(P_{x,y}) \\le 0.5^2$, we stop and do not proceed to $r>2$.  Otherwise, we also include the higher-order moments. Then we predict the moments with order $r$ determined by the covariance after update. We note that in the update steps, we always maintain the full observation model to incorporate the multimodality of sensors. \n\nWith this method, we show similar results compared to the full MEM-KF in the localization case. After a few iterations, we can do $r=2$ in the prediction, and the belief converges to a unimodal distribution that can be well represented by $r=2$. \n\nIn this case, the effort to recover the MED on $SE(2)$ reduces from evaluating $s(4, 4) = 70$ to $s(2, 4) = 15$ terms when recovering the distribution or computing the moments. Note that with equality $SO(2)$ constraints, there is still no closed form on $SE(2)$.\n\n**Q3:**\n\nWe thank the reviewer for these suggestions! The primary scenario we consider is in robotics, where the dynamics and measurements are exact polynomials. \n\nThe most relevant latent polynomial model is the deep Koopman operator that is linear:\nSee: Lusch, B., Kutz, J.N. & Brunton, S.L. Deep learning for universal linear embeddings of nonlinear dynamics. Nat Commun 9, 4950 (2018).\n\nHowever, such a system requires a very high-dimensional latent to make the dynamics linear systems. \n\nWe will make this limitation and the applicable scenario clear in the revised version."}, {"review_id": "b1ynDk3icQ", "reviewer": "Reviewer_VyDk", "summary": "The paper proposes a Max Entropy Moment Kalman Filter (MEM-KF), to tackle the failure of Kalman Filter (KF) in non-linear, non-Gaussian state space model settings.Unlike KF which assumes Gaussian distributions, MEM-KF uses a Moment-Constrained Max-Entropy Distribution (MED) to model state x and observation y. This paper provides quite solid mathematic deductions for 1) why MED can asymptoticly approximate any Moment-Determinate Distributions 2) how to use MED to solve the predict and update steps of a recursive Bayes filter. Although the computation complexity is relatively high, simulations show that MEM-KF outperforms traditional KF and KF variants in a non-linear, non-Gaussian localization problem setting.", "questions": "1. Significance of the Work\n\nCan MEM-KF handle partially known or imperfect state-space models? For example, what happens if the dynamics f, observation model h, or noise distributions w, v are only partially known or imperfect? This is often the case in real-world applications.\n\nAre moment-determinate distributions—which are neither heavy-tailed nor exhibit rapidly growing higher-order moments—appropriate for modeling practical systems such as robot localization? It would be valuable to justify whether this assumption aligns with common real-world settings.\n\n2. Benchmarking\n\nWhat is the runtime of MEM-KF compared to other Kalman filter variants (e.g., EKF, UKF, PF)? It would be helpful to include a runtime comparison table to evaluate computational efficiency.\n\nBeyond model-based algorithms, when sample size N grows, how is MEM-KF compared with neural network(NN)-based methods e.g. KalmanNet [1]? Including NN-based algorithms in the benchmark comparison would provide a more comprehensive evaluation.\n\n3. Side Questions\n\nHow should one choose the highest moment order r, considering the trade-off between accuracy and computational efficiency? For example, how much does the computational complexity increase per predict and update step as r increases, and what is the corresponding performance gain?\n\nHas the numerical stability of the algorithm been evaluated, particularly when dealing with high-order moments, which may take on large values and affect optimization?\n\n[1] Revach, Guy, et al. \"KalmanNet: Neural network aided Kalman filtering for partially known dynamics.\" IEEE Transactions on Signal Processing 70 (2022): 1532–1547.", "limitations": "See Questions.", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 3, "significance": 2, "originality": 3, "strengths_and_weaknesses": "Strengths:\nThe theoretical framework and claims in this paper are well-supported by rigorous mathematical proofs. The writing is generally clear and accessible, making the technical content relatively easy to follow and understand.\n\nWeaknesses:\nThe practical applicability of the proposed MEM-KF is somewhat questionable. The strength of MEM_KF is on improved prediction accuracy at the cost of increased computational overhead. However, the prediction accuracy of MEM-KF relies on complete and accurate model information for belief propagation—including precise knowledge of the system dynamics f, observation model h, and the noise distributions w and v. In many real-world scenarios, however, these models are imperfect or only partially known, which limits the utility of such a method in practice.\nAdditionally, the Moment-Constrained Max-Entropy Distribution (MED) is not a novel contribution of this paper. The idea and its solution via convex optimization were already developed in [1]. As such, the paper would benefit from a more thorough discussion of why MED is particularly suitable in this context. For example, MED approximates well the moment-determinate distributions, which are neither heavy-tailed nor exhibit rapidly growing higher-order moments. Is this assumption realistic in practical applications such as robot localization? A deeper justification and empirical support for this modeling choice would significantly strengthen the paper.\n\nReference:\n[1] Bandyopadhyay, K., et al. \"Maximum entropy and the problem of moments: A stable algorithm.\" Physical Review E, 71.5 (2005): 057701.", "was_revised": true, "final_justification": "From the rebuttal discussions, my main concerns about the followings are resolved.\n1. The suitability of this algorithm for practical filtering cases: it can address partially known state space model; the moment-determinancy is a reasonable assumption.\n2. The originality of applying MED to filtering, and justification of its suitability.\n3. Compared with SoTA machine learning based filtering methods, e.g. KalmanNet, it has the advantage of robustness against data distribution shift.\nWhile the low-dimensional limitation is still there, I believe flexible choice of r and GPU implementation could mitigate this issue.\nOverall, I decide to raise my score to weak accept 4.", "rebuttal": "We thank the reviewer for these great suggestions and appreciations of this work:\n\n**Originality of MED**\n\nWe agree that MED is not invented by us. We mainly refer to the multivariate version ``Lawrence R Mead and Nikos Papanicolaou. Maximum entropy in the problem of moments. Journal of Mathematical Physics, 1984.''\n\nHowever, to the best of our knowledge, there is no literature applying MED in filtering. The use of MED enables a convex optimization-based method with asymptotical performance guarantees, which has not been seen before.\n\nAdditionally, we considered a constrained feasible set by polynomial equality that needs to incorporate the quotient ring structure, which is novel in the polynomial systems. This part is shown in Appendix. F. Compared to [40], we are using a line-search method with BFGS instead of Newton's method, which is more efficient to avoid integration for the exact Hessian matrix that has far more terms. \n\nWe will have a remark about the difference between our method and the literature in the revised version.\n\n**Q1-1: Partially known model**\n\nYes, the MEM-KF is capable of handling the partially or unknown models. The philosophy is that we model the unknown part in a max-entropy manner, which avoids overconfidence.\n\nIn our example shown in Fig. 4, the actual landmark association is not as modeled for the update; thus, both $h$ and the associated noise $v$ remain unknown. The KF-type baselines need to choose the most likely $h$ from four possible landmark positions. In the lower row, the observation comes more from one landmark, while the model we use considers the measurement to be uniformly distributed from all landmarks. We note that the uniform distribution model has the largest entropy, which does not lead to overconfidence in a certain pattern. \n\nWe show that the max-enropy model we use for update leads to the best generalization to unseen environments, as it does not overfit. In the lower row, even though the distributions are completely different, our method can still track the trajectories. While the baseline using the most-likely association strategy has a clear bias.\n\nWe will highlight this part in the revision. \n\n**Q1-2: determinancy**\n\nWe thank the reviewer for raising this important point. We show that the moment-determinacy is a reasonable assumption in engineering applications.\n\nAs established in the classical Hausdorff moment problem, any multivariate distribution supported on a bounded domain such as $[0, 1]^n$ is moment-determinate. This result extends to general compact subsets of $\\mathbb{R}^n$ via standard scaling arguments.\n\nSee: Kleiber, C. and Stoyanov, J. (2013). Multivariate distributions and the moment problem. Journal of Multivariate Analysis, 113, pp.~7--18.\n\nIn engineering contexts such as robotics, it is typically reasonable to assume that process disturbances and sensor noise are bounded. Otherwise, it would imply unphysical behavior—such as injecting infinite energy into a mechanical system. As such, the assumption of compactly supported distributions is not only mathematically convenient but also physically justified.\n\nWhile uncertainty may grow under long-term propagation, especially in partially unobservable systems, over any \\emph{finite time horizon} and with \\emph{bounded noise}, the state distribution remains supported on a compact domain (barring pathological singularities). Therefore, the use of moment-determinate distributions is a valid and practically meaningful assumption in our setting.\n\nA more rigorous integration strategy is to manually choose a compact support when recovering the distribution. When the compact support in the integration is sufficiently large to include the actual support of noise or state, we can ensure the unique convergence. We also note that as long as the actual support is contained, there is no overconfidence in the state estimation. To see this, consider the optimization $P_i$ to get the MED:\n$$ E_i = \\max_{p(x)} \\int_{\\mathcal{K}_i}p(x)\\operatorname{ln}{p(x)}dx $$\nWe can see that entropy satisfy $E_2 \\ge E_1$ if $\\mathcal{K}_1 \\subset \\mathcal{K}_2$ by checking:\n$$p_2(x)=\n\\begin{cases}\np^*_1(x), & x \\in \\mathcal{K}_1, \\text{or}, \n0,      & x \\in \\mathcal{K}_2 \\setminus \\mathcal{K}_1.\n\\end{cases}$$ is also a feasible point to $P_2$. Thus, we can always choose a larger support to ensure no overconfidence in the estimation.\n\nIn conclusion, assuming a sufficiently large compact support is sufficient to ensure the moment determinacy assumption is realistic in real-world settings.\n\nIn our practice, choosing a support sufficiently large has little difference compared to the entire space due to the role of Gaussian quadrature.\n\nWe will add the compact regularization condition to the revised version.\n\n**Q2-1 Runtime**\n\nWe note that all the UKF, EKF, UKFM, and InEKF can be run in real-time in the cases considered in this work, both prediciton and update can be completed less than $1e-4s$. For the PF with $1e5$ particles, the prediction and update steps with resampling and CPU parallelization requires less than $0.2s$. \n\nFor the MEM-KF in the localization example, the mean (std.dev) time (secs) is $6.0165e^{-4}$ $(2.0275e^{-4})$ in the prediction, $6.8753$ ($8.7022$) in the update, and $34.5102 $ ($26.2063$) in recovering the MED from moments. The Laptop is equipped with Intel i7-11850H CPU. The convex optimization is run until the KKT violation is below $1e-4$. \n\nWe will add the table in the revision. \n\n**Q2-2 KalmanNet**\n\nWe thank the author for suggestions on using KalmanNet as a benchmark. The KalmanNet is a data-driven method that requires labeled ground-truth data to supervise the network that determines the Kalman gain. The loss function is designed to minimize the predicted loss over the entire trajectory. \n\nThe problem with KalmanNet is the overfitting to the training data set, and it has no guarantee that it will generalize.\n\nIn Fig. 4, we choose the uniform association strategy that has the largest entropy. We can find that such a strategy with minimized over-confidence leads to generalization to different associations. The proposed methods generalize to both scenarios with guarantees from the upper bound of entropy.\n\nIn our tests, we train the KalmanNet with ground truth from the simulator. We find that if the landmark has a certain pattern, such as coming from one of the landmarks, the NN overfits to it and has compatible results with our methods. When the environment changes, e.g., the measurement comes from another landmark, the results degrade dramatically.\n\nWe will highlight this part in the revision. \n\n**Q3-1**\n\nThe feasible set to problem ${P}_r$ will only gets more smaller as $r$ goes up to account for more moment constraints. With smaller support, the entropy of the recovered distribution will only get lower. Thus, the trade-off is that one should try the highest possible $r$ before the numerical computation becomes unstable or out of computational budget.\n\nAn empirical way to determine the sufficient $r$ is by the condition number of the Hessian matrix $H_r := \\int \\phi(x)\\phi(x)'p(x)dx$ of the recovered distribution by the moments $\\phi_{2r}(x)$. The insight is that if the distribution is around some lower-dimensional manifolds, when $r$ gets higher, there will be sufficiently many monomial basis to approximate such manifolds, say $c(x) = \\theta'\\phi(x) = 0$. Thus, a certain subspace of the Hessian will have very small eigenvalues. Say that $\\int c'\\phi(x)\\phi'(x)cp(x) dx = c'Hc = 0$. For example, the example in Fig. 2 needs to be terminated at $r = 12$ as a $12$-th order polynomial is sufficiently good to approximate the pattern. When $r>12$, the $H$ is very degenerate, that makes the gradient descent harder. For the binary distribution with four peaks, a fourth-order polynomial is also sufficient to have four local minimum to represent the pattern. Thus, the moment matrix with $r>4$ will have degeneracy to make the optimization ill-conditioned.\n\n**Q3-2:**\n\nWe show some theoretical analysis of the scalability of the proposed method. The complexity is mainly determined by the dimension of the system $n$ and the order of moments $r$. The number of the monomials for such system is $$s(r, n):= |\\mathbb{N}_r^n| = \\begin{pmatrix}\n    n+r \\\\ \n    n\n\\end{pmatrix} = \\frac{(n+r)!}{n!r!}.$$ \n\nThus, the terms to evaluate scales polynomially w.r.t the order and exponentially w.r.t the dimension. (See responses to other reviewers for details.)\n\nFor the numerical integration, the number of quadrature scales exponentially w.r.t. $n$. Say that the number of quadratures is $c^n$.  \n\nCompute the moments: Requires evaluating $s(r, n)$ terms of moments by numerical integration. These require $s(r, n) \\times c^n $ operations. \n\nReconstruction of distribution: We need to run the convex optimization, that each gradient descent to compute the moments of $\\phi(x)$, which involves numerical integrations of $s(r, n)$ terms of the $ n$-dimensional system.  \n\nUpdate with measurements: This step requires $s(r, n)$ operation of addition to update the coefficients given measurements and prior belief.\n\n**Q3-3:**\nThe performance monotonically increases with more moments. The marginal gain depends on the information contained in the moments. For the binary case, from r=2 to r=4 improves a lot, while higher does not help. For the case in Fig. 2, the approximation improves until r=12. \n\n**Q3-4:**\n\nThe moment problem is numerically challenging when the state is off the center, which results in large numbers. We admit a centralize-and-regularize strategy to avoid the problem of ill-posed moments when the order gets large. The idea is to apply a linear transformation $z = Ax + b$ to ensure the $z$ has zero mean and smaller covariance for integration. Obtaining $A$ and $b$ is the same as shifting a Gaussian distribution, given that the covariance and mean of $x$ are known. The moment of $z$ can be computed by another linear transformation."}, {"review_id": "xVJTh5ticP", "reviewer": "Reviewer_WTDJ", "summary": "This paper proposes a generalization of the Kalman filter called the MaxEnt Moment Kalman Filter (MEM-KF). MEM-KFs represent uncertainty through maximum entropy distributions with constrained moments of polynomials with pre-specified degree. Unlike the Kalman filter, which maintains two moments (mean and covariance), MEM-KF can represent more expressive distributions and is therefore well-suited to nonlinear filtering problems. The contributions of the work include a complete framework for modeling and filtering using MEM-KFs along with experiments to validate its characteristics on simulated data.", "questions": "1. How does MEM-KF compare numerically to particle filter and the various extensions of Kalman Filters? How do these methods compare when applying to problems complicated dynamics such as Lorenz-63 or those involving real-world robotics data?\n2. How do each of the components of MEM-KF scale with respect to $n$ and $r$ (memory and computation)? What is a realistic upper bound for the state-space dimensionality one could expect to perform on say a single GPU?\n3. With respect to numerical integration specifically, how does MED recovery behave as the dimensionality of the state space increases?", "limitations": "Additional discussion about scaling with respect to state-space dimensionality and polynomial degree would improve the paper.", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 3, "strengths_and_weaknesses": "## Strengths\n- **Meaningful Problem Domain**. The problem of nonlinear filtering is becoming increasingly important, with applications ranging from continual learning to computational neuroscience. \n- **Originality**. This is a creative extension of the Kalman filter that applies moment-constrained maximum entropy distributions to the problem of Bayesian filtering.\n- **Qualitative Exploration**. There are a good number of ablations and plots to demonstrate the qualitative performance of MEM-KFs.\n\n## Weaknesses\n- **Experiments**. Although the experiments contain a number of visualizations, the paper lacks a results table for the SE(2) experiment. It would also be informative to include experiments with more complicated state-space dynamics, such as state estimation in a Lorenz-63 system, or real-world data.\n- **Analysis of scaling.** Since MEM-KFs represent additional computation relative to the Kalman filter, there should be some analysis (ideally theoretical or at the least empirical) about how the required memory and computation scale with dimensionality n and degree r. This is important because of the variety of sub-procedures in the filtering algorithm (moment computation, propagation, point estimation, and so on). The small scale of the experiments indicates that this is indeed a significant limitation of the method. Though not a dealbreaker *per se*, it is important to discuss these limitations thoroughly in the paper.\n- **Practicality.** The proposed algorithm relies on numerical integration, which may be impractical for high-dimensional state spaces.", "was_revised": true, "final_justification": "The authors provided a theoretical analysis about scalability and quantitative results for the SE(2) experiments. They also explained why it is difficult to include additional experiments. The limitations on scalability are nontrivial. However, assuming these limitations are acknowledged, the paper meets the criteria for acceptance. I have increased my quality rating and overall rating accordingly.", "rebuttal": "We thank the reviewer for these great suggestions!\n\n**Q1-1:**\n\nThe results compared with KFs are shown in Fig. 4. The numerical difference is very obvious as these methods are inherently incapable of handling multimodality. They just drift for the problem we show.\n\nFor the upper row, the RMSE of the position drift for EKF, InEKF, UKF, UKFM and MEM-KF are $$\\text{Fig. 4, upper row: } 2.0893, 1.4062, 1.1619, 1.3638, 0.1155$$ and $$\\text{Fig. 4, lower row: } 1.2618, 1.4925, 1.4638, 1.4999, 0.1159 $$ for the lower row. We can see that the MEM-KF is at least 10 times more accurate in RMSE than the baseline in terms of the position error. \n\nThe comparison with PF using different particles is shown in Fig. 7 (trajectories) and Fig. 8 (numerical statistics) in the Appendix. G. \n\nWe will move the PF plot and the table about the runtime to the main script in the revised version for clarity! \n\n**Q1-2:**\n\nWe sincerely appreciate the reviewers’ suggestion to further include experiments on Lorenz systems or real-world robotics data. We fully agree that such additions can further demonstrate the general applicability of MEM-KF. However, due to the tight rebuttal timeline and the technical challenges in implementing stable numerical integration, configuring GPU-based solvers for real-world pipelines (usually $SE(3)$ pose data using $\\deg 12$ polynomial) or the Lorentz system ($n=4$ but hard to be solved by the multivariate recursion-based adaptive integration scheme), we were unable to include such experiments in this version.\n\nThat said, we emphasize that the core contribution of MEM-KF lies in the deterministic and convex formulation. Our current evaluations, including the challenging SE(2) localization task with unknown data association, already demonstrate both the theoretical soundness and practical viability of our method in structured and realistic settings.\n\nWe see extending to more chaotic or high-dimensional systems as promising future work and will highlight this direction in the revised version.\n\n**Q2-1, complexity:**\n\nWe show some theoretical analysis of the scalability of the proposed method. The complexity is mainly determined by the dimension of the system $n$ and the order of moments $r$. The number of the monomials for such system is $$s(r, n):= |\\mathbb{N}_r^n| = \\begin{pmatrix}\n    n+r \\\\ \n    n\n\\end{pmatrix} = \\frac{(n+r)!}{n!r!}.$$ \nFor fixed $r$, we have $s(r, n) \\sim \\mathcal{O}(n^r)$ by eliminating the $n!$ in the denominator: $$ s(r, n) = \\frac{(n+1)(n+2)\\cdots(n+r)}{r!}, \\frac{(n+r)^r}{r!} \\ge s(r, n) \\ge \\frac{n^r}{r!}.$$\nFor fixed $n$, we have $s(r, n) \\sim \\mathcal{O}(r^n) $ by eliminating the $r!$ in the denominator: $$ s(r, n) = \\frac{(r+1)(r+2)\\cdots(r+n)}{n!},  \\frac{(r+n)^n}{n!} \\ge s(r, n) \\ge \\frac{(r+1)^n}{n!}.$$\nThus, the terms to evaluate scales polynomially w.r.t the order and exponentially w.r.t the dimension. \n\nFor the numerical integration, the number of quadrature scales exponentially w.r.t. $n$. Say that the number of quadratures is $c^n$.  \n\nCompute the moments: Requires evaluating $s(r, n)$ terms of moments by numerical integration. These require $s(r, n) \\times c^n $ operations. \n\nPropagation of moments: We need to compute $x_{\\alpha, k+1}=Ax_{\\alpha, k}$,  which requires $\\mathcal{O}(s(r, n)^2)$ multiplications and $\\mathcal{O}(s(r, n)^2)$ additions. To see this, we find each entry of $\\bar{x}_{\\alpha, k+1}$ needs $s(r, n)$ of multiplication and add them together.\n\nReconstruction of distribution: We need to run the convex optimization, that each gradient descent to compute the moments of $\\phi(x)$, which involves numerical integrations of $s(r, n)$ terms of the $ n$-dimensional system.  \n\nUpdate with measurements: This step requires $s(r, n)$ operation of addition to update the coefficients given measurements and prior belief.\n\nNormalization: This step requires integrating the unnormalized belief once. \n\nIn summary, the overall complexity scale exponentially w.r.t the dimension $n$ and polynomially w.r.t the order $r$. We will include these discussions in the revised version. \n\n**Q2-2, GPU implementation:**\n\nWe agree that the integration is the main bottleneck of the proposed method, given the complexity of deterministic numerical integration. However, we note that the proposed algorithm is suitable for parallelization and has provable asymptotic properties of polynomial approximations.\n\nFor numerical integration that uses adaptive quadrature, the idea is to partition the domain into grids of different sizes based on the tolerance and shape of the functions to integrate. The region with more peaks requires more grids, namely $N_g$. In each region, the point to evaluate to integration requires $N_p$ points. Thus, the operation will require computing the $ N_g * N_p * s(r, n)$ number of operations to evaluate the integrations. Such operations only require matrix multiplications that are suitable for parallelization. Usually, the $N_g$ scales exponentially w.r.t the dimension. \n \nIn our original implementation, we used an adaptive quadrature algorithm ``Robert W Johnson. Algorithm 988: Amgkq: an efficient implementation of adaptive multivariate gauss-kronrod quadrature for simultaneous integrands in octave/matlab. ACM Transactions on Mathematical Software (TOMS), 44(3):1–19, 2018.'' that supports the CPU parallelization to compute integration of vector value functions on grids.\n\nFor future work, a possible GPU implementation is presented in the following dissertation and open-sourced : ``A Portable Numerical Library for the Calculation of Multi Dimensional Integrals, Ioannis Sakiotis''\n\n**GPU memory expectation**\n\nWe consider the sparse grid shown in ``Cuba library for multidimensional numerical integration''. Note that the thesis with GPU implementation considers denser tensor-product-based grids. We here consider a less conservative one. \n\nSay that the term to integrate is $s(r, n)$. The number of partitioned regions approximately satisfies $100 * 1.3^n$. We consider the table in Section. 5 of the above paper. Each point needs $n$-number to represent its value. With the degree-9 rule, the points in each region is $C=[153, 273, 453, 717, 1105, 1689, 2605, 4117, 6745]$ from dimension 4 to 12. Thus, we can approximate the memory in GB by $C_{n+r}^r * (100*1.3^n) * C(n) * n / (1024^3) * 4$. \n\nThus, by approximation, we find that for $n = 6, r = 10$, the memory consumption is 39GB for the input using Float32. The dimension of the moment is approximately $8000$. The $r=10$ can represent highly nonlinear distributions, which can represent highly complicated distributions. \n\nFor $n=8, r=4$, the memory consumption is $13.3$ GB. While for $n=8, r=6$, the memory becomes $80.7 GB$. \n\nThus, with full moment space, we may expect dimensions up to **n=8** with **r=4** to be deployed on an Nvidia Orin for real-time computation. Then the actual time depends on the FLOPS of the GPU. \n\nWe note that there are possible sparse polynomial representations that may help to further lower the memory consumption, which can be interesting future directions. See:\n\nVictor Magron and Jie Wang, Sparse Polynomial Optimization: Theory and Practice, World Scientific Press, 2023\n\n**Q3: Dimension:**\n\nThe ability to recover the distribution does not change w.r.t the dimensionality. The programming $P_r$ provides a monotonically decreasing upper bound of the entropy given moment constraints with larger $r$. Such a statement does not change w.r.t the system dimension. As shown in the previous statement, the price to integrate the function increases exponentially, which is expected according to the properties of numerical integration."}, {"review_id": "IGOk7Xesdi", "reviewer": "Reviewer_Ana8", "summary": "This paper proposes the Max Entropy Moment Kalman Filter (MEM-KF) for dynamical systems with nonlinear dynamics and non-Gaussian noise.  The MEM-KF is focused on polynomial systems and can asymptotically approximate almost any distribution through a countable set of moment constraints.  Analogous to a Kalman Filter, the paper proposes prediction and measurement update steps whereby moments are propogated through the process model to recover a maximum entropy distribution (MED).  All steps in the MEM-KF can be solved via convex optimization.", "questions": "The authors claim (or at least allude to) in several places that MEM-KF is an optimal filter.  For example, the abstract begins as \"Designing optimal Bayes filters...\" and later in the introduction \"Thus, a natural question is, how can we provide provably optimal state estimation...\"  It is unclear in what context MEM-KF is considered an optimal filter.  Optimal in what sense?  The discussion on line3s 119-121 suggest that the SDP relaxation is not always exact, and thus this would violate any notion of optimality, no?", "limitations": "I would like to see more honesty in limitations on runtime execution.  The authors state that gradient calculations scale exponentially w.r.t. the dimension of the state-space in the \"Limitations and Future Work\" section, but this point is not raised in the main text anywhere.  This seems like a significant limitation.", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 4, "significance": 3, "originality": 3, "strengths_and_weaknesses": "**Strengths**\nVery well-written paper and clear presentation that studies a fundamental problem in dynamical systems.  The ideas appear novel and are carefully-presented.  Contributions are adequately contextualized w.r.t. most existing work (modulo connections to the assumed density filter (ADF) and expectation propagation (EP)).  Overall, a strong submission.\n\n**Weaknesses**\nThe most significant weakness of this paper is in the numerical evaluations.  Filtering for nonlinear / non-Gaussian systems has enormous application potential, and I would have expected to see more interesting and compelling experiments than a few low-dimensional simulated models.  Moreover, as I understand it the authors only provide one experiment that involves the full filter.  Indeed, experiments in Sec. 7.2 and 7.3 appear only to explore the density estimation problem for a single timestep.\n\nI would argue that the primary comparison to these systems should be to the particle filter (PF), which the authors have provided but have removed to the appendix.  Compared to the PF the MEM-KF does about what one would expect, namely it is competitive up to a certain number of particles, at which point the PF provides superior estimates.  Critically, the authors do not provide any runtime statistics, so it is unclear how MEM-KF compares to PF with adequate number of particles in terms of runtime.  My intuition is that computing the higher-order moments required by MEM-KF will scale polynomially in the moment order, which could be prohibitive.  Worse yet the authors suggest that the gradient calculations in Eqs. (4) and (5) scale exponentially w.r.t. the dimension of the systems, likely why they only consider low-dimensional systems.\n\nSome detailed comments:\n* L101: \"that sums up to $r$\" --> \"that sums up to no more than $r$\"\n* L103: $\\gamma_i$ undefined (i.e. with the index)\n* L107: The definition of $M_r$ is recursive (it appears on both sides) and so it lacks an explicit definition\n* L107-108: The authors refer to $\\bar{x}_\\alpha$ as a moment, but the integral definition isn't a moment (it lacks a probability measure)\n* L119-121: This discussion is a bit confusing, the SDP relaxation is not exact but the authors argue that all optimization is convex?  Please clarify.\n* Eq. (1) : Shouldn't the measurement model be $y_k = h(x_k,v_k)$?\n* L136: \"Assume that state at time $k+1$ is solely dependent on...measurement $y_k$...\" The state process evolves independent of the measurement, perhaps the authors mean the posterior belief.\n* L154: $p_t$ undefined (with the subscript)\n* L162: Please provide the proof in the appendix for completeness.\n* Eq. (3): Shouldn't this be up to proportionality?  Where is the normalizer?\n* L210,212,218: All references to $p(x_{k+1})$ should be replaced by $p^-(x_{k+1})$. \n* First line in the for-loop, Algorithm 1: $x_{\\beta,k-1}$ should be $\\bar{x}_{\\beta,k-1}$ \n* L317: References Figure 8, I think this should reference Figure 4?", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the appreciation and great suggestions! \n\n**The optimality:**\n\nThe optimality is in the sense of asymptotic approximation and in the computation. For the former one, the MEM-KF is guaranteed to converge to any moment determinate distribution with the upper bound of entropy minimized (with additional compact support conditions in engineering problems as suggested by the reviewer VyDk to ensure the true distribution is indeed moment determinate). Thus, the MEM-KF is capable of asymptotically approximating any distribution on a compact domain without overconfidence. \n\nFor the main iterations, all steps of MEM-KF can be **exactly** solved by convex optimization with guaranteed optimality. \n\nThe extraction of point estimation by SDP relaxation is off the main loop and does not violate the belief update and propagation. The SDP is provably tight for estimation with noise in moderate magnitude and can always provide a lower bound approximation. Thus, the extraction of point estimation has guaranteed optimality either as an exact solution or a non-trivial lower bound. \n\nSee: Diego Cifuentes, Sameer Agarwal, Pablo A Parrilo, and Rekha R Thomas. On the local stability of semidefinite relaxations. Mathematical Programming, pages 1–35, 2020.\n\nWe will make this point clearer in the revised version and highlight the **asymptotical optimality** specifically. \n\n**Runtime**\n\nWe note that all the UKF, EKF, UKFM, and InEKF can be run in real-time in the cases considered in this work, both prediciton and update can be completed less than $1e-4s$. For the PF with $1e5$ particles, the prediction and update steps with resampling and CPU parallelization requires less than $0.2s$. \n\nFor the MEM-KF in the localization example, the mean (std.dev) time (secs) is $6.0165e^{-4}$ $(2.0275e^{-4})$ in the prediction, $6.8753$ ($8.7022$) in the update, and $34.5102 $ ($26.2063$) in recovering the MED from moments. The Laptop is equipped with Intel i7-11850H CPU. The convex optimization is run until the KKT violation is below $1e-4$. \n\nWe will add the table in the revision. \n\n\n**The computational limits:**\n\nWe show some theoretical analysis of the scalability of the proposed method. The complexity is mainly determined by the dimension of the system $n$ and the order of moments $r$. The number of the monomials for such system is $$s(r, n):= |\\mathbb{N}_r^n| = \\begin{pmatrix}\n    n+r \\\\ \n    n\n\\end{pmatrix} = \\frac{(n+r)!}{n!r!}.$$ \nFor fixed $r$, we have $s(r, n) \\sim \\mathcal{O}(n^r)$ by eliminating the $n!$ in the denominator: $$ s(r, n) = \\frac{(n+1)(n+2)\\cdots(n+r)}{r!}, \\frac{(n+r)^r}{r!} \\ge s(r, n) \\ge \\frac{n^r}{r!}.$$\nFor fixed $n$, we have $s(r, n) \\sim \\mathcal{O}(r^n) $ by eliminating the $r!$ in the denominator: $$ s(r, n) = \\frac{(r+1)(r+2)\\cdots(r+n)}{n!},  \\frac{(r+n)^n}{n!} \\ge s(r, n) \\ge \\frac{(r+1)^n}{n!}.$$\nThus, the terms to evaluate scales polynomially w.r.t the order and exponentially w.r.t the dimension. \n\nFor the numerical integration, the number of quadrature scales exponentially w.r.t. $n$. Say that the number of quadratures is $c^n$.  \n\nCompute the moments: Requires evaluating $s(r, n)$ terms of moments by numerical integration. These require $s(r, n) \\times c^n $ operations. \n\nPropagation of moments: We need to compute $x_{\\alpha, k+1}=Ax_{\\alpha, k}$,  which requires $\\mathcal{O}(s(r, n)^2)$ multiplications and $\\mathcal{O}(s(r, n)^2)$ additions. To see this, we find each entry of $\\bar{x}_{\\alpha, k+1}$ needs $s(r, n)$ of multiplication and add them together.\n\nReconstruction of distribution: We need to run the convex optimization, that each gradient descent to compute the moments of $\\phi(x)$, which involves numerical integrations of $s(r, n)$ terms of the $ n$-dimensional system.  \n\nUpdate with measurements: This step requires $s(r, n)$ operation of addition to update the coefficients given measurements and prior belief.\n\nNormalization: This step requires integrating the unnormalized belief once. \n\nIn summary, the overall complexity scale exponentially w.r.t the dimension $n$ and polynomially w.r.t the order $r$. \n\nWe will include these discussions in the revised version. \n\n**minor issues:**\n\nL101: We will make this change! \n\nL103: It's the coefficients, we will make this change! \n\nL107: The LHS takes the probability moments as the arguments, while the RHS is the polynomial variables.\n\nL107-108: We will correct this point!\n\nL119-121: The (POP) is not convex. For estimation problems, the SDP relaxation usually give a tight relaxation, thus exact. The rare cases is when the outliers are too big. Such extraction is not required by MEM-KF iterations. To model the observation as MED, we need to separate the noise to one side of the measurement model. \n\nL136: We will correct this point! The dynamics do not depend on the measurement.\n\nL154: We will make this point clear in the revised version. \n\nL162: We will try to adapt the proof in the revised version. The proof relying on analysis is very long, as shown in [39]. The normalizer is the constant lumped in side of the $\\exp(\\cdot)$ and optimized in the convex optimization. \n\nL210, 212, 218: We will make this change! \n\nL317: We will make this change!\n\n**Format**\n\nWe will make sure the format is correct in the camera-ready version!"}]

Author rebuttal:

Meta-review: Submission 11180 proposes a novel Max Entropy Moment Kalman Filter (MEM-KF), to tackle the failure of Kalman Filter (KF) in non-linear, non-Gaussian state space model settings. The MEM-KF uses a Moment-Constrained Max-Entropy Distribution (MED) to model states and observations. The authors claim that the MED can asymptotically approximate almost any distribution given an increasing number of moment constraints. They also show how to use MED to solve the predict- and update- steps of a recursive Bayes filter.

In particular, it is possible to highlight:

- Theoretical contributions seem to be grounded and well-supported by rigorous mathematical proofs.
- This work brings novel contributions that are of interest to the ML community.
- The paper is well-motivated and well-structured. The writing is generally clear and accessible, making the technical content relatively easy to follow and understand.

Despite the positive aspects mentioned above, there remain reservations regarding the trade-off between computational complexity, memory usage, and accuracy. While it is understood that submissions may not address every aspect expected by the ML community, these limitations need to be further discussed in the paper. The authors have already discussed these points during the reviewer-author phase, emphasizing that their model could potentially benefit from GPU acceleration and other power-up computational techniques. The reviewers have globally welcomed the authors' responses and further discussions.

For the aforementioned reasons, the decision for this submission is an acceptance. However, the authors need to properly prepare the camera-ready version while discussing the computational limitations and other remarks acknowledged in the "Author Final Remarks by Authors".
Decision: Accept (poster)