PSI-SINDy: Post-Selection Inference for Sparse Identification of Nonlinear Dynamics

Ashraful Islam; Shuichi Nishino; Tomohiro Shiraishi; Ichiro Takeuchi

stat.MLOct 8, 2026

PSI-SINDy: Post-Selection Inference for Sparse Identification of Nonlinear Dynamics

Authors: Ashraful Islam, Shuichi Nishino, Tomohiro Shiraishi, Ichiro Takeuchi

Organizations: RIKEN · Nagoya University

Abstract

Sparse identification of nonlinear dynamics (SINDy) is a data-driven framework for discovering governing dynamics from time-series data by identifying a sparse subset of candidate dynamical terms from a prespecified library. In this work, we develop a statistical inference framework for quantifying the reliability of dynamical terms selected by SINDy through hypothesis tests and confidence intervals. A key difficulty is that using the same noisy trajectory for both selecting dynamical terms and assessing their statistical significance can introduce selection bias. Post-selection inference provides a principled framework for addressing such bias, and we propose PSI-SINDy, a post-selection inference method tailored to SINDy. Direct application of existing post-selection inference techniques is challenging because SINDy involves measurement error in the candidate terms and shared noise between the response and design. To address these challenges, PSI-SINDy uses data thinning to decompose a single observed trajectory into four mutually independent views with distinct roles in selection and inference. This construction enables inference for selected dynamical terms while accounting not only for selection bias but also for measurement-error and shared noise effects. We establish the theoretical validity of PSI-SINDy under stated conditions and evaluate its performance through numerical experiments on simulated and experimental dynamical-system data.

Figures & tables

Figure 1: Primary selected-BLA inference conditional on selecting

x^{2}

(

\sigma_{X}=0.06

\delta=0

). (a) Rejection at each support’s true BLA,

\gamma=0

; dotted line at 0.05. (b) Power against the reference null ( 18 ),

M=800

. (c) Median width of PSI-SINDy’s 95% BLA intervals,

\gamma=0

; counts and MC intervals in Tables 2 – 3 .

$\sigma_{X}$	Term	$N_{\mathcal{S}}$	Full	No gap	No feature correction
0.06	$1$	812	5.42	7.51	6.28
0.06	$x$	812	5.30	7.51	$5.30^{\dagger}$
0.06	$x^{2}$	812	5.42	6.28	$5.42^{\dagger}$
0.12	$1$	452	5.53	7.30	11.06
0.12	$x$	452	6.86	8.85	$6.86^{\dagger}$
0.12	$x^{2}$	452	5.53	6.86	$5.53^{\dagger}$

Table 1: Rejection at the true selected BLA (%; nominal 5%) in the matched polynomial study:

\gamma=0.20

\delta=0

M=800

\mathcal{S}=\{1,2,3\}

N_{\mathcal{S}}

counts datasets selecting this support. Full uses both adjustments. Paired differences and MC intervals are in Table 8 .

Appendix figures & tables17 assets

Supplementary material from the paper’s appendix.

Appendix

		Selected-BLA target
$M$	$N_{A}$	Naive	Bonferroni	PSI-SINDy
200	3425	7.36 [6.53,8.28]	1.84 [1.44,2.35]	4.41 [3.77,5.15]
800	3360	7.86 [6.99,8.82]	2.02 [1.60,2.56]	5.39 [4.67,6.20]
1600	3423	7.36 [6.53,8.29]	1.61 [1.24,2.09]	4.91 [4.23,5.68]
3200	3452	8.69 [7.80,9.68]	1.74 [1.35,2.23]	5.27 [4.58,6.07]
		Realized-IV target

Appendix

Table 2: Primary quadratic calibration, conditional on selection. Both panels report rejection percentages with pointwise 95% Wilson intervals at the stated true target. Widths are median coefficient-interval widths. Bonferroni uses multiplier four.

