Organizations: School of Mathematical Sciences, Tongji University, Shanghai 200092, China · Department of Information and Computational Sciences, Wuhan University, Wuhan 430072, China · Department of Mathematics, University of South Carolina, Columbia, South Carolina 29208, USA
Bayesian state estimation for high-dimensional nonlinear dynamical systems entails a fundamental tension between statistical fidelity and computational tractability, as particle weights can collapse, while Gaussian ensemble updates can miss non-Gaussian posterior structure. Score-based diffusion filters offer a sampling-based alternative, but existing training-free score filters often rely on heuristic likelihood corrections, which can compromise posterior accuracy by neglecting uncertainty about the system state associated with each noisy reverse particle. To address these issues, we propose AECSF, a training-free adaptive ensemble conditional score filter. AECSF constructs an analytically tractable score estimator from the conditional Tweedie identity, which recasts noisy posterior score estimation as estimating the conditional mean of the system state given a noisy reverse particle and the observation. To estimate these conditional means efficiently, AECSF employs a shared adaptive weighted proposal ensemble, while particle-specific conditional weights yield an estimate for each noisy reverse particle without separate proposal sampling. The proposal ensemble is updated using reverse-particle information within the same reverse-diffusion run to improve conditional-mean estimation. Theoretically, we characterize when a fixed weighted proposal measure yields the exact noisy posterior score. Under stated assumptions, we establish a bound relating conditional-mean estimation errors to reverse-sampling endpoint error. Numerical experiments demonstrate that AECSF improves the accuracy of posterior sampling and nonlinear filtering in high-dimensional problems with limited forecast ensembles.
Figures & tables
Figure 1: Gaussian endpoint energy discrepancies ( M=J=20 ). A pools four repetitions (80 particles/problem); B excludes same-run particle pairs and reference self-pairs. Colors denote fitted and data-generating retained-noise references. Large points/bars show means/sample SDs across ten problems; small points show individual problems. Definitions: Appendix D.1 .
Method
RMSE ↓
CRPS ↓
SSR →1
Norm. obs. misfit
Analysis time/cycle (s)
EnSF
0.2418±0.0006
0.1432±0.0002
1.7366±0.0042
1.8180±0.0028
0.0485±0.0010
LETKF
0.1263±0.0116
0.0578±0.0015
0.6119±0.0604
1.0393±0.0349
9.5530±0.1153
Binder et al.
5.1323±0.0082
4.1428±0.0076
(24.66±13.29)×10−5
28.4641±0.0390
0.1507±0.0015
AECSF
0.1545±0.0039
0.0776±0.0007
0.8998±0.0232
0.9209±0.0042
0.3052±0.0139
Table 1: Lorenz–96 results ( d=10,000 , 20 members; J=M=20 for AECSF). Entries are means ± sample SD across ten seeds, each summarized over cycles 50–199. Timings use serial runs on the same host with GPU analysis for all methods.
Figure 2: Lorenz–96 proposal-update ablations ( d=10,000 , M=20 ). Langevin variants share the reverse schedule and 294 steps. A–C summarize cycles 50–199; D gives timings for four variants. Large points/bars show means/sample SDs across ten seeds; small points show paired seed results. Fixed proposal positions retains weight updates; target-weight resampling uses ω .
Figure 3: Kuramoto–Sivashinsky accuracy versus analysis time (logarithmic horizontal axes). A: sampling steps at J=20 ; labels give T or Nt for IEnSF-L1. B: ensemble size at 20 steps; labels give J . AECSF varies forecast, reverse-particle and proposal counts together ( J=M ). IEnSF-L1 extends to J=80 and the other methods to J=160 . RMSE uses cycles 950–999 and timing all 1,000 cycles. Points/bars show means/sample SDs over five paired seeds; lines follow increasing steps or ensemble size.
Appendix figures & tables12 assets
Supplementary material from the paper’s appendix.
Appendix
Method
Fitted EV
True EV
Fitted Ecross
True Ecross
Seconds/run
Analytic conditional score (fitted prior)
0.008850
0.018799
0.000284
0.010230
0.069935
AECSF
0.019040
0.035886
0.009168
0.026011
0.300461
Movement-only
0.070636
0.090871
0.011818
0.032051
0.216578
Pointwise guidance
0.038129
0.044543
0.031485
0.037897
0.062395
Denoised guidance
0.098303
0.081641
0.091629
0.074965
0.066121
Appendix
Table 2: Gaussian endpoint summaries across ten problems. Each problem uses four repetitions of 20 particles. The two energy statistics are defined in Eqs. ( 70 )–( 71 ).
Method
Fitted mean MSE
Fitted variance MSE
True mean MSE
True variance MSE
Analytic conditional score (fitted prior)
0.011516
0.005614
0.018309
0.005799
AECSF
0.015528
0.008179
0.025495
0.009191
Movement-only
0.093729
0.028066
0.105807
0.029772
Pointwise guidance
0.020602
0.009073
0.024058
0.010208
Denoised guidance
0.058519
0.009065
0.047680
0.009776
Appendix
Table 3: Mean and diagonal-variance errors of individual 20-member Gaussian outputs. Each entry averages four repetitions within each of ten problems and then averages across problems. Both references include the retained endpoint noise.
Setting
Gaussian
Lorenz–96
Kuramoto–Sivashinsky
M/J
20/20
20/20
J/J
Reverse steps
100
100
T
Weight updates / resampling / move calls
99/50/98
99/50/98
(T−1)/(T/2)/(T−2)
Langevin steps per move
3
3
3
Langevin step size h
0.02
0.006
0.002
ϵa
0.001
0.0005
0.001
Appendix
Table 4: AECSF configurations. M counts proposal members and J counts reverse particles; the forecast ensemble has J members. The main benchmark configurations use M=J ; Table 8 varies M at fixed J . All configurations use stratified resampling and unit forecast and likelihood gains.
Method
RMSE ↓
CRPS ↓
SSR →1
Norm. obs. misfit
Sampling-update time (s/cycle)
AECSF
0.154482±0.003892
0.077631±0.000740
0.899823±0.023152
0.920890±0.004196
0.305238±0.013899
Direct posterior ULA
0.184311±0.007514
0.078474±0.001215
0.596951±0.023914
1.042469±0.013634
0.136624±0.003211
Appendix
Table 5: Matched-step Lorenz–96 comparison with 20 members at dimension 10,000. Both methods use 294 Langevin steps per analysis with step size 0.006. Entries are means and sample standard deviations across seeds 10–19, using cycles 50–199 within each seed. Timing brackets the sampling update and excludes model construction and final analysis processing.
Method
RMSE ↓
CRPS ↓
SSR →1
Norm. obs. misfit
Update time (s/cycle)
AECSF
0.154968±0.004083
0.077842±0.000623
0.896359±0.024766
0.921380±0.005849
0.289720±0.001397
Direct posterior ULA
0.171986±0.001801
0.094169±0.000456
1.116732±0.010192
1.149319±0.002589
0.296331±0.000446
Appendix
Table 6: Approximately time-matched posterior sampling on Lorenz–96 with 20 members in 10,000 dimensions. AECSF uses 294 Langevin steps; direct ULA uses 1,350 steps with step size 0.012. Entries are means and sample standard deviations across ten seeds, using cycles 50–199 within each seed.
Steps in first/last 49 stages
h
RMSE ↓
CRPS ↓
SSR →1
3/3
0.0060
0.154482
0.077631
0.899823
4/2
0.0060
0.156707
0.079617
0.939406
2/4
0.0060
0.154794
0.077192
0.872296
3/3
0.0045
0.154682
0.078449
0.931029
3/3
0.0075
0.155021
0.077493
0.879071
Appendix
Table 7: Local movement sensitivity on Lorenz–96. Entries are means across seeds 10–19 over cycles 50–199. Every setting uses 294 Langevin steps; varying the step size also changes the total Langevin integration time.
M
M/J
RMSE ↓
CRPS ↓
SSR →1
s/update
20
1
0.153809
0.079206
0.8935
0.2995
40
2
0.148537
0.077347
0.9335
0.3018
80
4
0.148571
0.077442
0.9395
0.2998
160
8
0.148918
0.077408
0.9428
0.2969
240
12
0.150529
0.077595
0.9352
0.2985
320
16
0.147748
0.077230
0.9571
0.2986
Appendix
Table 8: Proposal capacity in Lorenz–96 with d=1000 and J=20 . Entries are means over seeds 10–19, each summarized over cycles 50–199. All configurations run on one RTX 4090. The per-proposal Langevin budget is fixed, so total gradient work increases with M . Time measures the sampling update.
Figure 4: Lorenz–96 EnSF d=1000 tuning grid over seeds 0,…,9 . The rank-1 configuration transferred to Table 1 is ϵa=0.6,ϵb=0.025 .
Figure 5: Lorenz–96 LETKF d=1000 tuning grid over seeds 0,…,9 . The rank-1 configuration transferred to Table 1 uses localization radius 2 , neighbor size 5 , and inflation 1.1 .
Method
J
Steps
RMSE
CRPS
Energy
SSR
s/cycle
AECSF
20
20
0.0889±0.0018
0.0625±0.0005
0.0792±0.0007
2.1560±0.0417
0.0988±0.0040
AECSF
20
100
0.0601±0.0040
0.0431±0.0012
0.0545±0.0017
2.2261±0.1376
0.4681±0.0188
EnSF
20
20
0.4466±0.0047
0.2722±0.0027
0.3673±0.0038
1.8207±0.0106
0.0297±0.0006
EnSF
20
100
0.1206±0.0006
0.0684±0.0003
0.0876±0.0004
1.0333±0.0075
0.1480±0.0024
EnFF-F2P
20
20
0.2096±0.0014
0.1569±0.0013
0.2076±0.0014
0.0135±0.0006
0.0237±0.0004
EnFF-F2P
20
100
0.2077±0.0014
0.1557±0.0013
0.2058±0.0014
0.0138±0.0006
0.1147±0.0034
Appendix
Table 9: Representative Kuramoto–Sivashinsky results with 20 members. Entries are means ± sample standard deviations across seeds 30–34. Quality metrics use cycles 950–999; analysis time uses all 1,000 cycles. AECSF uses ϵb=0.015 . Steps denote T for AECSF, EnSF and EnFF, and Nt for IEnSF-L1.
Method
J
Steps
RMSE
CRPS
Energy
SSR
s/cycle
AECSF
20
20
0.0889±0.0018
0.0625±0.0005
0.0792±0.0007
2.1560±0.0417
0.0988±0.0040
AECSF
20
50
0.0638±0.0016
0.0467±0.0005
0.0591±0.0006
2.3116±0.0612
0.2402±0.0086
AECSF
20
100
0.0601±0.0040
0.0431±0.0012
0.0545±0.0017
2.2261±0.1376
0.4681±0.0188
AECSF
20
200
0.0598±0.0035
0.0418±0.0010
0.0530±0.0015
2.1292±0.1276
0.9723±0.0179
AECSF
20
500
0.0579±0.0027
0.0405±0.0008
0.0514±0.0011
2.1330±0.0982
2.3375±0.0765
EnSF
20
20
0.4466±0.0047
0.2722±0.0027
0.3673±0.0038
1.8207±0.0106
0.0297±0.0006
Appendix
Table 10: Complete Kuramoto–Sivashinsky sampling-budget comparison at J=20 . Parameter settings, statistical summaries and step notation follow Table 9 .
Method
J
Steps
RMSE
CRPS
Energy
SSR
s/cycle
AECSF
10
20
0.1155±0.0047
0.0738±0.0017
0.0942±0.0025
1.6032±0.0687
0.0992±0.0032
AECSF
40
20
0.0732±0.0012
0.0571±0.0003
0.0722±0.0004
2.6609±0.0429
0.0998±0.0033
AECSF
80
20
0.0650±0.0012
0.0546±0.0003
0.0691±0.0004
3.0259±0.0558
0.1034±0.0033
AECSF
160
20
0.0593±0.0010
0.0531±0.0002
0.0671±0.0003
3.3327±0.0530
0.1038±0.0031
EnSF
10
20
0.4867±0.0056
0.2943±0.0035
0.3982±0.0036
1.6259±0.0182
0.0297±0.0006
EnSF
40
20
0.4259±0.0042
0.2625±0.0030
0.3525±0.0032
1.9366±0.0116
0.0384±0.0007
Appendix
Table 11: Kuramoto–Sivashinsky ensemble-size comparison at 20 sampling steps. The shared J=20 settings are reported in Table 10 . Parameter settings, statistical summaries and step notation follow Table 9 .
Bayesian filtering of partially and noisily observed dynamical systems seeks to infer the evolving conditional distribution of the state of a dynamical system given observations, in an online fashion. This Bayesian filtering distribution is rarely available as a supervised learning target. However, one can often use the forecast model to generate synthetic trajectories, with corresponding synthetic observations. We introduce the proper scoring ensemble filter (PSEF), an ensemble data assimilation method trained using only synthetic trajectories. The analysis step is represented as a permutation-equivariant, transformer-based map. Training is based on strictly proper scoring rules---with the energy score used in our implementation---so that probabilistic accuracy is rewarded over the whole probability distribution. Under a realizability assumption, the population mean-field objective is minimized by the true Bayesian filtering distribution. Our methodology allows the same learned parameters to be shared between different ensemble sizes, subject to an ensemble-dependent fine-tuning. Numerical experiments show that the learned filter accurately approximates challenging filtering distributions, including highly non-Gaussian and multi-modal posteriors, and achieves stronger performance in data assimilation tasks than classical methods or learning-based methods with mean-squared-error objectives.
Eviatar Bach, Ricardo Baptista, Jochen Bröcker +2
Department of Meteorology, University of Reading Brian Hoskins Building, Reading, RG6 6ET, UK · Department of Mathematics and Statistics, University of Reading Pepper Lane, Reading, RG6 6AX, UK · National Centre for Earth Observation Brian Hoskins Building, University of Reading, Reading, RG6 6ET, UK +2
Data assimilation is the process of estimating the state of a dynamical system over time by combining model predictions with measurements. This task becomes challenging when the system is nonlinear and high-dimensional. To address this, score-based Bayesian filters have recently emerged. However, these methods still show unsatisfactory performance in certain cases, particularly under spatially sparse measurements. Such degradation stems from heuristic approximations of the likelihood score, whose errors can accumulate over time. This limitation arises because the methods simply adopt a classical forward process for generative modeling that transforms a data distribution toward a Gaussian distribution, which is independent of the measurement equation. Here, we propose a forward process tailored for filtering that transforms the system state toward the measurement space, enabling a theoretically sound formulation of the likelihood score. Based on this, we develop the Measurement-Aware Score-based Filter (MASF). We evaluate MASF on Kolmogorov flow, a high-dimensional fluid benchmark with up to O(105) dimensions, under diverse measurement operators, including nonlinear cases with a dimensional mismatch between the state and the measurements. MASF shows improved performance over existing score-based filters and ensemble-type Kalman filters. Notably, MASF achieves up to a 28.2× wall-clock speedup compared with the baselines when using amortized pretraining. Our implementation is available at https://github.com/tcnllab-oss/masf.
A central obstacle in nonlinear Bayesian filtering is representing the belief distribution. Moment-based filters address this by propagating polynomial moments and reconstructing a density from them. Recent work completes the predict-update loop via the maximum-entropy (MaxEnt) principle, but each step requires the partition function and its gradient, both n-dimensional integrals whose cost scales exponentially, restricting the demonstrated MaxEnt moment filtering to n≤4. We avoid the partition function entirely by combining score matching with Stein's identity. In our setting, score matching reduces the density fit to a single linear solve whose coefficients are assembled directly from the propagated moments. The same parameters then drive Stein's identity to close the moment hierarchy during prediction and to recover posterior moments after each Bayesian update, keeping the full predict-update loop free of partition function evaluation. The resulting Score Kalman Filter (SKF) reduces to the classical information-form Kalman filter as a special case and performs every step through linear algebra. On nonlinear coupled-oscillator networks, the SKF runs through n=20 and reports lower RMSE than the EKF, UKF, EnKF, and particle-filter baselines on the tested synthetic benchmarks.
Kaito Iwasaki, Anthony Bloch, Taeyoung Lee +1
Department of Mathematics University of Michigan · Department of Mechanical and Aerospace Engineering George Washington University · Department of Naval Architecture & Marine Engineering and Department of Robotics University of Michigan