Simon Lacoste-Julien

Fredrik Lindsten

Francis Bach

Sequential Kernel Herding: Frank-Wolfe Optimization for Particle Filtering

Recently, the Frank-Wolfe optimization algorithm was suggested as a procedure to obtain adaptive quadrature rules for integrals of functions in a reproducing kernel Hilbert space (RKHS) with a potentially faster rate of convergence than Monte Carlo integration (and "kernel herding" was shown to be a special case of this procedure). In this paper, we propose to replace the random sampling step in a particle filter by Frank-Wolfe optimization. By optimizing the position of the particles, we can obtain better accuracy than random or quasi-Monte Carlo sampling. In applications where the evaluation of the emission probabilities is expensive (such as in robot localization), the additional computational cost to generate the particles through optimization can be justified. Experiments on standard synthetic examples as well as on a robot localization task indicate indeed an improvement of accuracy over random and quasi-Monte Carlo sampling.

Introduction

In this paper, we explore a way to combine ideas from optimization with sampling to get better approximations in probabilistic models. We consider state-space models (SSMs, also referred to as general state-space hidden Markov models), as they constitute an important class of models in engineering, econometrics and other areas involving time series and dynamical systems. A discrete-time, nonlinear SSM can be written as

x t | x 1:(t-1) ∼ p(x t |x t-1 ); y t | x 1:t ∼ p(y t |x t ), (1) where x t ∈ X denotes the latent state variable and y t ∈ Y the observation at time t. Exact state inference in SSMs is possible, essentially, only when the model is linear and Gaussian or when the state-space X is a finite set. For solving the inference problem beyond these restricted model classes, sequential Monte Carlo methods, i.e. particle filters (PFs), have emerged as a key tool; see e.g., [START_REF] Doucet | A tutorial on particle filtering and smoothing: Fifteen years later[END_REF]; [START_REF] Cappé | Inference in Hidden Markov Models[END_REF]; [START_REF] Doucet | On sequential Monte Carlo sampling methods for Bayesian filtering[END_REF]. However, since these methods are based on Monte Carlo integration they are inherently affected by sampling variance, which can degrade the performance of the estimators.

Particular challenges arise in the case when the observation likelihood p(y t |x t ) is computationally expensive to evaluate. For instance, this is common in robotics applications where the observation model relates the sensory input of the robot, which can comprise vision-based systems, laser rangefinders, synthetic aperture radars, etc. For such systems, simply evaluating the observation function for a fixed value of x t can therefore involve computationally expensive operations, such as image processing, point-set registration, and related tasks. This poses difficulties for particle-filtering-based solutions for two reasons: (1) the computational bottleneck arising from the likelihood evaluation implies that we cannot simply increase the number of particles to improve the accuracy, and (2) this type of "complicated" observation models will typically not allow for adaptation of the proposal distribution used within the filter, in the spirit of [START_REF] Pitt | Filtering via simulation: Auxiliary particle filters[END_REF], leaving us with the standardbut inefficient-bootstrap proposal as the only viable option. On the contrary, for these systems, the dynamical model p(x t |x t-1 ) is often comparatively simple, e.g. being a linear and Gaussian "nearly constant acceleration" model [START_REF] Ristic | Beyond the Kalman filter: particle filters for tracking applications[END_REF].

The method developed in this paper is geared toward this class of filtering problems. The basic idea is that, in scenarios when the likelihood evaluation is the computational bottleneck, we can afford to spend additional computations to improve upon the sampling of the particles. By doing so, we can avoid excessive variance arising from simple Monte Carlo sampling from the bootstrap proposal.

Contributions. We build on the optimization view from [START_REF] Bach | On the equivalence between herding and conditional gradient algorithms[END_REF] of kernel herding [START_REF] Chen | Supersamples from kernel herding[END_REF] to approximate the integrals appearing in the Bayesian filtering recursions. We make use of the Frank-Wolfe (FW) quadrature to approximate, in particular, mixtures of Gaussians which often arise in a particle filtering context as the mixture over past particles in the distribution over the next state. We use this approach within a filtering framework and prove theoretical convergence results for the resulting method, denoted as sequential kernel herding (SKH), giving one of the first explicit better convergence rates than for a particle filter. Our preliminary experiments show that SKH can give better accuracy than a standard particle filter or a quasi-Monte Carlo particle filter.

2 Adaptive quadrature rules with Frank-Wolfe optimization 2.1 Approximating the mean element for integration in a RKHS

We consider the problem of approximating integrals of functions belonging to a reproducing kernel Hilbert space (RKHS) H with respect to a fixed distribution p over some set X . We can think of the elements of H as being real-valued functions on X , with pointwise evaluation given from the reproducing property by f (x) = f, Φ(x) , where Φ : X → H is the feature map from the state-space X to the RKHS. Let κ : X 2 → R be the associated positive definite kernel. We briefly review here the setup from [START_REF] Bach | On the equivalence between herding and conditional gradient algorithms[END_REF], which generalized the one from [START_REF] Chen | Supersamples from kernel herding[END_REF]. We want to approximate integrals E p [f ] for f ∈ H using a set of n points x (1) , . . . , x (n) ∈ X associated with positive weights w (1) , . . . , w (n) which sum to 1:

E p [f ] ≈ n i=1 w (i) f (x (i) ) = E p[f ], (2)
where p := n i=1 w (i) δ x (i) is the associated empirical distribution defined by these points and δ x (•) is a point mass distribution at x. If the points x (i) are independent samples from p, then this Monte Carlo estimate (using weights of 1/n) is unbiased with a variance of V p [f ]/n, where V p [f ] is the variance of f with respect to p. By using the fact that f belongs to the RKHS H, we can actually choose a better set of points with lower error. It turns out that the worst-case error of estimators of the form (2) can be analyzed in terms of their approximation distance to the mean element µ(p) := E p [Φ] ∈ H [START_REF] Smola | A Hilbert space embedding for distributions[END_REF][START_REF] Sriperumbudur | Hilbert space embeddings and metrics on probability measures[END_REF]. Essentially, by using Cauchy-Schwartz inequality and the linearity of the expectation operator, we can obtain:

sup f ∈H f H ≤1 |E p [f ] -E p[f ]| = µ(p) -µ(p) H =: MMD(p, p), (3)
and so by bounding MMD(p, p), we can bound the error of approximating the expectation for all f ∈ H, with f H as a proportionality constant. MMD(p, p) is thus a central quantity for developing good quadrature rules given by (2). In the context of RKHSs, MMD(p, q) can be called the maximum mean discrepancy [START_REF] Gretton | A kernel two-sample test[END_REF] between the distributions p and q, and acts a pseudo-metric on the space of distributions on X . If κ is a characteristic kernel (such as the standard RBF kernel), then MMD is in fact a metric, i.e. MMD(p, q) = 0 =⇒ p = q. We refer the reader to [START_REF] Sriperumbudur | Hilbert space embeddings and metrics on probability measures[END_REF] for the regularity conditions needed for the existence of these objects and for more details.

Frank-Wolfe optimization for adaptive quadrature

For getting a good quadrature rule p, our goal is thus to minimize µ(p) -µ(p) H . We note that µ(p) lies in the marginal polytope M ⊂ H, defined as the closure of the convex-hull of Φ(X ). We suppose that Φ(x) is uniformly bounded in the feature space, that is, there is a finite R such that Φ(x) H ≤ R ∀x ∈ X . This means that M is a closed bounded convex subset of H, and we could in theory optimize over it. This insight was used by [START_REF] Bach | On the equivalence between herding and conditional gradient algorithms[END_REF] who considered using the Frank-Wolfe optimization algorithm to optimize the convex function J(g) := 1 2 g -µ(p) 2 H over M to obtain adaptive quadrature rules. The Frank-Wolfe algorithm (also called conditional gradient) [START_REF] Frank | An algorithm for quadratic programming[END_REF]) is a simple first-order iterative constrained optimization algorithm for optimizing smooth functions over closed bounded convex sets like M (see [START_REF] Dunn | Convergence rates for conditional gradient sequences generated by implicit step length rules[END_REF] for its convergence analysis on general infinite dimensional Banach spaces). At every iteration, the algorithm finds a good feasible search vertex of M by minimizing the linearization of J at the current iterate g k : ḡk+1 = arg min g∈M J (g k ), g . The next iterate is then obtained by a suitable convex combination of the search vertex ḡk+1 and the previous iterate g k : g k+1 = (1 -γ k )g k + γ k ḡk+1 for a suitable step-size γ k from a fixed schedule (e.g. 1/(k + 1)) or by using linesearch. A crucial property of this algorithm is that the iterate g k is thus a convex combination of the vertices of M visited so far. This provides a sparse expansion for the iterate, and makes the algorithm suitable to high-dimensional optimization (or even infinite)this explains in part the regain of interest in machine learning in the last decade for this old optimization algorithm (see [START_REF] Jaggi | Revisiting Frank-Wolfe: Projection-free sparse convex optimization[END_REF] for a recent survey). In our setup where M is the convex hull of Φ(X ), the vertices of M are thus of the form ḡk+1 = Φ(x (k+1) ) for some x (k+1) ∈ X . Running Frank-Wolfe on M thus yields