$\gamma$	$N_{A}$	Bonferroni	PSI-SINDy	Bon. joint	PSI joint
0.02	3540	101 2.85 [2.35, 3.45]	239 6.75 [5.97, 7.63]	2.02	4.78
0.04	3901	210 5.38 [4.72, 6.14]	460 11.79 [10.82, 12.84]	4.20	9.20
0.08	4440	718 16.17 [15.12, 17.28]	1280 28.83 [27.52, 30.18]	14.36	25.60
0.16	4896	2306 47.10 [45.70, 48.50]	3175 64.85 [63.50, 66.17]	46.12	63.50

Appendix

Table 3: Selected-BLA power. The Bonferroni and PSI-SINDy cells give rejection counts, conditional power percentages and pointwise 95% Monte Carlo intervals, with denominator

N_{A}

. The last two columns give joint selection-and-rejection percentages with denominator 5,000.

Method	Type I error	95% MC interval for Type I error	Structural coverage	Median width
PSI-SINDy	4.70	[3.71, 5.95]	95.30	0.1914
Ensemble, SG5	0.00	[0.00, 0.28]	100.00	2.9573
Ensemble, longer SG	8.10	[6.78, 9.66]	91.90	0.1675
Residual, matched support	20.48	[18.43, 22.69]	79.52	0.1676

Appendix

Table 4: Focal structural-null comparison at

M=800

\sigma_{X}=0.06

, and

\gamma=\delta=0

. All methods are assessed on the same 1,382 datasets, out of 2,000, in which PSI selects the quadratic term. Type I error, structural-coefficient coverage, and pointwise 95% Wilson Monte Carlo intervals are expressed as percentages. Monte Carlo intervals refer to Type I error. Width is the median finite coefficient-interval width.

System	$M$	$n$	Slope BLA coverage	Width	Family events	FWER
Polynomial	200	452	92.92 [90.18,94.94]	0.11091	14	3.10 [1.85,5.13]
Polynomial	800	812	94.70 [92.94,96.04]	0.05536	32	3.94 [2.81,5.51]
Polynomial	1600	1046	94.55 [93.00,95.77]	0.03917	45	4.30 [3.23,5.71]
Polynomial	3200	1234	94.17 [92.72,95.34]	0.02769	40	3.24 [2.39,4.38]
Sine-augmented	200	673	94.50 [92.51,95.99]	0.07656	47	6.98 [5.29,9.16]
Sine-augmented	800	1096	94.62 [93.12,95.80]	0.03818	54	4.93 [3.80,6.37]

Appendix

Table 5: All prespecified support-conditioned groups. Coverage and FWER entries are percentages with pointwise 95% Wilson intervals.

n=N_{\mathcal{S}}

; median width is in coefficient units. Simultaneous Bonferroni coverage is

100-\mathrm{FWER}

Figure 2: Primary IV–BLA bridge accuracy conditional on selecting

x^{2}

, with

\gamma=\delta=0

. (a) Median absolute first-order approximation error divided by the oracle scale. (b) Median estimated-to-oracle gap-scale ratio; reference one. In ( D.3 ), the oracle scale is the standard deviation of the scaled first-order gap. Both

M

axes are logarithmic. Table 6 gives complementary support-conditioned diagnostics.

System	$M$	IV rejection (%)	Scale ratio	Remainder	SD
Polynomial	200	6.19 [4.32, 8.81]	0.986	0.0267	1.005
Polynomial	800	5.42 [4.06, 7.20]	0.996	0.0154	1.010
Polynomial	1600	6.12 [4.82, 7.74]	0.999	0.0103	1.003
Polynomial	3200	4.86 [3.80, 6.21]	0.998	0.0074	1.008
Sine-augmented	200	5.35 [3.89, 7.32]	0.998	0.0145	1.005
Sine-augmented	800	5.02 [3.88, 6.48]	0.999	0.0076	1.002

Appendix

Table 6: IV and bridge diagnostics for the same slope and exact-support groups as Table 5 , with the same assessment counts. Brackets give pointwise 95% Wilson Monte Carlo intervals for IV rejection. Here

\widehat{t}=\widehat{t}_{j,\mathcal{S}}

and

\tau=\tau_{M,j,\mathcal{S}}

