Pith. sign in

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 →

arxiv 2507.15805 v2 pith:656FGF7K submitted 2025-07-21 math.OC cs.ITcs.NAmath.ITmath.NA

classification math.OCcs.ITcs.NAmath.ITmath.NA MSC 34-0465-04
keywords ODEsystemsfunctionalrelationshipsconservedquantitiessparseregressionleastsquaresinvariantsnumericalsolutioninitial-valueproblems
verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

The paper proposes that hidden functional relationships among solution components of a first-order ODE initial-value problem can be discovered by sparse regression on a numerical trajectory. It constructs a large library of candidate functions of the solution components, evaluates them on a grid, and solves a least-squares system to find the few nonzero coefficients that make a linear combination vanish along the solution. If correct, the method exposes conserved quantities, invariants, and supporting manifolds that are not obvious from the equations, and it does so with only the assumption that the relation is sparse. The approach is demonstrated on an enzyme kinetics model and a seven-species yeast glycolysis oscillator.

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.

Watch

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

Editorial extensions of the paper, not claims the author makes directly.

  • 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.
Share X Bluesky LinkedIn Reddit HN

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Referee Report

4 major / 5 minor

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)
  1. [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).
  2. [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.
  3. [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.
  4. [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)
  1. [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.
  2. [Section 4.1] There is a typo in the sentence 'Figure 1 illustrates the preservation of contraint (14)'; 'contraint' should be 'constraint'.
  3. [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.
  4. [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.
  5. [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

1 steps flagged · score 6.0 of 10

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.

  1. 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 2 free parameters · 3 assumptions · 0 invented entities

The method rests on user choices (library, iteration count, tolerances) and on the assumption that a numerical trajectory lies exactly on the constraint manifold, so that the nullspace of the data matrix yields genuine invariants rather than empirical fits.

free parameters (2)
  • Number of iteration passes p
    The algorithm repeats the column-elimination step 'p' times in Step 3, but no criterion for p is given; it controls which columns survive and therefore what relations are found.
  • RREF numerical tolerance
    No tolerance is specified for deciding when a pivot is effectively zero; this determines whether Θ^T Θ is treated as the identity or as rank-deficient.
assumptions (3)
  • domain assumption The numerical solution x_grid is an exact representation of the true solution (η=0 assumption)
    Stated in Step 4: 'For high-order methods, we treat η as negligible, that is, η=0.' This ignores discretization error, which can create spurious near-null vectors.
  • domain assumption The candidate library Θ(x) is chosen so that any true conserved quantity lies in its span
    Section 2 states the choice of candidate functions 'may not be always clear'; the method's success depends on this choice, which is not justified.
  • domain assumption Null vectors of the data matrix correspond to genuine algebraic relations on the trajectory, not empirical fits
    The algorithm treats any nullspace vector as a discovered relationship, but for a 1D trajectory there are many interpolating functions; without independent validation the output may be an artifact.

how reviews work

0 comments
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 reproduced from arXiv: 2507.15805 by the authors.

Figure 1
Figure 1. Example 1: ξ1 = 12, ξ2 = −2, ξ3 = 3, ξ4 = −5, ξ5 = −5 Example 2: ξ1 = −14, ξ2 = 8, ξ3 = 5, ξ4 = 3, ξ5 = 3 Example 3: ξ1 = 3, ξ2 = −1, ξ3 = 0, ξ4 = −1, ξ5 = −1 Example 4: ξ1 = 7, ξ2 = −3, ξ3 = −1, ξ4 = −2, ξ5 = −2 If the initial data is changed to C1(0) = 2, E(0) = 2, S(0) = 1, P(0) = 1, (15) [PITH_FULL_IMAGE:figures/full_fig_p006_1.png] view at source ↗
Figure 2
Figure 2. Example 1: ξ1 = −12, ξ2 = 3, ξ3 = −2, ξ4 = 5, ξ5 = 5 Example 2: ξ1 = −20, ξ2 = 5, ξ3 = 4, ξ4 = 1, ξ5 = 1 Example 3: ξ1 = −8, ξ2 = 2, ξ3 = 1, ξ4 = 1, ξ5 = 1 Example 4: ξ1 = 20, ξ2 = −5, ξ3 = −3, ξ4 = −2, ξ5 = −2 Notice that for the two sets of initial data considered above, only ξ1 is different. It suggests that, for ξ2 = ξ3+ξ5, and ξ4 = ξ5, with ξ3 and ξ5 as free parameters, ξ2 ·C1+ξ3 ·E +ξ4 ·S +ξ5 ·P is constant wi… view at source ↗
Figure 3
Figure 3. Example 1: ξ1 = 27, ξ2 = −18, ξ3 = 2, ξ4 = −25, ξ5 = −2, ξ6 = 9, ξ7 = 5, ξ8 = −5, ξ9 = 4, ξ10 = −1, ξ11 = 6, ξ12 = −1, ξ13 = 3, ξ14 = 3, ξ15 = −3 Example 2: ξ1 = 31, ξ2 = −9, ξ3 = 0, ξ4 = −32, ξ5 = −1, ξ6 = 0, ξ7 = 2, ξ8 = 2, ξ9 = 9, ξ10 = 3, ξ11 = 6, ξ12 = −2, ξ13 = 2, ξ14 = 2, ξ15 = −4 Example 3: ξ1 = −3, ξ2 = 8, ξ3 = −4, ξ4 = 11, ξ5 = 0, ξ6 = 1, ξ7 = 5, ξ8 = −2, ξ9 = −8, ξ10 = 1, ξ11 = −3, ξ12 = −3, ξ13 = 3, ξ14 … view at source ↗
Figures from the paper (1 more)
Figure 4
Figure 4. Figure 4: Example 1: ξ1 = 3.7753, ξ2 = −4.5229, ξ3 = −9.2616, ξ4 = 15.6269, ξ5 = 17.2988, ξ6 = 10.0912, ξ7 = −1.0314, ξ8 = 6.5885, ξ9 = 1.0066, ξ10 = 3.9460, ξ11 = 2.9008, ξ12 = −5.3707, ξ13 = −0.4891, ξ14 = −19.2035, ξ15 = −5.5577, ξ16 = −6.1561, ξ17 = −26.0902, ξ18 = −1.3775, …

Discussion (0). Continue with ORCID to comment.

Reference graph

Works this paper leans on

21 extracted references · 19 canonical work pages

  1. [1]

    Time-Series Machine Learning Techniques for Modeling and Identification of Mecha- tronic Systems with Friction: A Review and Real Application,Electronics, Vol

    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

  2. [2]

    Formulating turbulence closures using sparse regression with embedded form invariance,Physical Review Fluids, Vol

    Beetham, S., Capecelatro, J. Formulating turbulence closures using sparse regression with embedded form invariance,Physical Review Fluids, Vol. 5, No. 8, 2020, 084611

  3. [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

  4. [4]

    Discovering governing equations from data by sparse identification of nonlinear dynamical systems,Proceedings of the National Academy of Sciences, Vol

    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

  5. [5]

    Data-driven discovery of coordinates and governing equations,Proceedings of the National Academy of Sciences, Vol

    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

  6. [6]

    Efficient Inference of Parsimonious Phenomenological Models of Cellular Dy- namics Using S-Systems and Alternating Regression,PloS ONE, Vol

    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

  7. [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

  8. [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

Show all 21 references
  1. [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

  2. [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

  3. [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

  4. [12]

    Mathematical Biology

    Friedman, A. Mathematical Biology. Modeling and Analysis,CBMS Reg. Conf. Ser. Math., No. 127, 2018, viii+100 pp

  5. [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

  6. [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

  7. [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

  8. [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

  9. [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

  10. [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

  11. [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

  12. [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

  13. [21]

    https://www.mathworks.com/help/matlab/ref/ode45.html

Pith tools

Reviewed August 6, 2026 · model on record in the stance chip above.