Guidance for twisted particle filter: a continuous-time perspective
Jianfeng Lu Yuliang Wang


a]Jianfeng LuE-mail:jianfeng@math.duke.edu
b]Yuliang WangEmail:YuliangWang_math@sjtu.edu.cn
[a]Department of Mathematics, Department of Physics, Department of Chemistry, Duke University, Durham, NC 27708, USA.
[b]School of Mathematical Sciences, Institute of Natural Sciences, Shanghai Jiao Tong University, Shanghai, 200240, P.R.China.


Guidance for twisted particle filter: a continuous-time perspective
    [
    
===================================================================






§ ABSTRACT 
The particle filter (PF), also known as the sequential Monte Carlo (SMC), is designed to approximate high-dimensional probability distributions and their normalizing constants in the discrete-time setting. To reduce the variance of the Monte Carlo approximation, several twisted particle filters (TPF) have been proposed by researchers, where one chooses or learns a twisting function that modifies the Markov transition kernel. In this paper, we study the TPF from a continuous-time perspective. Under suitable settings, we show that the discrete-time model converges to a continuous-time limit, which can be solved through a series of well-studied control-based importance sampling algorithms. This discrete-continuous connection allows the design of new TPF algorithms inspired by established continuous-time algorithms. As a concrete example, guided by existing importance sampling algorithms in the continuous-time setting, we propose a novel algorithm called “Twisted-Path Particle Filter" (TPPF), where the twist function, parameterized by neural networks, minimizes specific KL-divergence between path measures. Some numerical experiments are given to illustrate the capability of the proposed algorithm.













§ INTRODUCTION

The particle filter (PF), or the sequential Monte Carlo (SMC), has found a wide range of applications in various areas such as computational statistics and machine learning. It has received increasing popularity when dealing with tasks including statistical inference problems for state space models <cit.> and complex static models <cit.>, capability and safety techniques in large language models (LLM) <cit.>, to name a few. Generally, the particle filter involves the simulation of an artificial particle system over time, particularly suited for estimating statistical quantities of the following form via Monte Carlo:

    Z = 𝔼[∏_k=0^n g_k(X_k)],

where g_k (0≤ k≤ n) are deterministic functions and (X_k)_k=0^n is a discrete-time Markov chain in ℝ^d with transition kernel P(x,dy). The simplest and most classical algorithm to solve (<ref>) is the bootstrap particle algorithm (BPF) proposed by Gordon, Salmond, and Smith in 1993 <cit.> (see Algorithm <ref> below). However, to achieve a desired level of precision, the required sample number N in the Monte Carlo approximation can be rather demanding, since the variance of ∏_k=0^n g_k(X_k) can be large. In order to reduce the variance, various “twisted" particle filters have been developed <cit.>. The general idea of the twisted particle filter (TPF) is as follows: Find a suitable (positive) twisting function φ(k,x) (1≤ k ≤ n, x∈ℝ^d), run the twisted Markon chain X^φ with the transition kernel P^φ_k(x,·) ∼φ(k,·)P(x,·) at the k-th step, and modify the function sequence g_k according to φ (denoted by g_k^φ) (see more details in (<ref>)–(<ref>) below). The twisted model wants to achieve: 

    
  (1) The statistical quantity Z is preserved (i.e. 𝔼[∏_k=0^n g_k(X_k)] = 𝔼[∏_k=0^n g_k^φ(X_k^φ)])
    
  (2) The variance of ∏_k=0^n g_k^φ(X_k^φ) is significantly reduced with a suitable choice of φ.

Notably, this approach is exactly the methodology of the importance sampling <cit.> via a change of measure for discrete path measures, and the key step in constructing a twisted model is seeking the optimal twisting function φ to minimize the variance above. In fact, most existing TPF algorithms adopt the “look ahead" strategies when constructing or learning the twisting function φ (see <cit.> for instance), due to the expression of the optimal twisting function (see (<ref>) and (<ref>) below). 



Similar minimization problems are well-studied in the continuous-time setting, where importance sampling task can be recasted into a stochastic optimal control problem. A fruitful line of research has established several algorithms for the problem <cit.>. Such control-based algorithms have wide range of applications in various areas including molecular dynamics <cit.>, mathematical finance <cit.>, etc. In this paper, for some continuous-time Markov process (X_t)_t≥ 0 and deterministic functions h, g, we are particularly interested in the target quantity of the form

    Z = 𝔼[exp(∫_0^T h(s,X_s)ds ) g(X_T)],

and the corresponding control-based importance sampling algorithms to reduce the variance of exp(∫_0^T h(s,X_s)ds ) g(X_T). See more details in Section <ref>.


Remarkably, although the discrete-time particle filter and the continuous-time control-based importance sampling are typically studied separately by researchers from different communities, they share very similar mathematical nature.  However, to the best of our knowledge, no attempt has been made to connect the twisted particle filter with the control-based importance sampling algorithms rigorously. In this paper, we build a bridge between these two problems under suitable settings for the discrete/continuous-time models (see Section <ref> for more details). Moreover, based on this connection, we propose a novel algorithm called “Twisted-Path Particle Filter" (TPPF), which is directly guided by a family of control-based importance sampling algorithms, as well as parts of their theoretical foundations including the classical Donsker-Varadhan variational principle <cit.>. The proposed TPPF algorithm is expected to be more versatile, less problem-dependent, and behave better in high dimensions (see detailed discussions at the beginning of Section <ref>). Some numerical examples are given in Section <ref> to validate the TPPF and compare it with existing approaches. 






 §.§ Related works




  
Twisted particle filters (TPF). The idea of twisting the model in the discrete-time setting is well established. Various twisted particle filter algorithms have been proposed and studied by researchers from different areas of computational statistics or machine learning. In <cit.>, the authors proposed the so-called “fully-adapted auxiliary particle filter" (FA-APF), where the “look-ahead" strategy was applied to twist the model, and they chose g_k(·) as the twisting function at each discrete time k. In <cit.>, the authors proposed the “iterated auxiliary particle filter" (iAPF) algorithm, where they made use of the recursive relation of the optimal twisting function (see details in (<ref>) below), and learned the twisting function iteratively using Galerkin approximation with Gaussian basis functions. Similar approaches were adopted in <cit.>, where the authors studied the problem from the view of discrete-time optimal control, and made some improvements to the learning structure for better numerical stability. Note that this fixed-point iteration method for approximating the optimal twisting function in <cit.> shares a similar idea with the temporal difference (TD) method <cit.> widely used in the reinforcement learning (RL) community. This fixed-point iteration structure for learning the twisting function has been applied and improved in some other results for twisted particle filters such as <cit.>, where <cit.> made use of the rejection sampling to generate samples from the twisted Markov chain, and <cit.> adopted suitable neural networks to parameterize the twisting function instead of the Galerkin approximation. Besides the TD-type methods, the authors of <cit.> proposed an “optimized auxiliary particle filter" (OAPF) algorithm, where they solve adaptively an convex optimization problem at each discrete time k to update the weights and positions of particles when simulating the discrete Markov chains. 



  
Learning-based TPF algorithms. Most TPF algorithms above are based on optimization over a Galerkin function space. More recently, designed for more practical tasks like large language modeling (LLM) in the machine learning community, some other deep-learning-based structures for approximating the optimal twisting function have been proposed. These learning-based algorithms include the Contrastive Twist Learning (CLT) in <cit.>, the Future Discriminators for Generation (FUDGE) in <cit.>, the SI𝒳O method in <cit.>, to name a few. Basically, in these frameworks, they parameterized the twisting function via a suitable neural network, chose the summation of some specially designed quantities along the time axis as the loss function (for instance, in <cit.>, they selected the KL-divergence at each time k determined by the current twisting function), and trained the network to obtain a good approximation of the optimal twisting functions. See more discussion on the choice of loss functions in Section <ref> below. Compared with Galerkin-based methods, such neural network approximations are expected to have better behavior in high dimensions. 




  
Continuous-time importance sampling (IS) via stochastic optimal control. As mentioned above, this paper connects the (continuous-time) control-based importance sampling and the (discrete-time) particle filters. From the theoretical perspective, the control-based importance sampling algorithms are well-studied in the importance sampling community. Various (continuous-time) importance sampling algorithms have been proposed, such as <cit.>. We are in particular interested in importance sampling algorithms for the path-dependent target (<ref>), and focus on a sequence of control-based ones due to the zero-variance property related to the solution of the stochastic optimal control problem (see for instance <cit.> or Proposition <ref> below). A large amount of control-based importance sampling algorithm admits the structure of the continuous-time policy gradient (PG) method <cit.> (also called the iterative diffusion optimization (IDO) in some other literature <cit.>), and the corresponding details will be discussed later in Section <ref>. There are various classical choices for the loss function in the PG iteration in the existing literature. In the “cross-entropy" algorithm <cit.>, the loss is chosen as KL-divergence between path measures (P^u^* P^u), where P^u is the path measure induced by the controlled SDE adding a control u to the drift, and P^u^* is the path measure corresponding to the optimal control u^*. Note that different from some other statistical inference problems, under the current setting, (P^u^* P^u) has explicit expression and is even quadratic after suitable parameterization for the control function u. Other choices of the loss functions include the relative entropy loss (P^u P^u^*) <cit.>, the variance loss (or log-variance loss) _P^v(dP^u^*dP^u) (or _P^v(logdP^u^*dP^u)) for some suitable basis path measure P^v <cit.>, etc. Some theoretical bounds for the KL-type losses above were established in <cit.>. Besides the PG-based algorithms, other related importance sampling methods include the well-known forward-backward stochastic differential equation (FBSDE) approaches <cit.>, where one approximates the target value Z via the solution of some SDE with given terminal-time state and a forward filtration.



  
Learning-based IS algorithms Many of the conventional methods above may suffer from the curse of dimensionality. In order to solve the related stochastic optimal control problem discussed above, a wide range of learning-based algorithms are available in reinforcement learning for the continuous-time setting. In these frameworks, the learning objective is usually parameterized by neural networks, and over decades a variety of these algorithms have been proposed and studied, such as the classical (soft) policy gradient <cit.>, actor-critic with temporal difference learning <cit.> or with deep backward SDE <cit.>, soft Q-learning <cit.>, etc.
Other improvements to these learning-based frameworks include adding entropy regularization to the optimal control problem <cit.>, and combining the metadynamics algorithms when X_s in (<ref>) is trapped in some metastable region <cit.>.




 §.§ Main contributions

To end the introduction, we summarize the novelty and main contribution of our work here.


  * First, to our knowledge, this paper is the first to rigorously build a bridge between (discrete-time) twisted particle filters and (continuous-time) control-based importance sampling algorithms. This connection is of great significance because through it we can (1) explore the continuous limit of some twisted particle filter algorithm and have a better understanding of its behavior from a continuous-time perspective; (2) propose novel twisted particle filter algorithms based on continuous-time importance sampling algorithms, which are well-studied and proved to have their own advantages when applied to suitable models.

  * Second, our proposed algorithm TPPF treats the KL divergence between path measures as a loss function, and learns the optimal twisting function by neural networks. This approach can overcome the curse of dimensionality in many models compared with other Galerkin-based approximations for the optimal twisting functions. Moreover, TPPF is more robust and has a wider range of applications, because the learning procedure of the twisting function is problem-independent, which is different from other popular Galerkin-based methods <cit.>.




The rest of the paper is organized as follows. 
In Section <ref>, we introduce the basic settings of both discrete and continuous time cases. We also discuss the importance sampling for both models via the twisting function or control variate, respectively. The existence and the zero-variance property of the optimal twisting function / optimal control are also discussed in Section <ref>. In Section <ref>, under suitable assumptions, we prove the convergence from the discrete-time model to the continuous-time model. Based on the connection built in Section <ref>, and motivated by existing algorithms for the continuous-time model, we proposed a novel particle filter algorithm called “Twisted-Path Particle Filter" (TPPF) in Section <ref>. Several numerical examples are given in Section <ref> to compare TPPF with existing approaches. A brief conclusion and some possible future work are summarized in Section <ref>. In the appendix, we provide the proofs of some technical lemmas and propositions.




§ IMPORTANCE SAMPLING FOR DISCRETE-TIME AND CONTINUOUS-TIME MODELS

In this section, we introduce both discrete-time and continuous-time models. For each model, after presenting the basic structure and the expression for the quantity of interest mentioned in (<ref>), (<ref>) above, we introduce the basic idea of variance reduction via importance sampling: Applying Girsanov's transform, we provide an unbiased estimate for the target quantity, and meanwhile the variance could be reduced if we choose a suitable change of measure. As introduced in Section <ref>, this change of measure is in practice achieved by twisting the Markov chain. The twisting procedure is realized by multiplying some suitable function to the original transition kernel in the discrete setting, or adding a control variate to the drift in the continuous setting. Afterwards, for both models, we present the expressions of the optimal twisting function (for discrete setting) / optimal control variate (for continuous setting). The optimal twisting function and the optimal control share similar zero-variance properties, similar backward-time evolution equations, and similar Feynman-Kac representations.




 §.§ Discrete-time model and optimal twisting

Let us begin with the discrete-time model. 
Fix a positive integer n.   Given a (time-homogeneous) Markov transition kernel P̂(x,dy) in ℝ^d (for simplicity we assume the corresponding transition density exists and still denote it by P̂(·,·)), a sequence of bounded, continuous, nonnegative function ĝ_k(·) (0 ≤ k ≤ n), and a (deterministic) initial point x ∈ℝ^d, the general discrete-time model is determined by the triple (P̂(·,·), (ĝ_k(·))_k, x): For the discrete Markov chain X̂_0:n (:= (X̂_0, X̂_1,…,X̂_n)) with transition kernel P̂(·,·) and initial state X̂_0 = x, the statistical quantity of interest is given by

    Z_dis(x) := 𝔼_x[∏_k=0^n ĝ_k(X̂_k)],

where 𝔼_x[·] := 𝔼[ ·|X̂_0 = x].

Such a general model is often called the Feynman-Kac model <cit.>, and the quantity Z_dis(x) corresponds to the terminal marginal measure in it. Note that throughout this paper, we consider the Markov chain with a homogeneous transition kernel and a fixed deterministic initial state. As common in the literature, such settings are only designed to simplify the arguments and notations, and all results in this paper can be easily extended to the Markov chains with inhomogeneous transition kernels and with random initial conditions.


A special case of the general model above is the Hidden Markov Model (a state space model with finite states), where the functions ĝ_0:n(·) are determined on the random observation values y_1:n with an observation conditional density function ĝ_y|x^k(·|X̂_k) at each time k:

    X̂_i ∼P̂(X̂_i-1, ·),   Ŷ_i ∼ĝ_y|x^i(·|X̂_i),   0≤ i ≤ n.

Then, given the observations y_0:n, the functions ĝ_0:n(·) are given by

    ĝ_k(·) := ĝ^k_x|y(· | Ŷ_j = y_j,j≤ k).

Consequently, conditioning on the observations y_0:n, the distribution of X̂_1:n is proportional to

    ĝ_0(x_0)∏_k=1^n ĝ_k(x_k) P̂(x_k-1,x_k) dx_1:n,

and its normaling constant is Z_dis(x) defined in (<ref>). Moreover, in order to compare with the time-continuous model driven by the Brownian model, later on we will consider the Gaussian transition (denoted by P̂^η(·,·)) with a given drift function b(·) and time step η. Namely, the more special case where

    P̂^η(x,dy) := (4πη)^-d/2exp(-1/4η|y-x-η b(x)|^2 ) dy.




Now, with the triple (P̂(·,·), (ĝ_k(·))_k, x) of the discrete model, we aim to calculate Z_dis(x). A standard approach is the particle filter method (also known as sequential Monte Carlo) <cit.>, where the expectation in Z_dis(x) is simulated via Monte Carlo. A classical particle filter algorithm is the following bootstrap particle filter (BPF) <cit.>:


Note that for each k, in (<ref>), particles are resampled according to the weights W_k-1^1:N defined in (<ref>).
An effective way to implement (<ref>) is the ancestor-prediction method (see for instance <cit.>), namely, at k-th step,

    
  * Sample ancestors A_k-1^i ∼𝒞(W_k-1^1, …, W_k-1^N), 1≤ i ≤ N.
    
  * Sample predictions ζ_k^i ∼P̂(ζ^A_k-1^i_k-1,·), 1 ≤ i ≤ N.

Above, 𝒞(·,…,·) denotes the categorial distribution to sample the indexes. Moreover, the resampling step does not need to occur at each step, and one improvement is the κ-adapted resampling <cit.>, where the resampling step only occurs when the effective sample size is smaller than κ N for some fixed κ∈ (0,1). 

For any w^1:N with ∑_i=1^N w^i = 1, the effective sample size ESS(w^1:N) is defined by

    ESS(w^1:N) := 1 / ∑_i=1^N (w^i)^2,

and it is a good standard to evaluate a particle filter algorithm - higher ESS usually means better performance.



One of the shortcomings of BPF is the relatively high variance, especially when the dimension d or the length of the Markov chain n is large.
The twisted particle filter was introduced <cit.> to reduce the variance of Z_dis(x) defined in (<ref>) based on importance sampling.  The general idea is to use the (discrete-time) change of measure, or equivalently, a sequence of twisting functions φ(k,·) (1≤ k ≤ n). Then, we consider the twisted model with the triple (P̂^φ(·,·), (ĝ_k^φ(·))_k, x) defined by

    P̂^φ_k(x,dy) := φ(k,y) /P̂[φ](k,x)P̂(x,dy),  P̂[φ](k,x) = ∫φ(k,y') P̂(x,dy')   1≤ k≤ n,


    ĝ^φ_k(x) := ĝ_k(x) ℓ^φ_k(x),   0≤ k≤ n,

with

    ℓ_0^φ(x) := P̂[φ](1,x),  ℓ_n^φ(x) := 1/φ(n,x),  ℓ^φ_k(x) := P̂[φ](k+1,x)/φ(k,x),   1≤ k≤ n-1.

Then it is not difficult to check that the twisted model preserves the quantity (<ref>). We conclude this property in the following proposition (see Appendix <ref> for a detailed proof):

Consider the twisted model defined in (<ref>)–(<ref>). Recall the quantity Z_dis(x) defined in (<ref>). Then

    Z_dis(x) = 𝔼_x[∏_k=0^nĝ_k^φ(X̂^φ_k)],

where X̂^φ is the twisted Markov chain with the (time-inhomogeneous) transition kernel P̂_k^φ(·,·) , namely, X̂^φ_k∼P̂^φ_k(X̂^φ_k-1, ·) (1≤ k ≤ n). 



Proposition <ref> is in fact the Girsanov's transform in the discrete-time setting. Denote 𝒫̂ the law of the untwisted Markov chain X̂_0:n with the transition kernel P̂(·,·), and denote 𝒫̂^φ the law of the twisted Markov chain X̂^φ_0:n with the transition kernel P̂^φ(·,·). Then, the Girsanov transform gives

    Z_dis(x) = 𝔼_x^X∼𝒫̂[∏_k=0^n ĝ_k(X_k)] = 𝔼_x^X∼𝒫̂^φ[∏_k=0^n ĝ_k(X_k) d𝒫̂/d𝒫̂^φ(X)],

where the Radon-Nikodym derivative is given by

    d𝒫̂/d𝒫̂^φ(X) = ∏_k=0^n ℓ_k^φ(X_k) = ∏_k=1^n P̂[φ](k,X_k-1)/φ(k,X_k).

    


With the twisted model (<ref>) - (<ref>) and Proposition <ref>, it is natural to consider the following twisted particle filter (TPF) method <cit.>, whose output Z^N,φ_dis(x) defined below approximates our target quantity Z_dis(x) due to law of large number (see for instance, Proposition 3 in <cit.>, or Section 3.6 in <cit.>):







Now since the φ-twisted model (recall its definition (<ref>)–(<ref>)) obtains an unbiased estimation for the quantity Z_dis(x) by Proposition <ref>, it is natural to ask the following question: How do we choose φ to minimize the variance. 
In fact, it is possible to choose an optimal twisting function sequence φ^*(k,·) (1≤ k ≤ n), under which the random variable ∏_k=0^nĝ_k^φ^*(X̂^φ^*_k) is an unbiased estimate of Z_dis(x) with zero variance. Optimal twisting functions have been widely studied in the literature <cit.>. In particular, the optimal twisting function sequence φ^*(k,·) (1≤ k ≤ n) is defined recursively as follows:

    φ^*(k,x) = ĝ_k(x) ∫φ^* (k+1,y) P̂(x,dy),   0 ≤ k ≤ n-1,
       φ^*(n,x) = ĝ_n(x).

Then the optimal twisting function has the following useful properties: 

Consider the functions φ^*(k,·) (1≤ k ≤ n) defined by (<ref>). Recall X̂_0:n is a Markov chain with transition density P̂(·,·). Then

    
  * φ^*(k,·) satisfies
    
    φ^*(k,x) = 𝔼[∏_i=k^n ĝ_i(X̂_i) |X̂_k = x],    0 ≤ k ≤ n,    x ∈ℝ^d.

    In particular, Z_dis(x) = φ^*(0,x), ∀ x ∈ℝ^d.
    
  * (zero variance property) 
    Consider the Markov chain X̂^φ^*_0:n with the initial state X̂^φ^*_0 = x and the transition density
    
    P̂^φ^*_k(x,y) := φ^*(k,y) /∫φ^*(k,y') P̂(x,y')dy'P̂(x,y),   1≤ k≤ n,

    Define
    
    Ŵ(X̂^φ^*) := L̂(X̂^φ^*) ∏_k=0^n ĝ_k(X̂_k^φ^*)

    with L̂(X̂^φ^*) := ∏_k=0^nℓ_k^*(X̂^φ^*_k) and
    
    ℓ_0^φ^*(x) := ∫φ^*(1,y)P̂(x,dy),  ℓ_n^φ^*(x) := 1/φ^*(n,x),
       ℓ^φ^*_k(x) := ∫φ^*(k+1,y)P̂(x,dy)/φ^*(k,x),   1≤ k≤ n.

    Then, Ŵ(X̂^φ^*) is an unbiased estimate of Z_dis(x)  with zero variance. Namely,
    
    Ŵ(X̂^φ^*) = φ^*(0,x) = Z_dis(x),   ∀ x ∈ℝ^d.





Proposition <ref> has also been discussed in related literature such as <cit.>. We provide a proof in Appendix <ref> for completeness. Moreover, a direct corollary of the zero-variance property in Proposition <ref> is that, with the optimal twisting function, after Monte Carlo approximation with N samples, the output Z^N,φ^*_dis(x) of twisted particle filter (Algorithm <ref>) is a perfect approximation of the target Z_dis(x). We refer to Appendix <ref> for a detailed proof. 


Consider the optimal twisting function defined in (<ref>). Then for any positive integer N and x ∈ℝ^d,

    Z^N,φ^*_dis(x) = Z_dis(x).







Note that although the φ^*-twisted particle filter provides a perfect approximation for our target Z_dis(x), it is impossible to implement the twisted particle filter (Algorithm <ref>) associated with φ^* in practice. The reason is that we need the pointwise value of the twisting function so that X̂^φ^*_k+1 can be sampled from the distribution proportional to φ^*(k,·)P̂(X̂_k-1^φ^*,·). Therefore, it is crucial to find suitable approaches to approximate the optimal twisting functions φ^*(k,·) (1≤ k ≤ n) to lower the variance of the particle filter algorithm. Moreover, due to the law of large numbers <cit.>, the larger sample size N and the better approximation of φ^* would give a better approximation for Z_dis(x) (i.e. smaller variance of the numerical output Z_dis^N,φ(x)).

In Section <ref>, we will propose a new method to approximate φ^*, guided by similar method for the condinuous-time model described in the following subsection.









 §.§ Continuous-time model and optimal control


In this subsection, we first define the time-continuous model and the target quantity. Then we introduce the control-based importance sampling method based on the Girsanov's theorem. In the end, we introduce the optimal control and its zero-variance property.

For a fixed function b: ℝ^d →ℝ^d, we consider the following SDE

    X_t = X_0 + ∫_0^t b(X_s) ds + √(2)B_t,   X_0 = x∈ℝ^d.

where (B_t)_t≥ 0 is the Brownian motion in ℝ^d under the probability measure P.
The dynamics is thus a time-homogeneous Markov process with transition density (or Green's function) P_t(y|x) satisfying a Fokker-Planck equation

    ∂_t P_t(y|x) = -∇_y · (b(y)   P_t(y|x)) + Δ_y P_t(y|x),   P_0(y|x) = δ_x(y),

where δ_x(·) is the Dirac delta at x. Moreover, given functions h(t,x), g(x) and T>0, the general continuous model is determined by (P;h,g;x), and the corresponding statistical quantity of interest is

    Z_con(x) = 𝔼_x[e^∫_0^T h(s,X_s) dsg(X_T)],



Recall the example in Remark <ref>. Here, a corresponding special case is the latent SDE model, which is a state space model which has infinitely many states. In detail, the functions h(t,x) and g(x) are determined by the value of a random observation process Y_t conditional on X_t, i.e.

    X_t = x + ∫_0^t b(X_s) ds + √(2)B_t,    Y_t ∼ g_y|x^t(·| X_s, s≤ t),   0≤ t≤ T.

Denote the conditional density function (which is the density of Law(X_t | y_s,s≤ t) by

    g_t(·):=g_x|y^t(·|y_s:s≤ t),

and let

    h(t,x) = log g_t(x),    g(x) = g_T(x),   0≤ t≤ T,   x∈ℝ^d.

Then the target quantity Z_con(x) is

    Z_con(x) = 𝔼_x[e^∫_0^T log g_s(X_s) ds g_T(X_T)],

which is comparable with the discrete target Z_dis(x), since Z_dis(x) has the similar expression

    Z_dis(x) = 𝔼_x[∏_k=0^n ĝ_k(X̂_k)] = 𝔼_x[e^∑_k=0^n-1logĝ_k(X̂_k)ĝ_n(X̂_n)].

We will prove the convergence from Z_dis(x) to Z_con(x) in Section <ref> below for the discrete transition kernel of the form (<ref>).

We will give a concrete example in Example <ref> of Section <ref>, where the function g_t(·)=g_x|y^t(·|y_s:s≤ t) has an explicit expression. The example satisfies all the conditions required by the convergence analysis in Section <ref> below.


Now, considering the general continuous model determined by (P;h,g;x), we aim to calculate the quantity Z_con(x). Similar to the discrete model, we approximate it via Monte Carlo using N samples:

    Z_con(x) ≈ Z^N_con(x) := 1/N∑_i=1^N e^∫_0^T log h(s,X^i_s) ds g(X^i_T),

where X^1,… X^N are N independent realizations of the SDE (<ref>).

Next, we aim to reduce the variance of the random variable Z^N_con(x). In fact,

    (Z^N_con(x)) = 1/N (e^∫_0^T log h(s,X_s) ds g(X_T)).

Therefore, we aim to reduce the variance of the random variable e^∫_0^T log h(s,X_s) ds g(X_T) for X satisfying (<ref>). Here, we also consider variance reduction via importance sampling, which is based on a change of measure, as discussed in the discrete case. In detail, the change of measure is constructed by adding a control variate u to the drift and considering the following controlled SDE:

    X^u_t = X_0^u + ∫_0^t (b(X_s^u) + √(2)u(s,X_s^u) ) ds + √(2)W_t,    X_0^u = x ∈ℝ^d.

Then, denoting the path measures 𝒫 := Law(X_[0,T]) (X solving (<ref>)) and 𝒫^u := Law(X^u_[0,T]) (X^u solving (<ref>)), the change of measure gives

    Z_con(x)    = 𝔼_x^X∼𝒫[e^∫_0^T log h(s,X_s) ds g(X_T)]
       = 𝔼_x^X∼𝒫^u[e^∫_0^T log h(s,X_s) ds g(X_T) d𝒫/d𝒫^u(X)] = 𝔼_x[e^∫_0^T log h(s,X^u_s) ds g(X^v_T)d𝒫/d𝒫^u(X^u)],

where the Radon-Nikodym derivative has the following expression due to Girsanov's theorem:

    d𝒫/d𝒫^u(X^u) = exp(-∫_0^T u(s,X^u_s) · dB_s - 1/2∫_0^T |u(s,X^v_s)|^2 ds ).

Consequently, the following Monte-Carlo approximation is an unbiased estimate of the target Z_con(x):

    Z^N,u_con(x) := 1/N∑_i=1^N e^∫_0^T log h(s,X^i,u_s) ds g(X^i,u_T) d𝒫/d𝒫^u(X^i,u),

where X^1,u,…,X^N,u are N independent realizations of the controlled SDE (<ref>).

Our goal is to minimize the variance of e^∫_0^T log h(s,X^u_s) ds g(X^u_T)d𝒫/d𝒫^u(X^u) for u belonging to some suitable admissable set. Luckily, we can find an optimal control (similar to the discrete case) satisfying a zero-variance property, and then we only need to find suitable approaches to approximate the optimal control. In fact, the optimal control u^* is given by

    u^*(t,x) = √(2)∂_x log v^*(t,x),

where v^*(t,x) satisfies the following backward Kolmogorov (parabolic) equation:

    -∂_t v^*(t,x) = b(x) ·∇_x v^*(t,x) + Δ_x v^*(t,x) + h(t,x) v^*(t,x),    0 ≤ t ≤ T,
       v^*(T,x) = g(x).

Moreover, the following Feynman-Kac representation and zero-variance property hold:


Fix T>0 and x∈ℝ^d. Given functions h: ℝ_+×ℝ^d →ℝ, g: ℝ^d →ℝ, b: ℝ^d →ℝ^d,


  * (Feynman-Kac representation) The solution to the PDE (<ref>) has the following representation:

    v^*(t,x) = 𝔼[e^∫^T_t h(s,X_s) ds g(X_T) | X_t = x ],   0≤ t≤ T,   x∈ℝ^d,

where the process X_t solves the SDE (<ref>). In particular, Z_con(x) = v^*(0,x), ∀ x ∈ℝ^d.


  * (zero-variance property) Let P be the probability measure such that under which (B_t)_t≥ 0 is a Brownian motion in ℝ^d. Then there exists a probability measure P^* such that under P^*, the law of the controlled process with control u^* = √(2)∂_x log v^*

    X^u^*_t = X_0^u^* + ∫_0^t (b(X^u^*_s) + 2∂_x log v^*(s,X^u^*_s))ds + √(2)B_t,   X^u^*_0 = x

is the same of the law of X_t under P (recall that X_t satisfies (<ref>)). Moreover, define

    L(T) := exp(-∫_0^T ∂_x log v^*(s,X^u^*_s) dB_s -  1/2∫_0^T |∂_x log v^*(s,X^u^*_s)|^2 ds),

and

    W := e^∫_0^T h(s,X^u^*_s) dsg(X^u^*_T)L(T).

We have 

    W = v^*(0,x) = Z_con(x)    P^*-a.s.










§ CONVERGENCE ANALYSIS: FROM DISCRETE-TIME TO CONTINUOUS-TIME MODELS

In Section <ref> above, we have introduced the importance sampling for discrete/continuous-time models. Now we aim to build a bridge between the discrete-time results and the continuous-time results. The benefit of such connection is evident:
Given an existing importance sampling algorithm in the continuous-time setting (approximating the optimal control), we can design the corresponding particle filter algorithm (approximating the optimal twisting function) in the discrete-time setting. Obviously, this also works in the reverse direction. 
Indeed, we will provide a concrete example for algorithm design in Section <ref>.

Now, in order to study the connection between the discrete and continuous models, we need the following restriction for their transition kernels P, P̂ and the functions determining the targets Z_con(x), Z_dis(x), so that we can establish the convergence results rigorously below. In detail, fix the function b: ℝ^d →ℝ^d, time step η > 0, and T :=nη. We consider the Gaussian transition kernel for the discrete model:

    P̂^η(x,dy) := (4πη)^-d/2exp(-1/4η|y-x-η b(x)|^2 ) dy.

Recall that the transition kernel for the continuous model P_t is associated with the SDE (<ref>) with drift b(·) and volatility √(2), and it satisfies the Fokker-Planck equation (<ref>). Clearly, P̂^η is the one-step transition kernel of the corresponding Euler-Maruyama scheme. Moreover, for convergence analysis in this section, we consider the discrete model determined by (P̂^η; ĝ_0:n(·); x) and the continuous model determined by (P_t;log g(·,·),g_T(·);x). Correspondingly, the statistical quantities of interest are respectively

    Z_dis(x) := 𝔼_x[∏_k=0^n-1ĝ_k(X̂_k)ĝ_n(X̂_n)] = 𝔼_x[exp(∑_k=0^n-1∫_kη^(k+1)ηlog(ĝ_k(X̂_k))^η^-1ds)ĝ_n(X̂_n)],

and

    Z_con(x) = 𝔼_x[exp(∫_0^Tlog g(s,X_s) ds) g_T(X_T)].

In what follows, under suitable assumptions for the functions b(·), ĝ_k(·), g(·,·) and g_T(·), we will consider the convergence from discrete-time model to the continuous-time model. In detail, we will show that at the time step η→ 0, P̂^η converges to P_t, and Z_dis(x) converges to Z_con(x).

We need the following assumptions to ensure convergence:

We assume the following conditions for b:ℝ^d →ℝ^d, ĝ_k: ℝ^d →ℝ (0≤ k ≤ n), g:ℝ_+×ℝ^d →ℝ and g_T: ℝ^d →ℝ:

    
  (a) b(·) is L_b-Lipschitz (i.e. |b(x) - b(y)|≤ L_b |x-y|, ∀ x,y∈ℝ^d), and sup_0≤ i ≤ n𝔼|X̂_iη|^2 < ∞ uniform in η.
    
  (b) ĝ_n(x) → g_T(x), η^-1logĝ_k(x) →log g(kη, x) uniformly in x and k (0 ≤ k ≤ n-1) as η→ 0.
    
  (c) log g(t, x) is continuous in t ∈ [0,T] uniform in x. η^-1logĝ_kη(x) is Lipschitz in x uniformly for all k (0 ≤ k ≤ n-1) and η>0.
    
  (d)  ĝ_n(·) is bounded. For any t ∈ [0,T], the functional X_[t,T]↦ e^∫_t^T log g(τ, X_τ)dτ is continuous and bounded.


Above, condition (a) is used to ensure the existence and uniqueness of the strong solution for (<ref>) and the convergence of transition kernel in Proposition <ref> below. Other conditions are required for the convergence of Z_dis(x) to Z_con(x) in Proposition <ref>  below. Also, we will give a concrete example in Example <ref> below which satisfies all the conditions in Assumption <ref>.

We also remark here that the upper bound for the second moment sup_0≤ i ≤ n𝔼|X̂_iη|^2 can be proved if one assumes: (1) L_b-Lipschitz condition for b(·); (2) the following confining condition for b(·):

    x · b(x) ≤ -C_1 |x|^2 + C_2,   ∀ x ∈ℝ^d,

where C_1, C_2 are two positive constants. Moreover, under such assumptions, the upper bound sup_0≤ i ≤ n𝔼|X̂_iη|^2 can be shown to be independent of the time T (recall that T = nη). 

In the following proposition, we establish the convergence of transition kernels in terms of the KL-divergence or the total variation (TV) distance. Note that for two probability measures μ, ν on ℝ^d, the KL-divergence and TV distance are given by

    (μν) :=
    {   ∫_ℝ^dlogdμ/dν μ(dx),  μ≪ν,
       +∞,  otherwise.
    .


    TV(μ,ν) := sup_A∈ℬ(ℝ^d) |μ(A) - ν(A)|.



Consider the Markov transition kernels P̂^η and P_η (=P_t=η) defined in (<ref>), (<ref>), respectively. Suppose that the condition (a) in Assumption <ref> holds. Then for any x∈ℝ^d, as η→ 0, P̂^η(x,·) → P_η(x,·) in terms of KL-divergence or TV distance. In detail, for η < 1, there exists a positive constant C independent of η and d such that

    ( P̂^η(x,·)   P_η(x,·)   ) ≤ Cdη^2 → 0   as  η→ 0,

and

    TV( P̂^η(x,·)  , P_η(x,·)   ) ≤ C√(d)η→ 0   as  η→ 0.



Note that the Proposition considers the local truncation error of the Euler-Maruyama scheme in terms of TV distance and KL divergence, consistent with the existing results <cit.>: consider the solution X_t (t ≥ 0) to the continuous-time SDE (<ref>) and its Euler-Maruyama discretization X̂_nη (n ∈ℕ_+) with time step η, for any time interval length T such that η divides T, one can show that 

    D_KL(Law(X̂_T)Law(X_T)) ≲ Tη.

This then reduces to (<ref>) when T = η. Moreover, recent literature (for instance, <cit.>) has shown that, if one assumes stronger smoothness conditions for the drift b(·), it is possible to improve the upper bound for the time-discretization error in terms of the KL-divergence from O(Tη) to O(Tη^2). Consequently, the rate in (<ref>) is O(η^3). The same arguments also hold for TV distance via Pinsker's inequality. 


The next Proposition guarantees that under the settings of Assumption <ref>, the target Z_dis(x) converges to Z_con(x) as η→ 0.


For x ∈ℝ^d, recall the definitions of Z_dis(x), Z_con(x) in (<ref>), (<ref>), respectively. Then under Assumption <ref>,

    Z_dis(x) → Z_con(x)

pointwisely as η→ 0.






Under current assumptions, one cannot obtain an explicit convergence rate due to conditions (b) and (c) in Assumption <ref>. However, if we use a stronger version Assumptions <ref>, it is possible to obtain an explicit rate. In fact, if we replace conditions (b) and (c) in Assumption <ref> by the following stronger quantitative version 

    
  (b') g_T(·) is bounded. |ĝ_n(x) - g_T(x) | ≤ C_1η^α_1, |η^-1logĝ_kη(x) - log g_kη(x) | ≤ C_2η^α_2, 0 ≤ k ≤ n-1, where C_1, C_2, α_1, α_2 are positive constants independent of x, k.
    
  (c')log g_s(x) is α_3-Hölder continuous in t ∈ [0,T] uniform in x. η^-1logĝ_kη(x) is Lipschitz in x uniformly for all k.

then, following exactly the same derivation, one can easily obtain for small η,

    |Z_dis(x) - Z_con(x)| ≤η^α,

where α := min (1/2, α_1, α_2, α_3).


We end this section by giving a concrete example satisfying all the conditions in Assumption <ref>. In particular, it satisfies the time-continuity condition of log g(t,x), which might not be very direct at first glance.
[an example satisfying time-continuity and other assumptions]
Note that one important assumption is that the function log g(t,x) is time-continuous. Recall the state space models in Remark <ref> and Remark <ref>. In such model, if the observation Y_t at time t only depends on the value of X_t (for instance, Y_t ∼ N(X_t, Σ) with Σ being a constant covariance matrix), then the time-continuity for the function g^t_x|y would not hold, since the observation process Y itself is not continuous in time due to the randomness. However, we remark in the following that by considering time-dependency for g^t_y|x (recall its definition in (<ref>)), we can indeed find such time-continuous log g(t,x) 


In fact, we consider the following continuous-time state space model:

    X_t = x + ∫_0^t b(X_s) ds + √(2)B_t,    Y_t = X_t + B̃_t+t_0,    0≤ t≤ T

for two independent Brownian motions (B_t)_t≥ 0, (B̃_t)_t≥ 0. Here we introduce the positive constant t_0 only to avoid a potential singularity of 1/t at t=0: As will be discussed below, the convergence would require the time continuity of - |x - y_t|^2/2(t+t_0) in the closed interval [0,T], thus the choice of some t_0 > 0.

Correspondingly, for the discrete model, for 0≤ k≤ n, we define Ŷ_k by

    X̂_k+1 = X̂_k + η b(X̂_k)  + √(2)(B_(k+1)η - B_kη),  Ŷ_k= X̂_k + B̃_kη+t_0,   0≤ k ≤ n.


Now suppose that we have been given observation y_t (0≤ t≤ T) and ŷ_k (0≤ k ≤ n) (note that by our construction above, the observed process y_t of the continuous model is time-continuous). Clearly, the expressions for g(t,x) and ĝ_k(x) are explicit and are of the Gaussian form. Therefore in this case, we are able to convert the conditions in Assumption <ref> for g and ĝ to assumptions for y and ŷ. Basically, defining 𝒩(·;μ,Σ) to be the density of d-dimension Gaussian distribution with mean μ and covariance Σ, we need the followings:

    
  * 𝒩(x;ŷ_n,T+t_0) →𝒩(x;y_n,T+t_0), |x - ŷ_k|^2 → |x-y_k|^2 uniformly in x, k as η→ 0.
    
  * For any t ∈ [0,T], the functional X_[t,T]↦ e^∫_t^T  -|X_τ - y_τ|^2/2(t+t_0)dτ is bounded.
    
  * - |x - y_t|^2/2(t+t_0) is continuous in t ∈ [0,T] uniform in x. - |x - ŷ_k|^2/2(kη+t_0) is Lipschitz in x uniformly for all k.

It is then easy to see that these conditions generally require (1) time-continuity for the observation y_t, (2) the observation for the two models y_kη, ŷ_k are very close. Clearly, since X̂ in (<ref>) is the Euler-Maruyama discretization of X in (<ref>), the continuous observed processes y, ŷ stay close when the time step η is small. 

Furthermore, once these requirements are met, by Proposition <ref>, as η→ 0, the rescaled target

    Z̃_dis(x) := 𝔼_x[∏_k=0^n-1 (ĝ_k(X̂_k))^η^-1ĝ_n(X̂_n)]

converges pointwise to

    Z_con(x) = 𝔼_x[exp(∫_0^Tlog g(s,X_s) ds) g_T(X_T)].




As a final remark in this section, the time-continuity of the observation (in particular, the form (<ref>)) is only required when proving the convergence from the discrete model to the continuous model. This discrete-continuous bridge then allows us to construct new twisted particle filters motivated by existing continuous-time importance sampling algorithms (e.g., see Section <ref>). However, since we construct algorithms directly for the discrete setting, the new algorithm does not necessarily require such time-continuity condition of the observation. In fact, the proposed algorithms work for simple linear Gaussian models (see Section <ref>), whose observation can be defined by Y_k ∼ N(X_k, Σ) with Σ being a constant covariance matrix.


































§ TWISTED-PATH PARTICLE FILTER (TPPF)


In this section, we give a concrete example of how to explore new twisted particle filters guided by continuous-time importance sampling algorithms.
Recall that in Section <ref>, we have proved the convergence from the discrete model to the continuous model under suitable settings. As has been mentioned at the beginning of Section <ref>, this connection we found can guide us to explore novel particle filter algorithms based on existing methods in the continuous-time setting. In this section, motivated by mature algorithms from the continuous-time community, we propose a Twisted-Path Particle Filter (TPPF). In this new algorithm for the discrete-time model, we treat the KL-divergence between twisted path measures and the optimal path measures as the loss function, parameterize the twisting function via neural networks, and approximate the optimal twisting function via suitable optimization methods such as the stochastic gradient descent. We expect this proposed algorithm to have the following strengths:

    
  * By parameterizing the twisting function using neural networks, TPPF is expected to overcome the curse of dimensionality compared with algorithms based on Galerkin approximation.
    
  * Due to the universal approximation theorem of neural networks <cit.>, TPPF is more versatile and less problem-dependent, while in some other existing algorithms, the twisting function is learnt within a problem-dependent function class.
    
  * Compared with some other frameworks for learning the twisting function (see <cit.> for instance), which are based on the backward recursive relation (<ref>), our method has relatively stronger theoretical foundations (see Lemma <ref>, Proposition <ref>, Corollary <ref>  below).


In the rest of this section, we will first explain the motivation of TPPF by presenting one importance sampling algorithm mainly from a variational perspective. Then we will formulate the details of the TPPF algorithm under the guidance of the existing continuous-time algorithm. Some theoretical results as well as implementation details will also be provided.





 §.§ Motivation: importance sampling via variational characterization in continuous-time setting

In the continuous-time model, it has been relatively well-studied to find the optimal control variate from a variantional perspective. In order to seek the optimal control variate with the zero variance property, instead of solving the PDE we derived in (<ref>), people convert this to an optimal control problem from a variational perspective <cit.>. We first give some basis on the so-called Donsker-Varadhan variational principle <cit.>, which is independent of whether the model is time-continuous or time-discrete. 


Given some function W(·): Ω→ℝ and probability measure P on Ω, the following relation holds:

    -log𝔼^X ∼ P[exp(-W(X) )] = inf_Q ∈𝒫(Ω){𝔼^X∼ Q[W(X)] + (Q  P)},

where the minimum is over all probability measures on Ω.



Take Ω to be the path space corresponding to the continuous-time model and choose W(·) to be the path integral for some stochastic process (X_s)_0≤ s≤ T, as in <cit.>

    W(X) := -∫_0^T h(s,X_s) ds - log g(X_T)

for some h(·,·), g(·) defined in Section <ref>, the variational relation (<ref>) becomes

    -log𝔼^X_[0,T]∼ P[exp(∫_0^T h(s,X_s) ds + log g(X_T) )]
    
        =inf_Q∈𝒫(Ω),Q<<P𝔼^X_[0,T]∼ Q[∫_0^T -h(s,X_s) ds - log g(X_T) ] + (QP).

In particular, if the path measure P above is the law of solution to some X_[0,T] satisfies an SDE as introduced in Section <ref>

    dX_t = b(X_t)dt + √(2)dW,    0 ≤ t ≤ T,    X_0 = x,

then the probability Q above is characterized by the path measure of the controlled SDE after adding a control variate u to the drift. In detail, consider the controlled the SDE

    dX^u_t = b(X^u_t)dt + √(2)u(t,X_t^u)dt + √(2)dW,    0 ≤ t ≤ T,    X^u_0 = x.

By Girsanov's theorem, the dual relation (<ref>) becomes

    -log𝔼_x[exp(∫_0^T h_s(X_s) ds + log g_T(X_T) )]
    
        = inf_u∈𝒰𝔼_x[-∫_0^T h_s(X^u_s) ds + 1/2∫_0^T |u_s(X_s^u)|^2 ds - log g(X_T^u)],

and under the setting of Section <ref>, the set of admissible controls 𝒰 is usually chosen to be as follows (see for instance Section 1 in <cit.>)

    𝒰={u ∈ C^1(ℝ^d ×[0, T], ℝ^d):  u  grows at most linearly in  x }.

Denote

    J(u) := 𝔼_x[-∫_0^T h_s(X^u_s) ds + 1/2∫_0^T |u_s(X_s^u)|^2 ds - log g(X_T^u)].

It can be verified that the minimizer u^* of the functional J(u) yields a zero-variance importance sampler in Proposition <ref> for the continuous setting, namely,

    exp(∫_0^T h(s,X^u^*_s) ds + log g(X^u^*_T) ) dP/dP^u^*(X^u^*) = Z_con(x),   P^u^* - a.s.,

where P^u^* is the path measure induced by the controlled process X^u^* in (<ref>) associated with the optimal control u^*, and P is the path measure induced by the original process X in (<ref>). Therefore, it is reasonable to find the optimal control variates by solving the following stochastic optimal problem:

    min_u ∈𝒰 J(u) = 𝔼_x[-∫_0^T h_s(X^u_s) ds + 1/2∫_0^T |u_s(X_s^u)|^2 ds - log g(X_T^u)].

In practice, one may parameterize u(t,x) by u = u(θ;t,x) and iteratively solve the optimization problem. In details, given a suitable loss function L(u) = L(F(X^u_[0,T],u)) (e.g., L = J) in each iteration, one implement the followings:

    
  * With the current control u(t,x) = u(θ;t,x) simulate N realizations of the controlled SDE X^u_[0,T], and calculate the loss L and its derivative ∇_θ L using the N realizations of X^u_[0,T].
    
  * Update θ using ∇_θ L via some suitable method, for instance, the stochastic gradient descent.

As mentioned in Section <ref>, such an iterative framework is exactly the continuous-time policy gradient (PG) method <cit.>, and it is also called the iterative diffusion optimization (IDO) in some other literature <cit.>. As a remark, J(u) of the form (<ref>) offers a good choice of the loss L, which equals to (P^uP^u^*) up to a constant. Moreover, people also studied other forms of loss to approximate the optimal control u^*, for instance, in the so-called cross-entropy method <cit.>, people use the loss (P^u^*P^u) instead. More discussion on the choice of loss function and detailed derivations would be given in the next subsection for the TPPF algorithm, where we also treat the KL-divergence between path measures as the loss function.



 §.§ Approximating the optimal twisting function in discrete-time setting

Now we propose a novel framework to approximate the optimal twisting function φ^* defined in (<ref>), guided by the existing method for training u^* in the continuous-time setting. Different from the iteration method proposed in <cit.> (which is similar with the temporal difference (TD) learning in the reinforcement learning community <cit.>) based on the recursive formula in (<ref>), our method learns the optimal twisting function in the whole path via a neural network. Hence, we name the proposed algorithm the “Twisted-Path Particle Filter" (TPPF). As discussed at the beginning of Section <ref>, we expect the proposed method TPPF to (1) perform better in high dimensions; (2) be more robust and have wider applications; (3) have relatively stronger theoretical foundations.


In what follows, we first derive a similar variational principle and some relations between different choices of loss induced by the twisting function in the discrete-time setting. After that we would give our detailed algorithms corresponds to possible choices of loss function. Some detailed formulas and implement details will also be discussed at the end of this section.

First, note that the Donsker-Varadhan variational principle can also be applied to the discrete-time model in Section <ref>. In fact, after a change of measure, we have

    -log Z_dis(x) = inf_φ J(φ),

where φ = (φ_1,…,φ_n ), φ_i(·) (1 ≤ i ≤ n) are positive, continuous functions in ℝ^d, and as a direct result of the dual relation, by choosing W(X̂) = -∑_k=0^nlogĝ_k(X̂_k) for some discrete Markov chain (X̂_k)_k=0^n

    J(φ) := 𝔼_x[-∑_k=0^nlogĝ_k(X̂_k)] + (P^φP^1),

where P^φ is the path measure induced by the twisting function φ as in (<ref>), and the P^1 is the original path measure (or equivalently, with twisting function being 1). More precisely, under P^φ, the Markov chain X^φ evolves under the twisted transition density

    P̂^φ_k(x_k-1,x_k) = φ(k,x_k)P̂(x_k-1,x_k)/P̂[φ](k,x_k-1),   1 ≤ k ≤ n,

with the normalizing constant P̂[φ](k,x_k-1) defined by

    P̂[φ](k,x_k-1) = ∫φ(k,y)P̂(x_k-1,y)dy,

and under P the Markov chain X evolves under the untwisted transition density P̂(x_k-1,x_k). Consequently, the twisted path measure associated with the twisting function φ has the following explicit expression:

    P^φ(dx_1:n) = ∏_k=1^n φ(k,x_k)P̂(x_k-1,x_k)/P̂[φ](k,x_k-1) dx_1:n.


Similarly as in the continuous-time setting, it can be verified that the optimal twisting function φ^* with zero-variance property (discussed in Proposition <ref>) is the minimizer of the functional J(φ) above, and minimizing the functional J(φ) is equivalent to minimizing the KL-divergence (P^φP^φ^*). We summarize this result in the Proposition below:

Consider the functional J(φ) defined in (<ref>) and the optimal twisting function φ^* defined in (<ref>). Then it holds

    J(φ) = J(φ^*) + (P^φP^φ^*).

Moreover, one can derive the following explicit formula for J(φ) and (P^φP^φ^*):

    J(φ) = 𝔼_x[-∑_k=0^n logĝ_k(X̂_k^φ) + ∑_k=1^n logφ(k,X̂_k^φ)/P̂[φ](k,X^φ_k-1)]

and

    (P^φP^φ^*) = 𝔼_x[-∑_k=0^n logĝ_k(X̂_k^φ) + ∑_k=1^n logφ(k,X̂_k^φ)/P̂[φ](k,X^φ_k-1)] + logφ^*_0(x)





After a change of measure, it is easy to see that J(φ) or (P^φP^φ^*) has an alternative expression:

    (P^φP^φ^*) = 𝔼_x[exp(∑_k=1^nlogφ_k(X_k)/φ̃_k-1(X_k-1))
    (-∑_k=0^n log g_k(X_k) + ∑_k=1^n logφ_k(X_k)/φ̃_k-1(X_k-1))] + logφ^*_0(x).

This expression is in particularly useful in the cases where sampling from the original transition P̂(x_k-1,x_k) is much easier than sampling from the twisted transition P̂_k^φ(x_k-1,x_k). For instance, if φ is parameterized by some neural network, basic sampling methods like rejecting sampler (proposed in <cit.> for the twisted particle filter) would not be so efficient, suffering from low, unstable acceptance rate. Moreover, it is natural to doubt whether the Monte Carlo approximation for the loss (and its gradient) would deteriorate after a change of measure. So far we have got no theoretical guarantee for this, but empirically, the experiments show that (<ref>) can approximate the loss well. 


Guided by various variational-based methods in continuous-time setting (discussed in Section <ref>), we can also consider different loss functions other than J(φ) or (P^φP^φ^*). Before seeking blindly for other possible choices, let us first study the relationship between the relative variance and KL-divergence, since our final goal of approximating the optimal twisting function is just to reduce the variance. In fact, using a generalized Jensen's inequality, we are able to derive the following:

Given some function W(·): Ω→ℝ and probability measures ν, ν̃ on Ω which are absolutely continuous with each other, define

    Z := 𝔼^X ∼ν[e^-W(X)] = 𝔼^X̃∼ν̃[e^-W(X̃)dν/dν̃(X̃)]

and the relative variance with respect to ν̃ defined by

    r(ν̃) :=√(Var_ν̃(e^-Wdν/dν̃))/Z

Suppose there is an optimal probability measure ν^* with the zero-variance property

    dν^*/dν = e^-W/Z  ν-a.s.,

Then the following estimates hold:

    
  * r^2(ν̃) ≥ e^(ν^*ν̃)-1
    
  * If the constants m:=inf_E ν^*(E)/ν̃(E), M = sup_E ν^*(E)/ν̃(E) exist, then
    
    e^m(ν̃ν^*) + (ν^*ν̃)-1 ≤ r^2(ν̃) ≤ e^M(ν̃ν^*) + (ν^*ν̃)-1









Consider the discrete-time model with (P̂(·,·); ĝ_0:n(·);x), we have the following direct corollary of Proposition <ref>.

Take Ω = ℝ^n, W(X̂) = -∑_k=0^nlogĝ_k(X̂_k), ν = P^1, ν̃ = P^φ, and ν^* = P^φ^* (recall the definition of P^φ in (<ref>)) in Proposition <ref>. Then there exists 0< m ≤ M depending on φ, φ^* such that

    e^m(P^φ P^φ^*) + (P^φ^* P^φ)-1 ≤ r^2(ν̃) ≤ e^M(P^φ P^φ^*) + (P^φ^* P^φ)-1.




Motivated by Proposition <ref> and Corollary <ref> above, it is then reasonable to consider the loss function of the form 

    a(P^φ P^φ^*) + ( P^φ^* P^φ)

for some positive a. Moreover, when P^φ approximates the target P^φ^* well, m, M in Proposition <ref> are approximately 1, so it is reasonable to choose a=1, namely, treat 

    (P^φ P^φ^*) + (P^φ^* P^φ)

as a loss function for ν̃. Of course with these bounds, it is also reasonable to consider the loss

    (P^φ^* P^φ),

which corresponds to the loss in the well-studied cross-entropy method in the continuous-time settings <cit.> (and this is the reason why we denote it by L_CE). When implementing the TPPF algorithm, we need explicit expressions for the loss functions chosen above, so that we can calculate them via Monte Carlo approximation. In fact, denoting

    L_RE:=(P^φ P^φ^*),   L_CE := (P^φ^* P^φ),   L_RECE = L_RE+L_CE.

Simple calculations yield

    L_RE    =𝔼_x[-∑_k=0^n logĝ_k(X̂_k^φ) + ∑_k=1^n logφ(k,X̂_k^φ)/P̂[φ](k,X̂_k-1^φ)] + logφ^*(0,x)
       = 𝔼_x[exp(∑_k=1^nlogφ(k,X̂_k)/P̂[φ](k,X̂_k-1))(-∑_k=0^n logĝ_k(X̂_k) + ∑_k=1^n logφ(k,X̂_k)/P̂[φ](k,X̂_k-1))] 
            + logφ^*(0,x)  ,

where we can use either the first line (via the twisted Markov chain X̂^φ) or the second line (via untwisted Markov chain X̂) in the training. See Section <ref> for more details. Also, for L_CE, we have

    L_CE   = 1/φ_0^*(x)𝔼_x[exp(∑_k=0^n logĝ_k(X̂_k) ) (∑_k=1^n logĝ_k(X̂_k) - logφ^*(0,x) )]  (=constant)
          - 1/φ_0^*(x)𝔼_x[exp(∑_k=0^n logĝ_k(X̂_k) ) (∑_k=1^n logφ(k,X̂_k)/P̂[φ](k,X̂_k-1))],





Note that the coefficient 

    exp(∑_k=0^n logĝ_k(X̂_k) ) / φ_0^*(x)

in second line of (<ref>) may bring numerical instability to L_CE-training. First, it is possible that in some extreme cases, we have no idea of the value of φ^*(0,x) (which is our target Z_dis(x)), and sometimes we can't even get a good approximation of it. Second, this coefficient above is determined by ĝ_k and the (random) position of X̂_k, so it is numerically unstable especially when the length of the Markov Chain n is large. Moreover, this observation is also empirically true in some examples, see details in Section <ref>.

















Given the loss functions above, to approximate optimal twisting function φ^* in practice, it is remaining to parameterize φ and learn the parameters according to the loss functions. In our algorithm, we parameterize the twisting function by neural network φ_k(x) = φ(θ;k,x), where θ denotes all the parameters of the network, and k, x are the input. For numerical stability, it is better to treat logφ(θ;k,x)=NN(θ;k,x) as the output so that φ_k(x) = exp(NN(θ;k,x)). We refer to Section <ref> for more details and give the structure of the proposed TPPF in Algorithm <ref> below:





As a final remark, recall that when implementing Algorithm <ref>, we need to calculate the gradient ∇_θ L via Monte Carlo approximation. Note that when the loss function L is approximated via the untwisted Markov Chain X̂, the gradient can be calculated via auto-differentiation. However, if the loss L is approximated via the twisted Markov Chain X̂^φ (see the first line in (<ref>)), the computation of ∇_θ L is not that direct, since X̂^φ depends on the parameter θ. Here, we provide the explicit formula for ∂_θ_i L_RE, where L_RE is given in the first line of (<ref>). In fact, we first write L_RE into the form of the second line in (<ref>), where X̂ in it is independent of θ. Then, after a change of measure, we can write the gradient into the following:

    ∂_θ_iL_RE 
        =𝔼_x[(1-∑_k=1^n logĝ_k(X̂_k^φ) + ∑_k=1^n logφ(k,X̂_k^φ)/P̂[φ](k,X̂_k-1^φ))(∑_k=1^n log∂_θ_iφ(k,X̂_k^φ)/P̂[∂_θ_iφ](k,X̂_k-1^φ))]   .

For a further implementation detail, to obtain ∂_θ_iL_RE in (<ref>), there is no need to calculate ∂_θ_iφ. Instead, we can first calculate the value of

    𝔼_x[(1-∑_k=1^n logĝ_k(X̂_k^φ) + ∑_k=1^n logφ(k,X̂_k^φ)/P̂[φ](k,X̂_k-1^φ))(∑_k=1^n logφ(k,X̂_k^φ)/P̂[φ](k,X̂_k-1^φ))]   .

using Monte Carlo approximation, and then “detach" the first term above so that it would contain no gradient information. Consequently, auto-differentiation via (one-time) back-propagation gives the desired value of the gradient in (<ref>). 







§ NUMERICAL EXAMPLES

In this section, we test the proposed TPPF algorithm with loss function L_RE, L_CE, or L_RECE on different models including the linear Gaussian model, the Lorenz-96 model, and the stochastic volatility model. We compare our algorithm with well-known competitors including the bootstrap particle filter (BPF) <cit.>, the iterated auxiliary particle filter (iAPF) <cit.>, and the fully-adapted auxiliary particle filter(FA-APF) <cit.>. In most examples, TPPF performs better than other algorithms. 




 §.§ Implementation details


Before the main numerical experiment, let us give some details on how to parameterize the twisting function φ, how to sample from the twisted Markov transition kernel P̂^φ, and how to calculate the normalizing constant P̂[φ](k,x). In fact, in the following three experiments, We consider two ways of parameterization: the robust non-parametric way and the problem-dependent parametric way. In particular, for the linear Gaussian model, we use the problem-dependent implementation, while for the other two models we use the non-parametric implementation.



  * The non-parametric implementation. As discussed in Section <ref>, we approximate the twisting function by

    logφ_k(x) = NN(θ;k,x),

where NN(θ;k,x) is a neural network with parameters θ and inputs k, x (0≤ k ≤ n, x ∈ℝ^d). In our experiments, we set NN as DenseNet <cit.> with two hidden layers. Moreover, we add a tanh activation to the final layer so that the output NN(θ;k,x) takes value in (ϵ,1), where ϵ>0 is a hyperparameter. 

This bounded restriction is designed for the following reject sampling step when sampling from the twisted kernel P̂^φ. In fact, since we usually do not have much information of the current twisting function φ, so it is not easy to sample from the twisting kernel P̂^φ_k(x,·) ∼φ(k,·)P̂(x,·). Under the robust implementation setting, we make use of the reject sampling recently proposed in <cit.>: Using X^φ_k and the untwisted transition kernel P̂(X^φ_k,·), propose a new position X_pro, accept it with probability φ(X_pro). Repeat until first acceptance.

Moreover, in the non-parametric implementation setting, the normalizing constant is calculated via Monte Carlo approximation using Ñ samples (following <cit.>, we choose Ñ = 50 in both Lorenz-96 and NGM-78 models):

    P̂[φ](k,x) := ∫φ(k,y) P̂(x,dy) ≈1/Ñ∑_i=1^Ñφ(k,U_i),   U_i ∼P̂(x,·)   i.i.d.


Also, to make the training faster, we use the untwisted process X instead of X^φ when calculating the loss and the gradient (i.e. for L_RE, the loss is computed via the first line in (<ref>) and its gradient is computed using (<ref>)).


  * The parametric implementation. In the linear Gaussian model, we can calculate the analytical solution of Z_dis(x) using the Kalman filter <cit.>. Moreover, we can analytically calculate the optimal twisting function using the backward recursive relation (<ref>), and clearly the optimal twisting function is also Gaussian. Therefore, with so much knowledge of the solution, it is reasonable to consider a problem-dependent way to parameterize the twisting function to make the learning more efficient. In more details, via a mean-variance estimation framework, we set

    μ_k = NN_1(θ_1;k) ∈ℝ^d,  σ_k^2 = NN_2(θ_2;k) ∈ℝ_+,   0≤ k ≤ n.

And then set

    φ(k,x) = C_k N(x;μ_k,σ^2_k).

Note that the twisting function is invariant of the constant scaling, so we do not care about the constant that multiplies the Gaussian. This means there is no need to learn C_k in (<ref>), and in our experiment we just consider C_k = (2πσ_k^2)^-d/2 so that φ(k,x) = exp(-|x-μ_k|^2/2σ^2_k). Also, under such settings, we no longer need the reject sampling and the inner-loop Monte Carlo for calculating the normalizing constant, because everything can be calculated analytically. Consequently, compared with the robust way, the problem-dependent implementation is less time consuming and the neural network is easier to train.



Another important remark is about the relative fairness of the comparisons in the numerical experiments. We learn and run the particle filter using the same particle numbers. For the time complexity, we admit that as a training-based algorithm, our algorithm is more time-consuming, and in fact, there is not a completely fair comparison due to different non-optimal choice of network structures. In order to conduct a relatively fair comparison, we restrict the number of iterations for the training so that the total running time is comparable and similar to its competitors.












 §.§ Linear Gaussian model

As a first example, we consider the linear Gaussian model, which is also the discretization of an OU process X_t with linear Gaussian observations Y_k ∼ N(·;B_k X_k,Σ_OB)

    X_k+1 = X_k + Δ t A X_k  + √(Δ_t) N(0,Σ).

We choose Δ t = 0.01, T = 0.5 (so the length of the discrete Markov Chain is n = T/Δ t = 50), B_k = I_d, A = -I_d, σ = I_d, Σ_OB = I_d.

We consider the parametric implementation with configurations d ∈{2,5,15,20 }. The boxplots in Figure <ref> compare the TPPF (with L_RE, L_CE, or L_CERE) with competitors including BPF, iAPF and FA-APF using 1000 replicates. The red cross represents the mean and the red dash line represents the medium. We also report the empirical standard deviations in Table <ref>, and in Table <ref> we report the average relative effective sample size (ESS-r) defined by ESS(W_1:N) = 1 / (N∑_i^N W_i^2) with ∑_i=1^N W_i = 1.
Clearly, TPPF with L_RE and L_RECE can defeat BPF and FA-APF, especially when the dimension is high, and TPPF with L_CE suffers from the curse of dimensionality due to the numerical instability discussed in Remark <ref>. Moreover, the iAPF performs the best in this linear Gaussian model, because iAPF is learning the twisting function in a Gaussian function class, each time solving a standard restricted least square minimization problem. Therefore, our training-based method cannot perform as well as iAPF in this experiment. However, as we will see in the other experiments when the optimal twisting function is not in the Gaussian family, our algorithm performs better than iAPF.















 §.§ NGM-78 model

In contrast to the linear nature of the last model, here we consider an (artificial) nonlinear model frequently used when testing the performance of particle filters <cit.>. To our knowledge, this model was first used by Netto, Gimeno, and Mendes in 1978 <cit.>, so here we name it NGM-78 model for simplicity. The NGM-78 model describes the following discrete-time Markov chain in ℝ^d:


    X_n = a_0 X_n-1 + a_1 X_n-1/(1 + |X_n-1|^2) + f(n) + v_n,

and the observations

    Y_n = a_2  |Y_n|^2 + u_n,

where v_n ∼ N(0,σ^2_v I_d) and u_n ∼ N(0, σ^2_u I_d). In our experiments, we choose a_0 = 1/2, a_1 = 25, a_2 = 1/20, f ≡ 0, σ_u^2 = 1, σ_v^2 = 0.01, and n = 0,1,…, 50. We test our TPPF algorithm on the NGM-78 model with the configuration d ∈{1,2,5,10,15,20 } and compare TPPF with its competitors. The empirical standard deviation of log Z is reported in Table <ref>. As we can see, in constrast to the linear Gaussian model, in the current nonlinear settings, the linear structure (twisted functions learned in some Gaussian function family) of iAPF leads to its relatively worse behavior. On the other hand, our TPPF algorithms behave better than iAPF partly due to the stronger expressive ability of neural networks. Moreover, in this model the term P̂[g](k,x) (recall the definition in (<ref>)) used in FA-APF is calculated via Monte-Carlo approximation since there is not analytical solution for it. Consequently, we observe in Table <ref> that although the FA-APF outperforms other algorithms when d=1, it suffers from the curse of dimensionality and we cannot obtain reasonable estimates with the FA-APF in a feasible computational time.












 §.§ Lorenz-96 model


The Lorenz-96 is a nonlinear model with partial observation <cit.>. The discrete model can be viewed as the Euler-Maruyama scheme of the following interacting particle system consisting of d particles in ℝ^1:

    dX^i = (-X^i-1X^i-2 + X^i-1X^i+1 - X^i + α)dt + σ^2 dW^i,   1 ≤ i ≤ d,

where α∈ℝ, σ^2 ∈ℝ^+, W^i (1≤ i≤ d) are independent Brownian motions in ℝ^d, and the indices should be understood modulo d. The observation is through Y_t ∼𝒩(·; HX_t,Σ_OB), where H is a diagonal matrix with H_ii = 1 (1 ≤ i ≤ d-2), H_ii = 0 (d-1 ≤ i ≤ d). Note that this model does not have an analytical solution. We use the non-parametric implementation in this example.

Choose α = 3.0, Σ_OB = I_d, σ = 1. Consider d ∈{3,5,10 }. We report here the empirical standard deviations in Table <ref> below:



To test the sensitivity of the proposed method, we fix d = 3, Σ_OB=I_d, σ=1 and run the experiment for different α. See the results in Figure <ref>.



































As we can see from Table <ref> and Figure <ref>, the proposed TPPF algorithm (especially the ones using loss functions L_RE and L_RECE) behaves much better than BPF and FA-APF, and slightly better than iAPF. Moreover, we observe that for the current partial-observed nonlinear model, sometimes the result of iAPF is more biased (see (b) in Figure <ref>), despite its relatively low variance. 





























§ CONCLUSION

In this paper, we study the discrete-time twisted particle filters (TPF) from a continuous-time perspective. In detail, under suitable settings, we prove the convergence from discrete-time to continuous-time models. This then enables us to view existing control-based importance sampling algorithms as good guidance for constructing novel TPF algorithms in discrete-time settings. As a concrete example, we propose a novel TPF algorithm, Twisted-Path Particle Algorithm (TPPF), inspired by algorithms in continuous-time settings. In TPPF, we choose some specific KL-divergence between path measures as the loss function and learn the twisting function parameterized by a neural network. We also give some numerical examples to illustrate the capability of the proposed algorithm. Some possible future work may include: seeking other practical TPF algorithms guided by the continuous-time importance sampling algorithms, and understanding existing TPF algorithms rigorously by finding its continuous-time limit.





§ ACKNOWLEDGEMENTS

This work is supported in part by National Science Foundation via grant DMS-2309378. 
This work is done during Yuliang Wang's visit to Duke University. He thanks the Mathematics Department for their hospitality. 
The authors would thank Pieralberto Guarniero, Adam M. Johansen and Anthony Lee for help discussions on implementation details of iAPF method.



§ PROOFS OF SECTION <REF>



 §.§ Proof of Proposition <ref>


Recall the notation P̂[φ](k,x) = ∫φ(k,y') P̂(x,dy'). Then by definition,

    𝔼_x[∏_k=0^nĝ_k^φ(X̂^φ_k)]
    = ∫ĝ_0^φ(x) ∏_k=1^n  ĝ_k^φ(x_k) P̂_k^φ(x_k-1,x_k)dx_1:n
       =∫ĝ_0(x) (∏_k=1^n  ĝ_k(x_k))(P̂[φ](1,x)/φ(n,x_n)∏_k=1^n-1P̂[φ](k+1,x_k)/φ(k,x_k)∏_k=1^n φ(k,x_k)/P̂[φ](k,x_k-1)P̂_k(x_k-1,x_k))dx_1:n
       =∫ĝ_0(x) ∏_k=1^n  ĝ_k(x_k) P̂_k(x_k-1,x_k)dx_1:n = Z_dis(x).







 §.§ Proof Proposition <ref>


The proof is straightforward.

    
  * We first prove the first claim. Since φ^*(n,x) = ĝ_n(x), direct calculation yields

    φ^*(n-1,x) = ĝ_n-1(x) ∫ĝ_n(y) P̂(x,dy) = 𝔼[ĝ_n(X̂_n)ĝ_n-1(X̂_n-1)|X̂_(n-1) = x],
       φ^*(n-2,x) = ĝ_n-2(x) ∫𝔼[∏_i=n-1^n ĝ_i(X̂_i)|X̂_n-1 = y] P̂(x,dy)
       = 𝔼[∏_i=n-2^n ĝ_i(X̂_i)|X̂_n-2 = x],
       ⋯
       φ^*(k,x) = ĝ_k(x) ∫𝔼[∏_i=k+1^n ĝ_i(X̂_i)|X̂_k+1 = y] P̂(x,dy)
       = 𝔼[∏_i=k^n ĝ_i(X̂_i)|X̂_k = x],
       ⋯
       φ^*(0,x) = ĝ_0(x) ∫𝔼[∏_i=1^n ĝ_i(X̂_i)|X̂_1 = y] P̂(x,dy)
       = 𝔼[∏_i=0^n ĝ_i(X̂_i)|X̂_0 = x] = Z_dis(x).


  * Next, we prove the unbiased and zero-variance property of Ŵ. The unbiased property is a direct result of Proposition <ref> and Remark <ref>, and is valid for any twisting function sequence φ(k,·) (1≤ k ≤ n). To prove the zero variance property, we observe that under the probability measure P^*,

    Ŵ = ∏_k=0^nℓ_k^*(X̂_k) ∏_k=0^n ĝ_k (X̂_k) = ĝ_n(X̂_n)φ^*(0,X̂_0)/φ^*(n,X̂_n) = φ^*_η(0,X̂_0) = Z_dis(X̂_0).











 §.§ Proof of Corollary <ref>


By definition,

    Z^N.φ_dis(x) = ∏_k=0^n 1/N∑_i=1^N ĝ^φ^*_k(ζ_k^i),

and

    ĝ^φ^*_n(x)=ĝ_n(x) ℓ_n^φ^*(x) = ĝ_n(x)/φ^*(n,x)≡ 1,


    ĝ^φ^*_k(x) = ĝ_k(x) ℓ_k^φ^*(x) = ĝ_k(x) P̂[φ^*](k+1,x)/φ^*(k,x)≡ 1,   1≤ k ≤ n-1,


    ĝ^φ^*_0(x) = ĝ_0(x) ℓ_0^φ^*(x) = ĝ_0(x) P̂[φ^*](1,x) = φ^*(0,x).

Therefore, 

    Z^N,φ^*_dis(x) =φ^*(0,x)= Z_dis(x),  ∀ x ∈ℝ^d.









 §.§ Proof of Proposition <ref>




  * We first prove the Feynman-Kac representation formula. Fix 0≤ t≤ T, define the following process for s ∈ [T-t,T]

    Y(s) := e^∫_t^s h(τ,X_τ) dτ v^*(s,X_s).

Clearly,

    Y(t) = v^*(t,X_t),   Y(T) = e^∫_t^T h(τ ,X_τ) dτ g(X_T).

By Itô's formula,

    dY(s)    = e^∫_t^s h(τ,X_τ) dτ h(s,X_s) v^*(s,X_s) ds
         + e^∫_t^s h(τ,(X_τ) dτ∇_x v^*(s, X_s) ·(b(X_s) ds + √(2) dB_s )
         + e^∫_t^s h(τ,X_τ dτΔ_x v^*(s,X_s)ds
          + e^∫_t^s h(τ,X_τ) dτ∂_t v^*(s,X_s) ds

By (<ref>), we have

    dY(s) = √(2) e^∫_t^s h(τ,X_τ) dτ∇_x v^*(s, X_s) · dB_s.

Hence, Y(s) is a martingale, and consequently,

    𝔼[Y(T)| X_t = x] = 𝔼[Y(t)| X_t = x],

namely,

    v^*(t,x) = 𝔼[ e^∫_t^T h(τ, X_τ) dτ g(X_T)| X_t = x],  ∀ x ∈ℝ^d,  ∀ t∈[0,T].



  * The existence and expression of the Radon-Nikodym derivative L(T) is guaranteed by the classical Girsanov's theorem. We focus on the derivation of the zero-variance property (<ref>) here. 

For s ∈ [0,T], define the process

    ω_s := e^∫_0^s h(τ, X^u_τ)dτv^*(s,X^u_s)L(s).

Clearly, 

    ω_0 = v^*(0,X_0),  ω_T = W.

By Itô's formula,

    dω_s    = h(s, X_s^u^*)e^∫_0^s h(τ, X^u_τ)dτv^*(s,X^u_s)L(s) ds
         + e^∫_0^s h(τ, X^u^*_τ)dτ∂_s v^*(s,X^u^*_s)L(s)ds
         +e^∫_0^s h(τ, X^u^*_τ)dτ∂_x v^*(s,X^u^*_s) ·(b(X_s^u^*) ds + 2 ∂_x log v^*(s,X^u^*_t) ds + √(2)dB_s )L(s)
         +e^∫_0^s h(τ,X^u^*_τ)dτv^*(s,X^u^*_s)L(s)(-√(2)∂_x log v^*(s,X_s^u^*) dB_s )
         -e^∫_0^s h(τ, X^u_τ)dτL(s)(√(2)∂_x v^*(s,X_s^u^*) ·√(2)∂_x log v^*(s,X_s^u^*))ds.

Using the time evolution (<ref>) for u^*, we have

    dω_s ≡ 0.

Therefore

    v^*(0,X_0) = ω_0 = ω_T = W    P^*-a.s.




    





§ PROOFS OF SECTION <REF>



 §.§ Proof of Proposition <ref>


Consider the following SDEs with the same initial state:

    dX_t = b(X_t) dt + √(2)dB,   X_0 = x,   0≤ t ≤ T=nη,
       dX̅_t = b(X̅_kη)dt + √(2)dB,  X̅_0 = x,    t∈[kη,(k+1)η),   0≤ k≤ n-1.

Clearly, P_η(x,·) = Law(X_η) and P̂^η(x,·) = Law(X̅_η).
Denote by P^x_[0,t], P̂_[0,t]^x the path measures with the same initial x associated with the time interval [0,t] for any fixed t ≤ T. Then using the data processing inequality <cit.> and Girsanov's theorem <cit.>, it holds that

    (Law(X̅_t) Law(X_t) ) ≤(P̂_[0,t]^x  P_[0,t]^x ) ≤𝔼[∑_k=0^⌈ t/η⌉∫_kη^(k+1)η |b(X̅_t) - b(X̅_kη)|^2 dy].

By condition (a) in Assumption <ref>, b(·) is L_b-Lipschitz, and the second  moment for X̅_kη has uniform bound, then for any 0≤ k ≤⌈ t/η⌉, one has 

    𝔼|b(X̅_t) - b(X̅_kη)|^2   ≤ L_b^2𝔼|(t-kη)b(X̅_kη) + ∫_kη^t dB_s |^2
       ≤ 2L_b^2(2η^2(|b(0)|^2 + L_b^2 sup_0≤ i ≤ n𝔼|X̅_iη|^2) + η d ) ≤ Cdη,

where we need η < 1 and C = C(L_b, b(0), sup_0≤ i ≤ n𝔼|X̅_iη|^2 ) is a positive constant.
Combining (<ref>) and (<ref>), we know that

    (Law(X̅_t) Law(X_t) ) ≤ Cdtη.

And by Pinsker's inequality, we have

    TV(Law(X̅_t) , Law(X_t) ) ≤ C√(dtη).

Finally, taking t = η, we obtain the desired result.




 §.§ Proof of Proposition <ref>


Fix T>0, η>0, x∈ℝ^d. Recall that T = nη, and

    Z_dis(x) := 𝔼_x[∏_k=0^n-1ĝ_k(X̂_k)ĝ_n(X̂_n)] = 𝔼_x[exp(∑_k=0^n-1∫_kη^(k+1)ηlog(ĝ_k(X̂_k))^η^-1ds)ĝ_n(X̂_n)],


    Z_con(x) = 𝔼_x[e^∫^T_0log g(x,X_s) ds g_T(X_T)  ].

Then, for

    dX_t = b(X_t) dt + √(2)dB,   X_0 = x,   0≤ t ≤ T=nη,
       dX̅_t = b(X̅_kη)dt + √(2)dB,  X̅_0 = x,    t∈[kη,(k+1)η),   0≤ k≤ n-1,

we have

    Z_con(x) - Z_dis(x)
       = 𝔼_x[exp(∫^T_0log g(s,X_s) ds) g_T(X_T) ] - 𝔼_x[exp(∫^T_0log g(s, X̅_s) ds) g_T(X̅_T) ]
       +𝔼_x[exp(∫^T_0log g(s, X̂_s) ds) ( g_T(X̅_T) -ĝ_n(X̅_T) ) ]
       +𝔼_x[ĝ_n(X̅_T) (∏_i=0^n-1exp(∫_iη^(i+1)ηlog g(s,X̅_s)ds) - ∏_i=0^n-1exp(∫_iη^(i+1)ηη^-1logĝ_iη(X̅_iη)ds) ) ],


For the first term above, by condition (d) in Assumption <ref>, exp(∫^T_kηlog g(s, X_s) ds) g_T(X_T) is a continuous, bounded functional of X_[kη,T] with X_kη = x, denote by F_k(X_[kη,T]). Then, using the KL upper bound for path measures obtained in (<ref>), we have

    𝔼_x[exp(∫^T_0log g(s,X_s) ds) g_T(X_T) ] - 𝔼_x[exp(∫^T_0log g(s, X̅_s) ds) g_T(X̅_T) ]
       =𝔼_x[F_0(X_[0,T]) ] - 𝔼_x[F_0(X̅_[0,T]) ]
       =∫ F_0(y) (P_[0,T]^x(y) - P̂_[0,T]^x(y) )dy
       ≤ C  TV(P_[0,T]^x , P̂_[0,T]^x)
        ≤ C (P̂_[0,T]^x  P_[0,T]^x)^1/2≤ C√(Tdη)→ 0,

where we have used the Pinsker's inequality in the last line above.

For the second term, by conditions (b), (d) in Assumption <ref>, exp(∫^T_0log g(s, X̂_s) ds) is bounded, and ĝ_n(x) → g_T(x) uniformly in x. Therefore, as η→ 0,

    𝔼_x[exp(∫^T_0log g(s, X̂_s) ds) ( g_T(X̅_T) -ĝ_n(X̅_T) ) ] → 0.


For the third term, for s ∈ [iη,(i+1)η),

    |log g(s,X̅_s) - η^-1logĝ_iη(X̅_iη)| ≤ |log g(s,X̅_s) - log g(iη,X̅_s)|
     + |log g(iη, X̅_s) - η^-1logĝ_iη(X̅_s)| + |η^-1logĝ_iη(X̅_s) - η^-1logĝ_iη(X̅_iη)|.

By conditions(b), (c) in Assumption <ref>, as η→ 0, |log g(s,x) - log g(iη,x)|, |log g(iη,x) - η^-1logĝ_iη(x)| tend to 0 uniformly in s, i, x. Moreover, by condition (c) in Assumption <ref>, |logĝ_iη(X̅_s) - logĝ_iη(X̅_iη)| ≲ |X̅_s - X̅_iη | ≲η |b(X̂_iη)| + |W_s - W_iη|.
Hence, as η→ 0,

    𝔼_x[ĝ_n(X̅_T) (∏_i=0^n-1exp(∫_iη^(i+1)ηlog g(s,X̅_s)ds) - ∏_i=0^n-1exp(∫_iη^(i+1)ηη^-1logĝ_iη(X̅_iη)ds) )) ] → 0.

Combining all the above, we conclude that for any fixed x ∈ℝ^d,

    | Z_dis(x) - Z_con(x) |  → 0   as  η→ 0.












§ PROOF OF SECTION <REF>



 §.§ Proof of Lemma <ref>


We refer to <cit.> for a similar proof. The proof is based on the Jensen's inequality. In fact, by direct calculation and convexity of -log(·), we have

    -log𝔼^X ∼ P[exp(-W(X) )] = -log∫ e^-W(x) P(dx) = -log∫ e^-W(x)dP/dQ(x) Q(dx)
    ≤∫(W(x) - logdP/dQ(x)) Q(dx) =  𝔼^X∼ Q[W(X)] + (Q  P),

and the equality holds if and only if

    dQ/dP(X) = exp(-log𝔼^X ∼ P[exp(-W(X))] -W(X)) .





 §.§ Proof of Proposition <ref>


By definition (<ref>), we have

    J(φ)    = 𝔼_x[-∑_k=0^nlogĝ_k(X̂_k)]
         + ∫log(∏_k=1^nφ(k,x_k)P̂(x_k-1,x_k)/P̂[φ](k,x_k-1) / ∏_k=1^n P̂(x_k-1,x_k)) ∏_k=1^nP̂^φ_k(x_k-1,x_k) dx_1:n
       = 𝔼_x[-∑_k=0^n logĝ_k(X̂_k^φ) + ∑_k=1^n logφ(k,X̂_k^φ)/P̂[φ](k,X^φ_k-1)],

and

    (P^φP^φ^*) = ∫log(∏_k=1^nφ(k,x_k)/P̂[φ](k,x_k-1) / ∏_k=1^nφ^*(k,x_k)/P̂[φ^*](k,x_k-1))∏_k=1^nP̂_k^φ(x_k-1,x_k) dx_1:n
        =𝔼_x[∑_k=1^n logφ(k,x_k)/P̂[φ](k,X^φ_k-1)] - ∫log(φ^*(n,x_n)/φ^*(0,x)∏_k=0^n-1φ^*(k,x_k)/P̂[φ^*](k+1,x_k))∏_k=1^nP̂^φ_k(x_k-1,x_k) dx_1:n
       = 𝔼_x[-∑_k=0^n logĝ_k(X̂_k^φ) + ∑_k=1^n logφ(k,X̂_k^φ)/P̂[φ](k,X^φ_k-1)] + logφ^*_0(x),

where we have used the recursive relation (<ref>) in the last equality. Moreover, using the derived expression for J(φ) and (<ref>), we know that

    J(φ^*) = logφ_0^*(x).

Consequently,

    J(φ) = J(φ^*) + (P^φP^φ^*).










 §.§ Proof of Proposition <ref>


The proof relies on the following auxiliary results:

    
  * (Lemma <ref>, generalized Jensen's inequality) Given deterministic functions ϕ: Ω→ℝ and f: ℝ→ℝ. Assume that f is convex. Let λ, λ' be two probability measures on Ω that are absolutely continuous with each other. Define the functional
    
    𝒥(f,λ,ϕ) := 𝔼_λ f(ϕ) - f𝔼_λϕ.

    Then, 
    
    m𝒥(f,λ, ϕ) ≤𝒥(f,λ', ϕ) ≤ M𝒥(f,λ, ϕ),

    where m := inf_E∈ℬ(Ω)λ'(E)/λ(E), M := inf_E∈ℬ(Ω)λ(E)/λ'(E).
    
  * (Lemma <ref>, equivalence with the 𝒳^2-divergence)
    
    r^2(ν̃) = 𝒳^2(ν^* | ν̃),

    where the 𝒳^2-divergence is defined by
    
    𝒳^2(ν^* | ν̃):=𝔼_μ[|dν^*/dν̃|^2 - 1].


Now, by Jensen's inequality, we have

    (ν^*ν̃) = 𝔼_ν^*[logdν^*/dν̃] ≤log𝔼_ν^*[ dν^*/dν̃]

Combining (<ref>) and (<ref>), we have

    r^2(ν̃) = 𝒳^2(ν^* | ν̃) = 𝔼_μ[|dν^*/dν̃|^2 - 1] = 𝔼_ν^*[dν^*/dν̃] - 1 ≥ e^(ν^*ν̃) - 1.

For the second claim in Proposition <ref>, we choose λ' = ν^*, λ = ν̃, ϕ = dν^*/dν̃, and f = -log in (<ref>). Then,

    𝒥(f,λ,ϕ) = -𝔼_ν̃logdν^*/dν̃ + log𝔼_ν̃dν^*/dν̃ = (ν̃ν^*),

and

    𝒥(f,λ',ϕ)    = -𝔼_ν^*logdν^*/dν̃ + log𝔼_ν^*dν^*/dν̃
       = -(ν^*ν̃) + log (𝒳^2(ν^* | ν̃) + 1) = -(ν^*ν̃) + log (r^2(ν̃) + 1),

where we have used (<ref>) in the last equality. Combining (<ref>), (<ref>) and (<ref>), we obtain the second claim 

    e^m(ν̃ν^*) + (ν^*ν̃)-1 ≤ r^2(ν̃) ≤ e^M(ν̃ν^*) + (ν^*ν̃)-1,

where m:=inf_E ν^*(E)/ν̃(E) and M := sup_E ν^*(E)/ν̃(E).



The two auxiliary lemmas used in the proof of Proposition <ref> are given below:


Given deterministic functions ϕ: Ω→ℝ and f: ℝ→ℝ. Assume that f is convex. Let λ, λ' be two probability measures on Ω that are absolutely continuous with each other. Define the functional
    
    𝒥(f,λ,ϕ) := 𝔼_λ f(ϕ) - f(𝔼_λϕ).

    Then, 
    
    m𝒥(f,λ, ϕ) ≤𝒥(f,λ', ϕ) ≤ M𝒥(f,λ, ϕ),

    where m := inf_E∈ℬ(Ω)λ'(E)/λ(E), M := inf_E∈ℬ(Ω)λ(E)/λ'(E).    



We first show that

    m𝒥(f,λ, ϕ) ≤𝒥(f,λ',ϕ).

Taking E = Ω, we have 

    m := inf_E λ'(E)/λ(E)≤λ'(E)/λ(E) = 1.

Without loss of generosity, assume m<1. (If m=1, then λ≡λ', and the argument is trivial). (<ref>) is then equivalent to

    m𝔼_λ f(ϕ) - mf(𝔼_λϕ) ≤𝔼_λ'f(ϕ) - f(𝔼_λ'ϕ).

Since m ∈ (0,1), and m = inf_E λ'(E)/λ(E), the probability λ̅ := λ' - mλ/1-m is well-defined. Then, by Jensen's inequality, since f is convex, we have

    𝔼_λ'f(ϕ) - m𝔼_λ f(ϕ) = (1-m)𝔼_λ̅ f(ϕ) ≥ (1-m)f(𝔼_λ̅ϕ) = (1-m)f(𝔼_λ'ϕ - m𝔼_λϕ/1-m).

Using convexity of f again, we have

    (1-m)f(𝔼_λ'ϕ - m𝔼_λϕ/1-m) + mf(𝔼_λϕ) ≥ f(𝔼_λ'ϕ).

Therefore, (<ref>) holds, and thus

    m𝒥(f,λ, ϕ) ≤𝒥(f,λ',ϕ).

Assuming M>1 and using exactly the same arguments, we have

    M𝒥(f,λ, ϕ) ≥𝒥(f,λ', ϕ).





Recall the definitions of ν̃, ν^* and the relative variance r(ν̃) in Proposition <ref>. Then
    
    r^2(ν̃) = 𝒳^2(ν^* | ν̃),

    where the 𝒳^2-divergence is defined by
    
    𝒳^2(ν^* | ν̃):=𝔼_μ[|dν^*/dν̃|^2 - 1].




By definition, we have

    𝒳^2(ν^* | ν̃):=𝔼_μ[|dν^*/dν̃|^2 - 1] = 𝔼_ν̃| dν^*/dν̃|^2 - |𝔼_ν̃dν^*/dν̃|^2 = _ν̃(dν^*/dν̃) = _ν̃(dν^*/dνdν/dν̃).

Recall the zero-variance property of ν^*:

    dν^*/dν = e^-W/Z  ν-a.s. .

Consequently,

    _ν̃(dν^*/dνdν/dν̃) = _ν̃(e^-W/Zdν/dν̃) = _ν̃(e^-Wdν/dν̃)/Z^2 = r^2(ν̃).

Hence, 

    r^2(ν̃) = 𝒳^2(ν^* | ν̃).






plain


