Freeze, Then Select: Structured Field Adapters and Stability-Validated Weak Selection for PDE Discovery from Sparse Observations

Unknown

Conceptio › Archive › arXiv CS

arXiv CSopen access

Freeze, Then Select: Structured Field Adapters and Stability-Validated Weak Selection for PDE Discovery from Sparse Observations

Unknown · 2026 · arxiv_cs

arXiv CS · Papers · License: Open Access · 2026

Open Source ↗Direct PDF ↓

machine learning, deep learning, neural networks

Freeze, Then Select: Structured Field Adapters and Stability-Validated Weak Selection for PDE Discovery from Sparse Observations Juncheng Zhong1 , Chenghuang Shen3,4 , Jianfeng Liu4 , Zhengdong Xiao4 , Longjiu Luo4 , Qianrong Wang4 , Wenjun Xu5,6,* , Wenlian Lu1,2,3,* 1

School of Mathematical Sciences, Fudan University, Shanghai 200433, China Center for Applied Mathematics, Fudan University, Shanghai 200438, China 3 Shanghai Center for Mathematical Sciences, Fudan University, Shanghai 200438, China 4 Alibaba Group, Hangzhou, China 5 Shanghai Institute of Technical Physics, CAS, Shanghai 200083, China 6 National Key Laboratory of Infrared Detection Technologies, Shanghai Institute of Technical Physics, CAS, Shanghai 200083, China * Corresponding authors: [email protected], [email protected] 2

arXiv:2607.29665v1 [cs.LG] 31 Jul 2026

Abstract PDE discovery from sparse observations requires reconstructing a continuous field and selecting the correct differential terms. Our analysis of optimization paths in coupled neural PDE discovery reveals three behaviors: the exact support can persist to the end of training, appear only transiently, or fail to emerge. To decouple equation selection from neural optimization, we develop a freeze-then-select method combining a structured field adapter with Stability-Validated Weak Selection (SVWS). Trained from observations without a PDE residual, the adapter factorizes the field into learned spatial features and temporal coefficients represented by cubic splines. After freezing the field, SVWS identifies recurrent terms across independent weak-form systems, refits candidate supports, and selects the final equation on held-out weak-form systems. Beyond fixed libraries, we apply the same principle to expressions generated by genetic programming and recover the power-law form of an unknown nonlinear diffusion function from sparse, noisy observations. Across all six sparse MDBench regimes, our method attains the highest exact support recovery rate, with its clearest gains over classical and neural baselines on challenging Kuramoto–Sivashinsky dynamics.

Introduction

Data-driven discovery of governing equations seeks parsimonious dynamical laws from measurements of physical systems (Bongard and Lipson 2007; Schmidt and Lipson 2009; Brunton, Proctor, and Kutz 2016). For PDEs, sparseregression methods such as PDE-FIND select active differential terms from prescribed libraries (Rudy et al. 2017; Schaeffer 2017), while weak formulations reduce sensitivity to pointwise differentiation (Schaeffer and McCalla 2017; Messenger and Bortz 2021). These approaches require field values, derivatives, or weak integrals on a suitable space– time representation. When observations cover only a subset of spatial sensors or time frames, those quantities are not directly available. PDE discovery from sparse observations therefore entails two linked tasks: reconstructing a continuous field and selecting the governing terms from the quantities induced by that reconstruction.

Neural surrogates address reconstruction by providing continuous, differentiable fields from limited observations (Both et al. 2021; Chen, Liu, and Sun 2021; Stephany and Earls 2024a,b). Many neural PDE-discovery methods couple the surrogate to equation coefficients or operator selection. This coupling can be effective: the evolving equation can regularize reconstruction, while the surrogate supplies derivatives or weak integrals away from sensors. The reconstructed field and inferred support, however, evolve together along a single optimization path. Our analysis of these paths reveals three behaviors in coupled neural PDE discovery (Fig. 3): the exact support can persist to the end of training, appear only transiently, or fail to emerge. These observations motivate our freeze-then-select design, which separates equation selection from neural surrogate training (Fig. 1). We fit a continuous field from observations without a PDE residual, freeze it, and evaluate candidate equations on the same reconstruction. Earlier two-stage methods demonstrate the value of this separation: DL-PDE applies sparse regression after fitting a neural field, while an integral-form method performs genetic search on quantities obtained from a pretrained surrogate (Xu, Chang, and Zhang 2021; Xu, Zhang, and Wang 2021). For sparse observations, however, two questions remain: how to reconstruct a field that yields informative differential or weak-form quantities, and how to reduce sensitivity to the choice of weak-form system and sparsity threshold. We address these questions with a structured field adapter and Stability-Validated Weak Selection (SVWS). The structured adapter factorizes the field into learned spatial features and temporal coefficients represented in a cubic B-spline basis (Fig. 2). This design follows the familiar view of coherent dynamics as spatial patterns with amplitudes that evolve over time (Sirovich 1987; Holmes, Lumley, and Berkooz 1996). The learned features capture complex spatial structure, while the spline representation provides smooth temporal evolution and analytic time derivatives (Sun, Liu, and Sun 2021; Sun et al. 2022). After freezing the field, SVWS generates candidate supports across independently constructed weak-form systems and sparsity thresholds, re-

Sparse observations

Structured field adapter

: fitted to observations without a PDE residual

ral baselines on KS (Table 1). Adapter ablations and selector sensitivity analyses examine the two stages of the method (Table 2). Beyond fixed libraries, the SVWS extension recovers the power-law form of an unknown nonlinear diffusion function from sparse, noisy observations (Table 3). Contributions. (i) We develop a freeze-then-select PDEdiscovery method whose structured field adapter couples learned spatial features to explicit temporal spline coefficients and is trained from observations without a PDE residual. (ii) We introduce SVWS, which separates support generation, coefficient refitting, and validation across independently constructed weak-form systems. (iii) We analyze optimization paths in coupled neural PDE discovery, showing that the exact support can persist to the end of training, appear only transiently, or fail to emerge. (iv) We extend SVWS beyond fixed libraries to expressions generated by GP, demonstrating recovery of the power-law form of an unknown nonlinear diffusion function from sparse, noisy observations.

Related Work SVWS PDE-FIND

Candidate generation

Stability filtering

Coefficient refit

WSINDy

Held-out weak-form validation

GP-generated expressions

Selected equation

Figure 1: Freeze, then select. A structured adapter fits sparse observations without a PDE residual and is frozen before equation discovery. SVWS selects both fixed-library supports and GP-generated expressions; the frozen field can also be passed to classical selectors. tains recurrent terms, refits each support, and selects the final equation on held-out weak-form systems. In this selection stage, integration by parts avoids high-order pointwise derivatives of the reconstruction (Messenger and Bortz 2021; Tang et al. 2023), while shifted test-function grids provide multiple views of the same frozen field for stability assessment and held-out validation. Beyond fixed libraries, an extension of SVWS refits and validates symbolic expressions generated by genetic programming (GP) on independent weak-form systems from the frozen field. We evaluate the method on sparse MDBench regimes spanning KdV, Kuramoto–Sivashinsky (KS), and twodimensional advection–diffusion, with observations restricted to subsets of fixed spatial sensors or time frames (Ziaei Bideh, Georgievska, and Gryak 2026). Our method attains the highest rate of exact support recovery in each of the six regimes, with its clearest gains over classical and neu-

Sparse regression and model selection. SINDy identifies parsimonious governing equations from predefined libraries (Brunton, Proctor, and Kutz 2016), and PDE-FIND extends sparse identification to PDEs (Rudy et al. 2017). Later work develops implicit and relaxed formulations, ensemble methods, information criteria, and approaches to noisy or incomplete data (Mangan et al. 2017; Zheng et al. 2019; Kaheman, Kutz, and Brunton 2020; Reinbold, Gurevich, and Grigoriev 2020; Fasel et al. 2022; Xu, Zeng, and Zhang 2023). Stability selection retains terms that recur across resampled fits, while the one-standard-error rule selects the simplest model within one standard error of the minimum validation error (Meinshausen and Bühlmann 2010; Hastie, Tibshirani, and Friedman 2009). SVWS adapts these ideas to supports generated across independently constructed weak-form systems. Weak and integral formulations. Weak-form PDE discovery transfers derivatives from the field to smooth test functions, reducing reliance on pointwise differentiation. Variational system identification, integral sparse selection, WSINDy, and WeakIdent use this principle to improve robustness to noise (Wang, Huan, and Garikipati 2019; Schaeffer and McCalla 2017; Messenger and Bortz 2021; Tang et al. 2023). We apply weak-form selection after freezing the reconstructed field. Neural surrogates for PDE discovery. Physics-informed neural networks (PINNs) fit differentiable surrogates with data and physics objectives (Raissi, Perdikaris, and Karniadakis 2019; Karniadakis et al. 2021); balancing these objectives remains an important optimization challenge (Wang, Teng, and Perdikaris 2021; Krishnapriyan et al. 2021). Coupled neural discovery methods, including PDE-Net, DeepMoD, PINN-SR, PDE-LEARN, and Weak-PDE-LEARN, learn trainable representations together with differential operators or coefficients (Long et al. 2018; Long, Lu, and Dong 2019; Both et al. 2021; Chen, Liu, and Sun 2021; Stephany and Earls 2024a,b). Two-stage alternatives separate reconstruction from discovery: DL-PDE applies sparse regression

Spatial branch

Trainable Background

Fixed Fourier features

Spatial MLP

Spatial features

Temporal branch Centered cubic B-spline basis

Observation-only training objective (no PDE residual)

Observation fit

Feature diversity

Temporal smoothness