denote the estimated gap scale and oracle standard deviation of the scaled first-order gap for the stated coordinate and support. The scale ratio is the median

\sqrt{M}\widehat{t}/\tau

; the remainder is the median absolute scaled linearization error divided by

\tau

. SD denotes the empirical standard deviation of the studentized BLA statistic.

Study	Term	$M$	$n$	BLA coverage	Remainder	Scale ratio
Primary	$x^{2}$	200	3425	95.59	0.0607	1.0195
Primary	$x^{2}$	3200	3452	94.73	0.0156	1.0002
Sine drift	$x^{2}$	800	2989	87.82	0.8341	0.9232
Pendulum	$xy^{2}$	1600	275	85.82	0.9562	1.4048

Appendix

Table 7: Bridge diagnostics in the primary regime and two cases of undercoverage. Each row conditions on inclusion of the stated term and retains each support’s own target.

n

is the selected count; remainder and scale ratio are medians of ( D.3 ). Primary and sine rows use

\sigma_{X}=0.06

; the sine drift is

\dot{x}=-2\sin x

. The pendulum row uses its second equation and

\sigma_{X}=0.05

. Coverage is expressed as a percentage.

$\sigma_{X}$	$j$	$N_{\mathcal{S}}$	Full-method rejection (%)	Without gap scale	Without feature correction
0.06	1	812	5.42	$+2.09$ $[1.11,3.08]$	$+0.86$ $[-0.72,2.44]$
0.06	2	812	5.30	$+2.22$ $[1.20,3.23]$	$0^{\dagger}$
0.06	3	812	5.42	$+0.86$ $[0.23,1.50]$	$0^{\dagger}$
0.12	1	452	5.53	$+1.77$ $[0.55,2.99]$	$+5.53$ $[2.73,8.33]$
0.12	2	452	6.86	$+1.99$ $[0.70,3.28]$	$0^{\dagger}$
0.12	3	452	5.53	$+1.33$ $[0.27,2.38]$	$0^{\dagger}$

Appendix

Table 8: Matched component assessment at

M=800

and

\mathcal{S}=\{1,2,3\}

. Here

N_{\mathcal{S}}

is the number of datasets selecting exactly

\mathcal{S}

. The full-method column reports rejection at the true selected BLA (nominal 5%). The last columns give ablated-minus-full changes in percentage points with approximate paired 95% Monte Carlo intervals. Shared measurements, selection, views and bootstrap randomness give matched comparisons. A dagger marks unchanged decisions because, for the selected basis

\{1,x,x^{2}\}

, the intercept absorbs the quadratic feature correction.

	Rejection at the true selected BLA
$\sigma_{X}$	Full	Uncorrected	OLS
0.12	4.76 [4.20,5.39]	6.74 [6.08,7.47]	9.64 [8.85,10.49]
0.24	4.46 [3.92,5.07]	5.82 [5.20,6.50]	11.98 [11.11,12.91]

Appendix

Table 9: Additional intercept assessment at

M=800

, pooled over selected supports,

n=5000

per noise level. Entries are percentages with pointwise 95% Wilson intervals. OLS denotes ordinary least squares after whitening. The full and no-gap intervals target the same selected BLA.

Dynamics	$P$	Size	$n/500$	Coverage [MC]	Width	Size	$n/500$	Coverage [MC]	Width
		$M=800$				$M=3200$
Quadratic	4	3.010/4	473/500	95.77 [93.56,97.25]	0.1615	3.004/4	500/500	95.40 [93.19,96.92]	0.0806
	6	3.034/4	434/500	95.85 [93.54,97.36]	0.1614	3.002/4	495/500	95.56 [93.36,97.05]	0.0806
	8	3.052/5	426/500	95.54 [93.14,97.13]	0.1614	3.018/4	487/500	95.48 [93.26,97.00]	0.0806
	10	3.054/5	426/500	95.77 [93.42,97.31]	0.1614	3.014/4	489/500	95.50 [93.28,97.01]	0.0806
	15	3.056/5	424/500	95.75 [93.39,97.30]	0.1614	3.016/4	488/500	95.49 [93.27,97.00]	0.0806

