Differentiable Systematic Resampling for Variational Sequential Monte Carlo

Fredrik Cumlin; Saikat Chatterjee

stat.MLOct 8, 2026

Differentiable Systematic Resampling for Variational Sequential Monte Carlo

Authors: Fredrik Cumlin, Saikat Chatterjee

Organizations: Department of Information Science KTH Royal Institute of Technology Stockholm, Sweden

Abstract

Particle filters are a standard tool for nonlinear state estimation, but their resampling step is discrete, preventing gradient-based learning in variational sequential Monte Carlo. We introduce Differentiable Systematic Resampling (DSR), a temperature-controlled relaxation of systematic resampling, that preserves the CDF-ordered, banded structure of systematic resampling while enabling full gradient flow. DSR converges to exact systematic resampling as the temperature vanishes, and we prove a pointwise exponential convergence rate for the induced bias. Compared to optimal-transport-based differentiable resampling, DSR avoids iterative solvers and has substantially lower computational overhead. Experiments on stochastic dynamical systems and real-world handwriting data show that DSR achieves comparable or superior filtering and dynamics learning performance.

Figures & tables

Figure 1: Transportation maps for systematic resampling (left), DSR (middle), and EROT (right), for a sequential important resampling PF of the Lorenz-63 system. The transportation maps are calculated at step $150$ using $N=25$ particles, and DSR and EROT uses training relaxation parameters.

Method	KLD KF	NMSE (dB)	$\\|A-\hat{A}\\|_{F}$
Hard	$\mathbf{17.49}{\scriptstyle\pm 0.07}$	$\underline{-11.33}{\scriptstyle\pm 0.06}$	$0.035{\scriptstyle\pm 0.001}$
Soft	$\underline{17.52}{\scriptstyle\pm 0.09}$	$-10.33{\scriptstyle\pm 0.07}$	$0.025{\scriptstyle\pm 0.002}$
EROT	$22.09{\scriptstyle\pm 0.06}$	$-11.31{\scriptstyle\pm 0.02}$	$\underline{0.022}{\scriptstyle\pm 0.004}$
DSR	$21.77{\scriptstyle\pm 0.06}$	$\mathbf{-11.34}{\scriptstyle\pm 0.01}$	$\mathbf{0.019}{\scriptstyle\pm 0.001}$

Table 1: Performance on the linear Gaussian SSM (

d_{x}=4

, SMNR

=10

dB). KLD KF is the training-time divergence to the Kalman posterior; NMSE is evaluated with

N=1000

particles and hard systematic resampling;

\|A-\hat{A}\|_{F}

measures recovery of the true transition matrix. Results show mean

\pm

std over

3

runs.

Figure 2: Lorenz-63 system with nonlinear camera measurement model. (a) Sample state trajectory of Lorenz-63. Colored points indicate three observation times that are shown in the right-most figure. (b) Camera projection geometry: the 3D state $\mathbf{x}_{t}=(x_{1},x_{2},x_{3})$ projects onto an 8×8 pixel grid via a depth-dependent Gaussian point spread function (Eq. 10 ). The first two coordinates determine the 2D position, while the third coordinate (depth) controls the spread of the Gaussian. (c) Three sample camera observations at SMNR = 0 dB.

SMNR = $10$ dB
	State Estimation	Transition Model Quality
Method	NMSE (dB)	MSE	MSE-5-step	LL	KLD
Hard	$\underline{-28.86}{\scriptstyle\pm 0.41}$	$\underline{0.55}{\scriptstyle\pm 0.19}$	$5.50{\scriptstyle\pm 3.68}$	$\underline{-2.26}{\scriptstyle\pm 0.25}$	$\underline{5.41}{\scriptstyle\pm 1.47}$
Soft	$-27.82{\scriptstyle\pm 0.45}$	$0.86{\scriptstyle\pm 0.49}$	$6.45{\scriptstyle\pm 3.92}$	$-2.55{\scriptstyle\pm 0.48}$	$6.90{\scriptstyle\pm 2.43}$
EROT	$-28.76{\scriptstyle\pm 0.31}$	$0.56{\scriptstyle\pm 0.24}$	$\underline{3.39}{\scriptstyle\pm 2.59}$	$-2.32{\scriptstyle\pm 0.20}$	$5.74{\scriptstyle\pm 1.03}$
DSR	$\mathbf{-29.94}{\scriptstyle\pm 0.28}$	$\mathbf{0.47}{\scriptstyle\pm 0.24}$	$\mathbf{3.38}{\scriptstyle\pm 2.11}$	$\mathbf{-1.98}{\scriptstyle\pm 0.27}$	$\mathbf{3.60}{\scriptstyle\pm 1.04}$