Figure 2: Structured field adapter. A spatial network produces a background and spatial features with temporal coefficients represented in a centered cubic B-spline basis. The adapter is fit without a PDE residual and frozen before equation selection. to a fitted neural field, while an integral-form method performs genetic search on quantities obtained from a pretrained surrogate (Xu, Chang, and Zhang 2021; Xu, Zhang, and Wang 2021). We build on this separation with a structured adapter for reconstruction and SVWS after freezing. Structured field representations. Reduced-order models describe coherent dynamics through spatial modes and timedependent coefficients (Sirovich 1987; Holmes, Lumley, and Berkooz 1996). Coordinate networks provide continuous fields, with Fourier features and periodic activations improving high-frequency representation (Tancik et al. 2020; Sitzmann et al. 2020). Spline-based methods combine smooth reconstruction with sparse equation discovery (Sun, Liu, and Sun 2021; Sun et al. 2022). Our adapter follows this structure, learning spatial features with temporal coefficients represented in a spline basis. Symbolic PDE discovery. Symbolic regression searches compositional expressions beyond fixed candidate libraries (Koza 1992; Bongard and Lipson 2007; Schmidt and Lipson 2009; Cranmer 2023). PDE discovery beyond fixed libraries has used surrogate-assisted and symbolic genetic algorithms (Xu, Chang, and Zhang 2020; Chen et al. 2022), physicsinformed genetic programming (Cohen, Beykal, and Bollas 2024), and closed-form search (Kacprzyk, Qian, and van der Schaar 2023). Recent approaches also use reinforcement learning, abductive learning, differentiable weak-form networks, and generative models (Du, Chen, and Zhang 2024; Gao et al. 2025; Li et al. 2026; Xu et al. 2025). In our setting, genetic programming proposes expressions for an unknown nonlinear diffusion function, and the SVWS extension selects among them using weak-form systems from the frozen field.

Freeze-Then-Select PDE Discovery Our freeze-then-select method fits a continuous field to sparse observations without a PDE residual and then freezes it, so every candidate equation is evaluated on the same reconstruction. The structured field adapter provides this reconstruction, and Stability-Validated Weak Selection (SVWS) assigns support generation, coefficient refitting, and validation to independent weak-form systems. The frozen field can also be passed to classical PDE selectors. Beyond fixed libraries, the same refitting and validation principle applies to expressions generated by genetic programming. Problem setting. Let u : Ω × [0, T ] → R be an unknown scalar field on Ω ⊂ Rd , observed through sparse, possibly noisy samples Dobs = {(xi , ti , yi )}ni=1 ,

yi = u(xi , ti ) + ϵi ,

(1)

where ϵi denotes measurement noise. The sampling pattern may retain only a subset of spatial locations or time slices. In the fixed-library setting, the governing equation is assumed to be sparse in a prescribed library of J candidate terms: ∂t u =

J X

ξj⋆ Θj (u, ∇u, ∇2 u, . . .),

j=1

⋆

(2)

= {j : ξj⋆ ̸= 0}.

Given Dobs , the goal is to recover S ⋆ and its coefficients; support recovery is exact when Sb = S ⋆ . Sparse samples do not directly provide the continuous field, derivatives, or weak integrals required for equation selection. We therefore

reconstruct a continuous field adapter û from which these quantities can be evaluated. Structured field adapter. The adapter represents temporal evolution through smooth coefficients multiplying learned spatial features. After affine rescaling of x, a spatial network receives the fixed Fourier encoding γ(e x) = [sin(Be x), cos(Be x)] and outputs a background bη (x) and features Φη (x) ∈ RR . The rows of B are sampled once from a standard Gaussian distribution and then fixed. Let β(t) ∈ RM be a clamped cubic B-spline basis. We center it over the observed times, X 1 β̄(t) = β(t) − β(τ ), (3) |Tobs | τ ∈Tobs

so the dynamic component has zero empirical mean over the observed times and bη (x) represents the corresponding mean of the reconstructed field. With C ∈ RR×M , the reconstructed field is ûη,C (x, t) = bη (x) + Φη (x)⊤ C β̄(t).

(4)

Here cr (t) = Cr: β̄(t) is the coefficient of spatial feature r, while R and M specify the numbers of spatial features and temporal basis functions. We fit (η, C) from observations by (5)

The spatial regularizer RΦ penalizes deviations of the empirical feature Gram matrix from the identity, while Rt penalizes the integrated squared curvature of the temporal coefficients. The objective contains no PDE residual or PDE coefficients. Regularizer definitions and optimization settings are provided in Supplementary Sec. A.3. After fitting, we freeze ûη,C and construct every weak system from this shared reconstruction. Weak systems from the frozen field. Following weak-form PDE discovery (Messenger and Bortz 2021; Tang et al. 2023), we construct local weak systems using compactly supported space–time test functions ψℓ on interior patches Pℓ . Patch centers follow a regular grid whose phase is sampled independently for each weak system, with one row per center. The polynomial kernel, patch geometry, benchmark-specific settings, and formal error bounds are given in Supplementary Sec. B.1. R Define ⟨f, g⟩ℓ = Pℓ f g dx dt, and write each candidate as Θj (u) = Dαj qj (u), with the identity operator covering algebraic terms. Integration by parts gives bℓ (û) = −⟨û, ∂t ψℓ ⟩ℓ ,

Aℓj (û) = ⟨qj (û), Dα∗ j ψℓ ⟩ℓ . (6)

Here Dα∗ j denotes the weak adjoint. Because the test function and its required derivatives vanish at the patch boundary, integration by parts eliminates boundary terms and transfers derivatives from û to the test function. Stacking m patches yields A(û)ξ ≈ b(û),

A(û) ∈ Rm×J ,

b(û) ∈ Rm .

g {(Argen , brgen )}r=1 ,

(Afit , bfit ),

v {(Asval , bsval )}R s=1 . (8)

Each system uses a regular patch grid with an independently sampled phase. Generation systems propose supports, the fit system estimates their coefficients, and validation systems score the fitted equations; no system is reused across roles. After column normalization, sequential thresholded least squares (STLSQ) applied to generation system r at threshold λ ∈ Λ produces support Sr,λ (Brunton, Proctor, and Kutz 2016). Deduplicating these supports gives the candidate set C. Following stability selection and ensemble sparse discovery (Meinshausen and Bühlmann 2010; Fasel et al. 2022), we measure term recurrence by Rg

πj =

1 XX 1{j ∈ Sr,λ }, Rg |Λ| r=1 λ∈Λ

1X min |ûη,C (xi , ti ) − yi |2 + λΦ RΦ + λt Rt . η,C n i=1

If qj is locally Lipschitz on the reconstructed state range, the error in each weak entry is bounded by the local L2 reconstruction error and the norm of the corresponding testfunction derivative. Changing the grid phase changes the sampled local patches while keeping the reconstruction fixed, producing separate weak systems for support generation, coefficient refitting, and validation. Stability-Validated Weak Selection. SVWS assigns weak systems constructed from the same frozen field to three separate roles:

(7)

Jτ = {j : πj ≥ τ },

(9)

Cτ = {S ∈ C : S ⊆ Jτ }.

Thus πj is the selection frequency of term j, and Cτ retains supports whose active terms all recur with frequency at least τ . Using a threshold path exposes supports across multiple sparsity levels and reduces dependence on any single STLSQ threshold. The recurrence filter then excludes supports containing terms that do not recur across phases and thresholds. Let Dfit contain the fit-system column norms and define efit = Afit D−1 . Each eligible support is then refit on the A fit independent fit system by ridge regression: ˆe ξ S = arg

min e e ξ:supp( ξ)⊆S

efit ξe − bfit ∥22 + α∥ξeS ∥22 , ∥A

(10)

−1 ˆ ξˆS = Dfit ξeS .

The second line restores the coefficients to the original library scale. We then compute the normalized residual on each validation system, with ε = 10−12 for numerical stability: rs (S) =

∥Asval ξˆS − bsval ∥22 , ∥bsval ∥22 + ε R

v 1 X µ(S) = rs (S). Rv s=1

(11)

Let se(S) denote the standard error across these validation risks, and let Smin minimize their mean µ(S). The onestandard-error rule (Hastie, Tibshirani, and Friedman 2009) admits Cadm = {S ∈ Cτ : µ(S) ≤ µ(Smin ) + se(Smin )}, (12) Sb = arg min |S|, µ(S), −v(S) . S∈Cadm

The tuple is ordered lexicographically, and v(S) is the number of times S is generated. SVWS therefore selects the smallest support whose independently refitted equation lies within one standard error of the minimum validation risk; lower risk and higher generation count break ties. The reported coefficients are ξˆSb. Symbolic selection beyond fixed libraries. The proposal stage can generate symbolic expressions instead of fixedlibrary supports while retaining the frozen reconstruction, independent refitting, and validation on held-out weak systems. We consider K X ut ≈ ar Dαr qr (u), qr ∈ G, Dαr ∈ A, (13) r=1

where G is the expression class generated by a symbolic grammar and the operators in A have known weak adjoints. A proposed pair (qr , Dαr ) contributes the weak column ⟨qr (û), Dα∗ r ψℓ ⟩ℓ . Integration by parts therefore transfers the outer differential operator to the test function and avoids pointwise differentiation of qr (û). Independent GP searches on weak systems built from shifted patch grids propose a finite set of expressions (Koza 1992; Cranmer 2023). Because GP can produce many trees that share the same canonical structure but differ in their numerical parameters, each search groups them into expression families and retains the representative with the lowest proposal risk from each family. We merge these representatives and freeze the candidate pool before fitting their continuous parameters on an independent weak system. The fitted expressions are then compared on held-out weak systems using the mean validation risk and a fixed complexity penalty. Our experiment uses K = 1 and A = {∂xx }, giving ut = ∂xx q(u) after absorbing a1 into q. The outer operator is therefore fixed, while the nonlinear diffusion function q is discovered symbolically.

Experiments