Appendix

Table 10: Complete quadratic and quartic grids for the integrated ridge-STLSQ selector. Size is the selected-support mean/maximum over all 500 datasets;

n

counts focal selection (

x^{2}

for quadratic dynamics and

x^{4}

for quartic dynamics). Coverage is a percentage with its pointwise 95% Wilson interval, conditional on focal selection and pooled across supports. Width is the median selected-BLA interval width on the same assessed datasets. The quartic focal term is absent from the library at

P=4

, so its coverage and width are undefined. There are no unavailable intervals or unassessable selected targets in these cells.

Figure 3: Columns compare weak and integrated ridge-STLSQ selectors for quadratic and quartic dynamics; rows show exact support recovery, focal selected-BLA coverage, and median interval width. The focal term is

x^{2}

for quadratic dynamics and

x^{4}

for quartic dynamics. Both selectors and all library sizes reuse 500 datasets per dynamics/

M

cell. Recovery uses all datasets; coverage and width condition on focal selection and pool selected supports. Dotted lines mark nominal 0.95 coverage. Quartic

P=4

points are omitted: recovery is zero and focal inference is undefined because

x^{4}

is absent. Table 10 gives integrated-selector counts and Monte Carlo intervals.

Equation	Term	$n$	BLA	No gap	IV
VdP (2)	$1$	2000	94.85 [93.79,95.74]	13.85 [12.41,15.43]	4.25 [3.45,5.23]
VdP (2)	$x^{2}y$	2000	95.30 [94.28,96.14]	75.45 [73.52,77.29]	5.25 [4.36,6.32]
LV (1)	$xy$	2000	94.65 [93.58,95.55]	56.55 [54.37,58.71]	4.45 [3.63,5.44]
Lorenz (2)	$1$	2000	94.00 [92.87,94.96]	15.00 [13.50,16.63]	4.80 [3.95,5.83]
Lorenz (2)	$xz$	2000	98.00 [97.29,98.53]	55.50 [53.31,57.67]	4.90 [4.04,5.94]
Pendulum (2)	$x$	2000	90.50 [89.14,91.71]	78.75 [76.90,80.49]	4.50 [3.68,5.50]

Appendix

Table 11: Multistate assessment at

M=1600

\sigma_{X}=0.05

. BLA and No gap give selected-BLA coverage; IV gives rejection at the realized-IV target. Entries are percentages with pointwise 95% Wilson intervals.

n

counts selected datasets. VdP: Van der Pol; LV: Lotka–Volterra; parentheses identify the equation.

Study	Cells	$R_{\rm MC}$	Settings varied
Tuning	13	4000	Threshold $0.01,0.02,0.05,0.10,0.18$ ; $K=8,12,16,24$ ; radius $0.5,0.75,1$ ; allocations specified below. $M=800$ , $\sigma_{X}=0.06$ .
Noise scale	8	4000	$\sigma_{\rm work}/\sigma_{X}=0.5,0.75,0.9,1,1.1,1.25,1.5,2$ ; $M=800$ , $\sigma_{X}=0.06$ .
Misspecification	27	4000	Sine drift, exponential, logistic; $M=200,800,1600$ and $\sigma_{X}=0.02,0.06,0.12$ .
Misspecified signal	12	3000	Sine and exponential drifts, $c=0.25,0.5,0.75,1,1.5,2$ ; $M=800$ , $\sigma_{X}=0.06$ .
Components	6	5000	$M=200,800$ and $\sigma_{X}=0.06,0.12,0.24$ ; all arms retained.
Identification	12	3000	$M=200,800,3200$ and $\sigma_{X}=0.03,0.06,0.12,0.24$ .

Appendix

Table 12: Complete sensitivity grids.

R_{\rm MC}

is the number of independent datasets per cell,

\sigma_{\rm work}

is the supplied working noise standard deviation, and

c

multiplies the drift in the signal sweep. Except for tuning, the threshold is

0.02\sqrt{800/M}

. For

M=200,800,1600,3200

, bootstrap counts are 799,1199,1599,2399, respectively.