Table 2: Performance comparison on the Lorenz-63 system at SMNR = 10 and 0 dB. Results show mean

\pm

std over 10 runs. Best in bold , second-best underlined .

Figure 3: Gradient norm during training on Lorenz-63 (SMNR

=10

dB) for a representative run. The solid line shows the mean gradient norm per epoch across batches; the shaded region indicates the within-epoch standard deviation.

Method	Mean	Within-epoch std
Hard	$17399$	$17253$
Soft	$17013$	$8314$
EROT	$15771$	$6918$
DSR	$\mathbf{4504}$	$\mathbf{2720}$

Table 3: Gradient norm statistics on Lorenz-63 (SMNR

=10

dB), averaged over runs.

Figure 4: Handwriting dynamics from the CharacterTrajectories dataset. Sample trajectory of the letter ‘a’ captured at 200 Hz on a WACOM tablet. (a) state trajectory, with line thickness indicating pen pressure and colour gradient showing temporal progression from start (black dot) to finish. (b) noisy observations at SMNR = 10 dB.

	NMSE (dB)
Method	SMNR = 10 dB	SMNR = 0 dB
Hard	$-14.65{\scriptstyle\pm 0.20}$	$-7.66{\scriptstyle\pm 0.35}$
Soft	$-14.69{\scriptstyle\pm 0.29}$	$-7.08{\scriptstyle\pm 1.12}$
EROT	$\underline{-14.88}{\scriptstyle\pm 0.17}$	$\underline{-7.77}{\scriptstyle\pm 1.26}$
DSR	$\mathbf{-15.21}{\scriptstyle\pm 0.15}$	$\mathbf{-8.59}{\scriptstyle\pm 0.22}$

Table 4: Performance comparison of the estimated transition model on the CharacterTrajectories dataset. Results show mean

\pm

std over 10 runs with different seeds.

Method	Total (h)	Time/epoch (s)
Hard	$\mathbf{2.49}{\scriptstyle\pm 0.15}$	$\mathbf{17.94}{\scriptstyle\pm 1.11}$
Soft	$4.33{\scriptstyle\pm 0.05}$	$31.16{\scriptstyle\pm 0.39}$
EROT	$5.96{\scriptstyle\pm 0.11}$	$42.94{\scriptstyle\pm 0.77}$
DSR	$\underline{2.55}{\scriptstyle\pm 0.02}$	$\underline{18.39}{\scriptstyle\pm 0.16}$

Table 5: Training time on Lorenz-63.

Appendix figures & tables12 assets

Supplementary material from the paper’s appendix.

Appendix

Figure 5: Transportation maps for hard systematic resampling and DSR at varying temperatures $\tau\in\{0.01,0.05,0.1,0.5,1.0\}$ , alongside EROT ( $\varepsilon=0.5$ ), for a SIR particle filter of the Lorenz-63 system with $N=25$ particles.

Panel of Figure 6 — Figure 6: Empirical analysis of the DSR relaxation on the Lorenz-63 system with oracle dynamics at SMNR $=10$ dB. Error bars show $\pm 1$ standard deviation over 10 sequences.

	$f(\mathbf{x})=\mathbf{x}$			$f(\mathbf{x})=\mathbf{x}^{2}$