g k = k i=1 w (i) k Φ(x (i) ) = E p[Φ]
for some weighted set of points {w

(i) k , x (i) } k i=1 .
The iterate g k thus corresponds to a quadrature rule p of the form of (2) and

g k = E p[Φ],
and this is the relationship that was explored in [START_REF] Bach | On the equivalence between herding and conditional gradient algorithms[END_REF]. Running Frank-Wolfe optimization with the step-size of γ k = 1/(k + 1) reduces to the kernel herding algorithm proposed by [START_REF] Chen | Supersamples from kernel herding[END_REF]. See also [START_REF] Huszár | Optimally-weighted herding is Bayesian quadrature[END_REF] for an alternative approach with negative weights.

Algorithm 1 presents the Frank-Wolfe optimization algorithm to solve min g∈M J(g) in the context of getting quadrature rules (we also introduce the shorthand notation µ p := µ(p)). We note that to evaluate the quality MMD(p, p) of this adaptive quadrature rule, we need to be able to evaluate µ p (x) =

x ∈X p(x )κ(x , x)dx efficiently. This is true only for specific pairs of kernels and distributions, but fortunately this is the case when p is a mixture of Gaussians and κ is a Gaussian kernel. This insight is central to this paper; we explore this case more specifically in Section 2.3. To find the next quadrature point, we also need to (approximately) optimize µ p (x) over X (step 3 of Algorithm 1, called the FW vertex search). In general, this will yield a non-convex optimization problem, and thus cannot be solved with guarantees, even with gradient descent. In our current implementation, we approach step 3 by doing an exhaustive search over M random samples from p precomputed when FW-Quad is called. We thus follow the idea from the kernel herding paper [START_REF] Chen | Supersamples from kernel herding[END_REF] to choose the best N "super-samples" out of a large set of samples M . Thanks to the fact that convergence guarantees for Frank-Wolfe optimization can still be given when using an approximate FW vertex search, we show in Appendix B of the supplementary material that this procedure either adds a O(1/M = 1/k for all i ≤ k. A third alternative is to re-optimize J(g) over the convex hull of the previously visited vertices; this is called the fully corrective version [START_REF] Jaggi | Revisiting Frank-Wolfe: Projection-free sparse convex optimization[END_REF] of the Frank-Wolfe algorithm (hereafter referred as FCFW). In this case: (w (1) k+1 , . . . , w (k+1) k+1 ) = arg min w∈∆ k+1 w K k+1 w -2c k+1 w, where ∆ k+1 is the (k + 1)-dimensional probability simplex, K k+1 is the kernel matrix on the (k+1) vertices: (K k+1 ) ij = κ(x (i) , x (j) ) and (c k+1 ) i = µ p (x (i) ) for i = 1, . . . , (k + 1). This is a convex quadratic problem over the simplex. A slightly modified version of the FCFW is called the min-norm point algorithm and can be more efficiently optimized using specific purpose active-set algorithms -see Bach (2013, §9.2) for more details. We refer the reader to [START_REF] Bach | On the equivalence between herding and conditional gradient algorithms[END_REF] for more details on the rate of convergence of Frank-Wolfe quadrature assuming that the FW vertex is found with guarantees. We summarize them as follows: if H is infinite dimensional, then FW-Quad gives the same O(1/ √ N ) rate for the MMD error as standard random sampling, for all FW methods. On the other hand, if a ball of non-zero radius centered at µ p lies within M, then faster rates than random sampling are possible: FW gives a O(1/N ) rate whereas FW-LS and FCFW gives exponential convergence rates (though in practice, we often see differences not explained by the theory between these methods).

Example: mixture of Gaussians

We describe here in more details the Frank-Wolfe quadrature when p is a mixture of Gaussians p(x)

= K i=1 π i N (x|µ i , Σ i ) for X = R d and κ is the Gaus- sian kernel κ σ (x, x ) := exp(-1 2σ 2 x -x 2 ). In this case, µ p (x) = K i=1 π i ( √ 2πσ) d N (x|µ i , Σ i + σ 2 I d ).
We thus need to optimize a difference of mixture of Gaussian bumps in step 3 of Algorithm 1, a non-convex optimization problem that we approximately solve by exhaustive search over M random samples from p.

Sequential kernel herding

Sequential Monte Carlo

Consider again the SSM in (1). The joint probability density function for a sequence of latent states x 1:T := (x 1 , . . . , x T ) and observations y 1:T factorizes as p(x 1:T , y 1:T ) = T t=1 p(x t |x t-1 )p(y t |x t ), with p(x 1 |x 0 ) := p(x 1 ) denoting the prior density on the initial state. We would like to do approximate inference in this SSM. In particular, we could be interested in computing the joint filtering distribution r t (x 1:t ) := p(x 1:t |y 1:t ) or the joint predictive distribution p t+1 (x t+1 , x 1:t ) := p(x t+1 , x 1:t |y 1:t ). In parti-Algorithm 1 FW-Quad(p, H, N ): Frank-Wolfe adaptive quadrature Input: distribution p, RKHS H which defines kernel κ(•, •) and state-space X , number of samples N 1: Let g 0 = 0. 2: for k = 0 . . . N -1 do

3: Solve x (k+1) = arg min x∈X g k -µ p , Φ(x)
That is:

x (k+1) = arg min x∈X k i=1 w (i) k (κ(x (i) , x) -µ p (x)).

4:

Option (1): Let γ k = 1 k+1 .

5:

Option (2): Let γ k = g k -µp,g k -Φ(x (k+1) ) g k -Φ(x (k+1) ) 2
(LS)

6:

Update

g k+1 = (1 -γ k )g k + γ k Φ(x (k+1) ) i.e. w (k+1) k+1 = γ k ; and w (i) k+1 = (1 -γ k )w (i) k for i = 1 . . . k 7: end for 8: Return: p = N i=1 w (i) N δ x (i)
(i) t , x (i) 1:t } N
i=1 as in (2). We note that it is easy to marginalize p with a simple weight summation, and so we will present the algorithm as getting an approximation for the joint distributions r t and p t defined above, with the understanding that the marginal ones are easy to obtain afterwards. In the terminology of particle filtering, x (i) t is the particle at time t, whereas x (i) 1:t is the particle trajectory. While principally the PF provides an approximation of the full joint distribution r t (x 1:t ), it is well known that this approximation deteriorates for any marginal of x s for s t [START_REF] Doucet | A tutorial on particle filtering and smoothing: Fifteen years later[END_REF]. Hence, the PF is typically only used to approximate marginals of x s for s t (fixed-lag smoothing) or s = t (filtering), or for prediction.

