Sample Complexity of Scientific Discovery: PAC Learnability of Compositional Function Trees

Technical Summary: Sample Complexity of Scientific Discovery

Problem Statement
Scientific discovery via symbolic regression (SR) is frequently characterized as statistically and computationally intractable. The prevailing view suggests that the hypothesis space of symbolic expressions grows combinatorially with depth, rendering the search for optimal expressions NP-hard and the statistical generalization of such models unreliable. This paper challenges the conflation of computational hardness with statistical learnability. It posits that while finding the optimal symbolic structure is computationally difficult, the statistical generalization of compositional function trees may be well-behaved if the underlying operators possess stability properties. The core question addressed is whether the sample complexity of learning compositional operator trees explodes exponentially with depth, or if it is instead controlled by the stability (Lipschitzness) of the operators and the depth of the tree.

Methodology
The authors approach the problem through the lens of Probably Approximately Correct (PAC) learning theory, specifically utilizing Rademacher complexity to bound the generalization error.

1. Hypothesis Space Definition: The study defines a compositional hypothesis space, $H^d_{comp}$, consisting of rooted operator trees of depth $d$. The internal nodes are drawn from a finite vocabulary of smooth, differentiable operators ($H_{base}$), such as $\{+, \times, \sin, \cos, \exp\}$ and affine maps. The leaves are either input coordinates or bounded affine functions.
2. Theoretical Framework:
   • Lipschitz Assumptions: The analysis assumes that base operators are locally Lipschitz continuous with constants $L_h$ on the relevant data ranges. A global Lipschitz constant $L = \max L_h$ is used for bounding.
   • Rademacher Complexity: The authors derive bounds on the empirical Rademacher complexity, $R_n(H^d_{comp})$. They extend classical contraction arguments (Ledoux & Talagrand, 1991) and vector contraction inequalities (Maurer, 2016) to tree-structured compositions.
   • Recursive Bounds: By applying a union bound over the finite operator vocabulary ($K = |H_{base}|$) and utilizing vector contraction for the composition of $b$ children (where $b$ is the arity), the authors establish a recursive relationship between the complexity of depth $d$ and depth $d-1$.
3. Experimental Validation:
   • Synthetic Data: The authors constructed a modular PyTorch codebase to generate synthetic "physics-like" targets with controlled depths ($d \in \{1, 2, 3, 4\}$) and known ground-truth structures.
   • Differentiable Trees: Instead of using genetic programming or discrete search, they trained differentiable operator trees (not MLPs) using gradient descent (Adam) to fit the coefficients of the operators.
   • Metrics: They measured the generalization gap (test MSE minus train MSE) and compared it against the theoretical complexity term $(bL^d)/\sqrt{n}$, where $\hat{L}$ is an empirical estimator of the Lipschitz constant derived from gradient norms.

Key Contributions

1. Statistical Learnability of SR: The paper proves that the generalization capacity of compositional trees does not necessarily blow up exponentially with the number of distinct symbolic structures. Instead, it is governed by the depth $d$ and the Lipschitz constants of the operators.
2. Explicit Complexity Bounds: The authors derive a concrete upper bound for the Rademacher complexity:
   $$R_n(H^d_{comp}) \leq (K \cdot b \sqrt{2} L)^{d-1} R_n(H^1_{comp})$$
   Consequently, the excess risk scales as $O(L^d / \sqrt{n})$ when the vocabulary size $K$ and arity $b$ are treated as constants. This contrasts with bounds that depend on the cardinality of the discrete syntax space.
3. Separation of Search and Estimation: The work formally separates the discrete optimization problem (searching for the tree structure) from the continuous regression problem (fitting coefficients). The derived bounds apply to the latter, clarifying the statistical scaling one should expect once a structure is proposed.
4. Empirical Correlation: The experiments demonstrate a positive correlation between the observed generalization gap and the predicted complexity term $(bL^d)/\sqrt{n}$. The data shows that while the gap increases with depth (due to the accumulation of Lipschitz constants) and decreases with sample size ($1/\sqrt{n}$), the scaling aligns with the theoretical prediction.

Results

• Theoretical: Theorem 4.9 establishes that for bounded outputs and Lipschitz operators, the excess risk decays as $1/\sqrt{n}$, with the constant factor depending exponentially on depth via the product of Lipschitz constants. Corollary 4.11 translates this into a PAC sample complexity bound, showing that sample complexity scales with $L^{2d}$ rather than the combinatorial count of trees.
• Empirical: Across 400 experimental runs varying depth ($d$) and sample size ($n$), the generalization gap was found to track the theoretical term.
  • For shallow depths ($d=1, 2$), the Lipschitz proxy $\hat{L}$ remained stable, and the gap scaled primarily with $1/\sqrt{n}$.
  • For deeper compositions ($d=3, 4$) involving non-linearities like $\exp$, $\hat{L}$ increased, leading to larger gaps, consistent with the multiplicative accumulation of stability constants.
  • A power-law fit on log-log scales yielded a correlation of 0.48, indicating that the compositional complexity term captures a meaningful fraction of the variance in generalization gaps, though other factors (optimization error, noise) also play a role.

Significance and Claims
The paper claims to provide a "learnability certificate" for scientific discovery using symbolic regression. Its primary significance lies in shifting the perspective from the combinatorial explosion of symbolic structures to the compositional stability of the operators.

• Feasibility of Small-Data Discovery: The authors argue that for stable operator vocabularies and modest depths, scientific discovery from small datasets is statistically feasible. The statistical generalization is not the bottleneck; rather, the challenge remains the computational search over structures and the design of inductive biases.
• Clarification of Scaling: The work clarifies that the statistical scaling of nonlinear coefficient estimation is determined by the accumulation of Lipschitz constants along the computation graph, not by the enumeration of possible tree shapes.
• Limitations: The authors modestly note that their guarantees are uniform bounds and do not exploit structural priors (like Occam's razor) as effectively as PAC-Bayes approaches might. Furthermore, the empirical validation relies on controlled, in-distribution synthetic data with one-dimensional inputs; the behavior under distribution shift or misspecification is not guaranteed by these bounds. The estimator $\hat{L}$ measures stability on the fitted map rather than certifying a worst-case global constant.

In conclusion, the paper asserts that symbolic regression models, when viewed as compositional function trees with stable primitives, possess favorable statistical properties that allow for generalization even with limited data, provided the depth and operator stability are managed.