Method	bias	var	MSE	bias	var	MSE
Hard SR	$2.3\times 10^{-4}$	$1.2\times 10^{-3}$	$1.2\times 10^{-3}$	$1.2\times 10^{-4}$	$7.0\times 10^{-4}$	$7.0\times 10^{-4}$
Multinomial	$4.9\times 10^{-4}$	$5.5\times 10^{-3}$	$5.5\times 10^{-3}$	$2.8\times 10^{-4}$	$2.6\times 10^{-3}$	$2.6\times 10^{-3}$
DSR ( $\tau=0.01$ )	$2.4\times 10^{-3}$	$2.9\times 10^{-5}$	$3.4\times 10^{-5}$	$1.2\times 10^{-1}$	$3.8\times 10^{-7}$	$1.4\times 10^{-2}$
DSR ( $\tau=0.05$ )	$1.2\times 10^{-2}$	$8.6\times 10^{-6}$	$1.6\times 10^{-4}$	$2.3\times 10^{-1}$	$3.1\times 10^{-7}$	$5.3\times 10^{-2}$
DSR ( $\tau=0.1$ )	$2.0\times 10^{-2}$	$4.7\times 10^{-6}$	$3.9\times 10^{-4}$	$2.7\times 10^{-1}$	$2.0\times 10^{-7}$	$7.3\times 10^{-2}$

Appendix

Table 6: Bias, variance, and MSE of resampling estimates against the weighted estimate

I

, over

20000

seeds with

N=100

. Hard SR and multinomial are unbiased (their bias is Monte-Carlo error).

Temperature $\tau$	NMSE (dB)
0.01	$-28.31\pm 0.15$
0.05	$-29.56\pm 0.34$
0.10	$\underline{-29.94}\pm 0.28$
0.50	$\mathbf{-30.21\pm 0.29}$
1.00	$-29.91\pm 0.28$

Appendix

Table 7: Ablation study of the temperature parameter

\tau

on Lorenz-63 (SMNR = 10 dB).

$\tau$	NMSE (dB)	Median	Mean	Median	Mean	$L_{2}$ to hard SR
		Gradient norm		Within-epoch std
$0.01$	$-28.31{\scriptstyle\pm 0.15}$	$14842$	$2.39\times 10^{5}$	$8163$	$4.35\times 10^{5}$	$0.102{\scriptstyle\pm 0.009}$
$0.05$	$-29.56{\scriptstyle\pm 0.34}$	$6457$	$7228$	$3009$	$3676$	$0.277{\scriptstyle\pm 0.012}$
$0.10$	$-29.94{\scriptstyle\pm 0.28}$	$4174$	$4504$	$1885$	$2634$	$0.356{\scriptstyle\pm 0.017}$
$0.50$	$\mathbf{-30.21}{\scriptstyle\pm 0.29}$	$2546$	$2951$	$1155$	$1985$	$0.505{\scriptstyle\pm 0.026}$
$1.00$	$-29.91{\scriptstyle\pm 0.28}$	$2306$	$3041$	$1027$	$3453$	$0.479{\scriptstyle\pm 0.046}$

Appendix

Table 8: Gradient behaviour and approximation error as a function of

\tau

on Lorenz-63 (SMNR

=10

dB). The

L_{2}

column is the per-particle distance of the DSR transport to hard systematic resampling from the same particle set and offset (Appendix E ).

	State Est.	Transition Model Quality
Setting	NMSE (dB)	MSE	MSE-5-step	LL	KLD
Fixed ( $\tau=0.1$ )	$-29.94{\scriptstyle\pm 0.28}$	$0.47{\scriptstyle\pm 0.24}$	$3.38{\scriptstyle\pm 2.11}$	$-1.98{\scriptstyle\pm 0.27}$	$3.60{\scriptstyle\pm 1.04}$
Annealed	$-30.33{\scriptstyle\pm 0.10}$	$0.40{\scriptstyle\pm 0.11}$	$3.01{\scriptstyle\pm 0.63}$	$-1.86{\scriptstyle\pm 0.13}$	$3.15{\scriptstyle\pm 0.32}$

Appendix

Table 9: Fixed versus annealed

\tau

on Lorenz-63 (SMNR

=10

dB). Results show mean

\pm

std over 10 runs.

		State Est.	Transition Model Quality