Algorithm 2 presents the bootstrap particle filtering algorithm [START_REF] Gordon | Novel approach to nonlinear/non-Gaussian Bayesian state estimation[END_REF] from the point of view of propagating an approximate posterior distribution forward in time (see e.g. [START_REF] Fearnhead | Using random quasi-Monte-Carlo within particle filters, with application to financial time series[END_REF]. We describe it as propagating an approximation pt (x 1:t ) of the joint predictive distribution one time step forward with the model dynamics to obtain pt+1 (x t+1 , x 1:t ) (step 5), and then randomly sampling from it (step 3) to get the new predictive approximation pt+1 (x t+1 , x 1:t ). As pt is an empirical distribution, pt+1 is a mixture distribution (the mixture components are coming from the particles at time t):

pt+1 (x t+1 , x 1:t ) = 1 Ŵt N i=1 p(y t |x (i) t )w (i) t mixture weight p(x t+1 |x (i) t ) mixture component δ x (i) 1:t (x 1:t ). (4)
We denote the conditional normalization constant at time t by W t := p(y t |y 1:(t-1) ) and the global normalization constant by Z t := p(y 1:t ) = t u=1 W u . Ŵt is the particle filter approximation to W t and is obtained by summing the un-normalized mixture weights in (4); see step 4 in Algorithm 2. Randomly sampling from ( 4) is equivalent to first sampling a mixture component according to the mixture weight (i.e., choosing a past particle x (i) 1:t to propagate), and then sampling its next extension state x (i) t+1 with probability p(x t+1 |x (i) t ). The standard bootstrap particle filter is thus obtained by maintaining uniform weight for the predictive distribution (w

(i) t = 1
N ) and randomly sampling from (4) to obtain the particles at time t+1. This gives an unbiased estimate of pt+1 : E pt+1 [p t+1 ] = pt+1 . Lower variance estimators can be obtained by using a different resampling mechanism for the particles than this multinomial sampling scheme, such as stratified resampling [START_REF] Carpenter | Improved particle filter for nonlinear problems[END_REF] and are usually used in practice instead.

One way to improve the particle filter is thus to replace the random sampling stage of step 3 with different sampling mechanisms with lower variance or better approximation properties of the distribution pt+1 that we are trying to approximate. As we obtain the normalization constants W t by integrating the observation probability, it seems natural to look for particle point sets with better integration properties. By replacing random sampling with a quasi-random number sequence, we obtain the already proposed sequential quasi-Monte Carlo scheme [START_REF] Philomin | Quasirandom sampling for condensation[END_REF][START_REF] Ormoneit | Lattice particle filters[END_REF][START_REF] Gerber | Sequential quasi-Monte Carlo[END_REF]. The main contribution of our work is to instead propose to use Frank-Wolfe quadrature in step 3 of the particle filter to obtain better (adapted) point sets.

Sequential kernel herding

In the sequential kernel herding (SKH) algorithm, we simply replace step 3 of Algorithm 2 with pt = FW-Quad(p t , H t , N ). As mentioned in the introduction, many dynamical models used in practice assume Gaussian transitions. Therefore, we will put particular emphasis on the case when (more generally) p(x t |x 1:(t-1) , y 1:(t-1) ) is a mixture of Gaussians, with parameters for the mixture components that can be arbitrary functions of the state history x 1:(t-1) , y 1:(t-1) , and is thus still fairly general. We thus consider the Gaussian kernel for the FW-Quad procedure as then we can compute the required quantities analytically. An important subtle point is which Hilbert space H t to consider. In this paper, we focus on the marginalized filtering case, i.e. we are interested in p(x t |y 1:t ) only. Thus we are only interested in functions of x t , which is why we define our kernel at time t to only depend on x t and not the past histories. For simplicity, we also assume that H t = H for all t (we use the same kernel for each time step). Even though the algorithm can maintain the distribution on the whole history pt (x 1:t ), the past histories x 1:(t-1) are marginalized out when computing the mean map, for example µ(p t ) = E pt(x1:t) [Φ(x t )]. During the SKH algorithm, we can still track the particle histories by keeping track from which mixture component in (4) x t was coming from, but the past history is not used in the computation of the kernel and thus does not appear as a repulsion term in step 3 of Algorithm 1. We leave it as future work to analyze what kind of high-dimensional kernel on past histories would make sense in this context, and to analyze its convergence properties. The particle histories are useful in the Rao-Blackwellized extension that we present in Appendix A and use in the robot localization experiment of Section 4.3.

Convergence theory

In this section, we give sufficient conditions to guarantee that SKH is consistent as N goes to infinity. Let p t here denote the marginalized predictive instead of the joint. Let F t be the forward transformation operator on signed measures that takes the predictive distribution p t on x t and yields the unnormalized marginalized predictive distribution F t p t on x t+1 in the SSM. Thus for a measure ν, we get (F t ν)(•) := Xt p(•|x t )p(y t |x t )dν(x t ). We also have that p t+1 =1 Wt F t p t . For the following theorem, F t is a function space on X t+1 defined (depending on H t+1 ) as all functions for which the following semi-norm is finite:

1 f Ft := sup h H t+1 =1 Xt+1 f (x t+1 )h(x t+1 )dx t+1 .
Theorem 1 (Bounded growth of the mean map).

Suppose that the function f t : (x t+1 , x t ) → p(y t |x t )p(x t+1 |x t ) is in the tensor product function space F t ⊗ H t with the following defined nuclear norm:

f t Ft⊗Ht := inf i α i Ft β i Ht
, where the infimum is taken over all the possible expansions such that f t (x t+1 , x t ) = i α i (x t+1 )β i (x t ) for all x t , x t+1 . Then for any finite signed Borel measure ν on X t , we have:

µ(F t ν) Ht+1 ≤ f t Ft⊗Ht µ(ν) Ht .
Theorem 2 (Consistency of SKH). Suppose that for all 1 ≤ t ≤ T , f t is in F t ⊗ H t as defined in Theorem 1 and o t is in H t . Then we have:

2 µ(p T ) -µ(p T ) H T ≤ ˆ T + R o T -1 H T -1 W T -1 + ρ T -1 T -1 t=1 χ t ˆ t T -2 k=t ρ k ,
where

ρ t := ft F t ⊗H t Wt , χ t := t-1 k=1
Ŵk W k and ˆ t is the FW error reported at time t by the algorithm:

ˆ t := µ(p t ) -µ(p t ) Ht .
We note that χ t ≈ 1 as we expect the errors on W k to go in either direction, and thus to cancel each other over time (though in the worst case it could grow exponentially in t). If ˆ t ≤ and ρ t ≤ ρ, we basically have

µ(p T ) -µ(p T ) = O(ρ T ) if ρ > 1; O(T ) if ρ = 1;
and O( ) if ρ < 1 (a contraction). The exponential dependence in T is similar as for a standard particle filter for general distributions; see [START_REF] Douc | Long-term stability of sequential Monte Carlo methods under verifiable conditions[END_REF] though for conditions to get a contraction for the PF.

Importantly, for a fixed T it follows that the rates of convergence for Frank-Wolfe in N translates to rates of errors for integrals of functions in H with respect to the predictive distribution p T . Thus if we suppose that H is finite dimensional, that p t has full support on X for all t and that the kernel κ is continuous, then by Proposition 1 in [START_REF] Bach | On the equivalence between herding and conditional gradient algorithms[END_REF], we have that the faster rates for Frank-Wolfe hold and in particular we could obtain an error bound of O(1/N ) with N particles. As far as we know, this is the first explicit faster rates of convergence as a function of the number Note that the theorem does not depend on how the error of is obtained on the mean maps of the distribution; and so if one could show that a QMC point set could also achieve a faster rate for the error on the mean maps (rather than on the distributions itself as is usually given), then their rates would translate also to the global rate by Theorem 2.4 4 Experiments

Sampling from a mixture of Gaussians

We start by investigating the merits of different sampling schemes for approximating mixtures of Gaussians, since this is an intrinsic step to the SKH al-gorithm. In Figure 1, we give the MMD error as well as the error on the mean function in term of the number of particles N for the different sampling schemes on a randomly chosen mixture of Gaussians with K = 100 components in d = 2 dimensions. Additional results as well as the details of the model are given in Appendix C.1 of the supplementary material. In our experiments, the number of FW search points is M = 50,000. We note that even though in theory all methods should have the same rate of convergence O(1/ √ N ) for the MMD (as H is infinite dimensional), FCFW empirically improves significantly over the other methods. As d increases, the difference between the methods tapers off for a fixed kernel bandwidth σ 2 , but increasing σ 2 gives better results for FW and FCFW than the other schemes.

In the remaining sections, we evaluate empirically the application of kernel herding in a filtering context using the proposed SKH algorithm.

Particle filtering using SKH on synthetic examples

We consider first several synthetic data sets in order to assess the improvements offered by Frank-Wolfe quadrature over standard Monte Carlo and quasi-Monte-Carlo techniques. We generate data from four different systems (further details on the experimental setup can be found in Appendix C.2):

Two linear Gaussian state-space (LGSS) models of dimensions d = 3 and d = 15, respectively.

A jump Markov linear system (JMLS), consisting of 2 interacting LGSS models of dimension d = 2. The switching between the models is governed by a hidden 2-state Markov chain.

A nonlinear benchmark time-series model used by, among others, [START_REF] Doucet | On sequential Monte Carlo sampling methods for Bayesian filtering[END_REF]; [START_REF] Gordon | Novel approach to nonlinear/non-Gaussian Bayesian state estimation[END_REF]. The model is of dimension d = 1 and is given by:

x t+1 = 0.5x t + 25 x t 1 + x 2 t + 8 cos(1.2t) + v t , y t = 0.05x 2 t + e t ,
with v t and e t mutually independent standard Gaussian.

These models are ordered in increasing levels of difficulty for inference. For the LGSS models, the exact filtering distributions can be computed by a Kalman filter. For the JMLS, this is also possible by running a mixture of Kalman filters, albeit at a computational cost of 2 T (where T is the total number of time steps). For the nonlinear system, no closed form expressions are available for the filtering densities; instead we run a PF with N = 100,000 particles as a reference.

We generate 30 batches of observations for T = 100 time steps from all systems, except for the JMLS where we use T = 10 (to allow exact filtering). We run the proposed SKH filter, using both FW and FCFW optimization and compare against a bootstrap PF (using stratified resampling [START_REF] Carpenter | Improved particle filter for nonlinear problems[END_REF]) and a quasi-Monte-Carlo PF based on a Sobol-sequence point-set. All methods are run with N varying from 20 to 200 particles. We deliberately use rather few particles since, as discussed above, we believe that this is the setting when the proposed method can be particularly useful.

To assess the performances of the different methods, we first compute the root-mean-squared errors (RMSE) for the filtered mean-state-estimates over the T time steps, w.r.t. the reference filters. We report the median RMSEs over the 30 different data batches, along with the 25% and 75% quantiles, and the minimum and maximum values in Figure 2. The SKH algorithms were run for three different values of σ 2 ∈ {0.01, 0.1, 1}. Here, we report the results for σ 2 = 1 for the LGSS models and the JMLS, and for σ 2 = 0.1 for the nonlinear benchmark model. The results for the other values are given in Appendix C.2. The improvements are somewhat robust to the value of σ 2 , but in some cases significant differences were observed. As can be seen, both SKH methods improve significantly upon both QMC and the bootstrap PF.

For the two LGSS models, we also compute the MMD (reported in the rightmost column in Figure 2).

Vision-based UAV Localization

In this section, we apply the proposed SKH algorithm to solve a filtering problem in field robotics.

We use the data and the experimental setup described by [START_REF] Törnqvist | Particle filter SLAM with high dimensional vehicle model[END_REF]. The problem consists of estimating the full six-dimensional pose of an unmanned aerial vehicle (UAV).

Törnqvist et al. ( 2009) proposed a vision-based solution, essentially tracking interest points in the camera images over consecutive frames to estimate the egomotion. This information is then fused with the inertial and barometer sensors to estimate the pose of the UAV. The system is modelled on state-space form, with a state vector comprising the position, velocity, acceleration, as well as the orientation and the angular velocity of the UAV. The state is also augmented with sensor biases, resulting in a state dimension of 22. Furthermore, the state is augmented with the threedimensional positions of the interest points that are currently tracked by the vision system; this is a vary-ing number but typically around ten.

To deal with the high-dimensional state-vector, [START_REF] Törnqvist | Particle filter SLAM with high dimensional vehicle model[END_REF] used a Rao-Blackwellized PF (see Appendix A) to solve the filtering problem, marginalizing all but 6 state components (being the pose, i.e., the position and orientation) using a combination of Kalman filters and extended Kalman filters.

The remaining 6 state-variables were tracked using a bootstrap particle filter with N = 200 particles; the strikingly small number of particles owing to the computational complexity of the likelihood evaluation.

For the current experiment, we obtained the code and the flight-test data from [START_REF] Törnqvist | Particle filter SLAM with high dimensional vehicle model[END_REF]. The modularity of our approach allowed us to simply replace the Monte Carlo simulation step within their setup with FW-Quad. We ran SKH-FW with σ 2 = 10 and SKH-FCFW with σ 2 = 0.1, as well as the bootstrap PF used in [START_REF] Törnqvist | Particle filter SLAM with high dimensional vehicle model[END_REF], and a QMC-PF; all methods using N = 50, 100, and 200 particles. We ran all methods 10 times on the same data; the variation in SKH coming from the random search points for the FW procedure, and in QMC for starting the Sobol sequence at different points. For comparison, we ran 10 times a reference PF with N = 100,000 particles and averaged the results. The median position errors for 100 seconds of robot time (there are 20 SSM time steps per second of robot time) are given in Figure 3. The UAV is assumed to start at a known location at time zero, hence, all the errors are zero initially. Note that all methods accumulate errors over time. This is natural, since there is no absolute position reference available (i.e., the filter is unstable) and the objective is basically to keep the error as small as possible for as long time as possible. SKH-FW here gives the overall best results, with significant improvements over the bootstrap PF and the QMC methods for small number of particles. SKH-FW even gives similar errors for the last time step with only N = 200 particles as one of the reference PFs (using N = 100,000 particles). See Appendix C.2.1 for a discussion of the role of σ 2 for FCFW.

Runtimes.

In these experiments, we focused on investigating how optimization could improve the error per particle, as the gain in runtime depends on the exact implementation as well as the likelihood evaluation cost. We note that the FW-Quad algorithm scales as O(N M ) for N samples and M search points when using FW, by updating the objective on the M search points in an online fashion (we also empirically observed this linear scaling in N ). On the other hand, FCFW scales as O(N 2 M ) as the weights on the particles possibly change at each iteration, preventing the same online trick. SKH scales linearly with the number of time steps T (as a standard PF). For the UAV 2009) spent an average of 0.2 s per time step for N = 50 particles (linear in the number of particles as the likelihood evaluation is the bottleneck) on a XEON E5-2620 2.10 GHz PC. The overhead of using our Matlab implementation of FW-Quad with N = 50 is about 0.15 s per time step for FW and 0.3 s for FCFW; and 0.3 s for FW and 1.0 s for FCFW for N = 100 (we used M = 10,000 search points in this experiment). In practice, this means that SKH-FW can be run here with 50 particles in the same time as the standard PF is run with about 90 particles. But as Figure 3 shows, the error for SKH-FW with 50 particles is still much lower than the PF with 200 particles.