Quantity	Setting
Processed observations	Ten-point averages assigned their mean times; 10 Hz
Fitting window	$[20,180)$ s; $M=1600$ per device
Library and weak equations	$(1,u,u^{2},u^{3})$ ; $K=20$ prespecified weights
Thinning allocation	Four equal fractions, using the estimated working scale
Selection	STLSQ on view 1, mandatory constant, at most 10 iterations
Threshold candidates	$0,0.005,0.01,0.02,0.05,0.1,0.2,0.5$

Appendix

Table 13: Common protocol for the averaged-data coefficient analysis. All intervals have a working-model interpretation.

Feature	Estimate	Response	Gap	Total	Working interval
$1$	-0.520512	0.7403	0.0197	0.7406	$[-0.520657,-0.520367]$
$u$	0.289392	1.8457	0.0454	1.8463	$[0.289030,0.289754]$
$u^{3}$	-0.139550	1.5289	0.0339	1.5292	$[-0.139850,-0.139250]$

Appendix

Table 14: Device 1 coefficients in scaled coordinates. Scale columns give

10^{4}e_{j}

10^{4}\widehat{t}_{j}

and

10^{4}\widehat{s}_{j}

, suppressing the common selected-support subscript

\mathcal{S}

; estimates and intervals are not multiplied by

10^{4}

. Intervals have nominal level 95% under the estimated working noise model, without established coverage for these observations.

	Independent noise		Correlated noise
Device	Supplied scale	Estimated scale	Supplied scale	Estimated scale
1	94.87 (0.66)	94.80 (0.67)	79.90 (1.20)	61.10 (1.47)
2	95.33 (0.65)	95.33 (0.64)	83.40 (1.12)	67.43 (1.40)
3	95.67 (0.60)	95.67 (0.60)	82.37 (1.16)	67.87 (1.41)
4	94.57 (0.68)	94.57 (0.68)	80.63 (1.19)	62.57 (1.46)
5	95.85 (0.60)	95.83 (0.61)	84.30 (1.11)	74.20 (1.34)

Appendix

Table 15: Coverage of nominal 95% intervals in the discharge-template experiment. All four conditions use iid working covariance weighting. Supplied scale means the true marginal standard deviation, without supplying temporal correlation. Entries are mean percentages of selected coefficient intervals containing their own finite-grid BLA targets over 1,000 replications per device, with Monte Carlo standard errors in parentheses, in percentage points. The equal-device row uses eight fixed templates. Replication-level proportions account for dependence among coefficients; this is not a simultaneous coverage summary.

$M$	$\gamma$	$B_{\mathrm{boot}}$	$N_{A}$	Pipeline	Inference
200	0	799	21	0.0082 [0.0001, 0.0090]	0.0088 [0.0076, 0.0091]
800	0	1199	20	0.0382 [0.0001, 0.0396]	0.0391 [0.0382, 0.0398]
1600	0	1599	24	0.1010 [0.0987, 0.1038]	0.1013 [0.1000, 0.1045]
3200	0	2399	19	0.2673 [0.0002, 0.2708]	0.2701 [0.2680, 0.2709]
800	0.16	1199	30	0.0379 [0.0374, 0.0382]	0.0378 [0.0373, 0.0381]

Appendix

Table 16: Wall-clock runtime in seconds, reported as median [25th, 75th percentile]. Each setting uses 30 generated datasets;

N_{A}

counts those selecting the quadratic term. Pipeline summaries include all datasets, including those that stop after selection; inference summaries use only the

N_{A}

selected datasets. These early exits can produce small pipeline lower quartiles, and the different populations allow a pipeline median below the inference median.

Explore similar work

Apr 20, 2026cs.LG

AC-SINDy: Compositional Sparse Identification of Nonlinear Dynamics

