Additive Atomic Forests for Symbolic Function and Antiderivative Discovery

# Additive Atomic Forests for Symbolic Function and Antiderivative Discovery Generating Primitives, Product-Rule Closure, and Derivative-Matching Optimisation
Source: [https://arxiv.org/html/2605.08130](https://arxiv.org/html/2605.08130)
\(May 2026\)

###### Abstract

We present a framework for the simultaneous symbolic recovery of a function and its antiderivative from data\. The framework rests on three ideas\. First, a*derivative algebra*: the observation that the product ruledd​x​\[f⋅g\]=f′​g\+f​g′\\frac\{d\}\{dx\}\[f\\cdot g\]=f^\{\\prime\}g\+fg^\{\\prime\}and the chain rule, applied to a seed set of elementary functions, generate a self\-expanding system of function–derivative pairs—a living library that grows each time a new function is discovered\. Second, two complementary primitives—EML\(eu−ln⁡v\)\\,\(e^\{u\}\-\\ln v\), which is theoretically complete for all elementary functions, and SOL\(sin⁡u−cos⁡v\)\\,\(\\sin u\-\\cos v\), introduced here, which makes trigonometric atoms available at depth 1 instead of depth∼\\sim8—that seed the library with core atoms cheaply\. Third,*additive atomic forests*: finite sums of primitive trees, optionally composed via multiplicative nodes, whose derivatives are fitted to data by continuous optimisation or by exhaustive search over the library\. Because differentiation of each atom is determined by construction, the forest simultaneously encodes a symbolic expressionFFand its derivativeF′F^\{\\prime\}; no symbolic integration step is required\.

The library is not a fixed object: it self\-constructs from a small seed set by recursive application of the product rule, chain rule, and the two primitives, and it can grow as newly discovered functions are folded back in\. The larger the library, the richer the expressible class of candidate functions\. We give conditional completeness, additive\-depth, and analytic simultaneous\-recovery results for the framework\. Empirically, in our reported runs on 17 classification benchmarks, sparse atom combinations match or exceed XGBoost on 13 datasets while producing interpretable formulas\. On the Feynman symbolic regression benchmark, a self\-built library to depth 3 gives exact recovery on 31% of equations and relative\-MSE below 0\.01 on a further 40%\. We demonstrate the method on real scientific data by proposing candidate radial\-acceleration relations from the SPARC galaxy database\.

###### Contents

1. [1Introduction](https://arxiv.org/html/2605.08130#S1)
2. [2Derivative algebra of elementary functions](https://arxiv.org/html/2605.08130#S2)1. [2\.1The elementary functions](https://arxiv.org/html/2605.08130#S2.SS1) 2. [2\.2Closure under differentiation via the product rule](https://arxiv.org/html/2605.08130#S2.SS2)
3. [3Generating primitives](https://arxiv.org/html/2605.08130#S3)1. [3\.1EML \(Exponential Minus Logarithm\)](https://arxiv.org/html/2605.08130#S3.SS1) 2. [3\.2SOL \(Sine Of, Less cosine\)](https://arxiv.org/html/2605.08130#S3.SS2)
4. [4Canonical forms and parameterised trees](https://arxiv.org/html/2605.08130#S4)
5. [5Additive atomic forests](https://arxiv.org/html/2605.08130#S5)
6. [6Multiplicative nodes and enhanced forests](https://arxiv.org/html/2605.08130#S6)
7. [7The derivative\-matching principle](https://arxiv.org/html/2605.08130#S7)
8. [8Self\-expanding atom libraries](https://arxiv.org/html/2605.08130#S8)1. [8\.1Construction procedure](https://arxiv.org/html/2605.08130#S8.SS1)
9. [9Search algorithms](https://arxiv.org/html/2605.08130#S9)1. [9\.1Precomputation](https://arxiv.org/html/2605.08130#S9.SS1) 2. [9\.2K=1K=1: vectorised scalar regression](https://arxiv.org/html/2605.08130#S9.SS2) 3. [9\.3K=2K=2: vectorised Cramer’s rule](https://arxiv.org/html/2605.08130#S9.SS3) 4. [9\.4K≥3K\\geq 3: beam search](https://arxiv.org/html/2605.08130#S9.SS4) 5. [9\.5Verification](https://arxiv.org/html/2605.08130#S9.SS5)
10. [10Training \(gradient\-based mode\)](https://arxiv.org/html/2605.08130#S10)
11. [11Structural complexity](https://arxiv.org/html/2605.08130#S11)1. [11\.1Connection to the Liouville–Risch decomposition](https://arxiv.org/html/2605.08130#S11.SS1)
12. [12Comparison with existing methods](https://arxiv.org/html/2605.08130#S12)
13. [13Empirical results](https://arxiv.org/html/2605.08130#S13)1. [13\.1Classification: atoms vs\. XGBoost](https://arxiv.org/html/2605.08130#S13.SS1) 2. [13\.2Regression](https://arxiv.org/html/2605.08130#S13.SS2) 3. [13\.3Interpretable formulas and their derivatives](https://arxiv.org/html/2605.08130#S13.SS3) 4. [13\.4Feynman representability test](https://arxiv.org/html/2605.08130#S13.SS4) 5. [13\.5Demonstration: radial acceleration relation](https://arxiv.org/html/2605.08130#S13.SS5) 6. [13\.6Failure modes](https://arxiv.org/html/2605.08130#S13.SS6)
14. [14Extensions](https://arxiv.org/html/2605.08130#S14)1. [14\.1Multivariate canonical forms](https://arxiv.org/html/2605.08130#S14.SS1) 2. [14\.2Conserved quantity discovery](https://arxiv.org/html/2605.08130#S14.SS2)
15. [15Conclusion](https://arxiv.org/html/2605.08130#S15)
16. [References](https://arxiv.org/html/2605.08130#bib)

## 1Introduction

The derivative of a product is

dd​x​\[f​\(x\)⋅g​\(x\)\]=f′​\(x\)​g​\(x\)\+f​\(x\)​g′​\(x\)\.\\frac\{d\}\{dx\}\[f\(x\)\\cdot g\(x\)\]=f^\{\\prime\}\(x\)\\,g\(x\)\+f\(x\)\\,g^\{\\prime\}\(x\)\.\(1\)This identity, together with the chain rule, is the engine of the present work\. If we know the derivatives of a finite collection of elementary functions, the product rule and chain rule let us compute the derivative of any composition, product, or power of those functions—and therefore of any elementary function\. Powers are a special case:fn=f⋅f​⋯​ff^\{n\}=f\\cdot f\\cdots f, sodd​x​\[fn\]=n​fn−1​f′\\frac\{d\}\{dx\}\[f^\{n\}\]=nf^\{n\-1\}f^\{\\prime\}follows from repeated application of \([1](https://arxiv.org/html/2605.08130#S1.E1)\)\. This closure property means that a small set of seed functions, equipped with their derivatives, can bootstrap an arbitrarily large library of function–derivative pairs\.

This observation suggests a strategy for symbolic regression that differs fundamentally from existing approaches\(Cranmer,[2023](https://arxiv.org/html/2605.08130#bib.bib4); Udrescu & Tegmark,[2020](https://arxiv.org/html/2605.08130#bib.bib19); Schmidt & Lipson,[2009](https://arxiv.org/html/2605.08130#bib.bib16)\)\. Rather than searching for a function that fits data, we build a*self\-expanding library*of atoms—each a named elementary function paired with its analytic derivative—and search for a sparse linear combination of*derivatives*that matches the data\. When the match is exact, the corresponding linear combination of atoms is the antiderivative, available immediately in closed form\. No symbolic integration algorithm\(Risch,[1969](https://arxiv.org/html/2605.08130#bib.bib13),[1970](https://arxiv.org/html/2605.08130#bib.bib14)\)is invoked; no numerical integration of data is performed\.

The library is not a static catalogue\. It begins from a seed set \(rational powers, depth\-1 trees of the generating primitives\) and self\-constructs by applying the product rule, chain rule, and compositions at increasing depth\. Each newly discovered function—whether found by the search or supplied by the user—is folded back into the library together with its derivative, expanding the search space for subsequent problems\. The library thus functions as a persistent*knowledge base*: each verified discovery can enlarge the candidate set available for later problems\. Its size is bounded by available memory, and the larger it grows, the richer the expressible class of candidate functions\.

The framework has three layers:

1. \(i\)Primitives\.Two binary operators—EML \(eu−ln⁡ve^\{u\}\-\\ln v, introduced byOdrzywołek \([2026](https://arxiv.org/html/2605.08130#bib.bib11)\)\) and SOL \(sin⁡u−cos⁡v\\sin u\-\\cos v, introduced here\)—used with terminals11andxx—that seed the library with core atoms at depth 1\. EML produces exponential and logarithmic atoms natively and is theoretically complete \(conditional on the cited result\)\. SOL produces trigonometric atoms natively; it is not individually complete, but without it, every trigonometric atom in the library would require an EML tree of depth∼\\sim8, inflating the library by orders of magnitude for the same coverage\. The mixed grammarS→1​∣x∣​eml⁡\(S,S\)∣sol⁡\(S,S\)S\\to 1\\mid x\\mid\\operatorname\{eml\}\(S,S\)\\mid\\operatorname\{sol\}\(S,S\)retains completeness while keeping both exponential and trigonometric atoms shallow\.
2. \(ii\)Self\-expanding atom library\.A layered construction procedure that, given a depth budgetdmaxd\_\{\\max\}and a data grid, builds a library of function–derivative pairs by recursive application of the generating primitives, the product rule, and the chain rule\. The library grows with depth: a shallow build \(dmax=1d\_\{\\max\}=1\) yields hundreds of atoms; a deep build \(dmax=3d\_\{\\max\}=3\) yields tens of thousands\. Previously discovered functions are retained across problems, so the library accumulates knowledge over time\.
3. \(iii\)Additive atomic forests\.A finite sumF=∑kck​TkF=\\sum\_\{k\}c\_\{k\}T\_\{k\}of atoms \(or, more generally, of parameterised trees composed via multiplicative nodes\), whose derivativeF′=∑kck​Tk′F^\{\\prime\}=\\sum\_\{k\}c\_\{k\}T\_\{k\}^\{\\prime\}is fitted to data\. The additive structure represents top\-level sums outside the primitive trees, rather than forcing addition to be encoded inside a single generating\-primitive tree\.

The intended contribution is the conjunction of these features: symbolic output for both a derivative\-matched function and its antiderivative, obtained without a separate symbolic\-integration step, together with a library that can feed each verified discovery back into a growing knowledge base\.

## 2Derivative algebra of elementary functions

### 2\.1The elementary functions

Letℰ\\mathcal\{E\}denote the field of elementary functions in one variable\(Liouville,[1833](https://arxiv.org/html/2605.08130#bib.bib8); Ritt,[1948](https://arxiv.org/html/2605.08130#bib.bib15)\): the closure of the rational functionsℂ​\(x\)\\mathbb\{C\}\(x\)underexp\\exp,log\\log, and composition\. This includes all algebraic, trigonometric, hyperbolic, and inverse trigonometric functions\.

### 2\.2Closure under differentiation via the product rule

###### Proposition 2\.1\(Product\-rule closure\)\.

Let𝒜=\{\(φk,φk′\)\}k=1M\\mathcal\{A\}=\\\{\(\\varphi\_\{k\},\\varphi\_\{k\}^\{\\prime\}\)\\\}\_\{k=1\}^\{M\}be a finite set of function–derivative pairs, with eachφk∈ℰ\\varphi\_\{k\}\\in\\mathcal\{E\}\. Define the*product\-rule closure*𝒜¯\\overline\{\\mathcal\{A\}\}as the smallest set containing𝒜\\mathcal\{A\}and closed under:

1. \(a\)Products:if\(φ,φ′\),\(ψ,ψ′\)∈𝒜¯\(\\varphi,\\varphi^\{\\prime\}\),\(\\psi,\\psi^\{\\prime\}\)\\in\\overline\{\\mathcal\{A\}\}, then\(φ⋅ψ,φ′​ψ\+φ​ψ′\)∈𝒜¯\(\\varphi\\cdot\\psi,\\;\\varphi^\{\\prime\}\\psi\+\\varphi\\psi^\{\\prime\}\)\\in\\overline\{\\mathcal\{A\}\}\.
2. \(b\)Compositions:if\(φ,φ′\)∈𝒜¯\(\\varphi,\\varphi^\{\\prime\}\)\\in\\overline\{\\mathcal\{A\}\}andffis an elementary function with known derivativef′f^\{\\prime\}, then\(f∘φ,\(f′∘φ\)⋅φ′\)∈𝒜¯\(f\\circ\\varphi,\\;\(f^\{\\prime\}\\circ\\varphi\)\\cdot\\varphi^\{\\prime\}\)\\in\\overline\{\\mathcal\{A\}\}\.
3. \(c\)Linear combinations:if\(φ,φ′\),\(ψ,ψ′\)∈𝒜¯\(\\varphi,\\varphi^\{\\prime\}\),\(\\psi,\\psi^\{\\prime\}\)\\in\\overline\{\\mathcal\{A\}\}andc∈ℝc\\in\\mathbb\{R\}, then\(φ\+c​ψ,φ′\+c​ψ′\)∈𝒜¯\(\\varphi\+c\\psi,\\;\\varphi^\{\\prime\}\+c\\psi^\{\\prime\}\)\\in\\overline\{\\mathcal\{A\}\}\.

Then𝒜¯\\overline\{\\mathcal\{A\}\}contains, for everyf∈ℰf\\in\\mathcal\{E\}built from elements of𝒜\\mathcal\{A\}by products, compositions, and sums, the pair\(f,f′\)\(f,f^\{\\prime\}\)\.

###### Proof\.

Immediate from the product rule, chain rule, and linearity of differentiation\. ∎

## 3Generating primitives

###### Definition 3\.1\(Generating primitive\)\.

A*generating primitive*is a pair\(P,c\)\(P,c\)whereP:ℂ2→ℂP\\colon\\mathbb\{C\}^\{2\}\\to\\mathbb\{C\}is a binary operation andc∈ℂc\\in\\mathbb\{C\}is a distinguished constant\. The variablexxis also available as a terminal\. The pair satisfies:

1. \(a\)Generation\.Everyf∈ℰf\\in\\mathcal\{E\}can be expressed as the evaluation of a finite expression generated by the grammarS→c​∣x∣​P​\(S,S\)S\\to c\\mid x\\mid P\(S,S\)\.
2. \(b\)Differentiability\.There exist functionsA,BA,Bsuch that for differentiableu​\(x\),v​\(x\)u\(x\),v\(x\), dd​x​P​\(u,v\)=A​\(u,v\)⋅u′\+B​\(u,v\)⋅v′\.\\frac\{d\}\{dx\}\\,P\(u,v\)=A\(u,v\)\\cdot u^\{\\prime\}\+B\(u,v\)\\cdot v^\{\\prime\}\.\(2\)

A pair\(P,c\)\(P,c\)satisfying only condition \(b\) but not \(a\) is called a*supplementary primitive*; it contributes atoms to the library but does not generate all ofℰ\\mathcal\{E\}on its own\.

We use two complementary primitives\. The first, EML, was introduced byOdrzywołek \([2026](https://arxiv.org/html/2605.08130#bib.bib11)\)and is a generating primitive in the sense of[Definition˜3\.1](https://arxiv.org/html/2605.08130#S3.Thmtheorem1)\(satisfying both conditions \(a\) and \(b\)\); the second, SOL, is new to this work and is a supplementary primitive \(satisfying \(b\) but not \(a\)\)\.

### 3\.1EML \(Exponential Minus Logarithm\)

eml⁡\(u,v\)=eu−ln⁡v,c=1\.\\operatorname\{eml\}\(u,v\)=e^\{u\}\-\\ln v,\\qquad c=1\.\(3\)Differentiation coefficients:A​\(u,v\)=euA\(u,v\)=e^\{u\},B​\(u,v\)=−1/vB\(u,v\)=\-1/v, giving

dd​x​eml⁡\(u,v\)=eu⋅u′−v′v\.\\frac\{d\}\{dx\}\\,\\operatorname\{eml\}\(u,v\)=e^\{u\}\\cdot u^\{\\prime\}\-\\frac\{v^\{\\prime\}\}\{v\}\.\(4\)Generation:eml⁡\(x,1\)=ex\\operatorname\{eml\}\(x,1\)=e^\{x\}andeml⁡\(1,eml⁡\(eml⁡\(1,x\),1\)\)=ln⁡x\\operatorname\{eml\}\(1,\\operatorname\{eml\}\(\\operatorname\{eml\}\(1,x\),1\)\)=\\ln x\. Sinceexp\\expandln\\lngenerateℰ\\mathcal\{E\}overℂ\\mathbb\{C\}, the cited EML result implies that the grammarS→1​∣x∣​eml⁡\(S,S\)S\\to 1\\mid x\\mid\\operatorname\{eml\}\(S,S\)is complete, provided the cited completeness theorem is applicable to the present terminal set\.

### 3\.2SOL \(Sine Of, Less cosine\)

sol⁡\(u,v\)=sin⁡\(u\)−cos⁡\(v\),c=1\.\\operatorname\{sol\}\(u,v\)=\\sin\(u\)\-\\cos\(v\),\\qquad c=1\.\(5\)Differentiation coefficients:A​\(u,v\)=cos⁡\(u\)A\(u,v\)=\\cos\(u\),B​\(u,v\)=sin⁡\(v\)B\(u,v\)=\\sin\(v\), giving

dd​x​sol⁡\(u,v\)=cos⁡\(u\)⋅u′\+sin⁡\(v\)⋅v′\.\\frac\{d\}\{dx\}\\,\\operatorname\{sol\}\(u,v\)=\\cos\(u\)\\cdot u^\{\\prime\}\+\\sin\(v\)\\cdot v^\{\\prime\}\.\(6\)At depth 1 with leaves from\{1,x\}\\\{1,x\\\}:

sol⁡\(x,1\)\\displaystyle\\operatorname\{sol\}\(x,1\)=sin⁡x−cos⁡1,\\displaystyle=\\sin x\-\\cos 1,dd​x\\displaystyle\\qquad\\tfrac\{d\}\{dx\}=cos⁡x,\\displaystyle=\\cos x,\(7\)sol⁡\(1,x\)\\displaystyle\\operatorname\{sol\}\(1,x\)=sin⁡1−cos⁡x,\\displaystyle=\\sin 1\-\\cos x,dd​x\\displaystyle\\qquad\\tfrac\{d\}\{dx\}=sin⁡x,\\displaystyle=\\sin x,\(8\)sol⁡\(x,x\)\\displaystyle\\operatorname\{sol\}\(x,x\)=sin⁡x−cos⁡x,\\displaystyle=\\sin x\-\\cos x,dd​x\\displaystyle\\qquad\\tfrac\{d\}\{dx\}=cos⁡x\+sin⁡x\.\\displaystyle=\\cos x\+\\sin x\.\(9\)

## 4Canonical forms and parameterised trees

Fix a primitive\(P,c\)\(P,c\)\(generating or supplementary\)\.

###### Definition 4\.1\(Canonical form of depthdd\)\.

The canonical form of depthddis a full binary tree with2d2^\{d\}leaves and2d−12^\{d\}\-1internal nodes, all applyingPP\. Each leafℓ\\ellcomputes

Lℓ​\(x;αℓ,βℓ\)=σ​\(αℓ\)⋅1\+σ​\(βℓ\)⋅x,L\_\{\\ell\}\(x;\\,\\alpha\_\{\\ell\},\\beta\_\{\\ell\}\)=\\sigma\(\\alpha\_\{\\ell\}\)\\cdot 1\+\\sigma\(\\beta\_\{\\ell\}\)\\cdot x,\(11\)whereσ\\sigmais softmax over the pair\(αℓ,βℓ\)\(\\alpha\_\{\\ell\},\\beta\_\{\\ell\}\)\. After snapping, each leaf selects either the constant11or the variablexx\.

###### Proposition 4\.2\(Recursive derivative\)\.

The derivative of a canonical form is computable bottom\-up:

Leaf:Lℓ′​\(x\)=σ​\(βℓ\),\\displaystyle\\quad L^\{\\prime\}\_\{\\ell\}\(x\)=\\sigma\(\\beta\_\{\\ell\}\),\(12\)Node:dd​x​P​\(u,v\)=A​\(u,v\)⋅u′\+B​\(u,v\)⋅v′,\\displaystyle\\quad\\tfrac\{d\}\{dx\}P\(u,v\)=A\(u,v\)\\cdot u^\{\\prime\}\+B\(u,v\)\\cdot v^\{\\prime\},\(13\)inO​\(2d\)O\(2^\{d\}\)operations, the same cost as forward evaluation\.

## 5Additive atomic forests

###### Definition 5\.1\(Atomic forest\)\.

An*atomic forest*of widthKKis

F​\(x;𝚯\)=∑k=1KTk​\(x;𝜽k\),F\(x;\\,\\boldsymbol\{\\Theta\}\)=\\sum\_\{k=1\}^\{K\}T\_\{k\}\(x;\\,\\boldsymbol\{\\theta\}\_\{k\}\),\(14\)where each atomTkT\_\{k\}is a canonical form \(EML or SOL\) of bounded depth, and the derivative isF′=∑kTk′F^\{\\prime\}=\\sum\_\{k\}T\_\{k\}^\{\\prime\}by linearity\.

###### Theorem 5\.2\(Conditional completeness\)\.

Assume the EML completeness theorem cited in[Section˜3](https://arxiv.org/html/2605.08130#S3)holds for the terminal set\{1,x\}\\\{1,x\\\}\. Then, for everyf∈ℰf\\in\\mathcal\{E\}, there existK,d∈ℕK,d\\in\\mathbb\{N\}and parameters𝚯∗\\boldsymbol\{\\Theta\}^\{\*\}such thatF​\(x;𝚯∗\)=f​\(x\)F\(x;\\boldsymbol\{\\Theta\}^\{\*\}\)=f\(x\)on its domain\.

###### Proof\.

Under the stated assumption, everyf∈ℰf\\in\\mathcal\{E\}is representable as a single EML tree of some finite depthd0d\_\{0\}with leaves in\{1,x\}\\\{1,x\\\}\. SettingK=1K=1recovers this single\-tree representation as a forest\. ∎

###### Theorem 5\.3\(Additive depth bound\)\.

Iff1,…,fK∈ℰf\_\{1\},\\ldots,f\_\{K\}\\in\\mathcal\{E\}are representable by atoms of depthsd1,…,dKd\_\{1\},\\ldots,d\_\{K\}, thenf=f1\+⋯\+fKf=f\_\{1\}\+\\cdots\+f\_\{K\}is representable by an additive forest of widthKKand per\-atom depthmaxk⁡dk\\max\_\{k\}d\_\{k\}\.

###### Proof\.

Use one atom for each summandfkf\_\{k\}and sum the resulting atoms in the forest\. The forest depth is the maximum of the depths of the summand atoms, because the top\-level addition is represented by the forest structure rather than by extra primitive nodes inside any single tree\. ∎

## 6Multiplicative nodes and enhanced forests

The additive forest avoids encoding top\-level sums inside a single primitive tree\. We similarly avoid encoding many products inside a single primitive tree by introducing a structural product node\.

###### Definition 6\.1\(Multiplicative node\)\.

Given two child nodesAAandBB\(each of any type\), the multiplicative node computes

M​\(x\)=A​\(x\)⋅B​\(x\),M′​\(x\)=A′​\(x\)​B​\(x\)\+A​\(x\)​B′​\(x\)\.M\(x\)=A\(x\)\\cdot B\(x\),\\qquad M^\{\\prime\}\(x\)=A^\{\\prime\}\(x\)\\,B\(x\)\+A\(x\)\\,B^\{\\prime\}\(x\)\.\(15\)This is the product rule \([1](https://arxiv.org/html/2605.08130#S1.E1)\) applied as a structural element\. It enables expressions such asx⋅exx\\cdot e^\{x\}to be represented as a shallow product of already available children, without requiring the product to be encoded inside a single EML tree\.

###### Definition 6\.2\(Enhanced forest\)\.

An*enhanced forest*is an additive forestF​\(x;𝚯\)=∑k=1KTk​\(x;𝜽k\)F\(x;\\boldsymbol\{\\Theta\}\)=\\sum\_\{k=1\}^\{K\}T\_\{k\}\(x;\\boldsymbol\{\\theta\}\_\{k\}\)in which each atomTkT\_\{k\}may be any of: aLeafNode\(constant11or identityxx\), anEMLMasterFormula\(canonical form witheml\\operatorname\{eml\}\), aSOLMasterFormula\(canonical form withsol\\operatorname\{sol\}\), or aMultiplicativeNode\(product of two children of any type\)\.

###### Example 6\.3\.

To recover∫\(ex\+x​cos⁡x\+sin⁡x\)​dx\\int\(e^\{x\}\+x\\cos x\+\\sin x\)\\,\\mathrm\{d\}x, the enhanced forest can use two atoms:T1T\_\{1\}: EMLMasterFormula\(depth=1\)→\\todiscoversexe^\{x\};T2T\_\{2\}: MultiplicativeNode\(LeafNode, SOLMasterFormula\(depth=1\)\)→\\todiscoversx​sin⁡xx\\sin x\. Indeed,

dd​x​\(ex\+x​sin⁡x\)=ex\+x​cos⁡x\+sin⁡x\.\\frac\{d\}\{dx\}\\bigl\(e^\{x\}\+x\\sin x\\bigr\)=e^\{x\}\+x\\cos x\+\\sin x\.Each atom operates at depth≤1\\leq 1\.

## 7The derivative\-matching principle

###### Definition 7\.1\(Derivative\-matching problem\)\.

Given data\{\(xi,yi\)\}i=1N\\\{\(x\_\{i\},y\_\{i\}\)\\\}\_\{i=1\}^\{N\}from a functionff, find a forestF​\(x;𝚯\)F\(x;\\boldsymbol\{\\Theta\}\)minimising

ℒ​\(𝚯\)=1N​∑i=1N\|F′​\(xi;𝚯\)−yi\|2\.\\mathcal\{L\}\(\\boldsymbol\{\\Theta\}\)=\\frac\{1\}\{N\}\\sum\_\{i=1\}^\{N\}\\bigl\|\\,F^\{\\prime\}\(x\_\{i\};\\boldsymbol\{\\Theta\}\)\-y\_\{i\}\\bigr\|^\{2\}\.\(16\)

###### Theorem 7\.2\(Analytic simultaneous recovery\)\.

LetI⊂ℝI\\subset\\mathbb\{R\}be a connected open interval on which the expressions are defined, and letf=g′f=g^\{\\prime\}for someg∈ℰg\\in\\mathcal\{E\}\. Suppose that, after snapping, the atomic forestF​\(x;𝚯∗\)F\(x;\\boldsymbol\{\\Theta\}^\{\*\}\)satisfiesF′​\(x;𝚯∗\)=f​\(x\)F^\{\\prime\}\(x;\\boldsymbol\{\\Theta\}^\{\*\}\)=f\(x\)on an infinite subset ofIIwith an accumulation point inII\. Then: \(i\)F′​\(x;𝚯∗\)=f​\(x\)F^\{\\prime\}\(x;\\boldsymbol\{\\Theta\}^\{\*\}\)=f\(x\)onII; \(ii\)F​\(x;𝚯∗\)=g​\(x\)\+CF\(x;\\boldsymbol\{\\Theta\}^\{\*\}\)=g\(x\)\+Cfor some constantCC; and \(iii\) both the matched function and its antiderivative are available as explicit symbolic expressions from𝚯∗\\boldsymbol\{\\Theta\}^\{\*\}alone\. No symbolic integration algorithm is invoked\.

###### Proof\.

After snapping, each atom is a concrete elementary function, soF′F^\{\\prime\}is elementary onII\. Two elementary functions that agree on a set with an accumulation point in a common analytic domain agree identically on that domain\. HenceF′=f=g′F^\{\\prime\}=f=g^\{\\prime\}, soF=g\+CF=g\+C\. ∎

## 8Self\-expanding atom libraries

The enhanced forest \([Section˜6](https://arxiv.org/html/2605.08130#S6)\) is trained by gradient descent—a non\-convex problem requiring multi\-restart\. An alternative is to*materialise*the library: precompute a large collection of atoms with their derivatives, and search for a sparse linear combination\. This reduces the problem toℓ1\\ell\_\{1\}\-regularised regression\(Tibshirani,[1996](https://arxiv.org/html/2605.08130#bib.bib18)\), which is convex\.

The library is not a fixed catalogue\. It is constructed by a*procedure*that, given a depth budget and a data grid, builds atoms from a seed set by recursive application of the generating primitives, the product rule, and the chain rule\. The procedure can be run to any depth the hardware allows: a shallow build yields hundreds of atoms; a deep build yields tens of thousands\. Crucially,*the library grows over time*: every function discovered by the search \(or supplied externally\) is folded back in as a new seed, and its derivative is computed automatically by the product and chain rules\. The library thus functions as a persistent knowledge base that becomes more capable with each problem it solves\.

### 8\.1Construction procedure

The library self\-constructs in layers, each adding atoms of greater compositional depth\.

#### Layer 0: rational powers\.

The seed:xp/qx^\{p/q\}forp∈\{−4,…,15\}p\\in\\\{\-4,\\ldots,15\\\}andq∈\{1,2,3,4\}q\\in\\\{1,2,3,4\\\}, with derivative\(p/q\)​xp/q−1\(p/q\)\\,x^\{p/q\-1\}\. Candidates producing non\-finite values on the data grid are discarded\.

#### Layer 1: depth\-1 generating\-primitive trees\.

For each pair of linear inner functionsℓ1=a1​x\+b1\\ell\_\{1\}=a\_\{1\}x\+b\_\{1\},ℓ2=a2​x\+b2\\ell\_\{2\}=a\_\{2\}x\+b\_\{2\}with integer slopes and offsets:

eml⁡\(ℓ1,ℓ2\)=ea1​x\+b1−ln⁡\(a2​x\+b2\),\\displaystyle\\operatorname\{eml\}\(\\ell\_\{1\},\\ell\_\{2\}\)=e^\{a\_\{1\}x\+b\_\{1\}\}\-\\ln\(a\_\{2\}x\+b\_\{2\}\),derivative by \([4](https://arxiv.org/html/2605.08130#S3.E4)\),\\displaystyle\\text\{derivative by \\eqref\{eq:eml\-diff\}\},\(17\)sol⁡\(ℓ1,ℓ2\)=sin⁡\(a1​x\+b1\)−cos⁡\(a2​x\+b2\),\\displaystyle\\operatorname\{sol\}\(\\ell\_\{1\},\\ell\_\{2\}\)=\\sin\(a\_\{1\}x\+b\_\{1\}\)\-\\cos\(a\_\{2\}x\+b\_\{2\}\),derivative by \([6](https://arxiv.org/html/2605.08130#S3.E6)\)\.\\displaystyle\\text\{derivative by \\eqref\{eq:sol\-diff\}\}\.\(18\)EML atoms requireℓ2\>0\\ell\_\{2\}\>0on the grid; SOL atoms are always well\-defined\. The two operators contribute roughly equal numbers of atoms at this layer: EML seeds the exponential and logarithmic families, SOL seeds the trigonometric family\. Together they ensure thatea​x\+be^\{ax\+b\},ln⁡\(a​x\+b\)\\ln\(ax\+b\),sin⁡\(a​x\+b\)\\sin\(ax\+b\), andcos⁡\(a​x\+b\)\\cos\(ax\+b\)are all available at depth 1 for any integeraain the specified range\. The range of slopes and offsets is a hyperparameter; wider ranges produce more atoms\.

#### Layer 2: compositions with quadratic arguments\.

For quadratic inner functionsg​\(x\)=q​x2\+a​x\+cg\(x\)=qx^\{2\}\+ax\+c: the atomsexp⁡\(g\)\\exp\(g\),−ln⁡\(g\)\-\\ln\(g\),sin⁡\(g\)\\sin\(g\),cos⁡\(g\)\\cos\(g\),1/g1/g, and selectedarctan\\arctanandarcsin\\arcsinatoms, each with chain\-rule derivative\.

#### Layer 3: products and cross\-terms \(product\-rule expansion\)\.

This is where the product\-rule closure \([Proposition˜2\.1](https://arxiv.org/html/2605.08130#S2.Thmtheorem1)\) is applied structurally\. For each pair of atoms\(φi,φj\)\(\\varphi\_\{i\},\\varphi\_\{j\}\)from previous layers: the productφi⋅φj\\varphi\_\{i\}\\cdot\\varphi\_\{j\}with derivativeφi′​φj\+φi​φj′\\varphi\_\{i\}^\{\\prime\}\\varphi\_\{j\}\+\\varphi\_\{i\}\\varphi\_\{j\}^\{\\prime\}\. The most distinctive atoms at this layer are the EML×\\timesSOL cross\-products—atoms such asex​sin⁡xe^\{x\}\\sin x,e−x​cos⁡2​xe^\{\-x\}\\cos 2x, andx​e−x​sin⁡xxe^\{\-x\}\\sin x—which combine exponential and trigonometric structure at shallow combined depth\. These mixed atoms would require depth∼\\sim9 or more in a pure EML library; having both primitives makes them available at depth 1 plus a product node\. Variable\-times\-atom productsx⋅φix\\cdot\\varphi\_\{i\}are also included\.

#### Layer 4: depth\-2 and depth\-3 nestings\.

For each atomφi\\varphi\_\{i\}from earlier layers, generate all nestings\{exp⁡\(φi\),sin⁡\(φi\),cos⁡\(φi\),−ln⁡\(φi\),1/φi,arctan⁡\(φi\)\}\\\{\\exp\(\\varphi\_\{i\}\),\\,\\sin\(\\varphi\_\{i\}\),\\,\\cos\(\\varphi\_\{i\}\),\\,\-\\ln\(\\varphi\_\{i\}\),\\,1/\\varphi\_\{i\},\\,\\arctan\(\\varphi\_\{i\}\)\\\}with chain\-rule derivatives\. The same operations are applied again for depth\-3 nestings\. All depth\-2 and depth\-3 atoms are also multiplied byxx\.

#### Deduplication\.

After each layer, atoms are deduplicated: if an existing atomψ\\psisatisfies\|corr​\(ψ,φ\)\|\>0\.999\|\\text\{corr\}\(\\psi,\\varphi\)\|\>0\.999on the data grid, the candidate is rejected\.

#### Knowledge\-base growth\.

When the search \([Section˜9](https://arxiv.org/html/2605.08130#S9)\) discovers a functionffthat is not already in the library, the pair\(f,f′\)\(f,f^\{\\prime\}\)is added\. On subsequent problems, this atom and all its compositions and products are available as candidates\. The library thus accumulates domain\-specific knowledge: a library trained on physics problems will contain Planck distributions, Yukawa potentials, and Lorentz factors; one trained on biomedical data will contain saturation curves and dose–response functions\.

#### Scaling\.

The library size grows roughly exponentially with the depth budgetdmaxd\_\{\\max\}: a depth\-1 build produces hundreds of atoms, depth\-2 produces thousands, and depth\-3 produces tens of thousands\. The search algorithms \([Section˜9](https://arxiv.org/html/2605.08130#S9)\) scale gracefully: theK=1K=1scan is linear in library size, and theK=2K=2Gram\-matrix search is quadratic but fully vectorised on GPU\. The library can grow as large as available memory permits\. Its expressive power grows with the number and diversity of atoms, while validation performance remains a statistical model\-selection question\.

## 9Search algorithms

Given a targetf​\(x\)f\(x\)\(symbolic or data\), the solver seeks coefficientsc1,…,cKc\_\{1\},\\ldots,c\_\{K\}and atom indicesi1,…,iKi\_\{1\},\\ldots,i\_\{K\}such that

f​\(x\)≈∑k=1Kck​φik′​\(x\)\.f\(x\)\\approx\\sum\_\{k=1\}^\{K\}c\_\{k\}\\,\\varphi^\{\\prime\}\_\{i\_\{k\}\}\(x\)\.\(19\)When the match is exact, the antiderivative isF​\(x\)=∑kck​φik​\(x\)F\(x\)=\\sum\_\{k\}c\_\{k\}\\,\\varphi\_\{i\_\{k\}\}\(x\)\.

### 9\.1Precomputation

LetD∈ℝN×MD\\in\\mathbb\{R\}^\{N\\times M\}be the derivative matrix withDi​k=φk′​\(xi\)D\_\{ik\}=\\varphi\_\{k\}^\{\\prime\}\(x\_\{i\}\), whereM=\|𝒜\|M=\|\\mathcal\{A\}\|is the current library size\. Precompute the Gram matrixG=D⊤​D∈ℝM×MG=D^\{\\top\}D\\in\\mathbb\{R\}^\{M\\times M\}and target projection𝐝=D⊤​𝐲∈ℝM\\mathbf\{d\}=D^\{\\top\}\\mathbf\{y\}\\in\\mathbb\{R\}^\{M\}in a single GPU matrix multiplication\. Both are reused across all candidate subsets\.

### 9\.2K=1K=1: vectorised scalar regression

For each atomkk\(no intercept\):ck∗=dk/Gk​kc\_\{k\}^\{\*\}=d\_\{k\}/G\_\{kk\}andMSEk=‖𝐲‖2/N−dk2/\(N​Gk​k\)\\operatorname\{MSE\}\_\{k\}=\\\|\\mathbf\{y\}\\\|^\{2\}/N\-d\_\{k\}^\{2\}/\(N\\,G\_\{kk\}\)\. Evaluated for allMMatoms simultaneously inO​\(M\)O\(M\)after theO​\(N​M\)O\(NM\)precomputation\. A constant offset in the derivative target is represented by the identity atomφ​\(x\)=x\\varphi\(x\)=x, whose derivative isφ′​\(x\)=1\\varphi^\{\\prime\}\(x\)=1\. The constant atomφ​\(x\)=1\\varphi\(x\)=1has derivative zero; it represents the arbitrary integration constant and does not improve derivative matching\.

### 9\.3K=2K=2: vectorised Cramer’s rule

For each pair\(i,j\)\(i,j\), the2×22\\times 2systemG\[i​j\],\[i​j\]​𝐜=𝐝\[i​j\]G\_\{\[ij\],\[ij\]\}\\,\\mathbf\{c\}=\\mathbf\{d\}\_\{\[ij\]\}is solved by Cramer’s rule\. The MSE for all\(M2\)\\binom\{M\}\{2\}pairs is computed as a single sequence of broadcasting operations onM×MM\\times Mmatrices—no Python\-level loop\. Pairs are filtered by determinant threshold \(\|Δ\|\>10−30\|\\Delta\|\>10^\{\-30\}\), coefficient cap \(\|ck\|≤20\|c\_\{k\}\|\\leq 20\), and finite MSE\.

### 9\.4K≥3K\\geq 3: beam search

Initialise with the top 200 results atK−1K\-1; for each, try adding every atom from the library and solve the\(K\+1\)×\(K\+1\)\(K\\\!\+\\\!1\)\\times\(K\\\!\+\\\!1\)Gram subproblem on GPU; keep the top 500 candidates\. Cost:O​\(beam×M\)O\(\\text\{beam\}\\times M\)per level\.

### 9\.5Verification

When a candidate achievesMSE<10−15\\operatorname\{MSE\}<10^\{\-15\}, it is verified on an independentxx\-grid to exclude numerical coincidence\. Verified discoveries are added to the library for future use\.

## 10Training \(gradient\-based mode\)

For the enhanced forest \([Definition˜6\.2](https://arxiv.org/html/2605.08130#S6.Thmtheorem2)\), training minimises the derivative\-matching loss \([16](https://arxiv.org/html/2605.08130#S7.E16)\) by Adam\(Kingma & Ba,[2015](https://arxiv.org/html/2605.08130#bib.bib6)\)with analytical backpropagation through the tree structure\.

#### Temperature annealing\.

The softmax temperature in each leaf decreases linearly from 1\.0 to 0\.1 over the first 80% of training, pushing weights toward one\-hot vectors\.

#### Log\-loss variant\.

For deep trees \(d≥3d\\geq 3\),log⁡\(1\+\|r\|2\)\\log\(1\+\|r\|^\{2\}\)replaces\|r\|2\|r\|^\{2\}to stabilise gradients when residuals are large\.

#### Gradient clipping\.

Gradients are clipped to norm≤5\\leq 5\.

#### Multi\-restart\.

Multiple random initialisations run in parallel\. Training stops when any restart achievesMSE<10−10\\operatorname\{MSE\}<10^\{\-10\}after snapping\.

#### Snapping\.

After training, each leaf’s softmax weights are replaced by the nearest one\-hot vector\. Constant\-valued trees are pruned\. The snapped formula is verified on the grid\.

## 11Structural complexity

###### Definition 11\.1\(Atomic complexity\)\.

The*atomic complexity*off∈ℰf\\in\\mathcal\{E\}is the pairκ​\(f\)=\(K∗,d∗\)\\kappa\(f\)=\(K^\{\*\},d^\{\*\}\): the lexicographically minimal width and per\-atom depth of an atomic forest representingff\. We callK∗K^\{\*\}the*additive complexity*andd∗d^\{\*\}the*compositional complexity*\.

###### Proposition 11\.2\(Derivative evaluation graph size\)\.

LetF=∑k=1KTkF=\\sum\_\{k=1\}^\{K\}T\_\{k\}be a forest whose atoms are binary trees of depth at mostdd\. ThenF′F^\{\\prime\}can be evaluated by a derivative computation graph of sizeO​\(K​2d\)O\(K2^\{d\}\), using one bottom\-up derivative pass through each atom\.

###### Proof\.

A full binary tree of depthddhasO​\(2d\)O\(2^\{d\}\)nodes\. At each leaf and internal node, the derivative value is computed from the stored child values and child derivatives using either the leaf rule or the corresponding chain/product rule\. Thus the derivative computation adds only constant work per node, and the full forest derivative is obtained by summing theKKatom derivatives\. ∎

### 11\.1Connection to the Liouville–Risch decomposition

The Liouville–Risch theorem\(Liouville,[1833](https://arxiv.org/html/2605.08130#bib.bib8); Risch,[1969](https://arxiv.org/html/2605.08130#bib.bib13)\)states that if∫f\\int fis elementary, then∫f=v0\+∑jcj​ln⁡\(vj\)\\int f=v\_\{0\}\+\\sum\_\{j\}c\_\{j\}\\ln\(v\_\{j\}\)for elementaryvjv\_\{j\}\. When the termsv0v\_\{0\}andln⁡\(vj\)\\ln\(v\_\{j\}\)are available as atoms, this decomposition corresponds to an additive forest with one atom forv0v\_\{0\}and one atom per logarithmic term\. Thus additive complexity is related to, but not identical with, the number of extensions appearing in a Risch\-style decomposition\.

## 12Comparison with existing methods

∗SINDy requires a predefined,*fixed*candidate library and does not recover integrals or conserved quantities\.‡In gradient\-based mode, the derivative match and antiderivative are obtained in a single optimisation\. In materialised\-library mode, library construction precedes the search; once the library is built, each query is a single pass\.

The unique property of atomic forests is the conjunction of all four columns: symbolic output for both the function and its antiderivative, obtained in a single pass, with a library that grows with each discovery\.

## 13Empirical results

Different experiments use libraries of different depth, reflecting the self\-constructing nature of the framework\. Classification benchmarks \([Section˜13\.1](https://arxiv.org/html/2605.08130#S13.SS1)\) use a depth\-1 library built by the standard construction procedure; the Feynman test \([Section˜13\.4](https://arxiv.org/html/2605.08130#S13.SS4)\) uses a self\-built library expanded greedily to depth 3; the demonstration application \([Section˜13\.5](https://arxiv.org/html/2605.08130#S13.SS5)\) uses a domain\-adapted library\.

### 13\.1Classification: atoms vs\. XGBoost

We expand each input variable through the atom library and fitℓ1\\ell\_\{1\}/ℓ2\\ell\_\{2\}\-regularised logistic regression on the expanded feature matrix\. All results are 5\-fold cross\-validated\. XGBoost\(Chen & Guestrin,[2016](https://arxiv.org/html/2605.08130#bib.bib2)\)serves as the black\-box baseline\.

Table 1:Classification accuracy \(5\-fold CV\)\. Atoms win on 13/17 datasets\.†The sklearn diabetes dataset is natively a regression task; it was binarised here by thresholding the target at the median\.Limitations\.All entries report raw accuracy only\. For datasets with imbalanced class ratios \(notably Haberman,∼\\sim73:27\), accuracy can be misleading; balanced accuracy, AUC, orF1F\_\{1\}should be reported in a full empirical study\.

### 13\.2Regression

On the synthetic regression task, the model reachesR2=0\.999R^\{2\}=0\.999, outperforming both XGBoost and ridge regression in this run\. On the sklearn diabetes dataset, atoms outperform XGBoost \(R2=0\.336R^\{2\}=0\.336\) while matching ridge \(R2=0\.479R^\{2\}=0\.479\), despite using a linear combination of nonlinear atoms rather than a polynomial expansion\.

### 13\.3Interpretable formulas and their derivatives

A distinctive property of the framework is that every fitted model is a named formula with an analytic gradient\. We highlight three examples\.

#### Haberman survival \(atoms: 75\.8%, XGBoost: 66\.0%\)\.

Withℓ1\\ell\_\{1\}regularisation, the model selects a single atom:F​\(x2\)=4\.404​arctan⁡\(x2\)−0\.859F\(x\_\{2\}\)=4\.404\\,\\operatorname\{arctan\}\(x\_\{2\}\)\-0\.859, with derivativeF′​\(x2\)=4\.404/\(1\+x22\)F^\{\\prime\}\(x\_\{2\}\)=4\.404/\(1\+x\_\{2\}^\{2\}\)\. The derivative quantifies the diminishing marginal effect of positive lymph nodes on mortality risk\.

#### Heart disease \(atoms: 83\.5%, XGBoost: 79\.8%\)\.

The dominant variable isx11x\_\{11\}\(number of major vessels coloured by fluoroscopy\)\. Thetanh\\operatorname\{tanh\}saturation in the fitted model indicates that going from 0 to 1 diseased vessel is the largest risk jump; the analytic derivative confirms:∂F/∂x11\|x11=0=21\.93\\partial F/\\partial x\_\{11\}\|\_\{x\_\{11\}=0\}=21\.93versus∂F/∂x11\|x11=2=1\.07\\partial F/\\partial x\_\{11\}\|\_\{x\_\{11\}=2\}=1\.07\.

#### Synthetic regression: high\-R2R^\{2\}recovery\.

Ground truth:y=3​e−x0\+2​sin⁡\(x1\)\+12​x22y=3e^\{\-x\_\{0\}\}\+2\\sin\(x\_\{1\}\)\+\\tfrac\{1\}\{2\}x\_\{2\}^\{2\}\. One fitted model wasy^=0\.97​e−x0\+1\.46​sin⁡\(x1\)\+1\.36​cos⁡\(x1\)\+0\.49​x22\+⋯\\hat\{y\}=0\.97\\,e^\{\-x\_\{0\}\}\+1\.46\\sin\(x\_\{1\}\)\+1\.36\\cos\(x\_\{1\}\)\+0\.49\\,x\_\{2\}^\{2\}\+\\cdots\. The trigonometric part has amplitude1\.462\+1\.362≈1\.99\\sqrt\{1\.46^\{2\}\+1\.36^\{2\}\}\\approx 1\.99and can be written as a phase\-shifted sinusoid, but it is not term\-by\-term identical to2​sin⁡\(x1\)2\\sin\(x\_\{1\}\)unless the phase shift is explained by preprocessing or an equivalent coordinate transformation\.

### 13\.4Feynman representability test

We test whether a self\-built library can*represent*45 equations from the Feynman Symbolic Regression Database\(Udrescu & Tegmark,[2020](https://arxiv.org/html/2605.08130#bib.bib19)\)\. The library is constructed from scratch using the procedure of[Section˜8\.1](https://arxiv.org/html/2605.08130#S8.SS1)with greedy expansion to depth 3: no physics\-specific atoms are provided; the library must discover the necessary building blocks itself\.

For each equation, we generate exact data from the known formula, construct the multi\-variable library, and run exhaustive search atK=1,2,3,4K=1,2,3,4\.

All 14 exact recoveries occur atK≤2K\\leq 2: 11 atK=1K=1\(Coulomb, Biot–Savart, ideal gas, etc\.\) and 3 atK=2K=2\. The 13 failures are structural: 8 require composed argumentsf​\(xi⋅xj/xk\)f\(x\_\{i\}\\cdot x\_\{j\}/x\_\{k\}\)not yet in the self\-built library, 3 are non\-separable, and 2 are rational\-of\-sums\. A library seeded with physics\-specific atoms \(Boltzmann factors, Lorentz factors\) or one that has accumulated knowledge from prior problems may cover more of these cases—one motivation for the growing knowledge\-base design\. Total time for this run was 76 s on an NVIDIA T4\.

### 13\.5Demonstration: radial acceleration relation

As a demonstration of the method on real scientific data, we apply the atom search to the radial acceleration relation \(RAR\) in galaxies\. This is an observational correlation, not a test of the method’s physical content: we do not account for per\-galaxy error bars, distance uncertainties, or the physical constraints that a serious RAR analysis would require\(McGaugh et al\.,[2016](https://arxiv.org/html/2605.08130#bib.bib9)\)\. The purpose is to show that the atom search can produce compact candidate fits on a well\-studied astrophysical dataset and to compare the resulting wMSE against standard reference curves\.

Using the SPARC database\(Lelli et al\.,[2017](https://arxiv.org/html/2605.08130#bib.bib7)\)111SPARC data publicly available at[http://astroweb\.cwru\.edu/SPARC/](http://astroweb.cwru.edu/SPARC/)\. We use the standard mass\-to\-light ratiosΥdisk=0\.5\\Upsilon\_\{\\rm disk\}=0\.5,Υbul=0\.7\\Upsilon\_\{\\rm bul\}=0\.7and the acceleration scalea0=1\.2×10−10a\_\{0\}=1\.2\\times 10^\{\-10\}m/s2\.\(175 galaxies, 3,389 data points\), we search forgobs/a0=∑kck​φk​\(gbar/a0\)g\_\{\\text\{obs\}\}/a\_\{0\}=\\sum\_\{k\}c\_\{k\}\\,\\varphi\_\{k\}\(g\_\{\\text\{bar\}\}/a\_\{0\}\)and compare against standard MOND interpolation functions\(McGaugh et al\.,[2016](https://arxiv.org/html/2605.08130#bib.bib9); Milgrom,[1983](https://arxiv.org/html/2605.08130#bib.bib10)\)\.

TheK=1K=1atom with lowest wMSE for the fitted responsey=gobs/a0y=g\_\{\\rm obs\}/a\_\{0\}isy=x2/3y=x^\{2/3\}, whose wMSE is slightly below that of the three reference MOND curves in this unconstrained search\. The margin is small \(0\.074 vs\. 0\.077–0\.079\) and should not be over\-interpreted, since the comparison does not account for per\-galaxy uncertainties or physical priors\. In MOND notation one usually writesgobs=ν​\(x\)​gbarg\_\{\\rm obs\}=\\nu\(x\)g\_\{\\rm bar\}withx=gbar/a0x=g\_\{\\rm bar\}/a\_\{0\}, so a fity=x2/3y=x^\{2/3\}corresponds toν​\(x\)=y/x=x−1/3\\nu\(x\)=y/x=x^\{\-1/3\}\. This does not satisfy the standard MOND asymptotic limits \(ν→1\\nu\\to 1for largexxandν→x−1/2\\nu\\to x^\{\-1/2\}for smallxx\), so it should be interpreted as the best unconstrained empirical atom in this run, not as a physically valid MOND interpolation function\.

### 13\.6Failure modes

The method fails predictably when: \(i\) high\-dimensional interactions exceed the cross\-product budget; \(ii\) the library is shallow and the target function requires deep compositions not yet discovered; \(iii\) inputs are categorical with no continuous nonlinear structure\. Failures of type \(ii\) are mitigated by deeper builds or by seeding the library with domain\-specific atoms—the self\-expanding design directly addresses this limitation\.

## 14Extensions

### 14\.1Multivariate canonical forms

Leaves are enlarged toLℓ=∑j=0nσ​\(αℓ,j\)​ϕjL\_\{\\ell\}=\\sum\_\{j=0\}^\{n\}\\sigma\(\\alpha\_\{\\ell,j\}\)\\,\\phi\_\{j\}withϕ0=1\\phi\_\{0\}=1,ϕj=xj\\phi\_\{j\}=x\_\{j\}\. Partial derivatives∂F/∂xj\\partial F/\\partial x\_\{j\}are computed by the same recursive rule\.

### 14\.2Conserved quantity discovery

Given𝐱˙=𝐆​\(𝐱\)\\dot\{\\mathbf\{x\}\}=\\mathbf\{G\}\(\\mathbf\{x\}\)and trajectory data, representHHas an atomic forest and minimise∑i\|∇F​\(𝐱i\)⋅𝐆​\(𝐱i\)\|2\\sum\_\{i\}\|\\nabla F\(\\mathbf\{x\}\_\{i\}\)\\cdot\\mathbf\{G\}\(\\mathbf\{x\}\_\{i\}\)\|^\{2\}\. The additive forest aligns with the typicalH=T\+VH=T\+Vdecomposition\(Greydanus et al\.,[2019](https://arxiv.org/html/2605.08130#bib.bib5)\)\.

## 15Conclusion

We have presented a framework for symbolic function and antiderivative discovery built on three ideas: the product\-rule closure of function–derivative pairs, two complementary primitives—EML providing theoretical completeness over the elementary functions, SOL providing practical depth reduction for the trigonometric family—and additive atomic forests that avoid encoding top\-level sums inside one tree while multiplicative nodes provide shallow structured products\.

The framework operates in two modes\. In the gradient\-based mode, enhanced forests with mixed EML, SOL, and multiplicative nodes are trained by derivative\-matching loss, producing a candidate derivative match together with the corresponding symbolic primitive in a single optimisation\. In the materialised\-library mode, a self\-expanding collection of atoms is searched by GPU\-accelerated exhaustive and beam methods\. For a fixed candidate subset the coefficient fit is a convex least squares or regularised\-regression problem; the sparse subset search itself remains combinatorial\.

The self\-expanding library is a distinctive feature: it accumulates knowledge across problems, growing more capable with each discovery\. A library trained on physics problems acquires physics atoms; one trained on biomedical data acquires saturation curves and dose–response functions\. The library’s expressive power grows with size, although validation performance still depends on regularisation, model selection, and data quality\.

Empirically, in the reported runs, sparse atom combinations match or exceed XGBoost on 13/17 classification benchmarks while producing interpretable formulas with analytic derivatives\. A self\-built library to depth 3 gives exact recovery on 31% of the tested Feynman equations and relative\-MSE below 0\.01 on a further 40%\. A demonstration on the galactic radial acceleration relation illustrates how the method can propose compact functional forms from real data\.

## References

- Brunton et al\. \(2016\)S\. L\. Brunton, J\. L\. Proctor, and J\. N\. Kutz\.Discovering governing equations from data by sparse identification of nonlinear dynamical systems\.*Proceedings of the National Academy of Sciences*, 113\(15\):3932–3937, 2016\.
- Chen & Guestrin \(2016\)T\. Chen and C\. Guestrin\.XGBoost: A scalable tree boosting system\.In*Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining*, pp\. 785–794, 2016\.
- Chen et al\. \(2018\)R\. T\. Q\. Chen, Y\. Rubanova, J\. Bettencourt, and D\. Duvenaud\.Neural ordinary differential equations\.In*Advances in Neural Information Processing Systems 31 \(NeurIPS\)*, 2018\.
- Cranmer \(2023\)M\. Cranmer\.Interpretable machine learning for science with PySR and SymbolicRegression\.jl\.*arXiv preprint arXiv:2305\.01582*, 2023\.
- Greydanus et al\. \(2019\)S\. Greydanus, M\. Dzamba, and J\. Yosinski\.Hamiltonian neural networks\.In*Advances in Neural Information Processing Systems 32 \(NeurIPS\)*, 2019\.
- Kingma & Ba \(2015\)D\. P\. Kingma and J\. Ba\.Adam: A method for stochastic optimization\.In*Proceedings of the 3rd International Conference on Learning Representations \(ICLR\)*, 2015\.
- Lelli et al\. \(2017\)F\. Lelli, S\. S\. McGaugh, and J\. M\. Schombert\.SPARC: Mass models for 175 disk galaxies with Spitzer photometry and accurate rotation curves\.*The Astronomical Journal*, 153\(6\):240, 2017\.
- Liouville \(1833\)J\. Liouville\.Premier mémoire sur la détermination des intégrales dont la valeur est algébrique\.*Journal de l’École Polytechnique*, 14:124–148, 1833\.
- McGaugh et al\. \(2016\)S\. S\. McGaugh, F\. Lelli, and J\. M\. Schombert\.Radial acceleration relation in rotationally supported galaxies\.*Physical Review Letters*, 117\(20\):201101, 2016\.
- Milgrom \(1983\)M\. Milgrom\.A modification of the Newtonian dynamics as a possible alternative to the hidden mass hypothesis\.*The Astrophysical Journal*, 270:365–370, 1983\.
- Odrzywołek \(2026\)A\. Odrzywołek\.All elementary functions from a single binary operator\.*arXiv preprint arXiv:2603\.21852*, 2026\.
- Raissi et al\. \(2019\)M\. Raissi, P\. Perdikaris, and G\. E\. Karniadakis\.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, 2019\.
- Risch \(1969\)R\. H\. Risch\.The problem of integration in finite terms\.*Transactions of the American Mathematical Society*, 139:167–189, 1969\.
- Risch \(1970\)R\. H\. Risch\.The solution of the problem of integration in finite terms\.*Bulletin of the American Mathematical Society*, 76\(3\):605–608, 1970\.
- Ritt \(1948\)J\. F\. Ritt\.*Integration in Finite Terms: Liouville’s Theory of Elementary Methods*\.Columbia University Press, 1948\.
- Schmidt & Lipson \(2009\)M\. Schmidt and H\. Lipson\.Distilling free\-form natural laws from experimental data\.*Science*, 324\(5923\):81–85, 2009\.
- Stachowiak \(2026\)T\. Stachowiak\.Algebraic structure behind Odrzywołek’s EML operator\.*arXiv preprint arXiv:2604\.23893*, 2026\.
- Tibshirani \(1996\)R\. Tibshirani\.Regression shrinkage and selection via the lasso\.*Journal of the Royal Statistical Society: Series B*, 58\(1\):267–288, 1996\.
- Udrescu & Tegmark \(2020\)S\.\-M\. Udrescu and M\. Tegmark\.AI Feynman: A physics\-inspired method for symbolic regression\.*Science Advances*, 6\(16\):eaay2631, 2020\.