Benchmarks and observation protocols. We evaluate on public MDBench trajectories for three PDEs (Ziaei Bideh, Georgievska, and Gryak 2026): KdV: ut = −6uux − uxxx , KS: ut = −uux − uxx − uxxxx , (14) 2D AD: ut = 0.25ux + 0.5uy + 0.5uxx + 0.5uyy . We use the released grids, initial conditions, and trajectories without resimulation. The main fixed-library comparison uses clean observations under two 20% sampling protocols. In S20, a random 20% subset of spatial locations is fixed and observed at every time point. In T20, the initial frame and randomly selected later frames form an irregular 20% temporal subset observed on the full spatial grid. For each PDE and protocol, we evaluate five test seeds and provide every method with the same observations for a given seed. Candidate libraries are shared across methods: KdV and KS use the eight-term library {1, u, u2 , u3 , uux , uxx , uxxx , uxxxx }, while 2D AD uses {1, u, u2 , ux , uy , uxx , uyy }. Baselines and implementation details. Ours combines the structured field adapter with SVWS. We compare with

PDE-FIND, WSINDy, DeepMoD, PINN-SR, Weak-PDELEARN, and DL-PDE (Rudy et al. 2017; Messenger and Bortz 2021; Both et al. 2021; Chen, Liu, and Sun 2021; Stephany and Earls 2024b; Xu, Chang, and Zhang 2021). PDE-FIND and WSINDy operate on a shared spectral reconstruction fitted by ridge regression to the sampled observations. Hyperparameters not specified by the original sources are selected on disjoint development seeds and fixed across all test settings. Implementation sources and reporting rules are documented in Supplementary Sec. C. In the fixed-library experiments, SVWS uses 12 generation systems, one fit system, and seven validation systems. All adapter and selector settings are fixed before testing; their complete configurations are provided in Supplementary Secs. A.3 and B.2. Metrics. Exact recovery requires Sb = S ⋆ and is reported as the number of successful runs out of five. We also report relative coefficient and field errors: ∥ξˆ − ξ ⋆ ∥2 ∥û − u⋆ ∥2,X ×T Eξ = ⋆ , Eu = , ∥ξ ∥2 + ε ∥u⋆ ∥2,X ×T where ε = 10−12 . For Eξ , each discovered equation is mapped to the shared physical library, with zero coefficients assigned to absent terms. Field error is evaluated on the full reference grid, and error summaries are medians over the five test runs. Fixed-library recovery. Table 1 reports the primary comparison. Our method recovers the exact support in all five runs across all six regimes. On KS, the best baseline reaches 3/5 exact recoveries in each observation protocol, whereas our method recovers all five and also attains the lowest median coefficient error. Supplementary Sec. G shows that our KS recovery remains exact across the evaluated sensor densities and that exact support is retained under 10% observation noise on KS and 2D AD. Adapter ablations and threshold robustness. Panel (a) of Table 2 focuses on KS, where adapter choice produces clear differences in field error and support recovery; Supplementary Sec. D.1 reports all six regimes. With matched parameter counts, the structured adapter is the only variant that recovers all five supports under both observation protocols and gives the lowest field error by a clear margin. Fourier features reduce reconstruction error for both the factorized adapter and the coordinate MLP, while the factorized representation remains more reliable for equation recovery. Panel (b) asks whether one STLSQ threshold transfers across equations and observation protocols, with λ0 denoting the benchmark-specific base threshold. Lower thresholds admit additional terms on KS, whereas the largest threshold removes governing terms in both 2D AD regimes. SVWS instead generates candidates over the full prespecified path and multiple weak systems, maintaining exact recovery across all four settings. Additional selector diagnostics are reported in Supplementary Sec. D.2. Support dynamics under joint optimization. Figure 3 traces support recovery along three coupled neural optimization paths. On 2D AD-S20, DeepMoD reaches the exact support and retains it to termination, whereas WeakPDE-LEARN visits it only transiently during training. On

KdV Method Ours PDE-FIND† WSINDy† DeepMoD PINN-SR‡ Weak-PDE-LEARN DL-PDE‡

2D AD

S20

T20

S20

T20

S20

T20

5/5 (0.036) 3/5 (0.044) 5/5 (0.002) 5/5 (0.001) 5/5 (0.015) 0/5 (0.326) 5/5 (0.049)

5/5 (0.008) 4/5 (0.019) 4/5 (0.011) 5/5 (0.001) 5/5 (0.018) 0/5 (0.323) 5/5 (0.094)

5/5 (0.035) 0/5 (0.136) 1/5 (0.782) 0/5 (1.010) 0/5 (0.998) 0/5 (1.238) 3/5 (0.927)

5/5 (0.055) 3/5 (0.512) 1/5 (0.818) 0/5 (0.993) 0/5 (0.997) 0/5 (1.752) 2/5 (0.905)

5/5 (0.010) 0/5 (0.894) 0/5 (0.248) 5/5 (0.008) 0/5 (0.326) 0/5 (0.675) 5/5 (0.031)

5/5 (0.009) 0/5 (0.885) 5/5 (0.099) 5/5 (0.012) 0/5 (0.308) 0/5 (0.682) 5/5 (0.018)

Applied after fitting a shared spectral reconstruction to the sampled observations by ridge regression. ‡ Reimplemented from the published method.

†

Table 1: Fixed-library discovery from sparse observations. Cells report exact-support recovery (/5), with median relative coefficient error Eξ in parentheses. S20 observes 20% fixed sensors; T20 observes 20% time frames. Bold marks the highest recovery; underlining marks the lowest Eξ among ties. (a) Adapter ablations on KS

KS-S20

KS-T20

Adapter

Exact

Structured adapter Structured adapter without Fourier features Coordinate MLP with Fourier features Coordinate MLP without Fourier features

5/5 3/5 3/5 0/5

0.082 0.659 0.440 0.905

5/5 5/5 1/5 0/5

0.082 0.650 0.437 0.867

(b) Robustness to the STLSQ threshold Selection

KS-S20

KS-T20

2D AD-S20

2D AD-T20

5/5 1/5 5/5 5/5

5/5 3/5 4/5 5/5

5/5 5/5 5/5 0/5

SVWS STLSQ (0.25λ0 ) STLSQ (λ0 ) STLSQ (2λ0 )

Table 2: Adapter ablations and STLSQ threshold robustness. Panel (a) compares adapters with matched parameter counts on KS; panel (b) compares SVWS with single-system STLSQ along the prespecified threshold path on the same frozen fields. Eu is median relative field error. (a) Persists to termination

(b) Appears only transiently

DeepMoD · 2D AD 20% fixed spatial sensors (S20)

Weak-PDE-LEARN · 2D AD 20% fixed spatial sensors (S20)

DeepMoD · KS full observation grid (S100)

Support F1

⋆

exact support

⋆ 1

First 12%

.67

⋆

• Terminal: exact support 0

100

0 ⋆ First exact: 1.8% • Terminal: {uy , uyy }

50 Training progress (%)

100

• Terminal: {u, uux , uxxxx } Missing: uxx 0

100

Figure 3: Representative support trajectories under joint neural optimization. DeepMoD retains the exact 2D AD support on S20, Weak-PDE-LEARN visits it only transiently, and DeepMoD never visits the exact KS support on S100. Stars mark first exact visits and circles terminal states. Each reported behavior occurs in all five runs for its corresponding setting; complete paths are reported in Supplementary Sec. E.

Noise (%)

Power-law recovery (/5)

Parameters on recovered forms κ̂

SVWS (ours) 20% fixed sensors 0 5/5 0.0978 ± 0.0019 5 5/5 0.0976 ± 0.0026 10 5/5 0.0982 ± 0.0041 PySR full-grid q(u) estimate 0 1/5 0.1000 5 3/5 0.0995 ± 0.0005 10 2/5 0.0988 ± 0.0014

m̂ 1.758 ± 0.027 1.682 ± 0.116 1.552 ± 0.092 1.730 1.761 ± 0.043 1.740 ± 0.039

Table 3: Symbolic recovery of q ⋆ (u) = 0.1u1.73 (five runs per noise level). SVWS uses 20% fixed sensors; PySR is applied to a full-grid weak-form estimate of q(u). Parameters are mean ± sample SD over successful recoveries. KS-S100, DeepMoD never reaches the exact support at a recorded checkpoint despite observing the full spatial grid. These behaviors have different implications for checkpoint reporting: a transient visit makes the reported equation checkpoint-dependent, whereas an unvisited support cannot be recovered by changing the checkpoint. Freezing the field decouples equation selection from this path: candidates are refit on the same reconstruction and compared on held-out weak systems rather than selected from a training checkpoint. Symbolic discovery beyond a fixed library. We next examine whether the same freeze-then-select design supports symbolic discovery beyond a fixed library. The structured adapter is trained using observations from 20% of the spatial sensors, after which genetic programming (GP) proposes candidate expressions and SVWS selects among them using independent weak systems. We consider the nonlinear diffusion equation ut = ∂xx q(u), q ⋆ (u) = 0.1u1.73 , under 0%, 5%, or 10% observation noise. Here q is the unknown nonlinear diffusion function. Since ∂xx eliminates additive constants, q and q + C define the same PDE. The target exponent 1.73 is absent from the initialization grid and must be reached through mutation and continuous refitting. We run two GP proposal searches on independently shifted weak systems. Each search groups candidates by canonical structure, and the retained representatives are merged into a frozen pool. A separate weak system fits their continuous constants and exponents, while three held-out systems evaluate the fitted expressions. Final selection combines the mean validation risk with a fixed complexity penalty. A run counts as power-law recovery when the selected expression canonicalizes to q̂(u) = κ̂um̂ + C. Table 3 reports recovery frequency together with the corresponding parameter estimates. Supplementary Sec. F documents the complete protocol. PySR (Cranmer 2023) does not reconstruct a field from sparse PDE observations, so we apply it to a weak-form estimate of q(u) from the full noisy grid. SVWS consistently recovers the target power-law structure from sparse sensor observations across all tested noise levels.

PySR yields accurate parameter estimates for its recovered power laws, but selects the target structure less consistently from the full-grid q(u) estimate. The comparison distinguishes accurate parameter fitting from reliable structural recovery and shows how the frozen field interface connects sparse observations to symbolic search.

