A zero-one law for one-shot system identification Nicolas Boulléa,1 , Diana Halikiasb , Samuel E. Ottoc , and Alex Townsendd a
Department of Mathematics, Imperial College London, London, SW7 2AZ, UK; b Courant Institute, New York University, New York, NY 10012, USA; c Sibley School of Mechanical and Aerospace Engineering, Cornell University, Ithaca, NY 14853, USA; d Department of Mathematics, Cornell University, Ithaca, NY 14853, USA
A
B
C
D
E
F
ODE
PDE
Can a model be identified from one experiment? We study analytic systems that are linearly parameterized by a combination of prescribed dictionary terms, such as partial differential operators and dynamical systems. For a single input-response pair, recovery is possible exactly when the evaluated dictionary terms are linearly independent. We prove a sharp zero-one law: either no input uniquely determines the coefficients, or almost every random input sampled from a nondegenerate Gaussian measure does. This dichotomy reduces one-shot system identification to a question about degenerate inputs and provides an a posteriori certificate for any recovered model. Numerical examples recover dynamical systems, nonlinear partial differential equations, and structured matrix families from single trajectory data, while also detecting when an extra probe is necessary. scientific machine learning | system identification | analytic systems | matrix recovery | experimental design
ow many experiments are needed to identify physical laws from data? A fundamental challenge in scientific machine learning (1) is to learn models without requiring an excessive amount of data, which may be costly or difficult to obtain. However, modern scientific machine learning can be strikingly data intensive. For example, recent machine learning weather models are trained on decades of reanalysis data (2), while geoscience models require large geo-datasets (3). This makes sample complexity a central issue: how many simulations, trajectories, or experiments are actually needed? In operator learning, this question has led to theory explaining when solution operators of partial differential equations (PDEs) can be learned efficiently from input-output data (4). Model identification concerns a different task. Rather than learn a black-box solution operator with a neural network (5), one seeks a human-interpretable formulation given by the coefficients of a governing law in a prescribed dictionary. While recent years have seen numerous breakthroughs in this area through the introduction of symbolic regression (6–8) and sparsity-promoting algorithms such as SINDy (9, 10), the potential identifiability of governing equations from single trajectory data has remained largely unexplained. Nearby fields address related but different problems. Operator inference learns reduced surrogate dynamics from simulation data (11). Structural identifiability and parameter estimation usually assume the model form and observation protocol are fixed, and then ask which parameters, or functions of parameters, are determined. Algebraic criteria have recently been developed for ordinary differential equations (12). Model recovery also appeared in specialized inverse problems with the establishment of global uniqueness for the Calderón problem to recover a conductivity coefficient from boundary measurements (13). Our work is motivated by recent developments in scientific machine learning (9, 10), but is genuinely different in its focus. We aim to characterize and quantify the set of
H
Matrix
arXiv:2607.15832v1 [math.NA] 17 Jul 2026
This manuscript was compiled on July 20, 2026
Fig. 1. Identifiability across analytic systems. The same rank criterion certifies oneshot recovery for nonlinear differential equations (A, Allen–Cahn; B, Navier–Stokes), dynamical systems (C, Lorenz; D, Duffing), and structured matrix families, while also detecting when an additional probe is needed (E, Hankel) and when one suffices (F, circulant).
inputs that guarantee recovery of the ground-truth model. We consider general analytic systems parameterized by a linear combination of dictionary terms and aim to identify the coefficients from a single input-response pair. The unknown object may be a differential equation, a continuous-time dynamical system, or a matrix drawn from a structured linear family (cf. Fig. 1). These settings all involve a possibly large known dictionary of analytic terms that can be evaluated at a single observed input-response pair, especially for applications where measurement alters the system irreversibly. The problem becomes a coefficient-identification problem where the unknown coefficients must be recovered from that observation, and selecting a good input is crucial to ensure recovery. Author contributions: N.B., D.H., S.E.O., and A.T. designed research; performed research; analyzed data; and wrote the paper. The authors declare no competing interest. 1
To whom correspondence should be addressed. E-mail: [email protected].
July 20, 2026
|
1
Results
A 2
For one-shot recovery, we assume that one can probe an unknown system with a forcing term f and observe the corresponding response u. Our starting point is the observation that coefficient identification reduces to a linear algebra task. We then aim to identify the coefficients in the following analytic system from (f, u):
B2 1
1
z 0 −1
y 0
−2 2 1
−1
y 0
L(u, f ) :=
N X ∗ i=1
ci Di (u, f ) = D0 (u, f ),
−1
[1]
where Di : U × F → V are prescribed analytic maps between Banach spaces and the coefficients c∗i ∈ R are unknown. Boundary conditions and initial conditions are built into the spaces U and F. Analyticity is a broad structural assumption here. It includes linear maps, polynomial and standard nonlinear dictionary terms, partial derivatives, weak forms of differential operators, and finite-dimensional matrix families. We assume that the solution map G : F → U associated with Eq. (1) is well-posed and analytic and that the input space F has a countable Schauder basis. This is the infinitedimensional analogue of a coordinate system, where each input can be expanded in a convergent series in countably many basis functions, and holds for many separable function spaces used in analysis and computation. Given one input f ∈ F to the system modeled by Eq. (1) with output u = G(f ), we define an auxiliary parameter-toPN recovery map Lf,u : RN → V as Lf,u (a) = a D (u, f ) i=1 i i for a ∈ RN . The coefficients c∗i in Eq. (1) are identifiable from the single observation (f, u) if and only if this map is injective. We then aim to characterize the set of degenerate inputs f ∈ F for which Lf,G(f ) is not injective, defined as
Fdeg (c∗ ) = f ∈ F | ∃a ∈ RN \ {0}, Lf,G(f ) (a) = 0 .
[2]
To express this question probabilistically, i.e., how likely it is that a randomly chosen input allows for identification, we consider sampling inputs from nondegenerate Gaussian measures on F, so that every nonzero continuous linear observation of the input is Gaussian with positive variance — a common practice in scientific machine learning (5). In finite dimensions, this is any Gaussian with nonsingular covariance. The following theorem establishes a zero-one law for one-shot system identification; either model identification is generically possible, or it fails for every input. We phrase this dichotomy in terms of a notion we call compactly slice-null sets, which are sets whose intersections with large enough finite-dimensional affine subspaces translated over any given compact set have zero Lebesgue measure on each slice. Theorem 1 (Zero-one law for system identification). Exactly one of the following is true: i) Fdeg (c∗ ) is compactly slice-null and γ(Fdeg (c∗ )) = 0 for every nondegenerate Gaussian measure γ on F. ii) Fdeg (c∗ ) = F.
Thus, for a fixed analytic system, if a single input succeeds in recovery, the set of degenerate inputs is negligible and recovery is almost surely possible for inputs sampled from a nondegenerate Gaussian measure. On the other hand, if no such input exists and Fdeg (c∗ ) = F, the obstruction is more fundamental, based on properties like a hidden symmetry or 2
|
−2 −2
−1
0
1
2
−2 −2
x
−1
0
1
2
x
Fig. 2. Zero level set of analytic functions. Lemma 1 ensures that the zero level sets of a two-dimensional (A) and a three-dimensional (B) analytic function, respectively illustrated as black lines and shaded surfaces, are negligible when intersected with a finite-dimensional slice.
a dictionary relationship that cannot be avoided by changing the forcing. The proof of Theorem 1 studies the prevalence of inputs that make Lf,G(f ) degenerate and relies on Lemma 1, which extends a standard result in geometric measure theory stating that zero level sets of finite-dimensional analytic functions are negligible (14). Lemma 1 (Zero level set of analytic functions). Let F be a Banach space with a countable Schauder basis and ϕ : F → R be a nonzero analytic function. Then, the zero level set of ϕ is compactly slice-null and ϕ−1 (0) has measure zero with respect to any nondegenerate Gaussian measure on Borel sets of F.
The proof of Lemma 1 combines a compactness argument with a restriction to finite dimensions to deduce that the zero level set of the analytic function is compactly slice-null, then exploits a disintegration property of Gaussian measures (15). To prove Theorem 1, we construct a nonzero analytic function ϕ : F → R whose zero level set contains Fdeg (c∗ ). Suppose that identification succeeds for one input f0 ∈ F. Then, one can construct linearly independent test functionals ψ1 , . . . , ψN ∈ V ∗ such that the following analytic function
ψ1 (Lf,G(f ) e1 ) .. ϕ : f 7→ det . ψN (Lf,G(f ) e1 )
··· .. . ···
ψ1 (Lf,G(f ) eN ) .. . ψN (Lf,G(f ) eN )
[3]
is nonzero at f0 and vanishes on the set of degenerate inputs, where {e1 , . . . , eN } is the canonical basis of RN . Applying Lemma 1 to ϕ in Eq. (3) yields the desired conclusion since the determinant is an analytic function of the input. We illustrate the zero-one law in Fig. 2 by plotting the analytic function ϕ for two-dimensional and three-dimensional analytic systems, along with their zero level sets. Here, we verify the conclusions of Lemma 1 by observing that these sets are compactly slice-null, in the sense that they have zero Lebesgue measure when intersected with a finite-dimensional slice. This implies that almost every input enables the identification of the coefficients in these models. Theorem 1 also provides us a practical identification algorithm that comes with a posteriori certificate for one-shot recovery. After observing the response u to a random Gaussian input f , one evaluates the dictionary terms Di (u, f ) through a finite number of test functionals ψ1 , . . . ψN : V → R and checks whether the resulting finite-dimensional recovery matrix is full rank, i.e., has nonzero determinant in Eq. (3). If this is the case, then the coefficients are uniquely determined by that input, as well as almost every other random input. Boullé et al.
Discussion The zero-one law applies to a wide range of analytic systems, including differential operators, continuous-time dynamical systems, and structured matrix families, PNas illustrated in Fig. 1. For differential operators L(u) := i=1 c∗i Di (u) = f , where the dictionary terms only depend on the state variable and partial derivatives of u, one can construct a smooth solution u0 such that the evaluated dictionary terms Di (u0 ) are linearly independent. Selecting f0 = L(u0 ) yields a nondegenerate input that identifies the coefficients. Analytic well-posedness then puts the system in alternative i) of Theorem 1: almost every random forcing identifies the coefficients. This covers weak and strong formulations and is illustrated by the Allen– Cahn example in Fig. 1A. PNFor ∗analytic continuous-time dynamical systems u̇(t) = c Di (t, u(t), f (t)), the input may be an initial condition, i=1 i a time-dependent forcing, or another experimentally chosen parameter. The zero-one alternative can go either way, but once one trajectory is found to identify the coefficients, almost every random input does. Thus even the chaotic Lorenz system (Fig. 1C) can be identified from a single trajectory in our dictionary model; the Duffing example (Fig. 1D) similarly recovers the governing second-order equation from one observed trajectory. The same principle also extends to matrix recovery. Suppose an unknown matrix belongs to the family A(c∗ ) = PN ∗ c A for matrices Ai ∈ Rm×n , and is observed only i=1 i i through matrix-vector products A(c∗ )x and A(c∗ )⊤ y. For any fixed collection of probes, identifiability is based on the injectivity of a finite-dimensional recovery map. Circulant matrices fall in the recovering alternative with one random probe, while Hankel matrices do not; stacking two probe-response pairs gives a new block recovery map that is full rank. Theorem 1 also applies to multi-shot recovery, where several input-response pairs are observed. In this case, the recovery map is a block map that stacks the individual maps for each input-response pair. The zero-one law then applies to this concatenated map, and the same dichotomy holds: either no collection of inputs identifies the coefficients, or almost every random collection does. Thus, we can study multishot recovery with the same theory by using the zero-one law applied to a concatenated observation space and increase the number of observations until the recovery map is full rank. The set of degenerate inputs can be precisely analyzed in specific cases. For single-variable polynomial ODEs, where the dictionary consists of polynomials of u and its derivatives, we find that the degenerate set Fdeg (c∗ ) does not contain rough functions that are non-differentiable at a dense set of points. On the other hand, the degenerate set is trivial for uniformly elliptic PDEs with known higher-order terms. In summary, the zero-one law provides a sharp characterization of one-shot system identification for analytic systems. It shows that the success of model recovery from a single experiment is determined by the existence of a single informative input, and that almost every random input will also be informative. This result has practical implications for experimental design and model discovery, as it shifts the focus from collecting many similar snapshots to designing experiments that break dictionary symmetries and using the observed linear system rank as a certificate of success. Our analysis is idealized in that it assumes exact model Boullé et al.
evaluation and noiseless observations. In practical settings with discretization and measurement error, near-degeneracy can still be detected through the smallest singular value of the recovery matrix, and multiple experiments can be combined to improve conditioning. Extending the zero-one framework to quantitative stability bounds under noise and to adaptive experiment design is a natural next step. Other open questions include the selection of test functionals and recovery from partial observations of the system. Materials and Methods A key technical result behind the zero-one law is Lemma 1, whose proof consists of applying the characterization of root sets of analytic functions (14) to a finite dimensional auxiliary map. To obtain the second part, we rely on a disintegration property of Gaussian measures on Banach spaces (15). The compactly slice-null property is illustrated in Fig. 2 by plotting the zero level set of two- and three-dimensional synthetic algebraic functions. For each experiment in Fig. 1, we observe one input-response pair (f, u), evaluate the prescribed dictionary terms, and assemble the finite-dimensional recovery matrix by applying selected test functionals to these evaluations. Coefficients are estimated by solving the resulting linear system using rank-revealing column-pivoted QR factorization, and are recovered with machine precision accuracy (see SI Appendix). Identifiability is certified a posteriori by full column rank of the recovery matrix. Data, Materials, and Software Availability. Code and datasets are publicly
available on GitHub at https://github.com/NBoulle/wias.
ACKNOWLEDGMENTS.
The work of A.T. was supported by the Defense Advanced Research Projects Agency (DARPA) through The Right Space (TRS) Disruption Opportunity (DARPA-PA-24-04-07) and National Science Foundation CAREER grant DMS-2045646. 1. GE Karniadakis, et al., Physics-informed machine learning. Nat. Rev. Phys. 3, 422–440 (2021). 2. I Price, et al., Probabilistic weather forecasting with machine learning. Nature 637, 84–90 (2025). 3. KJ Bergen, PA Johnson, MV de Hoop, GC Beroza, Machine learning for data-driven discovery in solid Earth geoscience. Science 363, eaau0323 (2019). 4. N Boullé, D Halikias, A Townsend, Elliptic PDE learning is provably data-efficient. Proc. Natl. Acad. Sci. USA 120, e2303904120 (2023). 5. Z Li, et al., Fourier neural operator for parametric partial differential equations in International Conference on Learning Representations. (2021). 6. J Bongard, H Lipson, Automated reverse engineering of nonlinear dynamical systems. Proc. Natl. Acad. Sci. USA 104, 9943–9948 (2007). 7. M Schmidt, H Lipson, Distilling free-form natural laws from experimental data. Science 324, 81–85 (2009). 8. SM Udrescu, M Tegmark, AI Feynman: A physics-inspired method for symbolic regression. Sci. Adv. 6, eaay2631 (2020). 9. SL Brunton, JL Proctor, JN Kutz, Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proc. Natl. Acad. Sci. USA 113, 3932–3937 (2016). 10. K Champion, B Lusch, JN Kutz, SL Brunton, Data-driven discovery of coordinates and governing equations. Proc. Natl. Acad. Sci. USA 116, 22445–22451 (2019). 11. B Kramer, B Peherstorfer, KE Willcox, Learning nonlinear reduced models from data with operator inference. Annu. Rev. Fluid Mech. 56, 521–548 (2024). 12. H Hong, A Ovchinnikov, G Pogudin, C Yap, Global identifiability of differential models. Commun. Pure Appl. Math. 73, 1831–1879 (2020). 13. J Sylvester, G Uhlmann, A global uniqueness theorem for an inverse boundary value problem. Ann. Math. 125, 153–169 (1987). 14. H Federer, Geometric Measure Theory. (Springer), Reprint of the 1969 edition, (1996). 15. VI Bogachev, Gaussian measures. (AMS), (1998).
July 20, 2026
|
3
Supporting Information Text Contents 1 Methods
1
2 Numerical examples A Allen–Cahn equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B Navier–Stokes equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C Lorenz system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . D Duffing oscillator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . E Hankel and circulant matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2 2 2 3 3 4
3 Zero level set of analytic functions
4
1. Methods The main contribution of this article is Theorem 1 from the main text which provides a zero-one law for the identifiability of linearly parameterized operators. In this supplementary text, we explain how this result leads to a weak identification procedure for linearly parameterized analytic operators, including settings in which the residuals are tested through functionals rather than inner products, and can be viewed as a generalization of Weak SINDy (1, 2) to arbitrary Banach spaces. We begin with a Hilbert-space formulation leading to a simple and computable recovery algorithm and consider an operator equation of the form: L(u, f ) =
N X ∗ i=1
ci Di (u, f ) = D0 (u, f ),
[1]
where c∗ = (c∗1 , . . . , c∗N ) ∈ RN is an unknown coefficient vector, u ∈ U and f ∈ F belong to Banach spaces. Here, D = {Di }N i=0 is a prescribed dictionary of analytic operators Di : U ×F → V, where V is a residual space. Given an observed input-output pair (f, u), we aim to determine whether the functions D1 (u, f ), . . . , DN (u, f ) ∈ V are linearly independent. When V is a real Hilbert space, this can be diagnosed through the associated Gram matrix defined as Gi,j = ⟨Dj (u, f ), Di (u, f )⟩V for 1 ≤ i, j ≤ N . When G is invertible, the dictionary terms are linearly independent and the coefficient vector c∗ is uniquely determined from the single pair (f, u) by one linear solve against the right-hand side b = [⟨D0 (u, f ), D1 (u, f )⟩V . . . ⟨D0 (u, f ), DN (u, f )⟩V ]⊤ as
⟨D1 (u, f ), D1 (u, f )⟩V .. ∗ c = . ⟨D1 (u, f ), DN (u, f )⟩V
··· .. . ···
−1
⟨DN (u, f ), D1 (u, f )⟩V .. . ⟨DN (u, f ), DN (u, f )⟩V
⟨D0 (u, f ), D1 (u, f )⟩V .. −1 = G b. . ⟨D0 (u, f ), DN (u, f )⟩V
If the Gram matrix is rank deficient, then the observed pair does not uniquely identify the coefficient vector; additional data or an additional normalization convention is required. This approach also extends to settings in which the residuals live in a general Banach space V, where there is no inner product. We choose linearly independent test functionals ψ1 , . . . , ψr0 ∈ V ∗ with r0 ≥ N and probe the dictionary terms through the pairings ψi (Dj (u, f )). We define the measurement map Ψ : V → Rr0 as
Ψ(g) = ψ1 (g), . . . , ψr0 (g) ,
g ∈ V.
For the observed pair (f, u), the tested dictionary matrix is Gi,j = ψi (Dj (u, f )) for 1 ≤ i ≤ r0 , and 1 ≤ j ≤ N . If G has full column rank, then the tested dictionary terms are linearly independent and the coefficients are uniquely determined from these tests. Conversely, full column rank is equivalent to independence in the residual space provided that Ψ is injective on span{D1 (u, f ), . . . , DN (u, f )}. In practice, we employ a column-pivoted QR factorization GP = QR for numerical stability to detect the numerical rank r of the tested dictionary matrix G, remove linearly dependent columns, and solve the resulting reduced triangular system by back substitution. This is important as the dictionary may contain redundant terms. If G is rank-deficient, the original coefficients are not uniquely identifiable as multiple coefficient vectors yield the same operator. In this case, the algorithm returns one representative coefficient vector for the reduced dictionary. Note that the formulation in Eq. (1) may also include spurious dictionary terms that are not present in the true operator L, and for which the associated coefficients are zero. The procedure is summarized in Algorithm 1. Unlike SINDy and Weak SINDy algorithms (1–4), we do not require any sparsity-promoting regularization to recover the coefficients, as we assume that we are in a noise-free regime to verify the applicability of Theorem 1. In the next section, we employ the weak identification algorithm to recover the coefficients of several analytic systems from a single input-output pair. In these examples, the weak formulation improves numerical stability by avoiding direct numerical differentiation of the observed data, as also observed in (1, 5). Nicolas Boullé, Diana Halikias, Samuel E. Otto, and Alex Townsend
1 of 5
Algorithm 1 Weak Identification of Analytic Systems r0 ∗ Require: Source term f , observed solution u, dictionary {Di }N i=1 , dual space family {ψi }i=1 ⊂ V 1: Compute the r0 × N matrix Gi,j = ψi (Dj (u, f )) ⊤ 2: Compute the right-hand side b = [ψ1 (D0 (u, f )) . . . ψr0 (D0 (u, f ))] 3: Compute a column-pivoted QR factorization GP = QR 4: Determine the numerical rank r from the diagonal of R 5: Keep the independent columns indexed by the first r pivots and remove the dependent columns ⊤ 6: Form the reduced upper-triangular system R1:r,1:r y = Q:,1:r b y 7: Recover the coefficient vector by undoing the permutation: c = P 0 8: return c
2. Numerical examples This section illustrates the recovery of analytic systems from a single input-output pair using Algorithm 1 on the six examples presented in Figure 1 of the main text. A. Allen–Cahn equation. We first consider the steady-state cubic-quintic Allen–Cahn equation modeling phase separation (6–9), which is given in strong form by 1 ∆u − 5u − u3 + u5 = f, u|∂Ω = 0, [2] 4
where Ω = [0, 1]2 , u ∈ U = H 2 (Ω) ∩ H01 (Ω), and f ∈ F = L2 (Ω). We discretize Eq. (2) using piecewise linear continuous Lagrangian finite elements in the Firedrake finite element software (10) on a 100 × 100 mesh. The weak formulation is obtained by multiplying Eq. (2) with a test function ϕ ∈ V = H01 (Ω) and integrating by parts, resulting in the following weak formulation: 1 − ⟨∇u, ∇ϕ⟩ − 5⟨u, ϕ⟩ − ⟨u3 , ϕ⟩ + ⟨u5 , ϕ⟩ = ⟨f, ϕ⟩. 4 where the boundary term vanishes because the test function has zero trace on the boundary of Ω. We solve the resulting equation using a random source term f ∼ GP(0, K), with squared-exponential covariance kernel K(x, y) = σ 2 exp(−∥x − y∥2 /2ℓ2 ), with variance σ 2 = 10 and length scale ℓ = 0.05. To sample such a forcing term in the finite-dimensional discretization, we evaluate the kernel at the finite element nodes, form the covariance matrix Kh , compute a Cholesky factorization Kh = LL⊤ , and set fh = Lξ with ξ ∼ N (0, I). The nonlinear system is then solved with Newton’s method using backtracking line search with absolute tolerance of 10−12 , while the linearized system is solved using a sparse LU direct solver (11). 0 To apply Algorithm 1, we choose test functionals from the finite-element discretization as follows: let {φi }ri=1 be the nodal basis of the finite element space, and define ψi : w 7→ ⟨w, φi ⟩. On a 100 × 100 element mesh there are 101 nodes per coordinate direction, hence r0 = 1012 nodal basis functions. We use the dictionary {uxx , uxy , uyx , uyy , 1, u, u2 , u3 , u4 , u5 } consisting of second-order derivatives and monomials up to degree five with right-hand side D0 (u, f ) = f . Note that this dictionary contains both spurious terms and the redundant mixed derivative terms uxy = uyx , which are removed by the QR factorization in Algorithm 1. We then assemble the weak forms against all basis functions, leading to a system of r0 = 1012 equations for the N = 10 coefficients in the dictionary, and solve the system using a column-pivoted QR factorization. The recovered coefficients match the ground truth ones with maximum absolute error of 1.5 × 10−12 . B. Navier–Stokes equations. In the second example, we consider the incompressible Navier–Stokes equations for the steady
lid-driven cavity flow on the unit square Ω = [0, 1]2 :
−∇ · (2νϵ(u)) + (u · ∇)u + ∇p = 0,
∇ · u = 0,
[3]
where ϵ(u) = 12 (∇u + ∇u⊤ ) is the strain rate tensor and ν = 10−3 is kinematic viscosity. The boundary conditions are given by g = (1, 0) on the top boundary and g = (0, 0) on the other three boundaries, and we seek a velocity-pressure pair (u, p) ∈ Vg × L20 (Ω), where Vg = {v ∈ H 1 (Ω; R2 ) | v|∂Ω = g}, as the pressure is only defined up to a constant. We solve Eq. (3) using a Taylor–Hood finite element discretization (12) on a 100 × 100 mesh, i.e., using vector-valued continuous piecewise quadratic Lagrangian elements for the velocity and piecewise linear Lagrangian elements for the pressure, implemented in the Firedrake finite element software (10) with the same solver as in Section A. Given a known pressure term p ∈ L20 (Ω), we aim to identify the momentum and continuity equations in Eq. (3) from the observed velocity u ∈ Vg . The weak formulation of Eq. (3) becomes ⟨2νϵ(u), ϵ(ϕu )⟩ + ⟨(u · ∇)u, ϕu ⟩ = ⟨p, ∇ · ϕu ⟩, ⟨∇ · u, ϕp ⟩ = 0,
2 of 5
[4a] [4b]
Nicolas Boullé, Diana Halikias, Samuel E. Otto, and Alex Townsend
for all test functions (ϕu , ϕp ) ∈ V0 × L20 (Ω). The Dirichlet conditions are imposed strongly through the choice of the finite element nodal basis. Although the boundary function g is discontinuous at the two top corners, the mismatch occurs only at isolated points (measure zero), so it does not alter the weak problem. We first focus on the momentum equation Eq. (4a), and define the dictionary of operators D = {Di }N i=1 , where NR = 16, as eightR second-order spatial derivatives terms capturing viscous stress and eight nonlinear advection terms: R { Ω ∂j ui ∂k ϕui dx, Ω uℓ ∂j ui ϕui dx} for i, j, k, ℓ ∈ {x, y}, with right-hand side D0 (u, p) = Ω p∇ · ϕu dx. We assemble the weak forms against all nodal basis functions of the finite element discretization, leading to a system of size 80802 × 16, which is solved by a column-pivoted QR factorization. Here, the system is rank-deficient due to the presence of the redundant mixed derivative terms in the dictionary: ⟨∂x ux ∂y ϕux ⟩ = ⟨∂y ux ∂x ϕux ⟩ and ⟨∂x uy ∂y ϕuy ⟩ = ⟨∂y uy ∂x ϕuy ⟩, which are removed by the QR factorization. The recovered coefficients match the ground truth ones with maximum absolute error of 3.8 × 10−15 . Note that this is expected since the test functions are chosen to be the basis functions of the finite element space used to discretize Eq. (4). Hence, the discrete solution satisfies the linear system corresponding to the finite element discretization of the partial differential equation exactly, up to the tolerance of the numerical solver. The continuity equation Eq. (4b) satisfies a homogeneous equation with zero right-hand side so its coefficients can only be recovered up to a multiplicative constant. We define the dictionary of operators D = {Di }N i=1 , where N = 4, as four first-order R spatial derivatives terms capturing the divergence of the velocity field: { Ω ∂j ui ϕp dx} for i, j ∈ {x, y}, with right-hand side D0 (u, p) = 0. We assemble the weak forms using piecewise linear continuous Lagrangian elements for the test function ϕp , corresponding to the dual space of our finite element discretization. This leads to a system Gc = 0, where G has size 10201 × 4 and has a non-trivial null space. We aim to find a basis for the one-dimensional null space of G and recover the coefficients up to a multiplicative constant by solving the following optimization problem: min ∥Gc∥22
c∈R4
s.t.
∥c∥2 = 1.
This is equivalent to finding the right singular vector of G corresponding to the smallest singular value, which can be computed efficiently using the singular value decomposition (SVD). After normalizing by the coefficient of ∂x ux , we recover the continuity equation ∂x ux + ∂y uy = 0 with maximum absolute error of 6.7 × 10−16 .
C. Lorenz system. As a first ordinary differential equation example, we consider the Lorenz system (13), which is a three-
dimensional system of ordinary differential equations given by ẋ = σ(y − x),
ẏ = x(ρ − z) − y,
ż = xy − βz,
[5]
where σ = 10, ρ = 28, β = 8/3 so that the system features chaotic dynamics. We numerically integrate Eq. (5) from the initial condition [1, 1, 1]⊤ using an explicit fifth order Runge–Kutta method implemented in the SciPy library (14, 15), and sample the trajectory at equispaced times in the interval [0, 40] with timestep ∆t = 10−3 . We then build a dictionary consisting of all monomials up to degree two in (x, y, z), along each component of the system, and a right-hand side D0 (x, y, z) = [ẋ, ẏ, ż]⊤ , so that the system decouples into three independent linear systems, one for each component of the Lorenz system. The weak formulation of Eq. (5) is obtained by multiplying the system with a test function ϕ ∈ V = H 1 ([0, 40]) and integrating by parts the right-hand side, resulting in the following weak formulation for the different components: ⟨c1 + c2 x + c3 y + c4 z + c5 x2 + c6 xy + c7 xz + c8 y 2 + c9 yz + c10 z 2 , ϕ⟩ = −⟨x, ϕ̇⟩, ⟨c1 + c2 x + c3 y + c4 z + c5 x2 + c6 xy + c7 xz + c8 y 2 + c9 yz + c10 z 2 , ϕ⟩ = −⟨y, ϕ̇⟩, ⟨c1 + c2 x + c3 y + c4 z + c5 x2 + c6 xy + c7 xz + c8 y 2 + c9 yz + c10 z 2 , ϕ⟩ = −⟨z, ϕ̇⟩.
Here, we consider sinusoidal test functions (16) of the form ϕ : t 7→ sin(kπt/40) for 1 ≤ k ≤ 50, and integrate the weak forms using the composite Simpson’s rule. This leads to a system of size 50 × 10 for each component of the Lorenz system, which is solved by a column-pivoted QR factorization. The recovered coefficients match the ground truth ones with maximum absolute error of 5.2 × 10−9 . D. Duffing oscillator. Next, we consider the forced Duffing equation modeling damped and driven oscillators, written in
first-order form as
u̇ = v, v̇ = −δv − αu − βu3 + γ cos(ωt), [6] with parameters (α, β, δ, γ, ω) = (1, 5, 0.02, 8, 0.5). We integrate Eq. (6) from the initial condition [u0 , 0]⊤ , where u0 ∼ N (0, 1), using an explicit fifth order Runge–Kutta method implemented in the SciPy library (14, 15), and store the solution (u, v) at equispaced time points on the interval [0, 100] with timestep ∆t = 10−3 . We aim to recover the second order equation ü + δ u̇ + αu + βu3 = γ cos(ωt) from the observed trajectory (u, v), where v = u̇. To do so, we select a dictionary containing all monomials up to degree three in (u, v), along with the high order term ü, and the right-hand side D0 = cos(ω·). We formulate the weak form of the system by multiplying the equation with sinusoidal test functions ϕk : t 7→ sin(kπt/100) for 1 ≤ k ≤ 50, similarly to the Lorenz system in Section C, as 1 − γ
Z 100 0
1 v, ϕ̇k dt + γ
Z 100 0
3
δv + αu + βu , ϕk dt =
Z 100 0
cos(ωt), ϕk dt
We discretize the integrals using Simpson’s composite rule, and obtain a linear system in the dictionary coefficients. After solving the resulting linear system and normalizing by the coefficient of ü, we recover the second order equation with maximum absolute error of 1.1 × 10−10 . Nicolas Boullé, Diana Halikias, Samuel E. Otto, and Alex Townsend
3 of 5
E. Hankel and circulant matrices. Finally, we consider the problem of recovering Hankel and circulant n × n matrices A(c∗ )
R L from right and left matrix-vector products: {(fi , ui = A(c∗ )fi )}qi=1 and {(gi , vi = A(c∗ )⊤ gi )}qi=1 . These matrices can be viewed PN as linearly parameterized operators A(c) = i=1 ci Ai (with N = 2n − 1 for Hankel matrices and N = n for circulant matrices), which we formulate as
Ai
X ci i=1 N
..
.
0 Ai A⊤ i
0
|
{z
..
Di (u,f )
.
Ai f1 u1 f1 .. .. .. . . . N X fqR Ai fq uq ci ⊤ R = R . g = 1 i=1 Ai g1 v1 . . . .. .. .. ⊤ vqL gqL Ai A⊤ i gqL
}
|
{z
Di (u,f )
}
[7]
| {z } D0 (u,f )
Note that this formulation gives us a principled way to generalize the recovery of linearly parameterized operators to multiple evaluations of the operator on different inputs. We generate random Hankel and circulant matrices of size n = 1000 by sampling their coefficients from a standard normal distribution. We then sample one random vector f1 ∼ N (0, In ) and compute the outputs u1 = A(c∗ )f1 , and apply a column-pivoted QR linear solver to the system resulting from Eq. (7).F or circulant matrices, the resulting system is full rank with probability one for a Gaussian input vector, leading to an error of 1.6 × 10−13 (measured in the ℓ2 -norm) in the recovered coefficients. For Hankel matrices, one right matrix-vector product gives at most n independent equations for 2n − 1 unknowns, so the system is rank deficient. Thus, we are in the case ii) of Theorem 1 in the main text, and we cannot recover the coefficients from a single input-output pair. To overcome this issue, we apply multi-shot recovery and sample two random vectors f1 , f2 ∼ N (0, In ) and compute the outputs u1 = A(c∗ )f1 and u2 = A(c∗ )f2 . We then assemble the system as in Eq. (7), noting that Theorem 1 applies with the operator A(c) replaced by the block A(c) 0 operator Ã(c) = , and the input-output pairs (f1 , u1 ) and (f2 , u2 ), which leads to a full rank system, with 0 A(c) coefficient recovery error of 3.5 × 10−11 . We then recovery the result reported in (17, Tab. 1), where Hankel matrices require two input-output pairs to be recovered, while circulant matrices only require one input-output pair.
3. Zero level set of analytic functions For finite dimensional input spaces, one can visualize the inputs for which the parameter-to-recovery map Lf,u : RN → V PN defined as Lf,u (a) = i=1 ai Di (u, f ), for a ∈ RN , is not injective. To do so, we select N linearly independent test functionals ψ1 , . . . , ψN ∈ V ∗ and consider the N × N matrix G with entries Gi,j = ψi (Dj (u, f )). The injectivity of Lf,u is equivalent to the invertibility of G, which is characterized by the non-vanishing of the determinant det(G). Hence, the degenerate set of inputs is the zero level set of the analytic function ϕ : F → R, defined as ϕ(f ) = det(G), which is a thin set in the input space by Lemma 1 of the main text. In this section, we illustrate this phenomenon with the examples plotted in Figure 2 of the main text. We begin with a two-dimensional example, where the input and residual spaces are R2 , with dictionary functions D1 (f ) = [cos(4f02 ), sin(2f1 )]⊤ and D2 (f ) = [sin(2f1 + 1), cos(2f0 )]⊤ for f = [f0 , f1 ]⊤ ∈ R2 . We select canonical test functionals in R2 ⊤ 2 (e⊤ 1 and e2 ) such that the resulting analytic function ϕ : R → R is given by ϕ(f ) = det
cos((2f0 )2 ) sin(2f1 )
sin(2f1 + 1) = cos((2f0 )2 ) cos(2f0 ) − sin(2f1 ) sin(2f1 + 1), cos(2f0 )
f = (f0 , f1 ) ∈ [−2, 2]2 .
In Figure 2A of the main text, we visualize the contour of ϕ for inputs f ∈ [−2, 2]2 , and highlight the zero level set ϕ(f ) = 0 in black. Lemma 1 states that the intersection of this set with a horizontal slice f0 = −0.3 (plotted as dark dots) has Lebesgue measure zero in the slice, which is confirmed by the numerical experiments. Note that one could also consider the intersection with the two-dimensional plane [−2, 2]2 , which yields the dark zero level set curves in the figure that are of measure zero in the two-dimensional input space Next, we consider a three-dimensional example with the dictionary functions D1 (f ) = [1 − 0.2f22 , f1 , 0]⊤ ,
D2 (f ) = [2f1 , sin(3f03 ) cos2 (f2 /2)(3 − f2 ), 0]⊤ ,
D3 (f ) = [0, 0, 1]⊤ ,
for f = [f0 , f1 , f2 ]⊤ ∈ R3 . Selecting canonical test functionals in R3 , we obtain the analytic function ϕ : R3 → R given by
"
1 − 0.2f22 f1 ϕ(f ) = det 0
2f1 sin(3f03 ) cos2 (f2 /2)(3 − f2 ) 0
#
0 0 = 1 − 0.2f22 sin(3f03 ) cos2 (f2 /2)(3 − f2 ) − 2f12 , 1
where we select f ∈ [−2, 2]3 . The zero level set of ϕ is visualized in Figure 2B of the main text as dark shaded surfaces, and we illustrate the intersection of this set with the plane f2 = 0 that has Lebesgue measure zero in this slice. 4 of 5
Nicolas Boullé, Diana Halikias, Samuel E. Otto, and Alex Townsend
References 1. DA Messenger, DM Bortz, Weak SINDy for partial differential equations. J. Comput. Phys. 443, 110525 (2021). 2. DA Messenger, DM Bortz, Weak SINDy: Galerkin-based data-driven model selection. Multiscale Model. & Simul. 19, 1474–1497 (2021). 3. SL Brunton, JL Proctor, JN Kutz, Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proc. Natl. Acad. Sci. USA 113, 3932–3937 (2016). 4. K Champion, B Lusch, JN Kutz, SL Brunton, Data-driven discovery of coordinates and governing equations. Proc. Natl. Acad. Sci. USA 116, 22445–22451 (2019). 5. H Schaeffer, SG McCalla, Sparse model selection via integral terms. Phys. Rev. E 96 (2017). 6. C Kuehn, Numerical Continuation and SPDE Stability for the 2D Cubic-Quintic Allen–Cahn Equation. SIAM/ASA J. Uncertain. Quantification 3, 762–789 (2015). 7. W van Saarloos, PC Hohenberg, Fronts, pulses, sources and sinks in generalized complex Ginzburg-Landau equations. Phys. D: Nonlinear Phenom. 56, 303–367 (1992). 8. T Kapitula, B Sandstede, Instability mechanism for bright solitary-wave solutions to the cubic–quintic Ginzburg–Landau equation. J. Opt. Soc. Am. B 15, 2757–2762 (1998). 9. RJ Deissler, HR Brand, Periodic, quasiperiodic, and chaotic localized solutions of the quintic complex Ginzburg-Landau equation. Phys. Rev. Lett. 72, 478 (1994). 10. DA Ham, et al., Firedrake User Manual (Imperial College London and University of Oxford and Baylor University and University of Washington), 1st edition (2023). 11. S Balay, et al., PETSc/TAO users manual, (Argonne National Laboratory), Technical Report ANL-21/39 - Revision 3.25 (2026). 12. C Taylor, P Hood, A numerical solution of the Navier-Stokes equations using the finite element technique. Comput. & Fluids 1, 73–100 (1973). 13. E Lorenz, Deterministic nonperiodic flow. J. Atmos. Sci. 20, 130–141 (1963). 14. P Virtanen, et al., SciPy 1.0: fundamental algorithms for scientific computing in Python. Nat. Methods 17, 261–272 (2020). 15. JR Dormand, PJ Prince, A family of embedded Runge-Kutta formulae. J. Comput. Appl. Math. 6, 19–26 (1980). 16. Z Chen, U Fasel, A Bizyaeva, Fourier Weak SINDy: Spectral Test Function Selection for Robust Model Identification. arXiv preprint arXiv:2604.20141 (2026). 17. D Halikias, A Townsend, Structured matrix recovery from matrix-vector products. Numer. Linear Algebr. Appl. 31, e2531 (2024).
Nicolas Boullé, Diana Halikias, Samuel E. Otto, and Alex Townsend
5 of 5