Conclusion

We have developed a method for Bayesian filtering problems using a combination of optimization and particle filtering. The method has been demonstrated to provide improved performance over both random sampling and quasi-Monte Carlo methods. The proposed method is modular and it can be used with different types of particle filtering techniques, such as the Rao-Blackwellized particle filter. Further investigating this possibility for other classes of particle filters is a topic for future work. Future work also includes a deeper analysis of the convergence theory for the method in order to develop practical guidelines for the choice of the kernel bandwidth. The next two columns are for the same mixture of Gaussians in higher dimension d = 5 with K = 100 components, but running FW-Quad with σ 2 = 1 (middle column) or σ 2 = 100 (last column). We see that using a higher σ 2 helps significantly in higher dimension. The dashed lines are linear fits with slopes reported next to the axes.

where N would be the error after N steps of a standard (non-approximate) Frank-Wolfe procedure (e.g. O(1/N ), though it could be O(1/N 2 ) if µ(p M ) is in the strict interior of the marginal polytope of X M as we show in Appendix G). Finally, we know that E[δ M ] ≤ R/ √ M , and we could also obtain a high probability bound for it as well using a concentration inequality with triangular arrays. This gives the guarantee for the MMD error of the SKH procedure with M random search points (with a term of O(1/M 1/4 )). Even though the rate is slow in M , the approach is motivated for problems where the bottleneck is the evaluation of the observation probability (which is only evaluated N times per time step) whereas M can be taken to be very large. We also note that if H is finite dimensional and the kernel κ is continuous, then an asymptotically faster rate of O(1/ √ M ) can be shown (see Appendix G), though with a worse constant that makes the comparison for smaller M less clear.

C Additional details on experiments

C.1 Mixture of Gaussians experiment

The parameters for the mixture of Gaussians p(x) = K i=1 π i N (x|µ i , Σ i ) were randomly sampled as follows:

• The means µ i 's are uniformly sampled on [-5, 5] -FCFW always performs best (this was observed similarly in [START_REF] Bach | On the equivalence between herding and conditional gradient algorithms[END_REF] but for other pairs of distribution / kernel).

-As d increases, the difference between the methods tapers off for a fixed σ 2 , but increasing σ 2 gives better results for FW and FCFW than the others (see for example the last column of Figure 4).

-The FW-LS results are identical to FW, and so we have excluded them from the plots for clarity.

-The improvement of QMC over MC decreases as the number of mixture components K increase. FW and FCFW are not affected by K as much.

QMC implementation. To generate quasi-random samples from the mixture of Gaussians, we generate a (d+1)-dimensional Sobol sequence using the Matlab qrandstream function. The last dimension is (naively) used to sample the mixture component by using the inverse transformation method for a discrete random variable.

The first d components are then used to sample from the corresponding multivariate Gaussian by transforming d independent standard normals. We note that Gerber and Chopin (2014) argued on the importance of sorting the discrete mixture components according to their location before choosing them with the standard inverse transformation method (in our naive implementation, the order is arbitrary and arising from how the mixture of Gaussians was stored). They propose a method for this that they called Hilbert sort, for which they could prove nice low-discrepancy properties. This approach might reduce the sensitivity to K of QMC. The worse results of QMC for the UAV experiment in Figure 3 might be explained by our naive implementation.

C.2 Synthetic data examples and additional results

In this section we provide additional details and results for the synthetic data examples. The LGSS models and the modes of the JMLS are generated randomly using the function drss from the Matlab Control Systems Toolbox. The four models that were considered are given by:

LGSS, d = 3 on the form

x t+1 = Ax t + v t , v t ∼ N (0, I),
y t+1 = Cx t + e t , e t ∼ N (0, 0.1) with (A, C) being an observable pair. The system has poles in -0.2825 and -0.3669 ± 0.0379i.

LGSS, d = 15 on the form

x t+1 = Ax t + v t , v t ∼ N (0, I), y t+1 = Cx t + e t , e t ∼ N (0, 0.1)
with (A, C) being an observable pair. The system has poles in 0.2456 ± 0.6594i, 0.4833, 0.3329, 0.0882 ± 0.2512i, -0.1485, -0.8045, -0.4848, -0.5252 ± 0.0368i, -0.6692 ± 0.0612i, -0.6604, and -0.6680.

JMLS on the form

P(r t+1 = |r t = k) = Π k , x t+1 = A rt x t + F rt v t , v t ∼ N (0, I), y t = C rt x t + G rt e t , e t ∼ N (0, 1),
with Π = 0.7 0.3 0.3 0.7 , and the two system modes corresponding to observable systems with poles in -0.4429, 0.0937, and -0.6576, 0.3109, respectively.

