Théophile Cantelobre

Carlo Ciliberto

Benjamin Guedj

Alessandro Rudi

Closed-form Filtering for Non-linear Systems

Sequential Bayesian Filtering aims to estimate the current state distribution of a Hidden Markov Model, given the past observations. The problem is well-known to be intractable for most application domains, except in notable cases such as the tabular setting or for linear dynamical systems with gaussian noise. In this work, we propose a new class of filters based on Gaussian PSD Models, which offer several advantages in terms of density approximation and computational efficiency. We show that filtering can be efficiently performed in closed form when transitions and observations are Gaussian PSD Models. When the transition and observations are approximated by Gaussian PSD Models, we show that our proposed estimator enjoys strong theoretical guarantees, with estimation error that depends on the quality of the approximation and is adaptive to the regularity of the transition probabilities. In particular, we identify regimes in which our proposed filter attains a TV ϵ-error with memory and computational complexity of O(ϵ -1 ) and O(ϵ -3/2 ) respectively, including the offline learning step, in contrast to the O(ϵ -2 ) complexity of sampling methods such as particle filtering.

Introduction

Sequential Bayesian Filtering is the task of inferring the distribution of unobserved variable X T from observations Y 1 , . . . , Y T where (X t , Y t ) n≥0 is a Hidden Markov Model. The distribution of X T given Y 1 , . . . , Y T is known as the filtering distribution (or optimal filter), denoted π ν T where ν ∈ P(E) is the distribution of X 0 (or prior knowledge on X 0 more generally). In a Hidden Markov Model, the conditional distributions of X n given X n-1 and Y n given X n are described by transition kernels Q and G respectively. In particular, Q(u, x) identifies the conditional probability of transitioning to the state X n = x at time n given the fact that the system was in the state X n-1 = u at time n -1, while G(x, y) is the probability of observing Y n = y given the fact that the system is in the state X n = x. If one has perfect knowledge of the initial distribution π 0 , the transition and observation kernels Q, G, then the filtering distribution can be computed recursively by applying Bayes rule:

π n (dx) = Q(u, dx)G(x, y n )π n-1 (du) Q(u, dx)G(x, y n )π n-1 (du) (1)
The recursive expression of the filtering distribution can be seen as a combination of two steps: the prediction of X n+1 given belief on X n given past observations, then the correction of this prediction based on the observation Y n+1 received. For complete introduction to Hidden Markov Models and inference, we refer the reader to [START_REF] Cappé | Inference in Hidden Markov Models[END_REF]. Two quintessential problems in filtering are studying the stability and robustness of the optimal filter, and computing an approximation of it in practice.

The stability of the optimal filter is related the robustness of the sequence π ν n with respect to the initial ν, since this distribution is unknown in most application settings and estimates/priors need to be used in practice. The optimal filter is stable when its dependence on the initial distribution decreases as more observations are provided (also known as "forgetting"). More formally, a filter is stable when the distance between π ν n and π µ n goes to zero as n → ∞ for any two initial distributions ν and µ. This problem has attracted considerable interest since the first contributions of [START_REF] Ocone | Asymptotic stability of the optimal filter with respect to its initial condition[END_REF] and [START_REF] Blankenship | Asymptotic behavior of nonlinear filters[END_REF]. Some reference relevant to this work include [START_REF] François | Stability and uniform approximation of nonlinear filters using the Hilbert metric and application to particle filters[END_REF]; [START_REF] Mcdonald | Exponential filter stability via Dobrushin's coefficient[END_REF]; [START_REF] Legland | Stability and approximation of nonlinear filters: an information theoretic approach[END_REF][START_REF] De Castro | Consistent estimation of the filtering and marginal smoothing distributions in nonparametric hidden markov models[END_REF]; [START_REF] Yu | Sensitivity of hidden markov models[END_REF]. A modern review of the literature and different approaches can be found in Kim (2022).

In general, the iteration in (1) is intractable. Two exceptions are when the state-space in finite or when the state-space is continuous but the transition kernels are Gaussian Linear Conditional Distributions. In the former case, the algorithm is known as the forward algorithm. In the latter, the algorithm is the well-known Kalman filter [START_REF] Kálmán | New results in linear filtering and prediction theory[END_REF]. The Kalman filter is known to compute the recursion (1) exactly where the Hidden Markov Model has linear dynamics and observations and independent Gaussian noise. If the state-space is finite, the Baum-Welch algorithm can compute filtering and smoothing distributions in closed-form using the forward-backward approach.

But many real-world systems do not have linear dynamics nor Gaussian uncertainty. On one hand, many algorithms have been devised to handle variations on the Kalman filter's assumptions, including the Extended Kalman Filter or the Unscented Filter [START_REF] Julier | Unscented filtering and nonlinear estimation[END_REF]) (see also [START_REF] Särkkä | Bayesian Filtering and Smoothing[END_REF]). These models approximate the state variable as a Gaussian, which excludes many systems where, for instance, multi-modality is present. On the other, Sequential Monte Carlo algorithms such as the Particle filter and variants were developed. These algorithms approximate the marginalization step above using sampling, and can handle multi-modality. This family of methods has strong theoretical guarantees (though under arguably stringent conditions), but do not give a closed-form expression of the approximate distribution and are known to be difficult to turn and prohibitively costly for online applications.

Gaussian PSD Models were introduced in Rudi and Ciliberto (2021) as models for probability distributions. Gaussian PSD Models are a special case of the family of models proposed in Marteau-Ferey et al. (2020). They generalize Gaussian Mixture Models by allowing for negative coefficients in the mixture. As originally highlighted in Rudi and Ciliberto (2021), Gaussian PSD Models enjoy appealing properties for applications involving Bayesian inference, with filtering as a special case: 1) They have optimal approximation guarantees with respect to a large family of probability densities 2) products and marginals of Gaussian PSD models can be efficiently computed in closed-form.

Approach & Contributions

In this paper, we study the problem of performing the iteration in (1) when knowledge of the transition Q and observation G probabilities are unknown and only an approximation in terms of two corresponding Gaussian PSD Models Q and Ĝ is available. We introduce a new algorithm to derive an estimator πn (dx) and study its relation with the true π n (dx). The proposed estimator extends previous filtering strategies, such as the Kalman filter and offers strong theoretical guarantees on a large family of application settings.

Our main contributions are:

1. A novel algorithm to tackle Sequential Bayesian Filtering, which recovers previously proposed estimators and can be applied to any filtering problem where the transition kernels admit a smooth density.

2. We show that the proposed estimator is both stable and robust with respect to a large family of application settings. These theoretical properties are adaptive to the regularity properties of the Hidden Markov Model.

3. The computational and space complexity of the proposed algorithm depends on the regularity of the transition kernels. For very regular kernels (e.g. infinitely differentiable) the algorithm has a computational complexity that is smaller than, for example, particle filtering.

Our paper is organized as follows: in Section 2, we describe Gaussian PSD Models and their properties, in particular their stability with respect to probabilistic operations. In Section 3, we devise an algorithm for learning Gaussian PSD Models from function evaluations and prove that optimal estimation rates are attained for smooth targets. In Section 4, we introduce PSDFilter, an approximate filtering algorithm which plugs Q and Ĝ in the iteration above. We prove this algorithm is robust to the choice of initial distribution and to the approximation error in Q and G. Sketches of the proofs of our main theorems Theorem 4 and Theorem 6 are presented in Section 5. Finally, in Section 6, we generalize Gaussian PSD Models to allow for a richer class of approximators, while retaining most of the desirable properties of Gaussian PSD Models.

Notation

We denote E = (-1, 1) d and F = (-1, 1) d ′ the state and observation space. P(E) is the set of probability measures on E and M + (E) the set of finite, positive measures on E. We assume that all measures admit a density with respect to the Lebesgue measure and use the abuse of notation µ(dx) = µ(x)dx. β > 0 is a smoothness parameter, σ is a mixing parameter for kernels and ϵ is the accuracy when doing function approximation. Denote R d + the set of vectors in R d with all positive components and S + (R M ) the set of positive definite matrices of size M .

Gaussian PSD Models

Gaussian PSD Models, introduced in Rudi and Ciliberto (2021), is a family of models for non-negative functions and, in particular, probability densities specializing the PSD Models from Marteau-Ferey et al. (2020). They are non-negative everywhere, admit a linear parametrization and can be learned from samples and function evaluations. They are characterized by a linear combination of kernels, with weights chosen such that the function is nonnegative. In this section, we recall the definition of Gaussian PSD Models of Rudi and Ciliberto (2021), show that they extend most well-established probability models and present how to perform operations such as multiplication or marginalization. Definition 1. A Gaussian PSD Model of order M is a function f : R d → R which can be written:

f (x) = M i=1 M j=1 A ij k η (x, x i )k η (x, x j ) (2)
where η ∈ R d + is the precision vector, X = (x i ) 1≤i≤M ∈ R M ×d are the anchor points and A ∈ S + (R d ) is the weight matrix. Such a function is denoted f (x; A, X, η) (or f (x; θ) for shorthand).