Discussion The fixed-library experiments show complementary contributions from field reconstruction and equation selection. On KS, the structured adapter improves both field accuracy and support recovery. Across KS and two-dimensional advection–diffusion, single-system STLSQ is sensitive to the sparsity threshold, while SVWS remains consistent on the same frozen fields (Table 2). In the external comparison, KdV and two-dimensional advection–diffusion are recovered by several baselines, whereas KS produces the clearest separation in support recovery and coefficient accuracy (Table 1). Together, these results show the value of combining structured reconstruction with selection across multiple weak systems and sparsity thresholds. Figure 3 provides further motivation for freezing. In coupled neural discovery, the reported equation depends on whether optimization reaches the correct support and whether that support survives to the selected checkpoint. A transient visit makes the result checkpoint-dependent, while a support absent from the recorded path cannot be recovered through checkpoint selection. The KS-S100 case further shows that complete spatial coverage does not ensure that the recorded path contains the correct support. The freeze-thenselect design avoids this checkpoint dependence by refitting candidate equations on the same reconstruction and comparing them on held-out weak systems. The symbolic experiment extends the same design from fixed-library supports to expressions proposed by GP. After the structured adapter converts sparse sensor observations into a frozen field, continuous parameters are fitted on one weak system and the candidate expressions are compared on held-out systems. SVWS consistently recovers the target power-law form, whereas PySR selects it less reliably from the full-grid q(u) estimate (Table 3). By fixing the outer operator ∂xx , the experiment isolates discovery of the nonlinear function and directly tests whether validation on the frozen field can extend beyond fixed-library support selection. The result supports this extension and motivates richer symbolic grammars that also search over multiple functional terms and differential operators.

Conclusion We presented a freeze-then-select method for PDE discovery from sparse observations. A structured field adapter reconstructs a continuous field without a PDE residual, and SVWS uses independent weak systems for support generation, coefficient refitting, and validation. The method recovers exact support across all six sparse MDBench regimes, with its clearest gains on KS; its symbolic extension recovers the power-law form of an unknown nonlinear diffusion function from sparse, noisy observations. Optimization-path analy-

sis further shows that coupled neural discovery may retain, lose, or never reach the exact support, motivating equation selection after reconstruction. By freezing the field before selection, the framework connects sparse measurements to both fixed-library and symbolic PDE discovery through a common validation principle.

Acknowledgments

This work was supported by the Science Fund for Excellent Research Groups of the National Natural Science Foundation of China (Grant No. 62588101) and by the Science and Education Integration Project of the Shanghai Institute of Technical Physics, CAS (Project No. SITPKJRH-2025-03). Code availability. Code and configurations will be made publicly available upon acceptance.

References

Bongard, J.; and Lipson, H. 2007. Automated reverse engineering of nonlinear dynamical systems. Proceedings of the National Academy of Sciences, 104(24): 9943–9948. Both, G.-J.; Choudhury, S.; Sens, P.; and Kusters, R. 2021. DeepMoD: Deep learning for model discovery in noisy data. Journal of Computational Physics, 428: 109985. Brunton, S. L.; Proctor, J. L.; and Kutz, J. N. 2016. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proceedings of the National Academy of Sciences, 113(15): 3932–3937. Chen, Y.; Luo, Y.; Liu, Q.; Xu, H.; and Zhang, D. 2022. Symbolic Genetic Algorithm for Discovering Open-Form Partial Differential Equations (SGA-PDE). Physical Review Research, 4(2): 023174. Chen, Z.; Liu, Y.; and Sun, H. 2021. Physics-informed learning of governing equations from scarce data. Nature Communications, 12: 6136. Cohen, B. G.; Beykal, B.; and Bollas, G. M. 2024. PhysicsInformed Genetic Programming for Discovery of Partial Differential Equations from Scarce and Noisy Data. Journal of Computational Physics, 514: 113261. Cranmer, M. 2023. Interpretable Machine Learning for Science with PySR and SymbolicRegression.jl. arXiv preprint arXiv:2305.01582. Du, M.; Chen, Y.; and Zhang, D. 2024. DISCOVER: Deep Identification of Symbolically Concise Open-Form Partial Differential Equations via Enhanced Reinforcement Learning. Physical Review Research, 6(1): 013182. Fasel, U.; Kutz, J. N.; Brunton, B. W.; and Brunton, S. L. 2022. Ensemble-SINDy: Robust sparse model discovery in the low-data, high-noise limit, with active learning and control. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 478(2260): 20210904. Gao, E.-H.; Ge, C.; Jiang, Y.; and Zhou, Z.-H. 2025. Discovering Symbolic Partial Differential Equation by Abductive Learning. In Advances in Neural Information Processing Systems, volume 38. Hastie, T.; Tibshirani, R.; and Friedman, J. 2009. The Elements of Statistical Learning: Data Mining, Inference, and Prediction. Springer, second edition.

Holmes, P.; Lumley, J. L.; and Berkooz, G. 1996. Turbulence, Coherent Structures, Dynamical Systems and Symmetry. Cambridge University Press. Kacprzyk, K.; Qian, Z.; and van der Schaar, M. 2023. D-CIPHER: Discovery of Closed-Form Partial Differential Equations. In Advances in Neural Information Processing Systems, volume 36. Kaheman, K.; Kutz, J. N.; and Brunton, S. L. 2020. SINDyPI: A Robust Algorithm for Parallel Implicit Sparse Identification of Nonlinear Dynamics. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 476(2242): 20200279. Karniadakis, G. E.; Kevrekidis, I. G.; Lu, L.; Perdikaris, P.; Wang, S.; and Yang, L. 2021. Physics-informed machine learning. Nature Reviews Physics, 3(6): 422–440. Koza, J. R. 1992. Genetic Programming: On the Programming of Computers by Means of Natural Selection. MIT Press. Krishnapriyan, A.; Gholami, A.; Zhe, S.; Kirby, R.; and Mahoney, M. W. 2021. Characterizing Possible Failure Modes in Physics-Informed Neural Networks. In Advances in Neural Information Processing Systems, volume 34, 26548–26560. Li, X.; Cui, X.; Qi, J.; Zhang, J.; Li, D.; and Yin, J. 2026. Weak-PDE-Net: Discovering Open-Form PDEs via Differentiable Symbolic Networks and Weak Formulation. arXiv preprint arXiv:2603.22951. Long, Z.; Lu, Y.; and Dong, B. 2019. PDE-Net 2.0: Learning PDEs from data with a numeric-symbolic hybrid deep network. Journal of Computational Physics, 399: 108925. Long, Z.; Lu, Y.; Ma, X.; and Dong, B. 2018. PDE-Net: Learning PDEs from Data. In Proceedings of the 35th International Conference on Machine Learning, volume 80 of Proceedings of Machine Learning Research, 3208–3216. Mangan, N. M.; Brunton, S. L.; Proctor, J. L.; and Kutz, J. N. 2017. Model selection for dynamical systems via sparse regression and information criteria. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 473(2204): 20170009. Meinshausen, N.; and Bühlmann, P. 2010. Stability Selection. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 72(4): 417–473. Messenger, D. A.; and Bortz, D. M. 2021. Weak SINDy for partial differential equations. Journal of Computational Physics, 443: 110525. Raissi, M.; Perdikaris, P.; and Karniadakis, G. E. 2019. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics, 378: 686–707. Reinbold, P. A. K.; Gurevich, D. R.; and Grigoriev, R. O. 2020. Using noisy or incomplete data to discover models of pattern-forming dynamics. Physical Review E, 101(1): 010203. Rudy, S. H.; Brunton, S. L.; Proctor, J. L.; and Kutz, J. N. 2017. Data-driven discovery of partial differential equations. Science Advances, 3(4): e1602614.

Schaeffer, H. 2017. Learning partial differential equations via data discovery and sparse optimization. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 473(2197): 20160446. Schaeffer, H.; and McCalla, S. G. 2017. Sparse Model Selection via Integral Terms. Physical Review E, 96(2): 023302. Schmidt, M.; and Lipson, H. 2009. Distilling free-form natural laws from experimental data. Science, 324(5923): 81–85. Sirovich, L. 1987. Turbulence and the dynamics of coherent structures. I. Coherent structures. Quarterly of Applied Mathematics, 45(3): 561–571. Sitzmann, V.; Martel, J. N. P.; Bergman, A. W.; Lindell, D. B.; and Wetzstein, G. 2020. Implicit Neural Representations with Periodic Activation Functions. In Advances in Neural Information Processing Systems, volume 33, 7462–7473. Stephany, R.; and Earls, C. J. 2024a. PDE-LEARN: Using deep learning to discover partial differential equations from noisy, limited data. Neural Networks, 174: 106242. Stephany, R.; and Earls, C. J. 2024b. Weak-PDE-LEARN: A weak form based approach to discovering PDEs from noisy, limited data. Journal of Computational Physics, 506: 112950. Sun, F.; Liu, Y.; and Sun, H. 2021. Physics-Informed Spline Learning for Nonlinear Dynamics Discovery. In Proceedings of the Thirtieth International Joint Conference on Artificial Intelligence, 2054–2061. Sun, L.; Huang, D. Z.; Sun, H.; and Wang, J.-X. 2022. Bayesian Spline Learning for Equation Discovery of Nonlinear Dynamics with Quantified Uncertainty. In Advances in Neural Information Processing Systems, volume 35, 6927– 6940. Tancik, M.; Srinivasan, P. P.; Mildenhall, B.; Fridovich-Keil, S.; Raghavan, N.; Singhal, U.; Ramamoorthi, R.; Barron, J. T.; and Ng, R. 2020. Fourier Features Let Networks Learn High Frequency Functions in Low Dimensional Domains. In Advances in Neural Information Processing Systems, volume 33, 7537–7547. Tang, M.; Liao, W.; Kuske, R.; and Kang, S. H. 2023. WeakIdent: Weak Formulation for Identifying Differential Equations Using Narrow-Fit and Trimming. Journal of Computational Physics, 483: 112069. Wang, S.; Teng, Y.; and Perdikaris, P. 2021. Understanding and Mitigating Gradient Flow Pathologies in PhysicsInformed Neural Networks. SIAM Journal on Scientific Computing, 43(5): A3055–A3081. Wang, Z.; Huan, X.; and Garikipati, K. 2019. Variational System Identification of the Partial Differential Equations Governing the Physics of Pattern-Formation: Inference under Varying Fidelity and Noise. Computer Methods in Applied Mechanics and Engineering, 356: 44–74. Xu, H.; Chang, H.; and Zhang, D. 2020. DLGA-PDE: Discovery of PDEs with Incomplete Candidate Library via Combination of Deep Learning and Genetic Algorithm. Journal of Computational Physics, 418: 109584.