Nonlinear benchmark model as described in the main text.

Additional results, for the different values of σ 2 ∈ {0.01, 0.1, 1} are reported in Figures 5678.

Simon Lacoste-Julien, Fredrik Lindsten, Francis Bach

C.2.1 Discussion of role of σ 2 for FCFW

The results in Figure 6 for the nonlinear benchmark show an interesting behavior for FCFW when the kernel bandwidth σ 2 is increasing. In particular, for σ 2 = 1 (rightmost plot), SKH-FCFW obtains the lowest error for all methods (and other σ's) at N = 20, but its error stays constant when increasing the number of particles while the other methods see their error decreasing. This phenomenon needs to be carefully studied further. Our current hypothesis is that when σ 2 is large, FCFW is too effective at myopically optimizing the MMD error for the mixture of Gaussians pt and yields a too small effective sample size (it sets many weights of particles to zero), thus hurting the particle filtering error. When d is small or when σ 2 is large, the Gaussian kernel matrix becomes rank-deficient due to numerical precision; we thus have numerically a finite dimensional H. In the case of the 1d nonlinear model, FCFW sometimes could optimize the MMD error within numerical precision (its square of the order of 1e-16) within 30 particles. FW-Quad would thus output only 30 particles even though we asked it to produce N > 30. This explains why increasing N did not translate in a reduction of filtering errors for SKH-FCFW with σ 2 = 1: the effective number of particles stayed much less than N . SKH-FW did not seem to suffer from this problem. This might partly explain why we were able to use the much bigger σ 2 = 10 for the UAV experiment with good results (Figure 3) whereas we used σ 2 = 0.1 for SKH-FCFW. LGSS (d = 3), using σ 2 = 0.01, σ 2 = 0.1, and σ 2 = 1 (left to right). Note that the MMD definition depends on σ 2 . This is why the MMD curves for PF and QMC are also changing with σ 2 (but their RMSE ones are not).

D Proof sketch for Theorem 1

Proof sketch. We assume that the function f t : (x t+1 , x t ) → p(y t |x t )p(x t+1 |x t ) is in the tensor product F t ⊗ H t , with F t defined as in the statement of the theorem. We consider the nuclear norm [START_REF] Jameson | Summing and nuclear norms in Banach space theory[END_REF]:

f t Ft⊗Ht = inf {αi,βi} ∞ i=1 i α i Ft β i Ht over all possible decompositions {α i , β i } ∞ i=1 of f t such that, for all x t , x t+1 f t (x t+1 , x t ) = i α i (x t+1 )β i (x t ).
In the following, let {α i , β i } ∞ i=1 be such a decomposition for f t . We have

(F t ν)(x t+1 ) = p(x t+1 |x t )p(y t |x t ) ft(xt+1,xt) dν(x t ) ∈ R µ(F t v) = (F t ν)(x t+1 )Φ(x t+1 )dx t+1 ∈ H t+1 . Now we have that µ(F t ν) Ht+1 = sup h H t+1 =1 | h, µ(F t ν) |, so we consider for some h ∈ H t+1 : h, µ(F t ν) = (F t ν)(x t+1 )h(x t+1
)dx t+1 (by linearity and reproducing property)

= f t (x t+1 , x t )dν(x t ) h(x t+1 )dx t+1 = f t (x t+1 , x t )h(x t+1 )dx t+1 dν(x t ) (Fubini's theorem) = i α i (x t+1 )β i (x t ) h(x t+1 )dx t+1 dν(x t ) = i β i (x t )dν(x t ) Eν [βi]
α i (x t+1 )h(x t+1 )dx t+1 (Fubini's theorem).

By (3), we have that |E ν [β i ]| ≤ β i Ht µ(ν) Ht . Thus we have:

µ(F t ν) Ht+1 = sup h H t+1 =1 | h, µ(F t ν) | = sup h H t+1 =1 i E ν [β i ] α i (x t+1 )h(x t+1 )dx t+1 ≤ i E ν [β i ] ≤ βi H t µ(ν) H t sup h H t+1 =1 α i (x t+1 )h(x t+1 )dx t+1 := αi F t ≤ i α i Ft β i Ht µ(ν) Ht .
This inequality was valid for any expansion {α i , β i } ∞ i=1 for f t , and thus we can take the infimum of the upper bound over all possible expansions to get: µ(F t ν) Ht+1 ≤ f t Ft⊗Ht µ(ν) Ht as we wanted to prove.

E Special case for the Gaussian kernel

In this section, we explore what form • F takes for the Gaussian kernel. We then show that f t F ⊗H is finite for a simple one-dimensional linear Gaussian SSM as long as σ is small enough (and thus H is big enough, as the size of H increases when σ decreases for the Gaussian kernel).

For the Gaussian kernel κ(x, y) = exp(-x -y 2 2 /2σ 2 ) =: q(x -y), its Fourier transform is q(ω) = e -ix ω q(x) = (2π) d/2 σ d e -ω 2 σ 2 /2 and

h 2 H = 1 (2π) d | ĥ(ω)| 2 q(ω) dω.
Moreover,

α, h L 2 = 1 (2π) d α, ĥ L 2 = 1 (2π) d ĥ(ω) q(ω) 1/2 q(ω) 1/2 α(ω)dω.
By applying Cauchy-Schwartz on L 2 on the RHS, this leads to

| α, h L 2 | h H 1 (2π) d |α(ω)| 2 q(ω)dω 1/2
. Thus, we can take the function space:

α 2 F = 1 (2π) d |α(ω)| 2 q(ω)dω = σ d (2π) d/2 |α(ω)| 2 e -ω 2 σ 2 /2 dω.
This allows for quite peaky distributions for the dynamics, as Diracs are authorized (with constant Fourier transform).

To compute an upper bound on f t F ⊗H , we simply need to find a decomposition of p(y t |x t )p(x t+1 |x t ) as a sum of terms α i (x t+1 )β i (x t ) and bound the appropriate norms of α i and β i .

We do this for a special case in the following section.

E.1 Bound for one-dimensional Gaussian distribution

We assume that x t+1 ∈ R d and y t ∈ R m , and that

p(x t+1 |x t ) = 1 ( √ 2πτ ) d e -1 2τ 2 xt+1-Axt 2 2 p(y t |x t ) = 1 ( √ 2πυ) m e -1 2υ 2 yt-Bxt 2 2
We only do the proof for d = m = 1 and y t = 0 (constant observations) to make the proof simpler. We conjecture that similar results hold more generally. We use the Mehler formula for w such that 2w 1-w 2 = 1 τ 2 [START_REF] Abramowitz | Handbook of mathematical functions: with formulas, graphs, and mathematical tables[END_REF], where H n is the n th Hermite polynomial (H n (x) := (-1) n e -x 2 d n dx n e -x 2 ):

e 2xyw/(1-w 2 ) = 1 -w 2 e (x 2 +y 2 )w 2 /(1-w 2 ) ∞ n=0 1 n! (w/2) n H n (x)H n (y). Thus p(x t+1 |x t )p(y t |x t ) = 1 ( √ 2πτ ) e -1 2τ 2 x 2 t+1 -1 2τ 2 A 2 x 2 t e (x 2 t+1 +A 2 x 2 t )w 2 /(1-w 2 ) 1 ( √ 2πυ) e -1 2υ 2 B 2 x 2 t 1 -w 2 ∞ n=0 1 n! (w/2) n H n (x t+1 )H n (Ax t ) = 1 -w 2 1 ( √ 2πτ ) 1 ( √ 2πυ) ∞ n=0 1 n! (w/2) n α n (x t+1 )β n (x t ) with, using -1 2τ 2 + w 2 1-w 2 = w 2 -w 1-w 2 = -w 1+w : α n (x t+1 ) = e -x 2 t+1 w/(1+w) H n (x t+1 ) β n (x t ) = e -A 2 x 2 t w/(1+w) e -1 2υ 2 B 2 x 2 t H n (Ax t ).
We thus now need to compute the norms of α n and β n , by first computing the Fourier transform. We use the representation:

H n (x) = n! 2iπ e -t 2 +2xt dt t n+1 ,
integrating over a contour around the origin, which leads to:

αn (ω) = R e -iωx H n (x)e -x 2 w/(1+w) dx = R e -iωx n! 2iπ e -t 2 +2xt e -x 2 w/(1+w) dt t n+1 dx = n! 2iπ e -t 2

R

e -iωx e +2xt e -x 2 w/(1+w) dx dt t n+1 .

We may now use to get, using -1

+ 4(1 + w)/4w = 1/w, αn (ω) = n! 2iπ e -t 2 π(1 + w) w exp (2t -iω) 2 (1 + w) 4w dt t n+1 = n! 2iπ exp - ω 2 (1 + w) 4w π(1 + w) w e t 2 /w exp -iωt(1 + w) w dt t n+1 = n! 2iπ exp - ω 2 (1 + w) 4w π(1 + w) w w -n/2 i n e -t2 exp ω t(1 + w) √ w d t tn+1
using the change of variable t = i √ w t. This leads to αn

(ω) = exp - ω 2 (1 + w) 4w π(1 + w) w w -n/2 i n H n ω(1 + w) 2 √ w |α n (ω)| exp - ω 2 (1 + w) 4w π(1 + w) w w -n/2 C exp ω 2 (1 + w) 2 8w (2 n n! √ π) 1/2 exp - ω 2 (1 -w 2 ) 8w π(1 + w) w w -n/2 C(2 n n! √ π) 1/2 using H n (x) C exp(x 2 /2)(2 n n! √ π) 1/2 , with C = π -1/4 , and -1+w 4w + (1+w) 2 8w = 1+w 4w -1 + 1+w 2 = -1-w 2 8w . We thus have α n 2 F = σ √ 2π |α n (ω)| 2 e -ω 2 σ 2 /2 dω w -n 2 n n! √ πC 2 π(1 + w) w σ √ 2π exp - ω 2 (1 -w 2 ) 4w e -ω 2 σ 2 /2 dω = w -n 2 n n! √ πC 2 π(1 + w) w σ √ 2π √ 2π σ 2 + (1 -w 2 )/2w = (w -n 2 n n!) × C(w, σ) Moreover, βn (ω) = R e -iωx H n (Ax)e -x 2 (A 2 w/(1+w)+B 2 /2υ 2 ) dx = n! 2iπ R e -iωx e -t 2 e +2Axt dt t n+1 e -x 2 (A 2 w/(1+w)+B 2 /2υ 2 ) dx = n! 2iπ e -t 2 R e -iωx e +2Axt e -x 2 (A 2 w/(1+w)+B 2 /2υ 2 ) dx dt t n+1 = n! 2iπ e -t 2 π A 2 w/(1 + w) + B 2 /2υ 2 exp [2At -iω] 2 4(A 2 w/(1 + w) + B 2 /2υ 2 ) dt t n+1 = n! 2iπ exp -ω 2 4(A 2 w/(1 + w) + B 2 /2υ 2 ) π A 2 w/(1 + w) + B 2 /2υ 2 × e -t 2 exp 4A 2 t 2 -4Atiω 4(A 2 w/(1 + w) + B 2 /2υ 2 ) dt t n+1 = n! 2iπ exp -ω 2 4(A 2 w/(1 + w) + B 2 /2υ 2 ) π A 2 w/(1 + w) + B 2 /2υ 2 × exp t 2 A 2 -(A 2 w/(1 + w) + B 2 /2υ 2 ) (A 2 w/(1 + w) + B 2 /2υ 2 ) exp -4Aitω 4(A 2 w/(1 + w) + B 2 /2υ 2 ) dt t n+1 = n! 2iπ exp -ω 2 4(A 2 w/(1 + w) + B 2 /2υ 2 ) π A 2 w/(1 + w) + B 2 /2υ 2 × exp t 2 A 2 /(1 + w) -B 2 /2υ 2 ) (A 2 w/(1 + w) + B 2 /2υ 2 ) exp -4Atiω 4(A 2 w/(1 + w) + B 2 /2υ 2 ) dt t n+1 = n! 2iπ exp -ω 2 4(A 2 w/(1 + w) + B 2 /2υ 2 ) π A 2 w/(1 + w) + B 2 /2υ 2 × exp t 2 1 -B 2 (1 + w)/2υ 2 A 2 ) w + B 2 (1 + w)/2υ 2 A 2 exp -4Atiω 4(A 2 w/(1 + w) + B 2 /2υ 2 ) dt t n+1 .
We can now perform the change of variable t = i t w+B 2 (1+w)/2υ 2 A 2 ) 1-B 2 (1+w)/2υ 2 A 2 = i t√ w, with w > w for B > 0.