$N$	Method	NMSE (dB)	MSE	MSE-5-step	LL	KLD
25	Hard	$\underline{-28.86}{\scriptstyle\pm 0.41}$	$\underline{0.55}{\scriptstyle\pm 0.19}$	$5.50{\scriptstyle\pm 3.68}$	$\underline{-2.26}{\scriptstyle\pm 0.25}$	$\underline{5.41}{\scriptstyle\pm 1.47}$
	Soft	$-27.82{\scriptstyle\pm 0.45}$	$0.86{\scriptstyle\pm 0.49}$	$6.45{\scriptstyle\pm 3.92}$	$-2.55{\scriptstyle\pm 0.48}$	$6.90{\scriptstyle\pm 2.43}$
	EROT	$-28.76{\scriptstyle\pm 0.31}$	$0.56{\scriptstyle\pm 0.24}$	$\underline{3.39}{\scriptstyle\pm 2.59}$	$-2.32{\scriptstyle\pm 0.20}$	$5.74{\scriptstyle\pm 1.03}$
	DSR	$\mathbf{-29.94}{\scriptstyle\pm 0.28}$	$\mathbf{0.47}{\scriptstyle\pm 0.24}$	$\mathbf{3.38}{\scriptstyle\pm 2.11}$	$\mathbf{-1.98}{\scriptstyle\pm 0.27}$	$\mathbf{3.60}{\scriptstyle\pm 1.04}$
50	Hard	$-28.30{\scriptstyle\pm 0.88}$	$\mathbf{0.39}{\scriptstyle\pm 0.08}$	$\mathbf{1.73}{\scriptstyle\pm 0.61}$	$\underline{-1.98}{\scriptstyle\pm 0.09}$	$\underline{4.30}{\scriptstyle\pm 0.66}$

Appendix

Table 10: Effect of particle count

N

on resampling methods for Lorenz-63 at SMNR

=10

dB. Results show mean

\pm

std over 10 runs. OOM indicates that training exceeded available GPU memory.

SMNR $=10$ dB
	State Est.	Transition Model Quality
Ordering	NMSE (dB)	MSE	MSE-5-step	LL	KLD
Morton	$-29.94{\scriptstyle\pm 0.28}$	$0.47{\scriptstyle\pm 0.24}$	$3.38{\scriptstyle\pm 2.11}$	$-1.98{\scriptstyle\pm 0.27}$	$3.60{\scriptstyle\pm 1.04}$
Random	$-29.96{\scriptstyle\pm 0.34}$	$0.44{\scriptstyle\pm 0.20}$	$2.31{\scriptstyle\pm 0.80}$	$-1.92{\scriptstyle\pm 0.21}$	$3.59{\scriptstyle\pm 0.80}$
SMNR $=0$ dB
Morton	$-18.27{\scriptstyle\pm 0.46}$	$0.54{\scriptstyle\pm 0.33}$	$3.70{\scriptstyle\pm 2.87}$	$-3.06{\scriptstyle\pm 0.24}$	$13.41{\scriptstyle\pm 1.94}$

Appendix

Table 11: Effect of particle ordering on DSR performance for Lorenz-63 at SMNR

=10

and

0

dB. Results show mean

\pm

std over 10 runs.

Method	NMSE (dB)	MSE	MSE-5-step	LL	KLD	Time/epoch (s)
	State Est.	Transition Model Quality
EROT ( $\varepsilon=0.5$ )	$-12.57{\scriptstyle\pm 1.64}$	$4.79{\scriptstyle\pm 5.65}$	$48.41{\scriptstyle\pm 60.72}$	$-4.28{\scriptstyle\pm 1.54}$	$33.82{\scriptstyle\pm 29.79}$	$42.94{\scriptstyle\pm 0.77}$
EROT ( $\varepsilon=1.0$ )	$-13.47{\scriptstyle\pm 0.95}$	$0.57{\scriptstyle\pm 0.43}$	$9.19{\scriptstyle\pm 7.58}$	$-3.49{\scriptstyle\pm 0.09}$	$20.20{\scriptstyle\pm 2.33}$	$40.35{\scriptstyle\pm 0.39}$
DSR ( $\tau=0.1$ )	$\mathbf{-18.28}{\scriptstyle\pm 0.46}$	$\mathbf{0.54}{\scriptstyle\pm 0.33}$	$\mathbf{3.70}{\scriptstyle\pm 2.87}$	$\mathbf{-3.06}{\scriptstyle\pm 0.24}$	$\mathbf{13.41}{\scriptstyle\pm 1.94}$	$\mathbf{18.39}{\scriptstyle\pm 0.16}$

