REVIEW 4 major objections 5 minor 21 references
Sparse Discovery of Functional Relationships in Solutions to Systems of Differential Equations
T0 review · 4 major / 5 minor · reviewed 2026-08-06 · deepseek-v4-flash
Pith's one-line read By solving a least-squares null-space problem on a library of candidate functions sampled along a numerical trajectory, the paper claims that sparse relations among ODE solution components can be discovered automatically.
desk verdict A nullspace computation presented as sparse discovery; the enzyme example works, but the central glycolysis claim is an unvalidated rank-threshold artifact. read the letter →
The pith
A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.
The reading
What carries the argument
The central object is the augmented library $\Theta(x)$ and its grid evaluation $\Theta_{\mathrm{grid}}$. The algorithm solves the normal equations $\Theta_{\mathrm{grid}}^T\Theta_{\mathrm{grid}}\xi=0$, repeatedly reduces $\Theta_{\mathrm{grid}}^T\Theta_{\mathrm{grid}}$ to reduced row echelon form, and removes columns whose rows have exactly one nonzero pivot, since those coefficients must vanish; the surviving columns span the space of candidate relations. This turns the search for conserved quantities into a null-space computation, with sparsity enforced by deleting provably zero coordinates rather than by a penalty term.
What would settle it
Take the enzyme system, add a small extra decay term $-\varepsilon P$ to $dP/dt$ so no exact invariant survives, and run the algorithm on a coarse grid with a loose row-reduction tolerance; if it still returns the unperturbed relation $-2C_1 - 3E + S + P = 0$, the method has reported numerical error as a conserved quantity.
Extended reading notes
Core claim
The central claim is that one can recover nontrivial constraints $F(x(t))=0$ from the numerical solution alone by building the augmented matrix $\Theta_{\mathrm{grid}}$ of candidate functions and solving $\Theta_{\mathrm{grid}}\xi=0$ in the least-squares sense. Iteratively computing the reduced row echelon form of $\Theta_{\mathrm{grid}}^T\Theta_{\mathrm{grid}}$ and discarding columns whose pivots isolate a single component eliminates all coefficients that must be zero; the remaining least-squares null vector gives the sparse relation. The enzyme example recovers invariant linear combinations such as $-2C_1 - 3E + S + P = 0$ and, with a richer library, a two-parameter family of linear constraints; the glycolysis example yields a relation among the 36 monomials of degree at most two in seven species. The author states that the results are least-squares in character, so the $\xi$-coefficients are the best available relation for the chosen library.
Load-bearing premise
The load-bearing premise is that the numerical trajectory is accurate enough that treating the residual $\eta$ as zero is valid, so the null-space computation of $\Theta_{\mathrm{grid}}^T\Theta_{\mathrm{grid}}$ with its implicit pivot threshold decides which coefficients are truly zero; if the grid is coarse or the threshold is poorly chosen, spurious relations will appear as genuine conserved quantities.
Editorial extensions
If this is right
- For the enzyme kinetics system, the method recovers invariant linear combinations that hold for all time, including $-2C_1 - 3E + S + P = 0$, and changing initial data changes only the constant term.
- Richer candidate libraries produce more relations: using $[1, x, x^2]$ returns a multi-parameter family of quadratic constraints, while $[1, x]$ returns a smaller family of linear ones.
- The same algorithm applied to the seven-species glycolytic oscillator finds a nontrivial degree-two monomial relation, while simpler libraries find none.
- Because the output is a least-squares null vector, the coefficients are the best available relation for the selected library and trajectory grid.
Reading between the lines
- A direct symbolic check of any discovered relation is to differentiate $F$ along the numerical solution; the method implicitly predicts the derivative is zero, so testing the residual $|dF/dt|$ on an independent fine grid is a cheap validation.
- The same machinery could be run in reverse: take a measured trajectory, treat it as data, and test whether a sparse relation exists; the main obstacle is that measurement noise makes the $\eta=0$ assumption fail.
- For systems with parameters, the discovered coefficient family may encode parameter-free invariants, as seen when different initial data change only $\xi_1$; this suggests a route to detect symmetries or conserved quantities shared across parameter values.
- A stringent test would be to add a small perturbation that provably destroys the invariant and see whether the algorithm's row-reduction threshold still reports it; this would reveal the sensitivity of the pivot-removal step to solver tolerance.
Editorial analysis
A structured set of objections, weighed in public.
Referee Report
Summary. The paper proposes an algorithm to discover functional relationships F(x(t))=0 among the components of solutions to systems of ODEs. The method constructs a library of candidate functions, samples it on a numerical solution grid, and then seeks a coefficient vector xi such that the sampled library matrix times xi vanishes. The algorithm iteratively discards columns whose coefficients are identified as zero by the reduced row echelon form of the normal matrix, and finally solves a homogeneous least-squares problem. Two examples are presented: an enzyme kinetics model and a glycolytic oscillator model, with MATLAB code provided on GitHub. The paper claims this is a novel sparse-identification framework for discovering conserved quantities and invariants from numerical solution data.
Significance. If the proposed method were correct, it would provide a useful tool for discovering conserved quantities and functional constraints from numerical solutions of ODE systems, with potential applications in model reduction and analysis of large nonlinear systems. The availability of a public MATLAB implementation and the enzyme kinetics example, where the method recovers known linear conservation laws, are positive aspects. However, the central claim that the algorithm discovers genuine functional relationships is not established: the least-squares formulation is mathematically misstated, the rank-reduction step lacks tolerances and a stopping criterion, the eta=0 assumption is inconsistent with the use of numerical solutions, and the reported relations, especially the quadratic invariant in Section 4.2, are not validated on independent data. As a result, the outputs are indistinguishable from numerical artifacts of rank-thresholding.
major comments (4)
- [Section 3, Eqs. (7)-(8)] The equivalence between the least-squares problem (7) and the homogeneous system (8) is mathematically incorrect. If eta is a fixed data vector, the normal equations for (7) are Theta^T Theta xi = Theta^T eta, not Theta^T Theta xi = 0. If eta is instead intended to represent the residual to be minimized in Theta xi = 0, then the problem is min ||Theta xi||, whose normal equations are indeed Theta^T Theta xi = 0, but this is not a least-squares solution of (7) with eta as given data. This conflation is load-bearing because the algorithm's output is defined as the solution set of (8).
- [Section 3, Step 3] The rank-reduction procedure relies on computing the reduced row echelon form of Theta^T Theta without a stated numerical tolerance or a stopping criterion; the number of passes p is an unexplained free parameter. For a generic trajectory of a nonlinear ODE, the library columns are linearly independent, so Theta^T Theta is positive definite and its RREF is the identity up to roundoff. The algorithm therefore implicitly decides rank by an unspecified threshold, and the output (e.g., the '25 basic and 11 free' solution reported in Section 4.2) is an artifact of that threshold rather than a demonstrated property of the trajectory. No analysis is provided of how the output depends on the grid size, the solver tolerance, or the rank tolerance.
- [Section 3, Step 4] The assumption eta=0 is inconsistent with Step 1's use of the numerical solver ode45: discretization and roundoff errors are not zero, and evaluating the candidate library on the approximate solution introduces additional error. For a true invariant F(x(t))=0, the sampled vector Theta_grid xi will not be exactly zero. The algorithm's exact homogeneous solve with the rank-reduced matrix must nevertheless produce a nontrivial null vector whenever the rank test declares a deficiency. No residual analysis, error bounds, or check on independent time points is provided, so the reported relation is not shown to hold beyond the fitted grid points.
- [Section 4.2] The claimed quadratic invariant for the glycolysis model is not validated. The paper does not demonstrate that the relation persists under changes of the time grid, the integration tolerance, or the initial condition; it does not compare the fitted relation against a separately computed trajectory; and it does not quantify the residual. In the absence of such validation, the reported relation is indistinguishable from a null vector of a numerically rank-deficient matrix, which is expected generically for a seven-dimensional trajectory and a 36-column quadratic library sampled at finitely many points.
minor comments (5)
- [Abstract and Section 2] The abstract and introduction emphasize that the method is 'sparse', but the algorithm in Section 3 does not use any sparsity-promoting regularization (such as L1 minimization or sequential thresholding); it computes exact null vectors of a homogeneous system. The role of sparsity in the algorithm should be clarified or the terminology revised.
- [Section 4.1] There is a typo in the sentence 'Figure 1 illustrates the preservation of contraint (14)'; 'contraint' should be 'constraint'.
- [Section 3, Step 3] The reuse of the notation Theta_grid for the reduced matrix after discarding columns is confusing; a different symbol or a superscript, such as Theta_grid^{(r)}, would improve readability.
- [Section 4.2] The sentence 'It is important to remember that the results are in the least-squares sense' is vague; the paper should specify what quantity is being minimized and with respect to which norm, especially because the algorithm solves a homogeneous system rather than a standard least-squares problem.
- [References] The cited MATLAB implementation on GitHub [18] should include a version identifier or commit hash, and any required MATLAB toolboxes, to support reproducibility.
Circularity Check
The central 'discovery' step reduces by construction: ξ is solved from Θ_grid ξ = 0 on the sampled trajectory, and the preservation figures evaluate the same rows; the glycolysis 25/11 relation is a fitted null vector, not an independently verified invariant.
-
fitted input called prediction
[Section 3, Step 4 (Eq. (8)); Section 4.2 and Figure 4]
"Preserving the same notation, ξ, for the unknown, the set of least-squares solutions of (7) coincides with the nonempty set of solutions of the linear system Θ T grid Θ grid ξ = 0, (8) ... The non-trivial solution set {ξ i1, ξ i2, . . . , ξ ik } of the system (8) is what we are looking for. ... Figure 4 illustrates the preservation of the above constraint corresponding to the initial data (23) for values of ξ specified in the caption."
The output coefficients are defined as solutions of Θ_grid ξ = 0 on the very same sampled rows used to build Θ_grid, so Figure 4's 'preservation' is the defining equation evaluated on its own input; the residual is zero by construction up to solver tolerance. For the glycolysis example no independent validation is given — no derivative check, no held-out trajectory, no residual on a different grid — so the reported 25-basic/11-free quadratic relation is the nullspace of the sampled library matrix after an unspecified rank threshold, i.e., a fitted vector presented as a discovered invariant.
full rationale
The paper's self-citation of [4] and [18] is not load-bearing: the algorithm is explicit and the method is not justified by an unverified theorem from the author. However, the central validation loop is circular in the internal sense: Step 4 defines the relation by solving Θ_grid ξ = 0 (with η set to zero) on the numerical trajectory, and the 'preservation' figures check the same trajectory used to construct the matrix. The enzyme linear invariant is genuinely supported because the paper verifies d/dt(ξ2 C1 + ξ3 E + ξ4 S + ξ5 P) = 0 analytically, so that part has independent content. The quadratic glycolysis relation, by contrast, is presented only as a null vector of the sampled 36-column library and its Figure checks the same data; its status as an ODE invariant is not established. This is partial circularity: the fitting step and the claimed verification coincide, though not every claimed result is tautological.
Assumptions & free parameters
free parameters (2)
- Number of iteration passes p
- RREF numerical tolerance
assumptions (3)
- domain assumption The numerical solution x_grid is an exact representation of the true solution (η=0 assumption)
- domain assumption The candidate library Θ(x) is chosen so that any true conserved quantity lies in its span
- domain assumption Null vectors of the data matrix correspond to genuine algebraic relations on the trajectory, not empirical fits
Cite this review
Pith. "Pith review of Sparse Discovery of Functional Relationships in Solutions to Systems of Differential Equations." pith.science (2026). https://pith.science/paper/656FGF7K
@misc{pith2026250715805,
author = {Pith},
title = {Pith review of: Sparse Discovery of Functional Relationships in Solutions to Systems of Differential Equations},
year = {2026},
howpublished = {\url{https://pith.science/paper/656FGF7K}},
note = {Machine review of arXiv:2507.15805}
}
read the original abstract
This work develops a framework to discover relations between the components of the solution to a given initial-value problem for a first-order system of ordinary differential equations. This is done by using sparse identification techniques on the data represented by the numerical solution of the initial-value problem at hand. The only assumption is that there are only a few terms that connects the components, so that the mathematical relations to be discovered are sparse in the set of possible functions. We illustrate the method through examples of applications.
Figures
Figures from the paper (1 more)
Reference graph
Works this paper leans on
-
[1]
Ayankoso, S., Olejnik, P. Time-Series Machine Learning Techniques for Modeling and Identification of Mecha- tronic Systems with Friction: A Review and Real Application,Electronics, Vol. 12, No. 17, 2023, 3669
work page 2023
-
[2]
Beetham, S., Capecelatro, J. Formulating turbulence closures using sparse regression with embedded form invariance,Physical Review Fluids, Vol. 5, No. 8, 2020, 084611
work page 2020
-
[3]
Sparse learning of stochastic dynamical equations,The Journal of Chemical Physics, Vol
Boninsegna, L., N¨ uske, F., Clementi, C. Sparse learning of stochastic dynamical equations,The Journal of Chemical Physics, Vol. 148, No. 24, 2018, 241723
work page 2018
-
[4]
Brunton, S.L., Proctor, J.L., Kutz, J.N. Discovering governing equations from data by sparse identification of nonlinear dynamical systems,Proceedings of the National Academy of Sciences, Vol. 113, No. 15, 2016, pp. 3932-3937
work page 2016
-
[5]
Champion, K., Lusch, B., Kutz, J.N., Brunton, S.L. Data-driven discovery of coordinates and governing equations,Proceedings of the National Academy of Sciences, Vol. 116, No. 45, 2019, pp. 22445-22451
work page 2019
-
[6]
Daniels, B.C., Nemenman, I. Efficient Inference of Parsimonious Phenomenological Models of Cellular Dy- namics Using S-Systems and Alternating Regression,PloS ONE, Vol. 10, No. 3, 2015, e0119821
work page 2015
-
[7]
A Robust SINDy Approach by Combining Neural Networks and an Integral Form, 2023, arXiv:2309.07193
Forootani, A., Goyal, P., Benner, P. A Robust SINDy Approach by Combining Neural Networks and an Integral Form, 2023, arXiv:2309.07193
arXiv 2023
-
[8]
Cluster-based network model identification and control of nonlinear systems.Nonlinear Dynamics, Vol
Hoffmann, M., Nageshrao, S.P., Haller, G. Cluster-based network model identification and control of nonlinear systems.Nonlinear Dynamics, Vol. 98, 2019, pp. 1537-1556
work page 2019
Show all 21 references
-
[9]
87, 1, 2021
Kang, S.H., Liao, W., Liu, Y., IDENT: Identifying Differential Equations with Numerical Time Evolution, Journal of Scientific ComputingVol. 87, 1, 2021
2021
-
[10]
Constrained sparse Galerkin regression,Journal of Fluid Mechanics, Vol
Loiseau, J.-C., Brunton, S.L. Constrained sparse Galerkin regression,Journal of Fluid Mechanics, Vol. 838, 2018, pp. 42-67
2018
-
[11]
Sparse reduced-order modeling: Sensor-based dynamics to full- state estimation.Journal of Fluid Mechanics, Vol
Loiseau, J.-C., Noack, B.R., Brunton, S.L. Sparse reduced-order modeling: Sensor-based dynamics to full- state estimation.Journal of Fluid Mechanics, Vol. 944, 2020, A36. NICOLAE TARFULEA: FUNCTIONAL RELATIONSHIPS IN SOLUTIONS TO ODE SYSTEMS 11
2020
-
[12]
Mathematical Biology
Friedman, A. Mathematical Biology. Modeling and Analysis,CBMS Reg. Conf. Ser. Math., No. 127, 2018, viii+100 pp
2018
-
[13]
Inferring biological networks by sparse identification of nonlinear dynamics,IEEE Transactions on Molecular, Biological, and Multi-Scale Communications, Vol
Mangan, N.M., Brunton, S.L., Proctor, J.L., Kutz, J.N. Inferring biological networks by sparse identification of nonlinear dynamics,IEEE Transactions on Molecular, Biological, and Multi-Scale Communications, Vol. 2, No. 1, 2017, pp. 52-63
2017
-
[14]
Data-driven discovery of partial differential equations, Science Advances, Vol
Rudy, S.H., Brunton, S.L., Proctor, J.L., Kutz, J.N. Data-driven discovery of partial differential equations, Science Advances, Vol. 3, No. 4, 2017, e1602614
2017
-
[15]
Temperature dependency and temperature compensation in a model of yeast glycolytic oscillations,Biophys Chem, Vol
Ruoff, P., Christensen, M., Wolf, J., Heinrich R. Temperature dependency and temperature compensation in a model of yeast glycolytic oscillations,Biophys Chem, Vol. 106, No. 179, 2003, pmid:14556906
2003
-
[16]
Learning dynamical systems and bifurcation via group sparsity,SIAM Journal on Applied Mathematics, Vol
Schaeffer, H., Tran, G., Ward, R. Learning dynamical systems and bifurcation via group sparsity,SIAM Journal on Applied Mathematics, Vol. 78, No. 1, 2017, pp. 327-347
2017
-
[17]
Sparse identification of truncation errors,Journal of Computational Physics, Vol
Schmelzer, M., Dwight, R.P., Noack, B.R. Sparse identification of truncation errors,Journal of Computational Physics, Vol. 416, 2020, 109517
2020
-
[18]
Sparse Discovery of Functional Relationships in ODE Systems
Tarfulea, N. Sparse Discovery of Functional Relationships in ODE Systems. GitHub, 2026, https://github.com/NicolaeTarfulea/Sparse-Discovery-of-Functional-Relationships-in-Solutions-to-Systems- of-Differential-Equations
2026
-
[19]
On the Convergence of the SINDy Algorithm,SIAM Journal on Scientific Computing, Vol
Zhang, S., Schaeffer, H. On the Convergence of the SINDy Algorithm,SIAM Journal on Scientific Computing, Vol. 41, No. 5, 2019, A3027-A3046
2019
-
[20]
SINDy-RL: Interpretable and Efficient Model-Based Rein- forcement Learning, 2024, arXiv:2403.09110
Zolman, N., Kutz, J.N., Fasel, U., Brunton S.L. SINDy-RL: Interpretable and Efficient Model-Based Rein- forcement Learning, 2024, arXiv:2403.09110
2024 arXiv
-
[21]
https://www.mathworks.com/help/matlab/ref/ode45.html
Reviewed August 6, 2026 · model on record in the stance chip above.
Discussion (0). Continue with ORCID to comment.