Xu, H.; Chang, H.; and Zhang, D. 2021. DL-PDE: DeepLearning Based Data-Driven Discovery of Partial Differential Equations from Discrete and Noisy Data. Communications in Computational Physics, 29(3): 698–728. Xu, H.; Chen, Y.; Cao, R.; Tang, T.; Du, M.; Li, J.; Callaghan, A. H.; and Zhang, D. 2025. Generative Discovery of Partial Differential Equations by Learning from Math Handbooks. Nature Communications, 16: 10255. Xu, H.; Zeng, J.; and Zhang, D. 2023. Discovery of partial differential equations from highly noisy and sparse data with physics-informed information criterion. Research, 6: 0147. Xu, H.; Zhang, D.; and Wang, N. 2021. Deep-Learning Based Discovery of Partial Differential Equations in Integral Form from Sparse and Noisy Data. Journal of Computational Physics, 445: 110592. Zheng, P.; Askham, T.; Brunton, S. L.; Kutz, J. N.; and Aravkin, A. Y. 2019. A Unified Framework for Sparse Relaxed Regularized Regression: SR3. IEEE Access, 7: 1404–1423. Ziaei Bideh, A.; Georgievska, A.; and Gryak, J. 2026. MDBench: Benchmarking Data-Driven Methods for Model Discovery. Proceedings of the AAAI Conference on Artificial Intelligence, 40(24): 19746–19754.

Supplementary Material

Freeze, Then Select: Structured Field Adapters and Stability-Validated Weak Selection for PDE Discovery from Sparse Observations This supplement provides implementation details and additional results for the experiments in the main paper. Sections A–C specify the experimental protocols, method configuration, and baseline implementations; Sections D–G present additional ablations, analyses of recorded optimization paths, symbolic discovery details, and analyses of observation density and noise. Section H summarizes the planned code release and computational resources.

Experimental Details

The fixed-library experiments use all five test seeds listed in Sec. A.1. Dense reference fields and ground-truth equations are used only for evaluation.

A.1

Benchmarks and Sparse Observation Protocols

We use the public MDBench trajectories for KdV, Kuramoto– Sivashinsky (KS), and two-dimensional advection–diffusion (2D AD). For each benchmark, we use the released coordinates, initial condition, and the first released trajectory (state index 0) without resimulation. Table S1 summarizes the complete grids and the resulting observation counts under the two sparse observation protocols. The primary fixed-library comparison uses clean observations under two complementary 20% sampling protocols. In S20, each seed samples spatial locations uniformly without replacement and observes the same sensor set at every time frame. In T20, the initial frame is retained and the remaining frames are sampled uniformly without replacement to form an irregular temporal subset observed on the complete spatial grid. We evaluate the five test seeds {1301, 1709, 2203, 2917, 3571} for every benchmark and protocol. For each benchmark, protocol, and seed, a shared observation record stores the sampled coordinates and values, the spatial masks, and the retained time indices. Every method receives the same record. Methods that fit continuous surrogates may evaluate those surrogates at additional coordinates for automatic differentiation or weak integration, but only the sampled values enter model fitting. Dense reference values are used exclusively to compute evaluation metrics. The shared physical candidate libraries are ΘKdV/KS = [1, u, u2 , u3 , uux , uxx , uxxx , uxxxx ], 2

Θ2D AD = [1, u, u , ux , uy , uxx , uyy ].

(S1) (S2)

These physical libraries are fixed before evaluation and shared by all methods. Algebraically equivalent conservative terms are mapped to the displayed convention before support and coefficient scoring; for example, ∂x (u2 ) = 2uux .

A.2

Metrics and Aggregation

We evaluate structural recovery separately from coefficient and field accuracy. Exact recovery requires equality between

the selected support Sb and the true support S ⋆ . We define support F1 as the harmonic mean of precision and recall over candidate terms. The relative coefficient and field errors are ∥ξˆ − ξ ⋆ ∥2 , ⋆ ∥ξ ∥2 + 10−12 ∥û − u⋆ ∥2 Eu = . ∥u⋆ ∥2 Eξ =

(S3)

Before computing Eξ , each selected equation is expressed in the shared physical basis defined in Sec. A.1, with zero coefficients assigned to absent terms. Field error is evaluated on the complete reference grid. The fixed library results report exact recoveries out of five runs and medians over the same five runs. Curves in the observation density study report medians and interquartile ranges. In the symbolic experiment, continuous parameters are summarized as mean ± sample standard deviation over runs that recover the target power-law structure.

A.3

Structured Field Adapter and Training

The structured field adapter reconstructs the continuous field as u bη,C (x, t) = bη (x) + Φη (x)⊤ C β̄(t) (S4)

using only the sampled observations. A spatial MLP jointly produces the background bη (x) and R spatial features Φη (x). Its input is the random Fourier encoding γ(x̃) = [sin(B x̃), cos(B x̃)], where the 64 rows of B are sampled once per run from a standard Gaussian distribution and then held fixed. The hidden layers use Swish activations. Here x̃ = x for KdV and KS, while each coordinate of 2D AD is divided by 5 before Fourier encoding. The temporal vector β̄(t) is a clamped uniform cubic B-spline basis centered by subtracting its empirical mean over the observed time locations, with K internal knots and M = K + 4 basis functions. The matrix C ∈ RR×M contains the trainable spline coefficients. S20 and T20 use the same adapter configuration; they differ only in the coordinates and values supplied to the reconstruction loss. Table S2 lists the benchmark-specific architecture sizes and maximum training epochs. The adapter is trained from observations without a PDE residual: min Lobs + λΦ RΦ + λt Rt , (S5) η,C

where Lobs =

1 X 2 |b uη,C (xi , ti ) − yi | , Nobs i

(S6)

1 Φη (Xg )⊤ Φη (Xg ) − I , Ng F R Z T X Rt = |c′′r (t)|2 dt,

RΦ =

r=1

(S7) (S8)

PDE

Spatial grid

Time frames

Space–time domain

S20 sensors

T20 frames

KdV KS 2D AD

512 1024 51 × 51

201 251 61

[−30, 30) × [0, 20] [0, 32π] × [0, 100] [−5, 5]2 × [0, 6]

102 204 520

40 50 12

Table S1: Benchmark grids and observation counts. S20 reports the number of fixed spatial sensors observed at every time frame; T20 reports the number of retained frames observed on the complete spatial grid. PDE

Rank R

KdV KS 2D AD

12 24 16

Hidden width × layers

Internal knots K

Basis size M

Trainable parameters

Maximum epochs

32 64 6

36 68 10

17,853 19,833 21,201

5,000 5,500 5,500

64×3 64×3 72×3

Table S2: Structured adapter settings. All models use Swish activations and the random Fourier encoding described above. Trainable parameter counts include the spatial MLP and C. where Xg contains Ng points sampled uniformly from the spatial domain and cr (t) = Cr: β̄(t) is the coefficient of spatial feature r. The Gram penalty RΦ uses 512 uniformly sampled spatial points per update, and the temporal curvature integral is approximated on a uniform 200-point time grid. Training uses AdamW with an initial learning rate of 10−3 , a cosine decay schedule over the maximum epoch budget, weight decay 10−2 , and global gradient-norm clipping at 1.0. During the first 400 epochs, C is fixed at zero to fit the background component. Afterward, η and C are optimized jointly. For KdV and KS, both regularization weights are zero for the first 800 epochs and increase linearly over the next 1,600 epochs; for 2D AD, the corresponding periods are 600 and 1,200 epochs. Their target weights are λΦ = 3 × 10−3 and λt = 5 × 10−4 . The weighted feature and temporal penalties are capped at 0.15 and 0.25 times the current observation loss, respectively. Each update samples at most 32 observed time frames and uses all spatial observations available at those frames. The MSE over all observations selects the lowest-loss reconstruction checkpoint among those tracked after the 400-epoch background phase. Early stopping begins at the same point and uses a patience of 800 epochs and a minimum decrease of 10−6 . The selected checkpoint is frozen before any weak system is constructed; equation selection does not feed back into adapter training.

B B.1

Stability-Validated Weak Selection

Weak System Construction and Continuity

Consider a candidate equation written as

ut =

X j

ξj Dαj qj (u),

(S9)

where αj ∈ Nd0 , Dαj = ∂x j , and qj is a scalar function of the state. The case αj = 0 covers algebraic terms, and α the formal adjoint is Dα∗ j = (−1)|αj | ∂x j . For a compactly supported test function ψℓ , integration by parts transfers the derivatives to ψℓ . The corresponding weak entries are bℓ (b u) = −⟨b u, ∂t ψℓ ⟩ℓ , (S10) Aℓj (b u) = ⟨qj (b u), Dα∗ j ψℓ ⟩ℓ , α

where ⟨·, ·⟩ℓ denotes integration over patch Pℓ . These entries form the regression system b(b u) ≈ A(b u)ξ. For a patch centered at cℓ , with spatial and temporal half-widths hℓ,a and hℓ,t , we use ρp (s) = (1 − s2 )p 1{|s| < 1}, d xa − cℓ,a t − cℓ,t Y ρp ψℓ (x, t) = ρp . hℓ,t hℓ,a a=1

(S11) (S12)