We present AC-SINDy, a compositional extension of the Sparse Identification of Nonlinear Dynamics (SINDy) framework that replaces explicit feature libraries with a structured representation based on arithmetic circuits. Rather than enumerating candidate basis functions, the proposed approach constructs nonlinear features through compositions of linear functions and multiplicative interactions, yielding a compact and scalable parameterization and enabling sparsity to be enforced directly over the computational graph. We also introduce a formulation that separates state estimation from dynamics identification by combining latent state inference with shared dynamics and multi-step supervision, improving robustness to noise while preserving interpretability. Experiments on nonlinear and chaotic systems demonstrate that the method recovers accurate and interpretable governing equations while scaling more favorably than standard SINDy.

Peter Racioppo

Independent Researcher Los Angeles, CA 90731

Latent Dynamics Modeling Sparse Identification of Nonlinear Dynamics Nonlinear System Identification Symbolic Regression

Ask PDF

Jun 14, 2026cs.LG

Multi-Fidelity SINDy: Sparse Discovery of Nonlinear Dynamical Systems with Fidelity-Weighted Measurements

Data from simulations and experiments are rarely noise-free and often exhibit heterogeneous levels of fidelity. Measurement uncertainty may vary across repeated observations, sensing devices, or even within a single experiment. This work addresses the problem of discovering nonlinear dynamical systems from such inhomogeneous data. We extend the Sparse Identification of Nonlinear Dynamical Systems (SINDy) framework to account for variable noise levels by combining Ensemble SINDy and Weak SINDy within a weighted regression formulation derived from generalized least squares. A statistical justification for the weighting strategy is also provided. The methodology is validated on several benchmark systems, including ordinary and partial differential equations. In addition, we show the benefit of multi-fidelity integration for forecasting the dynamics of a double pendulum system. The results confirm that the proposed approach mitigates the adverse effects of heteroscedastic noise and that repeated, low-cost, low-quality measurements can improve model recovery, in some cases matching or outperforming reconstructions obtained using only high-fidelity data.

Filippo Zacchei, Ana Larrañaga, Attilio Frangi +2

MOX – Department of Mathematics, Politecnico di Milano, Piazza Leonardo da Vinci 32, 20133 Milano, Italy · Department of Mechanical Engineering, University of Washington, Seattle, Washington 98195, United States · Department of Civil and Environmental Engineering, Politecnico di Milano, Piazza Leonardo da Vinci 32, 20133 Milano, Italy

Sparse Identification of Nonlinear Dynamics Nonlinear System Identification Multifidelity Modeling Heteroscedastic Regression

Ask PDF

Jul 16, 2026cs.LG

An Introduction to Sparse Identification of Nonlinear Dynamics for Engineering Applications

Many engineering problems involve phenomena whose governing equations are poorly characterized or only partially known. Surrogate modeling techniques such as neural networks can capture the behavior of these systems, but they typically demand large training datasets that are difficult to obtain in engineering contexts and yield models with limited physical interpretability. The Sparse Identification of Nonlinear Dynamics (SINDy) method addresses both limitations by performing sparse regression over libraries of candidate nonlinear terms, recovering interpretable governing equations from comparatively small datasets. Although SINDy has been demonstrated extensively on canonical benchmark systems, its application to practical engineering problems is less widely documented. This tutorial introduces the SINDy method and progressively builds toward its main extensions, from noise-robust weak-form and ensembling-based variants to constrained and parametrizable formulations. The paper and the accompanying tutorial (available at https://github.com/paullililili/SINDy4Engineers) is organized in three parts: the first introduces the standard SINDy algorithm and progressively extends it, inviting readers without prior knowledge to follow each step and adapt the methods to their own problems; the remaining two parts present detailed case studies on (1) the system identification of an unmanned aerial vehicle and (2) a chaotic thermosyphon heat exchanger. Through these examples, we aim to demonstrate that SINDy is simple to implement yet flexible enough to serve as a valuable identification tool for advanced engineering applications.

Yao Cheng Li, Ana Larrañaga, Steven L. Brunton +1

Department of Aeronautics, Imperial College London, SW7 2AZ, United Kingdom · Department of Mechanical Engineering, University of Washington, Seattle, WA 98195, USA · NSF AI Institute in Dynamic Systems, University of Washington, Seattle, WA 98195, USA

Sparse Identification of Nonlinear Dynamics Nonlinear System Identification

Ask PDF