Adopting the point of view in Marteau-Ferey et al. (2020), Gaussian PSD Models can equivalently be defined as functions of the form f

(x) = Φ η (x) ⊤ AΦ(x) where Φ η is defined as Φ η (x) = (k η (x, x 1 ) . . . k η (x, x M )) ⊤ ∈ R M . When defined on the product of Euclidean spaces R d × R d ′ with anchor points [X, Y ] (the
row-wise concatenation of X and Y ) and precision vector η = (η, η ′ ) (column-wise concatenation of η 1 and η 2 ), we denote denote the model f (x, y; A, [X, Y ], (η 1 , η 2 )). This split notation is justified by the fact that k η ((x, y)

, (u, v)) = k η 1 (x, u)k η 2 (y, v). Example 1 (Gaussian Mixture Model). Let p(x) = M k=1 α k p(x|µ k , η) where µ k ∈ R d , η k ∈ R d
+ , and α ∈ R d + with M k=1 α k = 1 and p(x|µ k , η) is the Gaussian density with mean µ and precision vector η. p is known as a Gaussian Mixture Model. p can be written as a Gaussian PSD Model of order M f (x; A, X, η/2) with A = diag(a) and X = (µ 1 . . . µ M ) ⊤ .

Example 2 (Squared linear Gaussian model). Let g(x) = w ⊤ Φ η (x) where w ∈ R d and Φ η (x) = (k η (x, x 1 ), . . . , k η (x, x M ) ⊤ . Then, f = g 2 can be written as a Gaussian PSD Model of order M with A = ww ⊤ and X = (x 1 , . . . , x M ) ⊤ . Indeed,

f (x) = (w ⊤ Φ η (x)) 2 = w ⊤ Φ η (x)w ⊤ Φ η (x) = Φ η (x) ⊤ ww ⊤ Φ η (x).
(3)

As pointed out in the introduction of this section, because

k η (x, u)k η (x, v) ∝ k 2η (x, u+v
2 ), a Gaussian PSD Model can be seen as a linear combination of Gaussians. It is important to note that the coefficients of the components can be non-negative, which makes them much more expressive then Mixture models. Consider for instance f (x) = (e -(x-2) 2 -e -(x-3) 2 ) 2 which is clearly non-negative and can be written as a Gaussian PSD Model but not as a Mixture model.

Operations on Gaussian PSD Models

Gaussian PSD Models are compatible with operations on probabilistic models such as integration, partial evaluation, product and marginalization. The operations are summarized in Proposition 1 and the algorithms, based on kernel evaluations and matrix-vector products are detailed in Rudi and Ciliberto (2021).

Proposition 1 (Closed form operations for Gaussian PSD Models). Let f (x, y; θ 1 ) and g(y, z; θ 2 ) be two Gaussian PSD Models of order M 1 and M 2 respectively, as in Definition 1. Then there exist some algorithms Integral, PartialEval, Product, Marginal such that • Integral over R d or over a hypercube f (x, y; θ 1 )dxdy can be computed exactly and in closed form by the algorithm Integral(θ 1 ) with a computational cost of O(M 2 1 d). • Partial evaluationf (x, y 0 ; θ 1 ) = h(x; θ ′ ) is a Gaussian PSD Model of order at most M 1 and θ ′ can be computed exactly and in closed form by the algorithm PartialEval(y 0 , θ 1 ) with a computational cost of O(M 2 1 d) • Product f (x, y; θ 1 )g(y, z; θ 2 ) = h(x, y, z; θ ′ ) is a Gaussian PSD Model of order at most M 1 × M 2 and θ ′ can be computed exactly and in closed form by the algorithm Product(θ 1 , θ 2 ) with a computational cost of O(M 2 1 M 2 2 d). • Marginalization f (x, y; θ 1 )dy = h(x; θ ′ ) is a Gaussian PSD Model of order at most M 1 and θ ′ can be computed exactly and in closed form by the algorithm Marginal(y, θ 1 ) with a computational cost of O(M 2 1 d). The proof of the proposition above can be found in (Rudi and Ciliberto, 2021, Appendix F). Note that a Markov transition g(x) = Q(u, x)f (u)du, when Q and f are Gaussian PSD models, can be decomposed in terms of product and marginalization and computed in closed form, with g again a Gaussian PSD Model. More importantly, g is of order M 1 (instead of the naïve M 1 × M 2 according to Proposition 1). This is summarized in the following proposition Proposition 2 (Constant order for Markov transition). If Q(u, x; θ Q and f (u; θ f ) are two Gaussian PSD Models of order M 1 and M 2 respectively, then

g(x) = Q(u, x; θ Q )f (u; θ f )du, computed via Product and Marginal is a Gaussian PSD Model of order M 1 .
The proof of Proposition 2 can be found in (Rudi and Ciliberto, 2021, Appendix F.5).

3 Learning transition and observation kernels with Gaussian PSD Models

In this section, we show that General Gaussian Models can be used to efficiently approximate smooth, non-negative functions in L ∞ using function evaluations. Let Ω = (-1, 1) d and f : Ω → R be the target function. We assume we can evaluate f (x) at any point x ∈ Ω. We assume that f is the sum of squared β-smooth functions. Formally, we introduce:

Assumption 1 (Smooth sum-of-squares assumption). There exist q ∈ N and β ≥ 0 and

f i ∈ W β 2 (Ω) ∩ L ∞ (Ω) such that f (x) = q i=1 f i (x) 2
. Assumption 1 is verified for most continuous dynamical models of interest. For instance, any transition kernel Q(x, y) ∝ e -∥Σ -1/2 (y-h(x)∥ 2 verifies the assumption. Rudi and Ciliberto (2021) provides a list of sufficient conditions, we recall in Proposition 3.

Proposition 3 (Generality of Assumption 1, Prop. 5 in Rudi and Ciliberto (2021)). A function f satisfies Assumption 1 on Ω = (-1, 1) d as soon as:

• f is a probability density and f ∈ W β 2 (Ω) ∩ L ∞ (Ω), and strictly positive on

[-1, 1] d ; • f is an exponential model f (x) = e -v(x) with v ∈ W β 2 (Ω) ∩ L ∞ (Ω) ;
• f is a mixture of models from (b) ;

• f is β + 2-times differentiable on [-1, 1] d ,
with a finite set of zeroes all in (-1, 1) d , and a positive definite Hessian in each zero.

In Section 3.1, we introduce the optimization problem we solve to learn f and present the learning algorithm. In Section 3.2, we prove that the obtained estimator f is

Learning algorithm

In Section 2, we showed that the square of any Gaussian Linear Model is a Gaussian PSD Model and its weight matrix is of rank 1 and given by A = aa ⊤ where a is the weight vector of the Gaussian Linear Model. We use this insight to efficiently approximate a smooth sum-of-squares function f with a Gaussian PSD Model f . This insight was first published by [START_REF] Marteau-Ferey | Sampling from arbitrary functions via psd models[END_REF] for probability densities. A full-rank estimator f can be also learned, by solving a Semi-Definite Programming problem using e.g. Newton's method.

In this work we propose to approximate f using f = ĝ2 where ĝ is a Gaussian Linear Model learned on g. Denoting the Linear Gaussian Model

ĝ(x; a, η, X) = M i=1 a i k η (x, xi ) where a ∈ R M , η ∈ R d
+ and X ∈ R M ×d , we introduce the optimization problem used to learn â from data points X ∈ R n×d :

min a∈R M 1 n n k=1 | f (x k ) -ĝ(x k ; a, η, X)| 2 + λ a ⊤ Ka (4)
where K is described by

K ij = k η (x i , xj ).
In Appendix B, we cast (4) as a kernel ridge regression problem, which can be efficiently solved for large values of n and M in O(n √ n) time using approximate kernel methods such as in Rudi et al. (2017).

Algorithm 1: LearnRankOne algorithm Data: f (x), M , N , η, λ X ← UniformSample(Ω, N ); Y ← ( f (x i )|x i ∈ X); X ← UniformSample(Ω, M ); â ← KernelRidgeRegression(k η , X, X, Y, λ) ; f (x) ← GaussianPSDModel(aa ⊤ , X, η); Result: f (x)

Learning rates

The Gaussian PSD Model f obtained from Algorithm 1 using function evaluations at uniformly sampled training points X approximates f in L ∞ (Ω) with optimal learning rates for L ∞ norm [START_REF] Wendland | Scattered data approximation[END_REF], if f is a β-smooth and bounded density, as formalized by Assumption 1. We build on the results in Rudi and Ciliberto (2021) and [START_REF] Marteau-Ferey | Sampling from arbitrary functions via psd models[END_REF]. The former studies the L 2 convergence of f to f when the training set is sampled from the target density, using the full-rank counterpart to Algorithm 1. The latter studies convergence in Hellinger distance using Algorithm 1. Both works use the insights of Rudi et al. (2015).

Algorithm 1 finds a solution to (4) in the reproducing kernel Hilbert space H η associated to k η where η is chosen as a function of the desired precision ϵ. Theorem 4 proves that f converges to f in L ∞ , with optimal rates [START_REF] Wendland | Scattered data approximation[END_REF]. In particular, to learn f uniformly to precision ϵ, N ≈ ϵ -2d 2β-d function evaluations and a model with M ≈ N anchor points suffices.

Theorem 4. Let β > d/2 and θ -1 < 1 + 2β/d. Let f : Ω → R such that f verifies Assumption 1. Set M ≥ C ′ (log( 1 ϵ )) d log( 1 δϵ )ϵ -d/β and n ≥ C ′ ϵ -2d/β log 1 δ .
Consider the set of anchor points X ∈ R M ×d and the set of training points X ∈ R n×d sampled independently and uniformly from Ω = (-1, 1) d . Let ϵ ≤ ϵ 0 . Let η = ϵ -2/β 1 d and H η the RKHS associated to k η . Let â solution to the kernel ridge regression problem defined in (4). We denote ĝ(•) = ĝ(• ; â, X, η) the estimator of g and f = ĝ2 . With probability at least

1 -3δ, ∥ f -f ∥ L ∞ (Ω) ≤ C∥ f ∥ 2 W β 2 (Ω) ϵ 1-d 2β (5)
where C, C ′ are constants depending only on β, d and independent of f and ϵ.

The proof of Theorem 4 can be found in Appendix B, and uses arguments from [START_REF] Marteau-Ferey | Sampling from arbitrary functions via psd models[END_REF]. A sketch is given in Section 5.

Gaussian PSD Models for filtering

Now that we know how to efficiently carry out closed-form filtering operations using Gaussian PSD Models and learn good approximations of non-negative functions using this family of models, we can compute an approximation of the filtering distribution.

Setting

Consider two discrete Markov Chains (X n ) ∈ E N -the hidden state chain -and (Y n ) ∈ F N -the observations. We assume that (X t , Y t ) has a Hidden Markov Model structure described by (ν, Q, G) where ν

∈ P(E) is an initial distribution, Q : E × E → R is a Markov kernel, and G : E × F → R is a transition kernel.
Formally, this can be summarized as: (1) the law of X n is fully-determined by the knowledge of X n-1 , i.e. P( 2) the law of Y n is fully-determined by the knowledge of X n ; i.e. 3) the law of X 0 is given by ν. In particular, we have the usual Markovian structure:

X n ∈ dx|X n-1 = x) = Q(x, dx) ; (

P(Y

n ∈ dy|X n = x) = G(x, dy) ; (
X n ⊥ ⊥ X l |X k and Y n ⊥ ⊥ Y l |X k for any n > k > l.
The goal of filtering is to compute the distribution of X n conditionally on past observation Y 1 , . . . , Y n . We denote this distribution π ν n (z 1:n , dx) where ν is initial distribution and z 1:n are the observations (not necessarily taken from the chain (Y n ). This distribution is known as the filtering distribution or optimal filter. When clear from context, we drop the dependence of z 1:n and ν.

Importantly, π n (dx) can be computed recursively using Q and G and beginning from ν:

π n = Rn (π n-1 ), where Rn (µ)(•) := R n µ (•) R n µ (E) and R n µ (•) := Q(u, •)G(x, y n )dµ(u) (6)
which recovers (1). Computing Rn (π n-1 ) is difficult in most circumstances since one must be able to compute products and marginals on probability distributions. Two notable exceptions include the Conditional Linear Gaussian Model (which corresponds to the Kalman filter) and when E is finite (which corresponds to the Baum-Welch algorithm).

PSD filter

To overcome these difficulties, we approximate Q and G from evaluations using Gaussian PSD Models then compute iteration (6) with these approximate kernels.

Given a sequence of observations (z k ) k≥1 , we define Rk (u, x) = Q(u, x) Ĝ(x, y k ) analogously to R n . The non-linear transformation Rk is defined for any positive, finite measure µ by

Rk (µ) = Rk µ Rk µ (E) . ( 7
)
and summarized in Algorithm 2. Note that at each step, πk is a valid, normalized density. However, Q and Ĝ are not properly normalized, i.e. Ĝ(u, F ) = 1 is not guaranteed. As long as the initial distribution π0 is a General PSD Model and Ĝ and Q are valid Gaussian PSD Models, πk is a Gaussian PSD Model for all k ≥ 0 and moreover is a valid density. Importantly, Corollary 5 shows that the order of πk is constant for k ≥ 1 and equal to the product of the orders of Q and Ĝ as shown in the following corollary which follows directly from Proposition 2.

Corollary 5 (π n has constant order for any n). Let π0 (x) a Gaussian PSD Model on

E of order M 0 , Q(u, x) a Gaussian PSD Model on E × E of order M Q and Ĝ(x, y) a Gaussian PSD Model on E × E of order M G . Let (z n ) ∈ F N . Let (π n
) n∈N the sequence of functions defined by the recursion (7). Then, for any n ≥ 0, πn is a Gaussian PSD Model of order at most M Q × M G and it is computed by Algorithm 2.

Algorithm 2: PSDFilter algorithm

Data: z 1 , . . . , z T , π0 , Q, Ĝ for k = 1, . . . , T do β ← Product(π k-1 , Q); Qπ k-1 (•) ← Marginal (β(u, •), [u]); Ĝk (•) ← PartialEval( Ĝ(•, y), y := y k ); πk ← Product Qπ k-1 , Ĝk ; Z ← Integral(π); π ← π/Z;
end Result: π1 , . . . , πT

Gaussian PSD Filter Stability and Robustness

In this section, we show that the algorithm introduced in Section 4.2 estimates a distribution which is close to the true filtering distribution in total variation distance. Theorem 6 combines stability properties of the optimal filter with robustness of the iteration (7) to the use of Q and Ĝ in place of Q and G. We introduce the following assumption.

Assumption 2 (R n is mixing). There exists 1 > σ > 0 and a probability density ξ ∈ P(E) such that for any n ∈ N, R n is σ-ξ-mixing, i.e. for any

u, x ∈ E × E, σξ(x) ≤ R n (u, x) ≤ 1 σ ξ(x). ( 8
)
Assumption 2 is a classical assumption for the study of filtering [START_REF] Cappé | Inference in Hidden Markov Models[END_REF]. We are ready to state the main theorem of this work.

Theorem 6 (PSD filter robustness and stability). Assume that Q and G verify Assumptions 1 and 2. Let ε > 0. When Ĝ, Q are learned using Algorithm 1 and M, N ∈ N are chosen such that

M ≥ C ′ (log (1/ε)) D log (1/εδ) ε -D β-D/2 , n ≥ C ′ ε -D β-D/2 log(1/δ)
then with probability at least 1 -6δ, the following holds: for any k ∈ N,

∥π k -πk ∥ T V ≤ C σ 2 1 -σ 2 1 + σ 2 k-1 ∥π 0 -π0 ∥ T V + ε σ , (9)
where π1 , . . . , πk are computed using Algorithm 2, moreover C = 2 log 3 and

C ′ depends only on ∥ √ Q∥ W β 2 , ∥ √ G∥ W β 2 , β, d.
Theorem 6 is proven in Appendix C and we include a sketch of the proof in Section 5.1.

The theorem above shows that the distance between the probability πk resulting from our algorithm and the optimal one π k corresponding to k steps of sequential Bayesian filtering with the true Q, G is bounded by two terms:

(1) the first term accounts for stability and goes to zero exponentially fast in the number of steps k and depends on how close we choose π0 with respect to π 0 ; (2) the second terms accounts for robustness, it does not increase in k and we can make it arbitrarily small by learning more precise Q, Ĝ, i.e. by increasing M, n. Assuming that π0 = π 0 the proposed algorithm Algorithm 2 for any k ∈ N achieves a an error ∥π k -πk ∥ T V ≤ ε/σ producing a Gaussian PSD model πk that satisfies

πk of order O(ε -2D β-D/2 ), and O(ε -6D β-D/2 ),
for Algorithm 2. The cost of learning Q, Ĝ via Algorithm 1 is paid only once at the beginning and is of O(ε -1.5D β-D/2 ), if we use fast algorithms as the ones recalled in Section 3.2. Note that the proposed approach is adaptive to the regularity of the kernels Q and G. In particular, if they are differentiable many times, i.e. β ≥ 4.5D, then the order of πk becomes only O(ε -1/2 ), i.e.

memory cost O(ε -1 ), computational cost O(ε -1.5 ).
This is remarkable since, for example, particle filter methods are bound to a computational complexity that cannot be smaller than ε -2 since they have to approximate an integral via sampling [START_REF] François | Stability and uniform approximation of nonlinear filters using the Hilbert metric and application to particle filters[END_REF].

5 Sketch of the proof of Theorem 4 and Theorem 6

In this section, we give the main arguments for the proof of Theorem 4. The proof can be found in Appendix B. We denote f the target function and g = √ g, and ĝ the solution to (4) and f = ĝ2 . As shown in Appendix B, by opening the proof of Proposition 11 in [START_REF] Marteau-Ferey | Sampling from arbitrary functions via psd models[END_REF], we can obtain bounds on ∥g -

ĝ∥ L 2 (Ω) is L 2 (Ω) and (essentially) ∥ĝ -g∥ W m 2 (Ω) as a function of ∥g∥ W β 2 (Ω)∩L ∞ (Ω)
and ϵ with optimal dependence on ϵ. We combine these results to bound ∥f -f ∥ L ∞ (Ω) with optimal rates in ϵ, which is needed to control the error terms in Theorem 6. Since

∥f -f ∥ L ∞ (Ω) ≤ 2∥g∥ L ∞ (Ω) + ∥g -ĝ∥ L ∞ (Ω) ∥g -ĝ∥ L ∞ (Ω) , we focus on controlling ∥g -ĝ∥ L ∞ (Ω) .
The first key argument is to separate g-ĝ as the sum of an approximation error g -g τ,ϵ and an estimation error g τ,ϵ -ĝ where g τ,ϵ ∈ H η and ∥gg τ,ϵ ∥ L ∞ (Ω) ≤ C∥g∥ϵ 1-ν . Using the triangle inequality, controlling ∥g-ĝ∥ L ∞ (Ω) reduces to controlling ∥g τ,ϵ -ĝ∥ L ∞ (Ω) . We apply the Gargliano-Niremberg inequality with well-chosen parameters to obtain a bound with the product of the Sobolev norm of the estimation error (which is controlled by the H η norm) and of the L 2 (Ω) norm, both of which are bounded by the learning approach. Combining all the terms yields the result in Theorem 4.

Sketch of the proof of Theorem 6

We start by proving a lemma which controls how errors committed at each step accumulate and are eliminated.

Lemma 7. Let T > 1. ∥π ν n -πn ∥ T V ≤ C σ τ n σ ∥π 0 -ν∥ T V An + δ n + C σ n k=1 τ n-k-1 σ δ k Bn ( 10
)
where

δ k = ∥ Rk (π k-1 ) -Rn (π k-1 )∥ T V .
The gist of Lemma 7 is that, using the triangle inequality, we can bound ∥π ν n -πn ∥ by A n = ∥π ν n -π π0 n ∥ and B n = ∥π π0 -πn ∥. Lemma 7 is proven by combining Proposition 21 and Proposition 20 in Appendix C.

The first source of discrepancy is the intialization error. Indeed, πk is initialized at π0 while π ν k is initialized at ν and we isolated this question above by considering the discrepancy between two optimal filters intiialized at ν and π0 . The behavior of A n is known as the stability or forgetting property, and is a property of the optimal filter, i.e. of the dynamical system considered, and not of the algorithms considered. In particuliar, since R n is mixing for any n ≥ 0, A n decreased exponentially as the application of the filter is contractive (under the mixing assumption).

The second source of discrepancy is the accumulation of errors committed at each step by applying the approximate iteration (7) in place of (6). A telescopic argument in the proof of Lemma 7 shows that the accumulation is limited because the forgetting property of the optimal kernel tends to make past errors δ j for j < k disappear thanks to the τ n-k-1 term. Our argument generalizes the argument in Gland and Oudjane ( 2004) and relies on the projective properties of the Hilbert metric [START_REF] Samuel | Hyperbolic contractivity and the hilbert metric on probability measures[END_REF]. We then combine Lemma 7 and Theorem 4 to obtain bound B n by a constant ε as small as we want.

6 Extension: filtering with Generalized Gaussian PSD Models

Throughout this work, we have focused on Gaussian PSD Models. In fact, many of the properties studied above remain valid for a more general family of models we introduce and study in this section.

Seen as a mixture (with potentially negatively weighted components), Gaussian Mixture Models have components which are aligned with the axes of the space. In cases where f (u, x) = Q(u, x) a transition kernel, we know that u and x are strongly "correlated" (indeed, Q(u, x) is a conditional density) and having non-diagonal precision matrices can be useful. When Q(u, x) is a Gaussian Linear Conditional Distribution, one Gaussian component with non-diagonal precision matrix is enough to approximate Q uniformly. We prove this in Theorem 10.

We introduce Generalized Gaussian PSD Models as the family of nonnegative combinations of such components. In this section, we show that Generalized Gaussian PSD Models share many of the properties of Gaussian PSD Models. In addition, we show in Theorem 10 that they generalize Kalman filters to more general initial distributions such as multi-modal models.

In this section, for P a positive semi-definite matrix, we denote for any x, y ∈ R d , k P (x, y) = e -(x-y) ⊤ P (x-y) and C(P ) = R d k P (x, 0)dx. Definition 2 (Generalized Gaussian PSD model of order M ). A Generalized Gaussian PSD model of order M is a function f : R d → R which can be written:

f (x) = T r(AB(x)) ( 11
)
where A is a positive semi-definite matrix of size M and for any x ∈ R d , B(x) is a positive semi-definite matrix with entries

B(x) ij = e C ij k P ij (x, x ij ) (12)
where

C ij ∈ R, x ij ∈ R d and P ij is a d × d positive semi-definite matrix, for any 1 ≤ i, j ≤ M .
We use the notation f (x; θ) where θ = (X, A, P, C) and X = (x ij ) 1≤i,j≤M , P = (P ij ) 1≤i,j≤M , and C = (C ij ) 1≤i,j≤M .

Like Gaussian PSD Models, Generalized Gaussian PSD Model are a generalization of Gaussian Mixture Models allowing for negative weights. Indeed, by developing ( 11) with ( 12), f can be written :

f (x) = M i=1 M j=1 A ij e C ij k P ij (x, x ij ).

Examples of Generalized Gaussian PSD Models

Below we show that Generalized Gaussian PSD Models generalize Gaussian Mixture Models, Gaussian PSD Models from Rudi and Ciliberto (2021) and indeed, any squared linear combination of Gaussian functions.

Example 3 (Gaussian Mixture Models are Generalized Gaussian PSD Models). Let a ∈ R d + such that d i=1 a i = 1. If f (x) = M i=1 a i p(x|µ i , P i ) with a i ≥ 0 and M i=1 a i = 1, then f verifies Definition 2 with A = diag(a), P ii = P i /2, x ii = x i /2 and C ii = -1 2 log(C(P i )). Example 4 (Gaussian PSD Models are Generalized Gaussian PSD Models).

Consider a Gaussian PSD Model

f (x) = Φ M (x) ⊤ AΦ M (x) with Φ M (x) = (k η (x, x 1 ), . . . , k η (x, x M ) ⊤ ∈ R M . Then, f a Generalized Gaussian PSD Model. Indeed, f (x) = T r(AΦ M (x)Φ M (x) ⊤ ). Each component of the positive semi-definite matrix Φ(x)Φ(x) ⊤ is a function of the form x → e C ij k 2η (x, x i +x j √ 2 ) with C ij = ∥x i ∥ 2 + ∥x j ∥ 2 - ∥x i +x j ∥ 2 2
, and f verifies Definition 2.

Example 5 (Squared Gaussian Linear Models are Generalized Gaussian PSD Models). Consider the function f

(x) = w ⊤ Φ M (x) 2 where Φ M (x) = (k P 1 (x, x 1 ), . . . , k P M (x, x M ) ⊤ and w ∈ R d . Then, f is a Generalized Gaussian PSD Model. Indeed, f (x) = T r Φ M (x) ⊤ ww ⊤ Φ M (x) = T r ww ⊤ Φ M (x)Φ M (x) ⊤ . Since ww ⊤ and Φ M (x)Φ M (x) ⊤ are both positive semi-definite matrices, f (x) ≥ 0 and Φ M (x)Φ M (x) ⊤ ij = e C ij k P ij (x, x ij ) with P ij = P i +P j , x ij = P -1/2 ij (P i x i + P j x j ).
6.2 Properties of Generalized Gaussian PSD Models 6.2.1 Closed-form stability with respect to probabilistic operations

Like Gaussian PSD Models, Generalized Gaussian PSD Models are closed under product, partial evaluation, and marginalization.

Proposition 8. Let f (x, y; θ 1 ) and g(y, z; θ 2 ) be two Generalized Gaussian PSD Models of order M 1 and M 2 respectively, where all precision matrices are positive definite.

• Integral over R d f (x, y; θ 1 )dxdy is given by the algorithm Integral(θ 1 ).

• Partial evaluationf (x, y 0 ; θ 1 ) = h(x; θ ′ ) is a Generalized Gaussian PSD Model of order M 1 and θ ′ is given by the algorithm PartialEval(y 0 , θ 1 ).

• Product f (x, y; θ 1 )g(y, z; θ 2 ) = h(x, y, z; θ ′ ) is a Generalized Gaussian PSD Model of order M 1 × M 2 and θ ′ is given by the algorithm Product(θ 1 , θ 2 ).

• Marginalization f (x, y; θ 1 )dy = h(x; θ ′ ) is a Generalized Gaussian PSD Model of order M 1 and θ ′ is given by the algorithm Marginal(y, θ 1 ).

Proof sketch for the product Using (11), f (x, y)g(y, z) can be written h(x, y, z) = T r(A 1 ⊗ A 2 × B 1 (x, y) ⊗ B 2 (y, z)). The entries of B 1 (x, y) ⊗ B 2 (y, z) are products of Gaussian functions which can simplified into the the form (12). The complete proof and description of each operation can be found in Appendix D.

Closed-form filtering iteration

An optimal filtering iteration (6) can be written using the four operations of Proposition 8.

Proposition 9. Let µ(x; θ µ ), q(x, x ′ ; θ q ) and g(x, y; θ g ) be three Generalized Gaussian PSD models with order M , M q and M g respectively. Let y ∈ R d such that q(u, x)g(x, y)µ(u)dudx > 0. The density µ ′ defined by µ ′ (x) = q(u,x)g(u,y)µ(u)du q(u,x)g(x,y)µ(u)dudx is a Generalized Gaussian PSD Model with order at most M × M q × M g whose parameters are given by FilterStep(y, θ µ , θ q , θ g ).

Applying Proposition 9 recursively as in Algorithm 2 to compute an approximate filter πn , the order of πn increases exponentially with n. A constant number of anchor points can be used by compression πn at each step, for example by learning a Gaussian PSD Model with a given number of anchor points (indeed πn is a smooth sum-of-squares), which is justified by Theorem 4.

Generalized Gaussian PSD Models generalize Kalman filters

Conditional Gaussian Linear Distributions are widely used in filtering and dynamical modeling since they cover transition or observation state-space equations such as X t+1 = F X t +b+ΣU t where U t is Gaussian noise considered in the Kalman filter and extensions.

Theorem 10 (Approximating a Conditional Gaussian Linear Distribution). Let p be a Conditional Gaussian Linear Distribution defined by p(y|x

) = N (y|F x + b, Σ) with F ∈ R d ′ ×d , b ∈ R d ′ and Σ ∈ S ++ d ′ (R).
Then for any ϵ > 0 and R > 0, there exists a Generalized Gaussian PSD Model of order 1

such that |p(y|x) -p(x, y)| ≤ ϵ,∀x, y ∈ R d × R d ′ such that ∥x∥ 2 2 + ∥y∥ 2 2 ≤ R 2
. This shows that a Generalized Gaussian PSD Model of order 1 can approximate a Conditional Gaussian Linear Model with arbitrary accuracy on a compact. Note that in this case, applying Algorithm 2 with Q and Ĝ such approximations (each of order 1) and π0 of order M yields an approximation of π n of constant order M . The proof of Theorem 10 can be found in Appendix D.5.

Generalized Gaussian PSD Models can be used to learn general transition kernels using non-convex optimization. We discuss this in Appendix D.6.

A Markov kernels & Hidden Markov Models

A.1 Tools and notation

Let (E, E) be a measurable space. We denote M(E) the set of finite signed measures on (E, E), M + (E) the set of finite positive measures, M 0 (E) the set of finite signed measures which sum to 0 and P(E) the set of probability distributions on (E, E). Let (F, F) be a second measurable space.

A.2 Total variation

Definition 3 (Total variation norm). Let ξ be a finite signed measure on (E, E). The total variation of ξ is ∥ξ∥ T V = ξ + (E) + ξ -(X) where ξ + , ξ -is the Jordan-Hahn decomposition of ξ. If ξ admits a density with respect to the Lebesgue measure, then ∥ξ∥ T V = |ξ(x)|dx.

Proposition 11. Let µ and ν be two finite measures on (E, E). Then,

∥μ -ν∥ = ∥µ -ν∥ µ(E) + |µ(E) -ν(E)| µ(E) (13)
In particular,

∥μ -ν∥ ≤ 2∥µ -ν∥ µ(E) (14)

A.3 Transition kernels

We report the essential results relative to transition kernels taken [START_REF] Cappé | Inference in Hidden Markov Models[END_REF].

Definition 4 (Transition kernel). A function Q : E × F → R + is an unnormalized transition kernel if:

• for all x ∈ E, Q(x, •) is a positive measure on (F, F); • for all A ∈ F, x → Q(x, A) is measurable. Q is normalized if for any x ∈ E, Q(x, F ) = 1. When E = F and Q is normalized, Q is said to be a Markov transition kernel.
By abuse of notation, when Q admits a density with respect to the Lebesgue measure, we denote it Q as well, i.e. Q(x, dy) = Q(x, y)dy.

Note that R n (u, x) = Q(u, x)G(x, z n ) is indeed an unnormalized transi- tion kernel.
Definition 5 (Effects of kernels). Let K be an unnformalized kernel on (E, E) × (F, F), µ ∈ M + (E) and f a bounded function f : E → R.

Then, Kµ ∈ M + (F ) with for any A ∈ F,

Kµ(A) = K(u, A)µ(du), (15)
and Kf : E → R is a bounded function with for any u ∈ E,

Kf (u) = K(u, A)µ(du). (16)

A.4 Hidden Markov Models

Definition 6 (Hidden Markov Model). Let (E, E) and (F, F) be two measurable spaces. Let Q and G denote a Markov transition kernel on (E, E) and G denote a transition kernel from (E, E) to (F, F). Let T be the Markov transition kernel defined on the product space

(E × F, E ⊗ F) by ∀(x, y) ∈ E×F, ∀C ∈ E⊗F, T [(x, y), C] = 1 C ((x ′ , y ′ ))Q(x, dx ′ )G(x ′ , dy ′ ) (17)
The Markov Chain {X k , Y k } k≥0 with Markov transition kernel T and initial distribution ν ⊗G, where ν is a probability distribution on (E, E) is called a Hidden Markov Model.

We denote P ν and E ν the probability measure and corresponding expectation associated with the process

{(X k , Y k )} over (E × F ) N , (E ⊗ F) ⊗N .
Throughout this work, we assume that for any x ∈ E, Q(x, •) ≪ λ(•) and G(x, •) ≪ λ(•) where λ is the Lebesgue measure over (E, E), and we denote Q(x, dx ′ ) = Q(x, x ′ )dx ′ and G(x, dy) = G(x, y)dy.

Definition 7 (Filtering distribution). Let ν be a probability distribution over (E, E) and n ≥ 0. We denote π ν n the conditional distribution of X n given Y 1:n , i.e.

• π ν n is a transition kernel from F n to E • π ν n satisfies for any bounded function f : E → R, E ν [f (X n )|Y 1:n ] = f (x)π ν n (Y 1:n , dx) (18)
A.5 Mixing kernels Definition 8 (Mixing kernel). We say a kernel K(x, dy) is mixing if there exists a positive constant σ > 0 and a non-negative measure ξ such that for any x ∈ E and A ∈ F,

σξ(A) ≤ K(x, A) ≤ 1 σ ξ(A). ( 19
)
If K is mixing for a measure ξ and a constant σ we write that K is σ-ξ-mixing.

Remark 1. Note that we can add the constraint that ξ be normalized. Indeed, if

K is σ-ξ-mixing with ξ(E) ̸ = 1, then K if σ-ξ-mixing with σ = σ × min ξ(E), 1 ξ(E) ≤ σ.
Proposition 12 (Sufficient condition for mixing when K admits a density).

If K(x, dy) = κ(x, y)dy and there exists σ > 0 and a measure density ξ such that for any x ∈ E and y ∈ F ,

σξ(y) ≤ κ(x, y) ≤ 1 σ ξ(y). ( 20
)
then K is mixing with constant σ and ξ(A) = 1 A (y)ξ(y)dy.

Proposition 13. If K is mixing with σ and ξ, then for any µ,

σξ(A) ≤ Kµ(A) ≤ 1 σ ξ(A). (21)
Proof.

σξ(A) ≤ Kµ(A) = 1 A (x) K(u, dx)µ(du) ≤ 1 A (x) 1 σ ξ(dx)µ(du) = 1 σ ξ(A) (22) A.6 Optimal kernel R n Definition 9. Let y n ∈ F . Let R n : E × E → R + a kernel defined by: for any bounded function f ∈ F b (E), and u ∈ E, R n (u, f ) = f (x)Q(u, x)G(x, y n )dx (23)
We call R n the optimal kernel.

For alternative definitions see for example Cappé et al., 2005, page 220.

A.7 Hilbert metric

Definition 10 (Comparable measures). Let µ and ν be two measures on (E, E). µ and ν are said to be comparable if there exists ∞ > a, b > 0 such that for any

A ∈ E, aν(A) ≤ µ(A) ≤ bν(A) (24)
Proposition 14. Let µ and ν be two comparable measures on (E, E) and let K be an unnormalized transition kernel on (E, E). Then, for any n ≥ 0, K n µ and K n ν are comparable.

Proof. By recursion,

bKν(A) ≤ Kµ(A) = 1 A (x) K(u, dx)µ(du) ≤ aKν(A) (25)
Definition 11 (Hilbert metric). Let µ and ν be two comparable probability distributions on (E, E). The Hilbert metric between µ and ν is defined as

h(µ, ν) = log sup A∈E,ν(A)>0 µ(A) ν(A) inf A∈E,ν(A)>0 µ(A) ν(A) (26)
Proposition 15. If µ and ν are two comparable probability distributions on (E, E) and furthermore they both admit densities then,

h(µ, ν) = log ess sup x µ(x) ν(x) ess inf x µ(x) ν(x) = log µ(x) ν(x) ∞ ν(x) µ(x) ∞ (27)
Definition 12 (Birkhoff contraction coefficient). Let K be an unnormalized transition kernel. Define τ (K) such that

τ (K) = sup 0<h(µ,ν)<∞ h(Kµ, Kν) h(µ, ν) (28)
where the supremum is taken over comparable, positive measures µ and ν.

Proposition 16 (Properties of the Birkhoff coefficient). Let K be an unnormalized transition kernel.

• τ is sub-multiplicative, i.e. τ (KL) ≤ τ (K)τ (L)

• τ ≤ 1 • if in addition K is σ-ξ-mixing, then τ (K) ≤ 1-σ 2 1+σ 2
These properties are proven in [START_REF] Samuel | Hyperbolic contractivity and the hilbert metric on probability measures[END_REF], which studies the Hilbert metric is detail.

Proof. We have τ (K) = tanh 1 4 ∆(K) where ∆(K) = sup µ,µ ′ h(Kµ, Kµ ′ ) (∆ for diameter). Since K is σ-ξ-mixing, if µ and ν are two finite measures on (E, E),

σ 2 Kµ(A) ≤ σξ(A) ≤ Kµ(A) ≤ 1 σ ξ(A) ≤ 1 σ 2 Kµ ′ (A) (29) So, h(K) ≤ log 1 σ 4 ≤ 1-σ 2 1+σ 2 .
Proposition 17 (Total variation -Hilbert comparaisons). Without any hypotheses on µ, ν finite measures, if μ = µ/µ(E) and ν = ν/ν(E) are normalized counterparts to µ and ν, then

∥μ -ν∥ ≤ 2 log(3) h(µ, ν) (30)
If in addition, K is an σ-mixing kernel,

h(Kµ, Kν) ≤ 1 σ 2 ∥µ -ν∥ (31)
These inequalities are proven in Cohen and Fausti (2023).

B Proof of Theorem 4

Notation and background results In this section, ∥ • ∥ without any subscript denotes ∥ • ∥ W β 2 (Ω) . Let f a target function defined on Ω = (-1, 1) d that verifies Assumption 1, and let g = √ f . Let ν = min(1, d/2β). Let ϵ > 0 and δ > 0. In this proof all constants C and exponents α, γ, ρ, . . . are independent of f, g, ĝ, f unless otherwise stated. Only β, ν, ñu and θ have importance. We recall also the Gagliardo-Nirenberg inequality, that will be used later.

Lemma 18 (Gagliardo-Nirenberg inequalityWendland (2004)). Let u ∈ L ∞ (R d )∩ W m,2 (R d ) then ∥u∥ L ∞ (Ω) ≤ C∥u∥ θ W m,2 (Ω) ∥u∥ 1-θ L 2 (Ω) (32)
where θ = d m and C is independent of u.

Setting the function space H η and g τ,ϵ Set τ = ϵ -2/β and λ = ϵ 2β+d β and H η the reproducing kernel Hilbert space associated to k η where η = τ 1 d .

Existence and properties of g τ,ϵ As a consequence of the Stein extension theorem (see Corollary A.3 of Rudi and Ciliberto (2021)), there exists a function g ∈ W β 2 (Ω) such that g|Ω = g and ∥g∥

W β 2 (Ω) ≤ ∥g∥ W β 2 (Ω)
and ∥g∥ L ∞ (Ω) ≤ C∥g∥ L ∞ (Ω) . According to Rudi and Ciliberto (2021) and Marteau-Ferey et al. ( 2022) (Proposition 7) there exists g τ,ϵ ∈ H η and some constants C 1 , C 2 depending only on β, d and independent of g such that:

-

g τ ∥ L ∞ (Ω) ≤ C 1 ϵ 1-ν ∥g∥ W β 2 (Ω) , (33)
∥g τ ∥ Hη ≤ C 2 ∥g∥ W β 2 (Ω) ϵ -d/2β . ( 34
)
Learning ĝ and f Let M, n ∈ N, draw X ∈ R n×d the training set and X ∈ R M ×d the set of anchor points. Define y ∈ R n such that y i = g(x i ) for any 1 ≤ i ≤ n and x i is the i-th row of X. We formalize (4) explicitly as a kernel ridge regression problem below:

min a∈R d 1 n |a T Φ η (X i ) -y i | 2 + λa T Ka (35)
where Φ η (x) = (k η (x, x1 ) . . . k η (x, xM )) T ∈ R M and K ∈ R M ×M is given by 35) is strongly convex and has a unique solution â. We denote ĝ the Gaussian Linear Model defined by â, η, and X, i.e. such that for any x ∈ R d , ĝ(x) = âT Φ η (x). Define f the Gaussian PSD Model f = ĝ2 where f (x) = Φ η (x) T ââ T Φ η (x). The analysis of [START_REF] Marteau-Ferey | Sampling from arbitrary functions via psd models[END_REF] shows that when M ≥ C ′ log d ( 1 ϵ ) log( 1 δϵ )ϵ -d/β , and n ≥ C ′ ϵ -d/β log 1 δ , then there exist two constants C 3 and C 4 independent of g and ĝ such that ĝ verifies the following inequalities each with probability at least 1 -δ,

K ij = k η (x i , xj ). Problem (
∥ĝ∥ H ≤ C 3 ∥g∥ W β 2 (Ω) ϵ -d/2β (36) ∥g -ĝ∥ L 2 (Ω) ≤ C 4 ∥g∥ W β 2 (Ω) ϵ. ( 37
)
Deriving the bound for g -ĝ in L ∞ (Ω) In particular, using the triangle inequality and combining (33) and (37), there exists a constant C 5 such that with probability at least 1 -δ,

∥g τ,ϵ -ĝ∥ L 2 (Ω) ≤ C 5 ∥g∥ W β 2 (Ω) ϵ. ( 38
)
We now have all the ingredients to bound ∥g -ĝ∥ L ∞ (Ω) in high probability.

First, notice that :

∥g -ĝ∥ L ∞ (Ω) ≤ ∥g -g τ,ϵ ∥ L ∞ (Ω) + ∥g τ,ϵ -ĝ∥ L ∞ (Ω) . (39)
We apply the Gagliardo-Nirenberg inequality (Lemma 18) to ∥g τ,ϵ -ĝ∥ L ∞ (Ω) :

∥ĝ -g τ,ϵ ∥ L ∞ (Ω) ≤ C∥ĝ -g τ,ϵ ∥ θ W m 2 (Ω) ∥ĝ -g τ,ϵ ∥ 1-θ L 2 (Ω) . (40)
with θ = d 2m and m ≥ d/2 (we fix m when we optimize the exponents below) and C is a constant independent of ĝ and g τ,ϵ .

The Sobolev norm ∥h∥ W m 2 (R d ) is upper bounded by the rkhs norm ∥h∥ Hη for h ∈ H η ⊂ W m 2 (R d ). Thus, applying the triangle inequality and bounds (34), ( 36) and ( 38) there exists a constant C 6 such that with probability at least 1 -2δ, ∥ĝ -

g τ,ϵ ∥ L ∞ (R d ) ≤ C 6 ∥g∥ϵ 1-θ-θd/2β . (41)
Combining ( 33) and ( 41), there exist two constants C, C ′ > 0 such that with probability at least 1 -2δ,

∥g -ĝ∥ L ∞ (Ω) ≤ C∥g∥ W β 2 (Ω) ϵ 1-ν + C ′ ∥g∥ W β 2 (Ω) ϵ 1-θ-θd/2β (42)
Choosing m = β + d 2 , there exists a constant C 7 > 0 such that for ϵ small enough with probability at least 1 -2δ,

∥g -ĝ∥ L ∞ (Ω) ≤ C∥g∥ W β 2 (Ω) ϵ 1-d/2β . ( 43
)
Bounding ∥g + ĝ∥ L ∞ (Ω) Using the triangle inequality and (36), there exists a constant C 8 > 0 and ρ > 0 such that with probability at least 1 -δ,

∥g + ĝ∥ L ∞ (Ω) ≤ 2∥g∥ L ∞ (Ω) + ∥g -ĝ∥ L ∞ (Ω) ≤ 2∥g∥ L ∞ (Ω) + C∥g∥ W β 2 (Ω) ϵ 1-d/2β . ( 44
) Bounding ∥f -f ∥ ∞ Notice that since a 2 -b 2 = (a -b)(a + b), ∥ f -f ∥ L ∞ (Ω) ≤ ∥g -ĝ∥ L ∞ (Ω) ∥g + ĝ∥ L ∞ (Ω) . (45)
Combining the above bounds:

∥ f -f ∥ L ∞ (Ω) ≤ C∥g∥ 2 W β 2 (Ω) ϵ 2-d/β + C ′ ∥g∥ 2 W β 2 (Ω) ϵ 1-d/2β
(46)

Since β > d/2, under the conditions on M, n, there exists a constant C ′′ = C + C ′ depending only on Ω, d, β and independent of f, g, ĝ, f such that with probability at least 1 -3δ,

∥ f -f ∥ L ∞ (Ω) ≤ C ′′ ∥g∥ 2 W β 2 (Ω) ϵ 1-d/2β
(47)

Recall that for any µ ∈ P(E), σξ(E) ≤ R n µ(E) ≤ 1 σ ξ(E), which allows us to control the denominator:

δ n ≤ 2 σξ(E) ∥ Rn (π n-1 ) -R n (π n-1 )∥ T V . (65)
Because we will be able to bound the quality of approximation between Q and Q (and between G and Ĝ), we split the above expression:

∥ Rn (π n-1 ) -R n (π n-1 )∥ T V ≤ πn-1 (u)|Q(u, x)G(x, z n ) -Q(u, x) Ĝ(x, z n )|dudx (66) ≤ Q(u, x)π n-1 (u)|G(x, z n ) -Ĝ(x, z n )|dudx An (67) + Ĝ(u, x)π n-1 (u)|Q(u, x) -Q(u, x)|dudx Bn . (68)
First, let us bound A n .

A n = Q(u, x)π n-1 (u)|G(x, y n ) -Ĝ(x, y n )|dudx (69) ≤ ∥G -Ĝ∥ L ∞ (E×E) πn-1 (u) Q(u, x)dxdu (70) = ∥G -Ĝ∥ L ∞ (E×E) . (71)
where we used that Q is a transition kernel (i.e. that Q(u, E) = 1 for all u ∈ E) and that πn-1 is a distribution. Second, let us bound B n :

B n = Ĝ(x, z n )π n-1 (u)|Q(u, x) -Q(u, x)|dudx (72) ≤ ∥Q -Q∥ L ∞ (E×E) πn-1 (u)du Ĝ(x, z n )dx (73) ≤ Vol(Ω)∥ Ĝ∥ L ∞ (E×E) ∥Q -Q∥ L ∞ (E×E) (74)
where we again used that πn-1 is a probability distribution.

C.3 Putting everything together

We assume that E = F without loss of generality (simply replace D = 2d for Q and D = d + d ′ for G). We choose β > D/2. Thus, here Ω = E × E ⊂ R D where D = 2d. Let γ > 0 and apply Theorem 4 with its parameter ϵ = γ, Then we have that there exist Ĝ and Q and two constants

C 1 , C 2 such that ∥ Ĝ-G∥ L ∞ (Ω) ≤ C 1 ∥ √ G∥ 2 W β 2 (Ω) γ 1-D 2β , ∥ Q-Q∥ L ∞ (Ω) ≤ C 2 ∥ Q∥ 2 W β 2 (Ω) γ 1-D 2β
where C 1 and C 2 are independent of Q, Q, Ĝ, G.

By combining these inequalities with Proposition 22, since ξ(E) = 1:

δ n ≤ 2 σ C∥ √ G∥ 2 W β 2 (Ω) γ 1-D 2β + C d ∥ Ĝ∥ 2 L ∞ (Ω) ∥ Q∥ 2 W β 2 (Ω) γ 1-D 2β . ( 75
) Since ∥ Ĝ∥ L ∞ (Ω) ≤ ∥G-Ĝ∥ L ∞ (Ω) +2∥ √ G∥ 2 L ∞ (Ω) ≤ ∥G-Ĝ∥ L ∞ (Ω) +2∥ √ G∥ 2 W β 2 (Ω) , δ n ≤ 2 σ C∥ √ G∥ 2 W β 2 (Ω) γ 1-D 2β + C + Cϵ 1-D 2β ∥ √ G∥ 2 W β 2 (Ω) ∥ Q∥ 2 W β 2 (Ω) γ 1-D 2β
Then, by upper-bounding the negligible terms, there exist two constants

C 1 , C 2 > 0 independent of Q, Q, G, Ĝ, such that : δ n ≤ C 1 C 2 ∥ √ G∥ 2 W β 2 (Ω) + C 2 ∥ √ G∥ 2 W β 2 (Ω) ∥ Q∥ 2 W β 2 (Ω) C ′ (Q,G) γ 1-D 2β σ (76)
Note that C ′ only depends on parameters on ∥

√ G∥ W β 2 (Ω) , ∥ √ Q∥ W β 2 (Ω) , d, β
, Ω and does not depend on Q, Ĝ, σ, γ.

We now apply Theorem 4. Let ε > 0 and δ > 0. As a consequence of the development above, if γ is chosen such that γ = ε C ′ 2β 2β-D /σ then (1) M, n correspond to the ones stated in the statement of the theorem (2) with probability at least 1 -6δ,

δ n ≤ ε σ (77)
and then (3),

D Computations on Generalized Gaussian PSD Models

The stability properties of Generalized Gaussian PSD Models under probabilistic operations rely at a high-level on the fact that T r(AB)T r(CD) = T r(A ⊗ CB ⊗ D) and that if A and B are positive semi-definite matrices then so is A ⊗ B.

D.1 Integral

Proposition 23 (Integration of a Generalized Gaussian PSD Model). Let f (x) = T r(AB(x)) with parameters {A, C, (P ij ), (µ ij )}. Then, Z = f (x)dx where

Z = T r(A • exp • (C) • C(P )) (78)
where exp • is the element-wise exponential map and C(P ) ∈ S(R M ) is decribed by C(P ) ij = C(P ij ).

We denote • the Hadamard product.

Proof. The proof is clear by linearity of the trace.

Remark 2 (Computational complexity). Because of the need to compute the determinant of P ij the computational complexity of the partial evaluation operation is O(M 2 d!).

D.2 PartialEval

Proposition 24 (Partial evaluation of a Generalized Gaussian PSD Model).

Let f (x, y) = T r(AB(x, y)) with parameters {A, C, (P ij ), (µ ij )} and y ∈ R d . Then, g(x) := f (x, y) = T r(A ′ B ′ (x)) with parameters A ′ , C ′ , (P ′ ij ), (µ ′ ij ) where

A ′ = A (79) P ′ ij = P ijxx (80) µ ′ ij = µ ijx + P -1 ijxx P ijxy (µ ijy -y) (81) C ′ ij = C ij + νP -1 xx ν -
Proof. We can compute B ′ (x) by expanding C -log(B(x, y) ij ) for any i, j.

Dropping the i, j dependence:

P 1/2 x y -µ 2 =x T P xx x + 2yP
Remark 3 (Computational complexity). Because of the need to compute the inverse of P xx , the computational complexity of the partial evaluation operation is O(M 2 d 3 x ).

D.3 Marginalization

Proposition 25 (Marginalization of a Generalized Gaussian PSD Model).

Let f (x, y) = T r(AB(x, y)) with parameters {A, C, (P ij ), (µ ij )}. Then, h(x) := f (x, y)dy = T r(A ′ B ′ (x)) with parameters A ′ , C ′ , (P ′ ij ), (µ ′ ij ) where

A ′ = A (88) P ′ ij = P -1 ij xx -1 (89) µ ′ ij = [µ ij ] x (90)
C ′ ij = C ij + log(C(P ij )) -log(C(P ′ ij )) (91)
Proof. We compute the integral component-wise, denoting Σ = (2P ) -1 :

B(x, y) ij dy = e C ij exp      - √ 2P 1/2 ij x y -µ ij 2 2      dy (92) = e C ij exp      - Σ -1/2 x y -µ ij 2 2      dy (93) = e C ij (2π) dx+dy |Σ| 1 (2π) dx+dy |Σ| exp      - Σ -1/2 x y -µ ij 2 2      dy (94) = e C ij (2π) dx+dy |Σ| 1 (2π) dx | [Σ] xx | exp    - [Σ] -1/2 xx x -[µ ij ] x 2 2    (95) = e C ij C P ij C P ′ ij exp -P ′ ij 1/2 x -[µ ij ] x 2 (96) = e C ′ ij exp -P ′ ij 1/2 x -µ ′ ij 2 (97)
where

C ′ ij = C ij + log(C P ij ) -log(C P ′ ij ), P ′ ij = P -1 ij xx -1 and µ ′ ij = [µ ij ] x .
Remark 4 (Computational complexity). Because of the need to compute the determinant of P as well as invert it, the computational complexity of the partial evaluation operation is O(M 2 max(d!, d 3 )). where µ ijkl = P -1 ijkl μijkl (111)

D.4 Product

C ijkl = C f ij + C g kl + μijkl P -1 ijkl μijkl -µ ij P ij µ ij -μkl Pkl μkl (112)
Remark 5 (Computational complexity). Because of the need to compute the inverse of P ′ the computational complexity of the product operation between models of order M and m is O(M 2 m 2 d 3 ).

D.6 Learning Generalized Gaussian PSD Models

From an approximation perspective, a Generalized Gaussian PSD Model is a Gaussian PSD Model in which one can optimize the anchor points x and precision matrices P of each kernel function. In the case of approximating transition kernels, this can yield significant improvements in model order. Indeed, a transition kernel Q(u, x) is a conditional probability distribution which depends in which the probability of the value x depends on the value u. This dependence is encoded in the combination of kernel evaluations but not in the kernel evaluations themselves.

To approximate a function f with a Generalized Gaussian PSD Model, we implicitly approximate the square-root of f using a Gaussian Linear Model:

min ĝ∈G M 1 n n i=1 f (x i ) -ĝ(x i ) 2 2 , (125)
where x i are sampled or chosen on a grid, and 125) is a smooth, nonconvex problem which can be solved approximately using off-the-shelf solver like L-BFGS [START_REF] Liu | On the limited memory bfgs method for large scale optimization[END_REF].

G M = { M j=1 α i k P i (x, µ i ) | µ i ∈ R d , P = R ⊤ i R i , R i ∈ R d×d , α i ∈ R for 1 ≤ j ≤ M }. (
In practice, we initialize the model by placing µ i is regions where f (µ i ) is large. In the case where f (u, x) is a transition kernel Q(u, x), one strategy is to chose [µ i ] u on a grid (or sampled uniformly) and then choose [µ i ] v such that f (µ i ) = sup x f ([µ i ] u , v). This is particularly interesting when Q is a non-linear Gaussian model Q(u, x) ∝ e -∥Σ -1/2 (x-h(u))∥ 2 for some non-linear transition model h.

Proposition 26 (C-

 26 Product of two Generalized Gaussian PSD Models). Let f (x, y) = T r(AB(x, y)) a generalized PSD model of order M with parameters {A, C, (P ij ), (µ ij )}. Let g(x, y) = T r( Ã B(x, y)) a generalized PSD model of order m with parameters Ã, C, ( Pkl ), (μ kl ) .Then, h(x) := f (x)g(x) = T r(A ′ B ′ (x)) is a generalized PSD model of order M m with parameters A ′ , C ′ , (P ′ ij ), (µ ′ ij )whereA ′ = A ⊗ Ã ijkl = C f ij + C g kl + μijkl P ijkl μijkl -µ ij P ij µ ij -μkl Pkl μkl (µ ij ] x [P ij µ ij ] y + [ Pkl μkl ] y [ Pkl μkl ] z Notice that: log (B f (x, y) ⊗ B g (y, z) ijkl ) = C f ij + C g klij P ij µ ij + μkl Pkl μkl (μijkl P -1 ijkl μijkl + µ ij P ij µ ij + μkl Pkl μkl µ ij ] x [P ij µ ij ] y + [ Pkl μkl ] y[ Pkl μkl ] z (x, y) ⊗ B g (y, z)) ijkl = exp C ijkl -P

D. 5 -

 5 Proof of Theorem 10Proof. Let P = L T Σ -1 L where L = (F -I) and P λ = P + λI.We have-log p(y|x) = -C Σ + ∥Σ -1/2 (F x + b -y)∥ 2 = -C Σ + ∥Σ -1/2 (Lu + b)∥ 2 (113) = -C Σ + uL T Σ -1 Lu -2b T Σ -1 Lu + b T Σ -1 b (114) If we define µ = P -1 λ β where β = L T Σ -1 b then, -log p(y|x) = -C + ∥P 1/2 λ (u -µ) ∥ 2 (115) = -C + u T P λ u -2µ T P λ u + µ T P λ µ (116) = -C + u T P u + λ∥u∥ 2 -2β T u + β T P -1 λ P λ P -1 λ β (117) = -C + u T P u + λ∥u∥ 2 -2β T u + β T P log(p(y|x)/p(y/x)) = -C Σ + uP u -2β T u + b T Σ -1 b + C -u T P u -λ∥u∥ 2 + 2β T u -β T P -1 λ β (120) = C -C Σ + b T Σ -1 b -λ∥u∥ 2 -β T P -1 λ β (121) Using Woodbury, b T Σ -1 b -β T P -1 λ β = b T Σ -1 -Σ -1 L L T Σ -1 L + λI -1 L T -Σ -1 b (122) = λb T λΣ + LL T -1 b (123) and -log(p(y|x)/p(y/x)) = C -C Σ + λb T λΣ + LL T -1 b -λ∥u∥ 2(124)With C = C Σ -λb T λΣ + LL T -1 b, p(y|x) p(x,y) = e λ∥u∥ 2 and the result follows.

  yP yy y + 2µ x P xy y + 2µ y P yy y -µP µ

			(82)
	where
	P =	P xx P xy P T xy P yy	.

  =x T P xx x + 2yP T xy x + yP yy y -2µ x P xx x -2µ x P xy y -2µ y P xy x -2µ y P yy y + µP (85) =x T P xx x -2(µ x P xx + µ y P xy -yP T xy ν )x + yP yy y -2µ x P xy y -2µ y P yy y + µP µ + yP yy y -2µ x P xy y -2µ y P yy y + µP µ.

		T xy x + yP yy y -2µP	x y	+ µP µ	(84)
					(86)
	= P 1/2 xx	x -P -1 xx ν	2 -νP -1 xx ν

Acknowledgments T.C. gratefully acknowledges support from the French National Agency for Research, grant ANR-18-CE40-0016-01. C.C. acknowledges the support of the Royal Society (grant SPREM RGS\R1\201149) and Amazon.com Inc. (Amazon Research Award -ARA). B.G. acknowledges partial support by the U.S. Army Research Laboratory and the U.S. Army Research Office, and by the U.K. Ministry of Defence and the U.K. Engineering and Physical Sciences Research Council (EPSRC) under grant number EP/R013616/1; B.G. also acknowledges partial support from the French National Agency for Research, grants ANR-18-CE40-0016-01 and ANR-18-CE23-0015-02. A.R. acknowledges partial support from the French government under management of Agence Nationale de la Recherche as part of the "Investissements d'avenir" program, reference ANR-19-P3IA-0001 (PRAIRIE 3IA Institute), and support from the European Research Council (grant REAL 947908).

C Proof of Theorem 6

In this section E = (-1, 1) d and Ω = E × F . Without loss of generality, we assume that F = E. Of course, any compact can be considered.

C.1 Propagation of one-step errors

We generalize the proof technique in [START_REF] François | Stability and uniform approximation of nonlinear filters using the Hilbert metric and application to particle filters[END_REF] to take into consideration general sequences of densities.

Proposition 19. Let π 0 ∈ P(E). Let (z k ) k≥1 and (π k ) k≥1 the optimal filter sequence computed on the (z k ) k≥1 and initialized at π 0 . Let (µ k ) ∈ P(E) N a sequence of distributions such that µ 0 = π 0 . Then, for any n ≥ 0:

with the notation that Rn:

Proof. Telescopic sum:

= µ n -Rn (µ n-1 ) + Rn (µ n-1 -Rn:n-1 (µ n-2 ) . . . + Rn:2 (µ 1 ) -Rn:1 (µ 0 )

Proposition 20 (Optimal filter stability). Let (π ν n ) and (π µ n ) two sequences of optimal filters initialized at ν and µ respectively, and computed on the same data sequence z 1 , . . . , z n . We assume that for all n ≥ 1, R n verifies Assumption 2. Then,

Proof. By (30),

Since R n , . . . , R 1 are mixing, R 1 µ and R 1 ν are comparable and we can apply Proposition 16 with the Hilbert contraction coefficient

And finally, applying (31),

Proposition 21. Let E = (-1, 1) d . Let σ > 0 and ξ ∈ P(E) such that the optimal kernel R n : E × E → R + is σ-ξ mixing for all n. Let µ n a sequence of positive, finite measures on E such that µ 0 = π 0 . Then,

where

Proof. Applying the triangle inequality with the total variation norm to (48),

We apply (30) from Proposition 17,

Since R k+1 and R k are mixing, R k+1 µ k and R k+1 R k µ k-1 are comparable and we can apply Proposition 16 with the Hilbert contraction coefficient

Since R k+1 is mixing, we apply (31) from Proposition 17:

where we use that the Hilbert contraction coefficient is sub-multiplicative and can be bounded away from 1 as a function of σ (we denote it τ σ this upper-bound given in Proposition 16). By denoting δ n = ∥μ n -Rn (μ n-1 )∥ T V we obtain the result.

We now prove Lemma 7.

Proof. We apply the triangle inequality to ∥π ν n -π0 ∥ T V :

and the result follows from Proposition 20 and Proposition 21.

C.2 Bound with smooth, bounded approximations

Let Q : E × E → R + and Ĝ : E × F → R + two bounded approximations of Q and G. Given a sequence (z k ) ∈ F N , define the approximate nonnegative kernel Rn (u, x) = Q(u, x) Ĝ(x, z n ) defined on E × E. We introduce (π n ) the sequence of probability distributions computed using the recursion πn = Rn πn-1 /R n πn-1 (E) = Rn (π n-1 ).

Proposition 22. Let δ n = ∥ Rn (π n-1 ) -Rn (π n-1 )∥ T V . Then,

where C d is a constant independent of Q, G, Ĝ, Q.

Proof. δ n is the total variation distance between two probability distributions. Recall that

for any µ, ν positive, finite measures on E and μ, ν their normalized counterparts. Thus,

∥ Rn (π n-1 ) -R n (π n-1 )∥ T V .

(64)