Because ρp and its derivatives through order p − 1 vanish at s = ±1, choosing p above the highest derivative order in each candidate library eliminates the boundary terms introduced by integration by parts. Spatial derivatives therefore act on the known test function rather than pointwise on the reconstructed field. This construction uses values from the frozen adapter while avoiding direct evaluation of its higher-order spatial derivatives. Patch centers form a regular grid inside the admissible domain. Each weak system applies an independently sampled phase shift to this grid. The field, observations, and candidate library remain unchanged; only the integration patches differ. This produces multiple views of the same frozen reconstruction and provides the variation used by SVWS to assess support stability. All integrals are evaluated using midpoint tensor grids. Table S3 lists the corresponding settings for each benchmark. For the exact weak entries, suppose that u, u b ∈ L2 (Pℓ ), the required derivatives of ψℓ lie in L2 (Pℓ ), and qj is Lj Lipschitz on an interval containing the essential ranges of u and u b over Pℓ . Then |bℓ (b u) − bℓ (u)| ≤ ∥b u − u∥L2 (Pℓ ) ∥∂t ψℓ ∥L2 (Pℓ ) , (S13) |Aℓj (b u) − Aℓj (u)| ≤ Lj ∥b u − u∥L2 (Pℓ ) × ∥Dα∗ j ψℓ ∥L2 (Pℓ ) .

(S14)

PDE

Patches per system

Nodes per patch

Spatial half-width

Temporal half-width

Kernel power p

KdV KS 2D AD

300 300 320

3,000 3,100 2,280

4.8 8.0 (1.6, 1.6)

1.0 1.0 0.48

8 8 4

Table S3: Weak system settings for each benchmark. Nodes per patch give the effective size of the midpoint tensor grid; all half-widths are in physical coordinates. These bounds follow from the Cauchy–Schwarz inequality and establish entrywise continuity of the exact weak system in the local reconstruction error under the stated range condition. The midpoint discretization satisfies analogous bounds in the corresponding quadrature-weighted norm. Comparing independently shifted systems then allows SVWS to favor supports that remain stable across weak discretizations. Support recovery also depends on library conditioning, trajectory excitation, and coefficient separation.

B.2

Selection Algorithm and Configuration

Algorithm S1 summarizes SVWS. It assigns support generation, coefficient fitting, and validation to disjoint weak systems built from the same frozen field. Independently shifted patch grids provide different weak views: generation systems produce candidate supports across sparsity levels, the fit system estimates their physical coefficients, and validation systems compare the fitted equations on held-out patches. All three stages operate after the adapter is frozen. We use Rg = 12 generation systems and the threshold multipliers {0.25, 0.5, 0.75, 1, 1.25, 1.5, 2} around the base STLSQ threshold. The resulting 12×7 = 84 support proposals are deduplicated, and term recurrence is computed over all generation-system and threshold-multiplier pairs. Terms with recurrence at least τ = 0.5 define the stable candidate pool. Each retained support is refitted by ridge regression on one independent fit system and evaluated on Rv = 7 independently shifted validation systems. The STLSQ solves and fixed-support refits both use ridge 10−4 , with at most 20 STLSQ iterations; the normalized validation risk uses ε = 10−12 . Columns are normalized during STLSQ and coefficient fitting, then mapped back to the physical library scale before validation. The base thresholds for KdV, KS, and 2D AD are 0.2, 0.4, and 0.4, respectively. These settings were selected using development seeds {42, 505, 606, 808, 909}, then shared across S20 and T20 and fixed for all test seeds.

External Baseline Implementations

All baselines receive the sampled observations and physical candidate libraries used in the main comparison. Each retains the optimization and support selection procedure specified by its source implementation or published algorithm. PDE-FIND and WSINDy operate on a shared spectral reconstruction of the samples, whereas the neural methods use the observations in their data objectives. Table S4 identifies the implementation, observation interface, and coefficient readout for each method. Settings not specified by the original source are chosen on development seeds and fixed before evaluation.

C.1

Implementation Sources and Reference Experiments

PDE-FIND and WSINDy are reimplemented following their published regression procedures. DeepMoD uses DeePyMoD v2.2.0, and Weak-PDE-LEARN uses the authors’ released Rational network. For the released neural methods, benchmark wrappers supply the MDBench observations and map native terms to the shared physical library; optimization, sparsification, and terminal reporting follow the released code. In particular, Dx (u2 ) in Weak-PDE-LEARN is mapped by Dx (u2 ) = 2uux . Source revisions and dependencies are included in the code. We reimplement PINN-SR and DL-PDE following the algorithms and training procedures described in their papers. PINN-SR uses alternating direction optimization. DL-PDE fits an observation-only neural field and applies STRidge to quantities evaluated by automatic differentiation. Their complete configurations are included in the code. We also evaluate the released neural pipelines on their reference configurations before introducing the MDBench interface and shared physical library. Weak-PDE-LEARN recovers the exact KdV-Sine support in all five runs using 4,000 measurements, 25% noise, and the paper-specified schedule for the Rational network. DeePyMoD recovers the exact support for its notebook KdV example in two of five runs under the public v2.2.0 configuration. These reference runs use each source’s own data, library, and reporting convention; they are distinct from the MDBench S20 and T20 experiments.

C.2

Configurations and Parameter Selection

Each external method uses one PDE-specific configuration, shared across S20, T20, and all five test seeds. PDE-FIND and WSINDy use ridge penalties of 10−3 for spectral reconstruction and 10−4 for coefficient refitting. DeepMoD uses a sparsification threshold of 0.1, four hidden layers of width 30 for KdV and KS and width 50 for 2D AD, and at most 100,000 iterations. Weak-PDE-LEARN uses five hidden layers of width 40 with Rational activations, 6,000 burn-in epochs, and 5,000 sparsification epochs. PINN-SR performs six alternating optimization cycles. DL-PDE applies STRidge to the fixed terminal surrogate checkpoint. Complete schedules are provided in the accompanying configurations. Table S5 reports the candidate values and selected settings. The PDE-FIND, WSINDy, and DL-PDE thresholds are selected on seeds {42, 505, 606, 808, 909} using aggregate performance across S20 and T20. We maximize exact recoveries, breaking ties by higher median support F1 , lower median fixed-support coefficient error, and the smaller threshold. The two KS adapter capacities are compared on

Algorithm S1: Stability-Validated Weak Selection on a frozen field

Input: frozen field û; candidate library; base threshold λ0 ; multipliers M; system counts Rg , Rv ; independent phase seeds; stability threshold τ ; ridge parameter α; numerical constant ε. Rg v 1 Build independently phased weak systems {(Argen , brgen )}r=1 , (Afit , bfit ), and {(Asval , bsval )}R s=1 from the same frozen field. 2 For each (r, m) ∈ {1, . . . , Rg } × M, run STLSQ at threshold mλ0 on the column-normalized generation system and record support Sr,m . P 3 Deduplicate the supports to form C. For each S ∈ C, record its generation count v(S) = r,m 1{Sr,m = S}. Compute term frequencies P πj = (Rg |M|)−1 r,m 1{j ∈ Sr,m } and stable terms Jτ = {j : πj ≥ τ }. 4 Retain Cτ = {S ∈ C : S ⊆ Jτ }. If this set is empty, set Cτ = C. efit = Afit D−1 . For every S ∈ Cτ , fit ξeS by ridge regression on (A efit , bfit ) and 5 Let Dfit contain the fit-system column norms and set A fit −1 e ˆ restore the original library scale with ξS = Dfit ξS . 6 On each validation system, compute rs (S) = ∥Asval ξˆS − bsval ∥22 /(∥bsval ∥22 + ε). Record the mean µ(S) and standard error se(S). 7 Let Smin ∈ arg minS∈Cτ µ(S) and admit every S satisfying µ(S) ≤ µ(Smin ) + se(Smin ). 8 Select the admissible support lexicographically by (|S|, µ(S), −v(S)). Output: selected support Sb and coefficients ξˆb fitted on the independent fit system and mapped to the original library scale. S

Method

Implementation

Observation interface

Reporting rule

PDE-FIND

Reimplementation following the published procedure Reimplementation following the published weak-form procedure Released DeePyMoD v2.2.0 implementation Reimplementation of the published ADO procedure Authors’ released implementation with the Rational network Reimplementation of the published two-stage procedure

Spectral reconstruction fitted to the sampled observations Shared spectral reconstruction evaluated through weak systems Sampled observations enter the neural data objective Sampled observations enter the alternating data objective Sampled observations are supplied through its data interface

STLSQ support; fixed-support ridge coefficients STLSQ support; coefficient refit on an independent weak system Terminal sparsity mask and unscaled constraint coefficients Last STRidge projection from alternating optimization Final sparsification checkpoint

A neural surrogate is fitted to the sampled observations

Terminal STRidge support; fixed-support ridge coefficients

WSINDy DeepMoD PINN-SR Weak-PDELEARN DL-PDE

Table S4: Implementations, observation interfaces, and reporting rules for the external baselines. independent development runs. Relative to (16, 32), the selected (24, 64) configuration preserves exact recovery and field accuracy under both protocols while reducing coefficient error. All selection runs are disjoint from the reported test seeds.

D.1

Additional Ablations

Adapter Representation Across Benchmarks

Table S6 extends the KS ablation in the main paper to all six fixed-library settings. Within each PDE, the four adapters have closely matched parameter counts and share the observations, optimization budget, weak-form systems, and SVWS configuration. The comparison isolates the effects of space– time factorization and fixed spatial Fourier features. KdV and 2D AD are recovered exactly by several adapters and therefore provide limited separation between representations. KS is more discriminating. The structured adapter recovers all five supports under both protocols, whereas the coordinate MLP with Fourier features recovers three under S20 and one under T20, and the coordinate MLP without Fourier features recovers none. Removing Fourier features from the structured adapter also substantially increases its KS field error. Fourier encoding alone therefore does not explain the recovery gain; its benefit is strongest when combined with the structured separation of space and time. Table S7 controls for coordinate MLP capacity. Increasing

the hidden width from 72 to 128 improves field accuracy and exact recovery, yet the 49,665-parameter MLP remains less reliable than the 19,833-parameter structured adapter.

Hyperparameter