Appendix

Table 12: EROT regularization sensitivity on Lorenz-63 at SMNR =

0

dB. Results show mean

\pm

std over

10

runs. DSR with

\tau=0.1

is included for reference. Best in bold .

Method	$d=1$	$d=2$	$d=3$	$d=5$	$d=10$
Hard SR	$0.000{\scriptstyle\pm 0.000}$	$0.000{\scriptstyle\pm 0.000}$	$0.000{\scriptstyle\pm 0.000}$	$0.000{\scriptstyle\pm 0.000}$	$0.000{\scriptstyle\pm 0.000}$
Soft	$0.000{\scriptstyle\pm 0.000}$	$0.000{\scriptstyle\pm 0.000}$	$0.000{\scriptstyle\pm 0.000}$	$0.000{\scriptstyle\pm 0.000}$	$0.000{\scriptstyle\pm 0.000}$
EROT	$0.000{\scriptstyle\pm 0.000}$	$0.000{\scriptstyle\pm 0.000}$	$0.747{\scriptstyle\pm 0.052}$	$1.000{\scriptstyle\pm 0.000}$	$1.000{\scriptstyle\pm 0.000}$
DSR (Morton)	$0.287{\scriptstyle\pm 0.005}$	$0.283{\scriptstyle\pm 0.005}$	$0.282{\scriptstyle\pm 0.005}$	$0.282{\scriptstyle\pm 0.005}$	$0.282{\scriptstyle\pm 0.005}$
DSR (w/o Morton)	$1.000{\scriptstyle\pm 0.000}$	$1.000{\scriptstyle\pm 0.000}$	$1.000{\scriptstyle\pm 0.000}$	$1.000{\scriptstyle\pm 0.000}$	$1.000{\scriptstyle\pm 0.000}$

Appendix

Table 13: Inter-modal mass (resampled mass in

|x_{0}|<3

) for a static two-component Gaussian as a function of dimension

d

. Lower is better;

0

means no mass placed between the modes. Results show mean

\pm

std over

10^{4}

simulations.

Method	$\log p_{\hat{\boldsymbol{\phi}}}$	$\mathbb{E}[\log\hat{p}_{N}]$	model error	Jensen’s gap
DSR	$-232.942{\scriptstyle\pm 0.087}$	$\mathbf{-235.138}{\scriptstyle\pm 0.385}$	$0.575{\scriptstyle\pm 0.087}$	$\mathbf{2.196}{\scriptstyle\pm 0.302}$
EROT	$-233.065{\scriptstyle\pm 0.104}$	$-235.851{\scriptstyle\pm 0.414}$	$0.699{\scriptstyle\pm 0.104}$	$2.786{\scriptstyle\pm 0.467}$
Hard	$\mathbf{-232.795}{\scriptstyle\pm 0.038}$	$-236.147{\scriptstyle\pm 0.924}$	$\mathbf{0.428}{\scriptstyle\pm 0.038}$	$3.352{\scriptstyle\pm 0.901}$
Soft	$-232.996{\scriptstyle\pm 0.138}$	$-237.215{\scriptstyle\pm 0.749}$	$0.630{\scriptstyle\pm 0.138}$	$4.218{\scriptstyle\pm 0.669}$

Appendix

Table 14: Evaluation-time evidence decomposition on the linear Gaussian SSM (

T=100

, SMNR

=10

dB). The true log-evidence is

-232.367

. Results show mean

\pm

std over model seeds.