This leads to βn

(ω) = n! 2iπ exp -ω 2 4(A 2 w/(1 + w) + B 2 /2υ 2 ) π A 2 w/(1 + w) + B 2 /2υ 2 × w-n/2 exp(-t2 ) exp 4A t√ wω 4(A 2 w/(1 + w) + B 2 /2υ 2 ) d t tn+1 = exp -ω 2 4(A 2 w/(1 + w) + B 2 /2υ 2 ) π A 2 w/(1 + w) + B 2 /2υ 2 × w-n/2 H 2A √ wω 4(A 2 w/(1 + w) + B 2 /2υ 2 ) | βn (ω)| exp -ω 2 4(A 2 w/(1 + w) + B 2 /2υ 2 ) π A 2 w/(1 + w) + B 2 /2υ 2 × w-n/2 (2 n n! √ π) 1/2 C exp 1 2 2A t√ wω 4(A 2 w/(1 + w) + B 2 /2υ 2 ) 2 cst × w-n/2 (2 n n! √ π) 1/2 exp(-(A, w, B)ω 2 ) with (A, w, B) = 1 4(A 2 w/(1 + w) + B 2 /2υ 2 ) - 1 8 w (A 2 w/(1 + w) + B 2 /2υ 2 ) 2 .
We have (A, w, 0) = 1-w 2 8A 2 w . Thus, by continuity if B if small enough, (A, w, B) > 0. Note that when we have B = 0 and A = 1, we recover previous results for α n .

We thus have

β n 2 H = 1 (2π) 3/2 σ | βn (ω)| 2 e ω 2 σ 2 /2 dω (2 n w-n n!)C(A, B, w, σ)
as long as σ 2 < 4 (A, w, B).

Thus, for σ small enough, we have the norm less than a constant times ∞ n=0

(w/ w) n/2 < ∞ since w > w. This shows that ∞ n=0 α n F β n H < ∞ and thus that C t is finite if the linear dependency parameter B and the kernel bandwidth σ 2 are small enough.

F Proof for Theorem 2

Proof We recall here that p t (x t ) = p(x t |y 1:(t-1) ) (we are in the marginalized setting). In the notation of the algorithm, we have pt+1 = 1 Ŵt F t pt . Let q t = p t Z t-1 be the un-normalized marginalized predictive distribution (and similarly, qt = pt Ẑt-1 ). We thus have pt+1 = 1 Ẑt F t qt . We use the metric inequality (as well as the linearity of the MMD in each of its argument as it is related to the RKHS norm; so a scalar multiplication of a distribution can be taken out of the MMD): MMD