Candidate values

Selected value(s)

KS adapter (R, K) (R, K) ∈ {(16, 32), (24, 64)} PDE-FIND threshold {0.05, 0.1, 0.2, 0.4, 0.8, 1, 2, 5} WSINDy threshold {0.05, 0.1, 0.2, 0.4, 0.8, 1, 2, 5} DL-PDE STRidge tolerance {0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1, 2, 5, 10, 20, 50, 100, 200, 500, 1000}

(24, 64) (0.2, 5, 0.05) (0.2, 2, 0.4) (1, 2, 2)

Table S5: Hyperparameter candidate sets and selected values. For the three baseline rows, selected values are ordered as KdV, KS, and 2D AD; each PDE-specific choice is shared across S20 and T20. KdV

2D AD

Adapter

S20

T20

S20

T20

S20

T20

(a) Exact support recovery Structured adapter (ours) Structured adapter without Fourier features Coordinate MLP with Fourier features Coordinate MLP without Fourier features

5/5 5/5 5/5 3/5

5/5 5/5 5/5 4/5

5/5 3/5 3/5 0/5

5/5 5/5 1/5 0/5

5/5 5/5 5/5 5/5

(b) Median field error Eu Structured adapter (ours) Structured adapter without Fourier features Coordinate MLP with Fourier features Coordinate MLP without Fourier features

0.1171 0.1336 0.0450 0.1699

0.0239 0.0979 0.0236 0.0716

0.0817 0.6588 0.4397 0.9046

0.0815 0.6503 0.4365 0.8674

0.0028 0.0055 0.0026 0.0139

0.0026 0.0046 0.0034 0.0124

Table S6: Adapter ablations across all six fixed-library settings. Within each PDE, the four variants have closely matched parameter counts and share the observations, optimization budget, weak-form systems, and SVWS configuration. Panel (a) reports exact-support recoveries over five seeds, and panel (b) reports median relative field error. Bold marks the best value in each column, including ties. KS-S20 Adapter Structured adapter (ours) Coordinate MLP with Fourier features (width 72) Coordinate MLP with Fourier features (width 128)

KS-T20

Parameters

Exact

19,833 19,873 49,665

5/5 3/5 4/5

0.082 0.440 0.340

5/5 1/5 3/5

0.082 0.437 0.315

Table S7: Capacity comparison on KS. The two coordinate MLPs differ only in hidden width; the structured adapter is included as a reference at a comparable parameter count. All rows use the same observations, training budget, and SVWS configuration. Eu is the median relative field error over five seeds.

D.2

Selector Components and Sensitivity

Table S8 evaluates selector sensitivity on KS and 2D AD over the same five test seeds. Panel (a) uses the clean frozen fields and configurations from the main comparison and extends the three-point comparison to the complete seven-point threshold path. Each STLSQ row uses the first generation weak system, so only the threshold changes. Low thresholds retain extra KS terms, whereas high thresholds remove governing terms from 2D AD. No single threshold recovers every run across all four settings. SVWS instead considers supports along the complete path and across multiple weak systems, recovering all five runs in every setting. Panel (b) compares selector rules on fields reconstructed from observations with 10% noise. Within each run, all variants share the frozen field, candidate pool, coefficient refits, and validation systems. Selecting the minimum mean validation risk over the full pool frequently retains extra terms, especially on KS. Removing stability filtering loses one 2D AD-S20 recovery. Term-wise majority voting matches

SVWS in these four settings, showing that recurrence alone can separate the correct terms once the field is reliably reconstructed. SVWS achieves the same exact recovery while combining recurrence filtering with held-out weak validation.

D.3

Adapter Objective

The adapter objective augments the observation loss with penalties that encourage feature diversity and temporal smoothness. To isolate their contribution, we compare the full objective with λΦ = λt = 0. The variants use the same architecture, optimizer, training budget, observations, and SVWS configuration, and both retain AdamW weight decay. Table S9 reports paired runs for all six S20 and T20 settings over the five test seeds. Both objectives recover every S20 support and all KdVT20 supports. Removing the explicit penalties causes one KS-T20 failure and increases the maximum field error from 0.037 to 0.085 on KdV-T20, from 0.108 to 0.185 on KS-T20, and from 0.004 to 0.087 on 2D AD-T20, while the S20 errors

(a) Sensitivity to a fixed STLSQ threshold on clean frozen fields Selection KS-S20 SVWS STLSQ (λ/λ0 = 0.25) STLSQ (λ/λ0 = 0.5) STLSQ (λ/λ0 = 0.75) STLSQ (λ/λ0 = 1) STLSQ (λ/λ0 = 1.25) STLSQ (λ/λ0 = 1.5) STLSQ (λ/λ0 = 2)

5/5 1/5 4/5 5/5 5/5 5/5 5/5 5/5

(b) Selector variants with 10% observation noise Selection KS-S20 SVWS Without stability filtering Minimum mean validation risk Term-wise majority vote

5/5 5/5 1/5 5/5

KS-T20

2D AD-S20

2D AD-T20

5/5 3/5 4/5 4/5 4/5 4/5 5/5 5/5

5/5 5/5 5/5 5/5 5/5 5/5 2/5 0/5

KS-T20

2D AD-S20

2D AD-T20

5/5 5/5 1/5 5/5

5/5 4/5 3/5 5/5

5/5 5/5 5/5 5/5

Table S8: Threshold and selector sensitivity on frozen fields. Panel (a) applies STLSQ to the first generation weak system along the prespecified threshold path on the clean fields from the main comparison. Panel (b) compares selector variants on fields reconstructed from observations with 10% noise. Entries are exact-support recoveries over five seeds. KdV

2D AD

Objective

Metric

S20

T20

S20

T20

S20

T20

Full objective

Exact recovery Median Eu Maximum Eu

5/5 0.117 0.143

5/5 0.024 0.037

5/5 0.082 0.127

5/5 0.082 0.108

5/5 0.003 0.003

5/5 0.003 0.004

No explicit penalties

Exact recovery Median Eu Maximum Eu

5/5 0.116 0.143

5/5 0.030 0.085

5/5 0.073 0.111

4/5 0.140 0.185

5/5 0.002 0.003

5/5 0.003 0.087

Table S9: Ablation of the adapter objective under the final benchmark configurations. The variants use identical architectures, observations, optimization budgets, and SVWS settings; the second sets λΦ = λt = 0 while retaining AdamW weight decay. Exact recovery and Eu are computed over the same five seeds. remain comparable. The explicit penalties therefore improve robustness primarily when temporal coverage is sparse.

Support Dynamics under Joint Optimization

To complement the representative trajectories in the main paper, we examine all five test seeds for the three settings visualized there. Table S10 also includes DeepMoD on KdVS20, a second setting in which the exact support persists to termination. The table summarizes whether the exact support is reached at any recorded checkpoint, is present at the checkpoint minimizing DeepMoD’s fixed held-out objective, and remains at termination. Figure S1 shows the complete support F1 trajectories for the three settings in the main paper. Because the two methods optimize different training objectives, we compare their paths using support F1 at each recorded checkpoint. DeepMoD support is given by its learned sparsity mask at threshold 0.1, whereas Weak-PDELEARN (WPL) support is defined by coefficients whose absolute value exceeds 10−3 . For DeepMoD, the fixed held-out objective is the sum of data MSE and PDE-residual MSE on its held-out split. Exact-support intervals are derived from

Method

Setting

DeepMoD KdV-S20 DeepMoD 2D AD-S20 DeepMoD KS-S100 WPL 2D AD-S20

Ever Held-out min. Terminal 5/5 5/5 0/5 5/5

5/5 5/5 0/5 –

5/5 5/5 0/5 0/5

Table S10: Exact support along recorded optimization paths, reported as runs out of five. Ever indicates an exact-support checkpoint; Held-out min. evaluates the DeepMoD checkpoint minimizing data MSE plus PDE-residual MSE on a fixed held-out split; Terminal denotes the final state. The held-out minimum is unavailable for Weak-PDE-LEARN because its test weak functions are resampled at each epoch. Ever is a retrospective path summary rather than a checkpoint selection rule. the resulting binary masks, and the recorded trajectories are plotted without smoothing. Across all five seeds, the trajectories exhibit the three behaviors highlighted in the main paper. DeepMoD reaches and retains the exact support in every KdV-S20 and 2D AD-S20 run. On KS-S100, it never reaches the exact support, despite complete spatial coverage. Weak-PDE-LEARN reaches the

(a) Persists to termination

(b) Appears only transiently

DeepMoD · 2D AD-S20

Weak-PDE-LEARN · 2D AD-S20

1.0

2.7%

0.5 0.0 1.0

seed: 1301

1.3%

0.5

Support F1

0.0 1.0

seed: 1709

1.8%

0.5 0.0 1.0

seed: 2203

3.1%

0.5 0.0 1.0

seed: 2917

1.9%

0.5 0.0

seed: 3571

100

Training progress (%)

Figure S1: Support F1 trajectories across all five seeds for the three settings shown in the main paper. Each row corresponds to one seed. The narrow subcolumn magnifies the first 12% of Weak-PDE-LEARN training, with labels marking the first exact-support visit for each seed. Curves show unsmoothed support F1 over training progress normalized within each recorded path. Teal segments mark exact-support intervals, stars mark the first exact visit, and terminal circles are teal when exact and vermilion otherwise. exact support in every 2D AD-S20 run but loses it before termination. Changing the checkpoint can therefore alter the reported equation when the exact support is transient, whereas checkpoint selection cannot recover a support that is absent throughout the recorded path. These paths motivate separating field fitting from equation selection: selection acts on one frozen reconstruction rather than inheriting the support at a particular joint-training checkpoint. In the representative Weak-PDE-LEARN trajectory shown in the main paper, the exact support first appears at epoch 92 with ût = 0.16ux + 0.52uy + 0.37uxx + 0.49uyy ; at epoch 5000, the terminal iterate retains only ût = 0.64uy + 0.77uyy .

Symbolic Discovery Beyond a Fixed Library