(a) $d=2$ , $N=4$ , $T=4$ ; $\\|\nabla\mathcal{L}\\|=9.64$
Method	Bias	Variance	MSE
Hard SR	$2.324{\scriptstyle\pm 0.014}$	$0.308{\scriptstyle\pm 0.016}$	$5.707{\scriptstyle\pm 0.067}$
DSR ( $\tau=1.0$ )	$1.256{\scriptstyle\pm 0.000}$	$0.241{\scriptstyle\pm 0.014}$	$1.818{\scriptstyle\pm 0.014}$
DSR ( $\tau=0.5$ )	$1.257{\scriptstyle\pm 0.001}$	$0.245{\scriptstyle\pm 0.014}$	$1.825{\scriptstyle\pm 0.014}$
DSR ( $\tau=0.1$ )	$0.764{\scriptstyle\pm 0.020}$	$0.389{\scriptstyle\pm 0.031}$	$0.973{\scriptstyle\pm 0.043}$
DSR ( $\tau=0.05$ )	$0.370{\scriptstyle\pm 0.050}$	$0.723{\scriptstyle\pm 0.037}$	$0.860{\scriptstyle\pm 0.052}$

Appendix

Table 15: Gradient bias, variance, and MSE against the unbiased VSMC gradient on the linear Gaussian SSM. Hard SR drops the score-function term. Results show mean

\pm

standard errors.

Explore similar work

Oct 29, 2025eess.SP

PyDPF: A Python Package for Differentiable Particle Filtering

State-space models (SSMs) are a widely used tool in time series analysis. In the complex systems that arise from real-world data, it is common to employ particle filtering (PF), an efficient Monte Carlo method for estimating the hidden state corresponding to a sequence of observations. Applying particle filtering requires specifying both the parametric form and the parameters of the system, which are often unknown and must be estimated. Gradient-based optimisation techniques cannot be applied directly to standard particle filters, as the filters themselves are not differentiable. However, several recently proposed methods modify the resampling step to make particle filtering differentiable. In this paper, we present an implementation of several such differentiable particle filters (DPFs) with a unified API built on the popular PyTorch framework. Our implementation makes these algorithms easily accessible to a broader research community and facilitates straightforward comparison between them. We validate our framework by reproducing experiments from several existing studies and demonstrate how DPFs can be applied to address several common challenges with state space modelling.

John-Joseph Brady, Benjamin Cox, Yunpeng Li +1

King’s College London · Max-Planck Institut für Physik · University of Edinburgh

Particle Filtering Bayesian Filtering Differentiable Programming State Space Models

Ask PDF

May 19, 2026cs.LG

Training-Free Bayesian Filtering with Generative Emulators

Bayesian filtering is a well-known problem that aims to estimate plausible states of a dynamical system from observations. Among existing approaches to solve this problem, particle filters are theoretically exact for non-linear dynamics and observations, but suffer from poor scalability in high dimensions. In this work, we show that diffusion-based emulators of dynamical systems can be used to implement, without additional training, an optimal variant of particle filters that has remained largely unexplored due to implementation challenges with classical numerical solvers. Experiments on nonlinear chaotic systems, including atmospheric dynamics, demonstrate that the proposed approach successfully scales particle filtering to high-dimensional settings.

Thomas Savary, François Rozet, Gilles Louppe

SAIL, Montefiore institute, University of Liège, Belgium.

Particle Filtering Bayesian Filtering Data Assimilation

Ask PDF

May 20, 2026cs.LG

Efficient Learning of Deep State Space Models via Importance Smoothing

Latent state space systems are ubiquitous in statistical modelling, arising naturally when time series are observed through noisy measurements. However, training deep state space models (DSSMs) at scale remains difficult. Two largely distinct strategies have emerged for training DSSMs. The first, auto-encoding DSSMs, trains generative models by optimising a variational lower bound. The second backpropagates through the outputs of classical sequential Monte Carlo (SMC) algorithms. Such approaches can train DSSMs for both discriminative and generative tasks, but their inherently sequential forward passes scale poorly on modern hardware. We propose \emph{parallel variational Monte Carlo} (PVMC), a new training method that bridges these paradigms and robustly trains DSSMs for both discriminative and generative tasks. Across a set of benchmark experiments, PVMC matches or exceeds state-of-the-art performance while training

10\times

faster than the fastest competing SMC-based approach.

John-Joseph Brady, Nikolas Nusken, Yunpeng Li

Centre for Oral, Clinical and Translational Sciences, King’s College London, London, United Kingdom · Department of Mathematics, King’s College London, London, United Kingdom

Variational Inference State Space Models

Ask PDF