MMD(p t+1 , p t+1 ) ≤ MMD(p t+1 , 1 Ẑt F t qt ) (I) FW error := ˆ t+1 + MMD( 1 Ẑt F t qt , 1 Z t F t qt ) (II) Normalization error + 1 Z t MMD(F t qt , F t q t ) (
( 1 Ẑt F t qt , 1 Z t F t qt ) = 1 Ẑt µ(F t qt ) H (A) ≤R 1 Z t |Z t -Ẑt | (B) ≤ ot H MMD(qt,qt) ≤ R o t H Z t MMD(q t , q t ).
For inequality (A), we note that Ft pt Ẑt = pt+1 which is a normalized distribution on x t+1 , this is why µ(p t+1 ) H ≤ R as Φ(x) H ≤ R ∀x ∈ X by assumption. For inequality (B), we have that

|Z t -Ẑt | = |E qt [o t ] -E qt [o t ]| ≤ ||o t || H MMD(q t , qt ) by (3) under the assumption that o t ∈ H.
Finally, the initialization error term (III) can be bounded by using Theorem 1 (with ν = qt -q t ):

MMD(F t qt , F t q t ) ≤ C t MMD(q t , q t ),
where C t := f t F ⊗H .

To control MMD(q t , q t ), we now only work on the un-normalized distributions: MMD(q t , q t ) ≤ MMD(q t , F t-1 qt-1 )

:= t + MMD(F t-1 qt-1 , F t-1 q t-1 ) ≤Ct-1 MMD(qt-1,qt-1) ≤ t u=1 u t-1 k=u C k ,
by repeating the arguments for smaller t's and unrolling the recursion (and recall that And so we can rewrite:

u Z t = ˆ u
Ẑu-1

Z t = ˆ u 1 W t t-1 k=u 1 W k u-1 k=1 Ŵk W k := χu .
We expect χ u = u-1 k=1

Ŵk W k ≈ 1 as the errors on the normalization constants could hopefully go in both direction and thus cancel each other, though in the worst case it could also grow with u. Substituting back in (6), we get what we wanted to prove:

MMD(p t+1 , p t+1 ) ≤ ˆ t+1 + (R o t H + C t ) W t t u=1 χ u ˆ u t-1 k=u C k W k . ( 7
)
Remark 1 (Bound for Ẑt ). For parameter estimation in a HMM, one would also be interested in the quality of approximation for Z t . We note that inequality (B) also gives us a bound on the relative error of our estimate Ẑt for the normalization constant:

|Z t -Ẑt | Z t ≤ o t H W t t u=1 χ u ˆ u t-1 k=u C k W k .
Remark 2 (Bound for joint predictive distribution p J t ). To be more precise, we could have used the notation H t and MMD t to be explicit that the RKHS considered was for functions of x t . For example Theorem 1 really says that MMD t+1 (F t qt , F t q t ) ≤ C t MMD t (q t , q t ). But since H t = H (in the isomorphism sense) for all t, we did not have to worry about this. On the other hand, as H t contains functions of x t only, we have that µ(p t ) is the same whether p t is the marginalized or the joint predictive distributions p J t (as for the joint, the expectation in the mean map definition will marginalize out the variables x 1:(t-1) as they do not appear in H t ). This means that if we consider the joint forward transformation F J t on a joint measures ν J on x 1:t : (F J t ν)(x t+1 , x 1:t ) := p(x t+1 |x t )p(y t |x t )ν J (x 1:t ), i.e. pJ t+1 = 1 Wt F J t p J t (now in the joint sense), then we have µ(F t p t ) = µ(F J t p J t ), and thus Theorem 1 also holds for the joint predictive distribution p J T .

Remark 3 (Bound without χ u ). The disadvantage of the bound (7) is the presence of the quantity χ u for which we did not provide an explicit upper bound (though we would expect it to be close to 1).

Thus if we let v k := k(g k -µ p ), then we get:

v k+1 2 = (ḡ k+1 -µ p ) + v k 2 = ḡk+1 -µ p || 2 + v k 2 + 2 v k , ḡk+1 -µ p ≤ 4R 2 + v k 2 + 2(kδ -r v k ) = v k 2 + 2 v k 2R 2 + kδ v k -r . ( 12
)
The last inequality used the crucial strict interior assumption that yielded (11). Now let C k := 1 r (2R 2 + kδ). Note that C k+1 ≥ C k . We will now proceed to show by induction that v k ≤ C k for k ≥ 1. Note that the bracket in ( 12) is negative if and only if v k ≥ C k (i.e. v k+1 ≤ v k in this case), giving the inspiration for the C k threshold.

First, we have that

v 1 = g 1 -µ p ≤ 2R ≤ 2R R r ≤ C 1 ,
by using the fact that R r ≥ 1 since a ball of radius r fitting in M implies that the maximum norm R of elements in M is at least r. Now suppose that v k ≤ C k , i.e. that v k = αC k for some α ∈ [0, 1]. Then (12) becomes:

v k+1 2 ≤ α 2 C 2 k + 2αC k [ rC k αC k -r] = α 2 C 2 k + 2C k r[1 -α].
The RHS is a convex function of α, and so it is maximized at the boundary of its domain. For α = 0, we get v k+1 2 ≤ 2C k r = 4R 2 + 2kδ. For α = 1, we get v k+1 2 ≤ C 2 k . And thus in general, supposing α ∈ [0, 1], we get that v k+1 2 ≤ max{2C k r, C 2 k } = C 2 k as:

C 2 k = 4R 2 R r 2 + 4kδ R r 2 + k 2 δ r 2 ≥ 2C k r = 4R 2 + 2kδ
using R r ≥ 1. This completes the induction step as this means that v k+1 ≤ C k ≤ C k+1 . Thus we conclude that v k ≤ C k for all k ≥ 1, i.e.

g k -µ p ≤ 1 k 2R 2 r + δ r .
This shows the faster O(1/k) rate ( 9). If we do not have µ p in the strict interior of M, i.e. r = 0, then we can unroll the inequality (12) to get:

v k+1 2 ≤ 4R 2 + 2kδ + v k 2 ≤ k l=1 4R 2 + 2lδ + v 1 2 ≤4R 2 ≤ (k + 1)4R 2 + k(k + 1) 2 2δ.
This thus shows (10):

g k -µ p 2 ≤ 1 k 4R 2 + δ.
This translates to a slower O(1/ √ k) rate on g k -µ p (with √ δ precision), but note that at least it does not have the log(k) factor from the rate by [START_REF] Freund | New analysis and results for the Frank-Wolfe method[END_REF].

Figure 1 :

 1 Figure 1: Top: MMD error for different sampling schemes where p is a mixture of 2d Gaussians with K = 100 components. Bottom: error on the mean estimate for the same mixture. The dashed lines are linear fits with slopes reported next to the axes. of particles than the standard O( 1 √ N ) for Monte Carlo particle filters. In contrast, Gerber and Chopin (2014, Theorem 7) showed a o( 1 √ N ) rate for the randomized version of their SQMC algorithm (note the little-o). 3Note that the theorem does not depend on how the error of is obtained on the mean maps of the distribution; and so if one could show that a QMC point set could also achieve a faster rate for the error on the mean maps (rather than on the distributions itself as is usually given), then their rates would translate also to the global rate by Theorem 2.4

Figure 2 :

 2 Figure 2: RMSEs (left and middle columns) for the four considered models and MMDs (right column) for the two LGSS models.

Figure 3 :

 3 Figure3: Median of position errors over 10 runs for each method. The errors are computed relative to the mean prediction over 10 runs of a PF with 100k particles (the variation of the reference PF is also shown for PF 100k). The error bars represent the [25%, 75%] quantile. The rightmost plot shows the error at the last time step as a function of N . 100 s of robot time represents 2,000 SSM time steps, it does not correspond to computation time.

Figure 4 :

 4 Figure4: Error on the mean function (top row) and MMD error (bottom row) for the mixture of Gaussians experiment. The first column is for K = 20 and d = 2. The next two columns are for the same mixture of Gaussians in higher dimension d = 5 with K = 100 components, but running FW-Quad with σ 2 = 1 (middle column) or σ 2 = 100 (last column). We see that using a higher σ 2 helps significantly in higher dimension. The dashed lines are linear fits with slopes reported next to the axes.

  d . • Σ i = σ 2 i I where σ 2 i is uniformly sampled on [0.1, 4.1].• π i are obtained by normalizing independent uniform random variables.

Figure 4

 4 Figure4present additional results for the mixture of Gaussians experiments. From our experiments, we make the following observations:

Figure 5 :

 5 Figure5: RMSE for JMLS, using σ 2 = 0.01, σ 2 = 0.1, and σ 2 = 1 (left to right).

Figure 6 :

 6 Figure6: RMSE for nonlinear benchmark model, using σ 2 = 0.01, σ 2 = 0.1, and σ 2 = 1 (left to right).

Figure 7 :

 7 Figure7: RMSE (top row) and MMD (bottom row) for LGSS (d = 3), using σ 2 = 0.01, σ 2 = 0.1, and σ 2 = 1 (left to right). Note that the MMD definition depends on σ 2 . This is why the MMD curves for PF and QMC are also changing with σ 2 (but their RMSE ones are not).

Figure 8 :

 8 Figure 8: RMSE (top row) and MMD (bottom row) for LGSS (d = 15), using σ 2 = 0.01, σ 2 = 0.1, and σ 2 = 1 (left to right).

Re

  -ax 2 +bx dx = π a e b 2 /4a

  III) Initialization error . The term (I) is the algorithmic Frank-Wolfe error ˆ t+1 := MMD(p t+1 , pt+1 ). The term (II) is the normalization error which can be bounded as follows:

  t+1 , p t+1 ) ≤ ˆ t+1 + (R o t H + C t )Finally, we transform back the t errors in the algorithmic quantities ˆ t that the FW algorithm measures:

	ˆ t+1 = MMD(p t+1 ,	1 Ŵt	F t pt ) = MMD(	1 Ẑt	qt+1 ,	1 Ŵt	F t	qt Ẑt-1	) =	1 Ẑt	MMD(q t+1 , F t qt ) =	1 Ẑt	t+1 .
													t-1 u=t (•) = 1 by convention).
	Combining the three terms, we thus get:
	MMD(p t u=1	u Z t	t-1 k=u	C k .	(6)

  To get an explicit upper bound for the error, we can repeat a similar argument but always working with the normalized quantities:MMD(p t+1 , p t+1 ) ≤ MMD(p t+1 , MMD(F t pt , F t p t )

							1 Ŵt	F t pt )	+ MMD(	1 Ŵt	F t pt ,	1 W t	F t pt )	+	1 W t	.
											(III) Initialization error
					(I) FW error := ˆ t+1		(II) Normalization error
	The term (II) is the normalization error which can bounded similarly as before as:
	MMD(	1 Ŵt	F t pt ,	1 W t	F t pt ) =	1 Ŵt	µ(F t pt ) H (A) ≤R	1 W t	|W t -Ŵt | (B) ≤ ot H MMD(pt, pt) MMD(p t+1 ≤ R o t H W t t
											ˆ u	ρk ,
											u=1	k=u

t , p t ).

Similarly as before, we also have for (III) that by Theorem 1 (with ν = pt -p t ) that MMD(F t pt , F t p t ) ≤ C t M M D(p t , p t ). Combining the three terms, we get:

MMD(p t+1 , p t+1 ) ≤ ˆ t+1 + R o t H + C t W t MMD(p t , p t ) ≤

Appearing in Proceedings of the 18 th International Conference on Artificial Intelligence and Statistics (AISTATS) 2015, San Diego, CA, USA. JMLR: W&CP volume 38. Copyright 2015 by the authors.

In general, the integral on Xt+1 should be with respect to the base measure for which the conditional density p(xt+1|xt) is defined. All proofs are in the supplementary material.

 2 We use the convention that the empty sum is 0 and the empty product is 1.

The rate holds on the approximation of integrals of continuous bounded functions.

We also note that a simple computation shows that for a Monte Carlo sample of sizeN , E µ(p) -µ(p) 2 H ≤ (R 2 -µ(p) 2 ) N.

We note that as XM is finite, we do not need to make the additional assumption that the kernel κ is continuous unlike in Proposition 1 of[START_REF] Bach | On the equivalence between herding and conditional gradient algorithms[END_REF].

Acknowledgements

We thank Eric Moulines for useful discussions. This work was partially supported by the MSR-Inria Joint Centre, a grant by the European Research Council (SIERRA project 239993) and by the Swedish Research Council (project Learning of complex dynamical systems number 637-2014-466).

Supplementary material A Extension for Rao-Blackwellization

A common strategy for improving the efficiency of the PF is to make use of Rao-Blackwellization-this idea can be used also with SKH. Rao-Blackwellization, here, refers to analytically marginalizing some conditionally tractable component of the state vector and thereby reducing the dimensionality of the space on which the PF operates. Assume that the state of the system is comprised of two components x t and z t , where the filtering density for z t is tractable conditionally on the history of x 1:t . The typical case is that of a conditionally linear Gaussian system, in which case the aforementioned conditional filtering density p(z t |x 1:t , y 1:t ) is Gaussian and computable using a Kalman filter (conditionally on x 1:t ). The Rao-Blackwellized PF (RBPF) exploits this property by factorizing:

where the conditional mean z t (x 1:t ) := E[z t |x 1:t , y 1:t ] and covariance matrix Σ t (x 1:t ) := V(z t |x 1:t , y 1:t ) can be computed (for a fixed trajectory x 1:t ) using a Kalman filter. The mixture approximation follows by plugging in a particle approximation of p(x 1:t |y 1:t ) computed using a standard PF. Hence, for a conditionally linear Gaussian model, the RBPF takes the form of a Mixture Kalman filter; see [START_REF] Chen | Mixture Kalman filters[END_REF]. Analogously to a standard PF, the SKH procedure allows us to to compute an empirical point-mass approximation of p(x 1:t |y 1:t ) by keeping track of the complete history of the state x 1:t . Consequently, by (5) it is straightforward to employ Rao-Blackwellization also for SKH; we use this approach in the numerical example in Section 4.3.

B Rates for SKH when using random search points

In this section, we show that we can get guarantees on the MMD error of the FW-Quad procedure when approximately finding the FW vertex in step 3 of Algorithm 1 using exhaustive search through M random samples from p. This means that despite not solving step 3 exactly, the SKH procedure with M random search points (under assumptions of Theorem 2) is still consistent as long as M grows to infinity.

The main idea is that the rates of convergence for the Frank-Wolfe optimization procedure still holds when the linear subproblem (step 3) is solved within accuracy of δ. More specifically, if we guarantee that the FW vertex ḡk+1 that we use satisfy J (g k ), ḡk+1 ≤ min g∈M J (g k ), g + δ during the algorithm, then the standard O(1/k) rate of convergence for FW carries through but within δ of the optimal objective (i.e. up to J(g * ) + δ). A simple modification of the argument by Jaggi (2013) (who used a shrinking δ during the FW algorithm) can show this for the step-size of γ k = 2 k+2 ; we give the proofs for the step-size of γ k = 1 k+1 as well as the potential faster rate O(1/k 2 ) for the MMD objective in Appendix G.

Let X M ⊆ X be the set of M search points, and p M be the empirical distribution for the M samples from p. Let δ M := µ(p M ) -µ(p) H which can be made small by increasing M . Consider the iteration k in FW-Quad where we do exhaustive search on X M in step 3. We thus have:

We can thus interpret step 3 as approximately solving (within δ M R M ) the linear subproblem for the Frank-Wolfe optimization of J M (g) := 1 2 g -µ(p M ) 2 H over the marginal polytope of X M . We thus get a rate of convergence to within δ M R M of min g J M (g) = 0. Finally, we have

where ρt := R ot H +Ct Wt , by unrolling the recursion for smaller t's.

The problem with bound (8) is that ρ > 1 usually due to the extra term R o t in its definition, which is why we preferred the tighter form (7).

Remark 4 (Removing the o t ∈ H condition in Theorem 2). We note that the condition o t ∈ H is not really necessary in Theorem 2. If o t / ∈ H, we can instead re-derive a similar argument as above but using

, where 1 is the constant unit function (here on x t+1 ). We then have |Z t -Ẑt | ≤ 1 H MMD (F t q t , F t qt ), where H is an augmented RKHS to ensure that it contains the constant function 1. We define H = H if 1 ∈ H already. If 1 / ∈ H, we define H to be the Hilbert sum of the RKHS H and the one generated by the constant kernel 1 (and thus H is a RKHS with kernel κ = 1 + κ where κ is the original kernel for H; see Berlinet and Thomas-Agnan (2004, Thm. 5)). We can show that running the Frank-Wolfe algorithm using the kernel κ yields exactly the same objective values and updates, and thus we can use the space H to analyze its behavior: Theorem 1 and an analog of Theorem 2 then hold, but with all norms defined with respect to H instead.

G Faster rates for FW with approximate vertex search for the MMD objective

In this section, we provide the proofs for the rate of convergence for FW on the MMD objective J(g) :=

H when an approximate vertex search is used (as mentioned in Appendix B) by extending the proofs from [START_REF] Chen | Supersamples from kernel herding[END_REF]. We consider the step-size γ k = 1 k+1 . 5 We note that the standard step-size for Frank-Wolfe optimization to get a

The best rate known for general objectives when using FW with γ k = 1 k+1 is actually O(log(k)/k) (Freund and Grigas, 2013, Bound 3.2). We make use of the specific form of the MMD objective here to prove the O(1/k) rate, as well as the faster O(1/k 2 ) rate under additional assumptions. For the rest of this section, we use • to mean • H .

Theorem G.1 (Rates for FW-Quad with approximate vertex search). Consider the FW-Quad Algorithm 1 where an approximate vertex search is used: g k -µ p , ḡk+1 ≤ min g∈M g k -µ p , g + δ, where ḡk+1 := Φ(x k+1 ) and δ ≥ 0. Suppose that µ p lies in the strict interior of M with a radius r > 0, i.e. a ball of radius r centered at µ p lies within M. Recall that max g∈M g ≤ R. Then we have the faster rate O(1/k 2 ) for the objective J:

If r = 0 (note that µ p ∈ M), then we can still get a standard O(1/k) FW rate:

Proof If a ball of a radius r centered at µ p lies within M, then we have that:

So the approximate vertex search yields ḡk+1 with the property:

By using the FW update g k+1 = γ k ḡk+1 + (1 -γ k )g k with the γ k = 1 k+1 step-size, we get:

5 We note that the rate extends to the line-search step-size as well as the improvement at each iteration can only be better in this case considering the proof technique that we use. 6 We also tried the γ k = 2 k+2 step-size in the mixture of Gaussians experiment of Section 2.3, but it gave similar results as the γ k = 1 k+1 step-size.

Consequence for SKH with random search points. Going back over the argument from Appendix B, we said that using M random search points in FW-Quad was similar to approximately solving (within δ M R M ) the linear subproblem for the Frank-Wolfe optimization of J M (g) := 1 2 g -µ(p M ) 2 H over the marginal polytope of X M . We note that the marginal polytope of X M is at most M -dimensional, and thus, by using a similar argument as in Proposition 1 of [START_REF] Bach | On the equivalence between herding and conditional gradient algorithms[END_REF], we could show that it contains a ball of radius r M > 0 centered at µ(p M ). 7 From (9) in Theorem G.1, we can thus conclude that:

This seems to give a faster rate, but the problem is that r M might shrink at an exponential rate with M if H is infinite dimensional. Thus even though δ M is O(1/ √ M ), the dependence of δ M R M r M might be worse than the previously quoted rate of O(1/M 1/4 ); the latter is thus the general worst-case.

On the other hand, under the additional assumption that H is finite dimensional and that there is a ball of radius r > 0 centered around µ p in M (by using Proposition 1 of [START_REF] Bach | On the equivalence between herding and conditional gradient algorithms[END_REF] for example), then for a sufficiently large M , we will have r M close to r. Thus for large M , we will have . This gives an asymptotically faster rate than the O R 1

rate given in Appendix B arising from (10), but the constant is worse by a factor of R r .

Supplementary References