The symbolic experiment tests whether the frozen reconstruction can also support discovery when candidate expressions are generated rather than drawn from a fixed library. The main paper reports aggregate recovery rates and parameter estimates. Here we describe the benchmark and symbolic search protocol, report example expressions from unsuccessful PySR runs, and list the principal configuration settings. We fix the outer operator ∂xx and search for the nonlinear diffusion function q(u) in ut = ∂xx q ⋆ (u),

q ⋆ (u) = 0.1u1.73

(S15)

Because ∂xx annihilates constants, q(u) is identifiable only up to an additive constant. We solve the equation on a periodic 192 × 121 grid over x ∈ [0, 1) and t ∈ [0, 0.3]. The initial condition is u(x, 0) = 0.75 + 0.30 sin(2πx) + 0.14 cos(4πx − 0.35) + 0.07 sin(6πx + 0.60). We generate the trajectory using second-order periodic finite differences and a BDF integrator with relative and absolute tolerances 2 × 10−7 and 10−9 , respectively. At each noise level, the structured adapter is trained on observations at 20% fixed spatial sensors over all time points, and five independent runs are evaluated. Independent zero-mean Gaussian noise is added with standard deviation equal to 0%, 5%, or 10% of the global standard deviation of the clean field. The target exponent 1.73 is excluded from the grid used to initialize the up terminals. Genetic programming (GP) constructs candidate expressions for q(u). We run two GP proposal searches, each using weak systems constructed from independently shifted patch grids. Within each search, expressions with the same canonical structure are grouped, and the expression with the smallest normalized weak-form residual in each family is retained. The representatives from the two searches are then merged and deduplicated. A separate weak system fits the continuous parameters of each retained expression, and three held-out weak systems evaluate the fitted candidates. SVWS selects from this fixed candidate pool by minimizing the mean heldout normalized residual plus 10−3 times expression complexity. The selected expression retains the parameters fitted

Noise (%) 0 5 10

Example PySR expression 0.100u

1.727

2.727

+ 0.0002u − 0.0001u + C 0.0667u2 + 0.0387u + C 0.0624u2 + 0.0437u + C

Table S11: Example PySR expressions from the unsuccessful run with the smallest prespecified seed at each noise level. Component

Configuration

Structured field adapter Adapter training

14 spatial features; 80 Fourier features; three 128-unit hidden layers; 40 internal knots; Softplus output AdamW with weight decay 10−2 ; initial learning rate 10−3 with cosine annealing; gradient norm clipped to 1.0; at most 3200 epochs; patience 600; at most 32 time frames per minibatch background warm start for 300 epochs; penalties enabled at epoch 500 and ramped over 900 epochs; terminal iterate (λΦ , λt ) = (3 × 10−3 , 5 × 10−4 ); weighted penalties capped at 0.12 and 0.20 times the observation loss leaves u, free constants, and up , p ∈ {0.5, 0.6, . . . , 3.0}; operators +, −, ×, and protected division with sampling probabilities (0.28, 0.18, 0.38, 0.16) denominator floor 10−4 ; intermediate outputs clipped to [−50, 50]; power input floor 10−6

Adapter schedule Regularization GP primitives

Protected operations

Table S12: Structured reconstruction and GP proposal settings for the symbolic experiment. before validation. A run counts as recovering the power-law family when the selected expression is algebraically equivalent to q̂(u) = κ̂um̂ + C for some additive constant C. The reported κ̂ and m̂ are summarized over runs that recover this family. SVWS operates on the field reconstructed from 20% fixed sensors. For PySR, we first estimate q(u) from the full noisy grid using weak equations and a fixed basis of seven centered cubic Bsplines in u, then apply PySR to the resulting 512 (u, q̂(u)) pairs. We use PySR 1.5.10 with model_selection=best; the principal settings are listed in Tables S12 and S13. Table S11 shows expressions selected in unsuccessful PySR runs, including forms with additional terms and polynomial alternatives. These examples complement the aggregate recovery rates in the main paper by illustrating the structures selected when the target power-law family is not recovered.

Component SVWS extension Proposal search

Evolution Parameter fit Continuous fit

Validation

PySR baseline PySR target construction

PySR search

Configuration 2 independent GP searches; per search, separate 128-patch fit and generation systems and one 192-patch proposal-validation system; (hx , ht ) = (0.105, 0.03); kernel power p = 8; population 72; 16 generations elite 8; mutation 0.35; crossover 0.55; maximum depth 6; complexity cap 7; complexity penalty 10−5 1 weak system; 128 patches; (hx , ht ) = (0.14, 0.045); kernel power p = 8 at most 32 candidates; exponents [0.2, 4]; constants [−5, 5]; scale [−10, 10]; at most six parameters, one division, and 35 function evaluations 3 systems; 192 patches each; (hx , ht ) = (0.105, 0.03); kernel power p = 8; 384 quadrature points; mean held-out risk plus 10−3 times expression complexity seven centered cubic B-spline basis functions in u; 512 weak patches; (hx , ht ) = (0.2, 0.03); kernel power p = 6; 384 quadrature points; 512 state values; penalties of 10−6 for second differences and the additive offset version 1.5.10; 100 iterations; 10 populations of 40; operators (+, ×, power); maximum size 15 and depth 6; power exponent complexity at most 1; parsimony 10−4

Table S13: Symbolic search, selection, and PySR settings.

G G.1

Sensitivity to Observation Density and Noise

Observation Density

Figure S2 summarizes the observation-density study for five methods spanning the proposed, classical, and neural approaches. KdV is evaluated with 5%, 10%, 15%, and 20% fixed spatial sensors, whereas KS and 2D AD are evaluated with 20%, 50%, 80%, and 100%. For each benchmark and seed, the sensor masks are nested across densities. The same five test seeds, method-specific configurations, and benchmark-specific candidate libraries defined in Sec. A.1 are used throughout each sweep. The 20% points are the runs reported in the main comparison. Curves show medians and interquartile ranges over the five seeds. Increasing sensor density generally reduces field error, but the corresponding gains in exact recovery depend on the method and PDE. Our method moves from 1/5 exact recoveries at 5% KdV sensors to 5/5 from 10% onward, and remains at 5/5 throughout the KS and 2D AD sweeps. PDE-FIND on KS reaches 5/5 from 50% coverage, while WSINDy on 2D AD reaches 3/5 at 50% and 5/5 from 80%. DeepMoD and PINN-SR do not recover the exact KS support at any tested density. Thus, denser observations improve field reconstruction but do not by themselves guarantee recovery of the correct equation terms.

G.2

Observation Noise

For the structured adapter with SVWS, we evaluate 10% observation noise on KS and 2D AD under both S20 and T20. Independent zero-mean Gaussian noise is added to the sampled values with standard deviation equal to 10% of the global standard deviation of the clean field. The sampled

Ours

PDE-FIND

KdV

WSINDy

DeepMoD

Kuramoto–Sivashinsky

PINN-SR

2D advection–diffusion

Support F1

1.0

0.5

0.0

Field error Eu (log scale)

Ours PDE-FIND WSINDy DeepMoD PINN-SR 10

Exact support (/5) 1 5 0 0 0 0 4 5 4 5

5 0 4 5 5

Exact support (/5) 5 5 0 5 1 1 0 0 0 0

5 3 5 5 5

5 5 1 0 0

Exact support (/5) 5 5 0 0 0 3 5 5 0 0

5 0 5 5 0

−1

−2

10% 15% 20% Observed spatial fraction

20%

50% 80% 100% Observed spatial fraction

20%

50% 80% 100% Observed spatial fraction

Figure S2: Effect of observation density under clean, nested fixed-sensor sampling. Curves and error bars show the median and interquartile range of support F1 and field error over five seeds. The middle panels report exact-support recoveries out of five. Filled markers and outlined columns identify the 20% setting used in the main comparison. coordinates and all adapter and selector settings are identical to those in the corresponding clean runs. Table S14 compares the noisy and clean results. Clean

10% noise

PDE

Obs. Exact Eξ

Eu Exact Eξ

KS KS

S20 T20

5/5 0.035 0.082 5/5 0.115 0.166 5/5 0.055 0.082 5/5 0.083 0.097

2D AD S20 2D AD T20

5/5 0.010 0.003 5/5 0.011 0.006 5/5 0.009 0.003 5/5 0.012 0.007

Table S14: Robustness to 10% observation noise. Exact reports recoveries out of five; errors are medians over the same seeds. The structured adapter with SVWS recovers the exact support in all 20 noisy runs. The largest increases in coefficient and field error occur on KS, whereas the 2D AD estimates remain close to their clean counterparts. Thus, at the tested noise level, support recovery is unchanged even though coefficient and field accuracy deteriorate.

Code, Data, and Computational Resources

Code and configurations will be made publicly available upon acceptance. The release will provide the implemen-

tations, configurations, and run-level results for the reported studies. It will include the structured adapter and SVWS, baseline wrappers, observation records for the primary fixedlibrary experiments, and scripts that verify the reported table values and rebuild the numerical figures for optimization paths and observation density. A README will provide execution commands, while a paper map will link reported results to their implementations, manifests, and included evidence. Environment files will list the required third-party dependencies. The release will specify the public MDBench trajectory and state index and include the script used to generate the symbolic diffusion benchmark. Experiments were run on two Ubuntu servers. One server had an Intel Xeon Silver 4309Y CPU, 125 GiB of memory, and two 24-GiB NVIDIA RTX 4090 GPUs; the other had an Intel Xeon w5-3525 CPU, 125 GiB of memory, and three 24-GiB NVIDIA RTX 4090 GPUs. Both environments ran Python 3.10. The first used PyTorch 2.5.1 with CUDA 12.1, and the second used PyTorch 2.7.1. Environment specifications and versions of the principal dependencies will accompany the release. Neural models were trained on one GPU per run using 32-bit floating-point model tensors; source trajectories were stored in 64-bit precision.

Record · ID 422234 · SHA-256 6b38282cf7f1c261
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.