S OLVE FOR THE H YPERPARAMETER , S KIP THE S EARCH : KOLMOGOROV-O PTIMAL S CALING L AWS FOR S PLINE R EGRESSION
Yong Yi Bay∗ Kathleen A. Yearick∗ PhD, University of Illinois at Urbana-Champaign
arXiv:2606.23575v1 [cs.LG] 22 Jun 2026
A BSTRACT Hyperparameter tuning almost always means search: fit the model at every value on a grid, score each by cross-validation, and keep the winner. For spline regression that search is unnecessary. The optimal resolution can be solved for in closed form, to the accuracy an exhaustive search would reach and at a fraction of the compute. Three facts make this possible. Classical approximation theory pins the squared bias to G−2β in the resolution G, which is exactly the Kolmogorov n-width of the smoothness class, so the spline family is approximation-optimal among linear methods of its size; the basis dimension is an explicit polynomial in G; and the leave-one-out (LOO) error of any linear smoother follows from a single fit through the PRESS identity. Balancing the two known curves gives the minimizer analytically. We carry the same calculus from one coordinate to many by replacing ambient input dimension with interaction order, the number of active low-order components in an ANOVA decomposition, and obtain a scaling law in which the optimal resolution and the optimal error are power functions of the effective density n/sr , with the input dimension absent from the exponent. The law becomes an algorithm. KORE (Kolmogorov-optimal Order-aware Resolution Estimation) fits two pilot resolutions, solves a leverage-calibrated 2 × 2 system for the bias and noise/variance scales, and evaluates the closed-form plug-in resolution with a tiny LOO certificate: a fixed dozen fits in place of a full grid sweep, with a consistency guarantee as n grows. Across additive and sparse pairwise targets up to 80 input dimensions, KORE matches exhaustive 3-fold cross-validation and the entire classical full-grid ladder of generalized cross-validation, Mallows’ Cp , AIC, and BIC in accuracy while fitting roughly 8× fewer models; on 36 real tabular datasets it ranks first among 21 methods in accuracy delivered per unit of compute, ahead of tuned gradient boosters and kernel machines. When a target’s complexity lives in low interaction order, solving for the resolution beats searching for it.
1
I NTRODUCTION
Hyperparameter selection is usually done by search: a model is trained at every candidate setting, validation error is recorded, and the best score wins. Cross-validation, the canonical realization of that workflow, is reliable but uninformative; it treats every hyperparameter as an opaque knob and offers no closed-form guidance for which values are worth trying or why one value beats its neighbors. The compute cost scales with the size of the grid multiplied by the number of folds and grows with every new model family. An attractive alternative is to solve for the optimal hyperparameter directly, in the same sense that the minimum of a known function is found by calculus rather than by tabulation. Why spline regression. Solving rather than searching requires a model class where prediction error has an analytically tractable dependence on the hyperparameter of interest. Three ingredients are needed: (i) an approximation theory that specifies how bias scales with the hyperparameter, (ii) an explicit formula for the number of free parameters, which controls variance, and (iii) a closed-form estimate of prediction error, so the bias-variance curve can be pinned down from a small number of fits rather than from an exhaustive ∗
Equal contribution. Correspondence: {yongyibay, kallie.a.yearick}@gmail.com.
1
grid. Spline regression is the natural home for this program. Classical B-spline approximation theory gives precise bias rates as a power of the knot resolution G (de Boor, 2001; Schumaker, 2007); these rates are exactly the Kolmogorov n-widths of the smoothness class, so among all linear methods of the same size the spline family is approximation-optimal (Kolmogorov, 1936; Melkman & Micchelli, 1978; Pinkus, 1985). The basis dimension is an explicit, countable polynomial in G. Because a spline fit is a linear smoother, exact leave-one-out error is available at O(1) per point through the PRESS identity (Allen, 1974) with no refitting. Together, these three properties express the entire bias-variance curve as a closed-form function of a single integer, and they make solving for the minimum a calculus problem rather than a search problem. Other popular model classes lack this combination. Neural-network hyperparameters (width, depth, learning rate) interact through a non-convex loss landscape, and no finite-sample power law links architecture to approximation error; the strength of deep models lies instead in extrapolating beyond the training support, a regime where they can outperform classical smoothers (Bay & Yearick, 2024). Tree ensembles (random forests, gradient boosting) have multiple interacting knobs whose joint effect on bias is data-dependent and resists closed-form treatment. Kernel methods come closest: bandwidth controls a bias-variance tradeoff, but the effective dimensionality of the smoothing problem is implicit in the kernel’s eigenspectrum and cannot be read off a design matrix. In each of these cases search is the only available option. Splines are the mature classical model family in which solving is tractable, and this paper studies how far that solvability extends into multiple coordinates. In spline regression, the key hyperparameter is the resolution G, the number of knot intervals per coordinate. It controls model capacity directly: too few intervals and the spline is too stiff to follow the true signal; too many and it starts fitting noise. For a grid of 20 candidate additive resolutions and 10 candidate sparse pairwise resolutions with 3-fold CV, exhaustive search costs (20 + 10) × 3 + 1 = 91 model fits in total, plus a final refit on the full data. When multiple families are compared the cost multiplies. KORE replaces that entire search with a constant number of fits by deriving a closed-form scaling law for the optimal resolution G⋆ and fitting at just two pilot resolutions to identify the law’s unknown constants. The question this paper answers is whether G⋆ can be predicted instead of searched. The scope is also explicit: the contribution is resolution selection once a structured spline family has been chosen. Penalized generalized additive models tune a continuous roughness penalty for a fixed basis by criteria such as GCV, AIC, or REML (Wood, 2017). KORE targets the complementary discrete question of selecting the right basis resolution when comparing low-order spline dictionaries. The closed-form selector applies inside the spline-ANOVA function class with bounded post-one-hot dimension; on signals dominated by high-order interactions or by deep categorical structure, tuned boosters retain the upper hand, and Section 4.9 delineates the boundary explicitly. The bias-variance tradeoff in resolution. Selecting the best resolution requires a quantitative description of how prediction error depends on G. Write Err(G) for the expected squared error of a spline predictor fitted at resolution G on a new test point. This error decomposes as p(G) Err(G) − σ 2 = A G−2β + τ , | {z } | {zn } variance bias2
(1)
where • σ 2 is the irreducible noise in the data, which no model can remove; • G is the resolution, the number of knot intervals per coordinate; • β is the smoothness exponent of the target function, set by the spline degree k; • p(G) is the basis size, the number of spline basis functions at resolution G; • n is the sample size; • A > 0 is the bias scale and τ > 0 is the noise/variance scale. The bias term A G−2β falls as G grows because a finer grid approximates the true function more closely. The variance term τ p(G)/n rises as G grows because more basis functions imply more parameters to estimate from n noisy observations. The sum traces a U-shaped curve. The resolution at the bottom of that U is the desired one. The entire paper is built on making this decomposition precise for structured multidimensional splines and then solving for its minimizer. 2
One coordinate: a clean classical answer. In one coordinate, G is the only knob and the answer is closedform. A spline with G knot intervals and degree k has G + k basis functions. The squared approximation bias decays as G−2β with β the smoothness of the target. The variance grows as (G + k)/n, since each basis function adds one degree of freedom to be estimated from n data points. Balancing the two forces gives the optimal resolution: G+k A G−2β = τ | {z } | {zn } bias2 (falling)
=⇒
G⋆ ∝ n1/(2β+1) .
variance (rising)
The exponent depends only on smoothness. The two constants A (bias scale) and τ (noise/variance scale) can be identified from fits at two different resolutions. For a single coordinate, the closed-form solution replaces search entirely. Multiple coordinates: the naive tensor product fails. The naive multivariate spline is the tensor product, the Cartesian product of d univariate bases. It produces (G + k)d basis functions, a count that explodes even in moderate dimension: d=1 d
(G + k) at G = 5, k = 3
d=2
8
64
d=5
d = 10
d = 20
32,768
9
∼ 1018
∼ 10
With 109 or more basis functions the variance term dominates, and any closed-form estimate of the optimum loses its practical edge over cross-validation. This is the curse of dimensionality. The escape is structural: most regression problems of practical interest do not need the full tensor product. Low-order structure: the escape.
Consider a function of ten coordinates that is purely additive,
f (x1 , . . . , x10 ) = f1 (x1 ) + f2 (x2 ) + · · · + f10 (x10 ). The natural spline basis is ten univariate blocks rather than a tensor product. The total basis size is 10(G + k) instead of (G + k)10 . At G = 5 with cubic splines that is 80 basis functions instead of about a billion. A function with sparse pairwise interactions has univariate terms plus a few bivariate terms fij (xi , xj ). If s pairs are active, the basis adds s(G + k)2 tensor-product columns for those pairs, giving a total of roughly d(G + k) + s(G + k)2 . While s is moderate, the basis remains tractable. In both cases the variance is governed by the number of active low-order components, not by the input dimension. Balancing bias against variance produces a different scaling law per structure family: Structure Full tensor product Additive Sparse pairwise
Variance scales with
Optimal G⋆ scales as
(G + k)d /n d(G + k)/n s(G + k)2 /n
(n)1/(2β+d) (n/d)1/(2β+1) (n/s)1/(2β+2)
For the full tensor product the exponent 1/(2β + d) shrinks with d and quickly becomes useless. For the additive family the exponent is 1/(2β + 1), the same as in one coordinate. For sparse pairwise families it is 1/(2β + 2). In neither case does the input dimension d appear in the exponent. The quantity that matters is the effective density: n/d for additive models, n/s for sparse pairwise models with s active pairs. From law to algorithm. The effective-density law fixes the shape of the error curve. Its two unknown constants A (bias scale) and τ (noise/variance scale) must still be identified from data. KORE fits at two pilot resolutions per family. Each fit yields an exact LOO error through the PRESS identity, with no refitting. Two LOO measurements at two resolutions form a 2 × 2 leverage-calibrated linear system in (A, τ ), solvable b † is then the unique positive root of the derivative in closed form. The continuous plug-in resolution G of the fitted excess-risk curve, and a small symmetric integer neighborhood around the rounded plug-in serves as a finite-sample certificate. The entire procedure costs about a dozen model fits. Exhaustive 3-fold cross-validation on the candidate grids (20 additive resolutions and 10 pairwise resolutions, three folds each, plus a final refit) costs (20 + 10) × 3 + 1 = 91. 3
1.2 Err(G) bias 2 variance KORE: closed-form G ⋆ CV: search all G
test error
1.0 0.8 0.6 0.4 0.2 2.5
5.0
7.5
10.0
12.5
15.0
17.5
20.0
resolution G
Figure 1. Law-driven versus search-driven resolution selection. Cross-validation evaluates every grid candidate (clay dots). KORE fits two pilot resolutions, identifies the bias and variance scale constants, and reports the closed-form optimum (navy star). The dashed curves show how squared bias (falling) and variance (rising) compose the U-shaped error curve.
Contributions.
The paper makes four contributions.
1. An intrinsic-order test-error law for structured spline regression (Proposition 1). The optimal resolution scales with the count of active highest-order interactions, not exponentially with input dimension (Theorem 1). The bias term of this law is the Kolmogorov n-width of the smoothness class, so the spline family is approximation-theoretically width-optimal and the selected resolution attains the Stone minimax rate (Remarks 4 and 5). 2. A search-free plug-in algorithm, KORE, that estimates (Af , τf ) from two pilot fits via a leverageb † , and certifies it with a calibrated 2 × 2 system, evaluates the closed-form plug-in resolution G f small LOO neighborhood (Algorithm 1). 3. A consistency theorem for the plug-in (Theorem 2): the estimated constants and the rounded plug-in converge to their population counterparts as n grows, and the integer selector matches exhaustive risk minimization once the oracle margin exceeds the plug-in error. An empirical companion (Section 4.3) confirms the convergence at d = 20. 4. An empirical accuracy-compute study up to d = 80. Across six controlled frontier tasks, KORE matches the entire classical full-grid ladder of 3-fold CV, GCV, Mallows’ Cp , AIC, and BIC at 8.1× fewer model fits, and the effective-density collapse holds at every dimension and density tested (Section 4.2, Section 4.4). On nine smooth named benchmarks the geometric-mean RMSE ratio against CV is 0.918 at 8.7× fewer fits (Section 4.5). Code is available at https://github.com/bay-yearick-lab/kore.
2
F ROM BIAS - VARIANCE TO SEARCH - FREE RESOLUTION LAWS
The title of the paper asks for more than a faster heuristic. It asks when a hyperparameter can be treated as a statistical estimand: a quantity determined by the training distribution, the model family, and the sample size, rather than by a blind validation sweep. For spline regression, the estimand is the resolution n o G•f := arg min Errf (G) − σ 2 , (2) G∈Gfstab
where f denotes the chosen structured spline family and Gfstab is the set of resolutions whose design matrices are numerically estimable. Cross-validation estimates (2) by evaluating every candidate. KORE estimates the constants in the risk curve and then solves for (2). This section gives the mathematical backbone for that claim. 4
2.1
S TRUCTURED SPLINE FAMILIES
The data are independent training samples (xi , yi ) with xi ∈ [0, 1]d and E[ε2i | xi ] = σ 2 .
E[εi | xi ] = 0,
yi = f (xi ) + εi ,
(3)
The analysis is stated for random design with a density bounded above and below on [0, 1]d ; the same formulas apply conditionally on a fixed design whenever the empirical Gram matrices of the spline bases have eigenvalues bounded away from zero and infinity on the stable range of G. This is the usual condition under which least-squares spline estimates behave like their population projections. Low-order ANOVA structure. The escape from the curse of dimensionality is structural rather than numerical. The signal f is assumed to decompose into centered low-order components, Z r X X f (x) = f0 + fu (xu ), fu (zu ) dzj = 0 for every j ∈ u, (4) t=1 u∈Ut
where u is a subset of coordinates, xu is the corresponding subvector, Ut is the set of active t-way components, st = |Ut |, and r is the highest active interaction order. The centering constraints make the decomposition identifiable: constants live in f0 , one-dimensional effects live in the main-effect blocks, and interactions are not allowed to re-create lower-order terms. The two families used throughout the paper are special cases of (4). Additive models have r = 1 and s1 = d, f (x) = f0 +
d X
fj (xj ),
j=1
while sparse pairwise models have r = 2 and an interaction graph E with s = |E| active edges, f (x) = f0 +
d X
fj (xj ) +
j=1
X
fij (xi , xj ).
(i,j)∈E
The resolution-indexed spaces. Fix spline degree k and let G be the number of knot intervals per coordinate. A univariate B-spline block has G + k columns before centering and m(G) = G + k − 1
(5)
centered columns after the constant direction is removed. An additive design therefore has one intercept plus d centered univariate blocks, Badd (G) = 1 | B1 | · · · | Bd , padd (G) = 1 + d m(G) = d(G + k) − (d − 1). A sparse pairwise design adds one row-wise tensor block for each active edge, Bpair (G) = 1 | B1 | · · · | Bd | (Bi ⊗ Bj )(i,j)∈E , ppair (G) = 1 + d m(G) + s m(G)2 . In general, if the highest active interaction order is r, the structured basis dimension is pr (G) = 1 +
r X
st m(G)t .
(6)
t=1
This count is the first key mathematical fact: the ambient dimension d enters only through the number of active components st . The dense tensor-product count (G + k)d has disappeared. 2.2
T HE RISK CURVE
Let Sr (G) be the structured spline space with dimension pr (G), and let fˆG be the ridge-stabilized leastsquares fit in that space, with ridge small enough that it only regularizes the linear solve. The expected test MSE is h i Err(G) = E(X,Y ),D (Y − fˆG (X))2 . (7) 5
Since Y = f (X) + ε and the test noise is independent of the training data, h i Err(G) = σ 2 + ED ∥fˆG − f ∥2L2 (PX ) .
(8)
The selection-relevant part is therefore g Err(G) := Err(G) − σ 2 .
(9)
The irreducible floor σ 2 matters for estimating the curve from data, but it does not affect the minimizer once the excess curve is known. Bias: approximation by a resolution-G space. Let ΠG f be the L2 (PX ) projection of f onto Sr (G). Assume each active component fu has smoothness β ≤ k + 1 on its |u|-dimensional domain. Tensor-product B-spline approximation theory gives ∥fu − ΠG,u fu ∥2L2 ≤ Cu,k,β G−2β ,
(10)
because the mesh width in every active coordinate is h = 1/G (de Boor, 2001; Schumaker, 2007; Stone, 1982). Summing over finitely many active components yields ∥f − ΠG f ∥2L2 (PX ) ≤ Ar G−2β ,
(11)
where Ar > 0 is a bias-scale constant determined by the component functions and the design distribution, but not by G or n. If the leading approximation term is nonzero, (11) is sharp up to lower-order terms. The bias rate is a Kolmogorov n-width. The exponent in (10) is not an artifact of the spline construction; it is the best rate any linear method of the same dimension could achieve. The Kolmogorov n-width dn (A, L2 ) =
sup inf ∥f − g∥L2
inf
L⊂L2 f ∈A g∈L dim L=n
(12)
is the smallest worst-case error attainable by any n-dimensional linear subspace over a function class A (Kolmogorov, 1936; Pinkus, 1985). For a univariate smoothness-β ball it decays as dn ≍ n−β , and univariate polynomial spline spaces of order k + 1 ≥ β are not merely rate-optimal but exactly optimal for this L2 width (Melkman & Micchelli, 1978; Pinkus, 1985). Each active component in (10) is approximated on a tensor product of |u| such width-optimal univariate blocks at per-coordinate resolution G, so the per-coordinate width is dm(G) ≍ m(G)−β ≍ G−β and the component’s squared error is ≍ G−2β . Summing over the sr active components, the bias scale Ar G−2β in (11) is, up to constants, the squared Kolmogorov width of the active-component family at resolution G. The resolution therefore indexes a Kolmogorov-width-optimal approximation family, and the standing choice β = k + 1 (Section 3) is precisely the spline order at which this width-optimality holds. This is the sense in which the closed-form selector is Kolmogorov-optimal: it never leaves the family of subspaces that realize the best possible linear approximation rate, and it tunes only where on that family the bias-variance balance sits. Variance: estimating a pr (G)-dimensional projection. On the stable range pr (G)/n ≤ κ < 1, least squares in a pr (G)-dimensional linear space has integrated estimation variance of order pr (G)/n. More precisely, under the Gram stability condition stated above, h i pr (G) pr (G) 2 ˆ ED ∥fG − ΠG f ∥L2 (PX ) = Br +o , (13) n n where Br > 0 is a variance-scale constant. In the ideal orthonormal homoskedastic case, Br = σ 2 ; for non-orthogonal but well-conditioned spline designs, Br absorbs the design correction. The important point is that the G-dependence is explicit because pr (G) is explicit. Proposition 1 (Intrinsic-order test-error law). Under the structured model (4), component smoothness β ≤ k + 1, and stable spline Gram matrices on the candidate range, the selection-relevant risk obeys Pr 1 + t=1 st m(G)t −2β g Err(G) = Ar G + Br + remn (G), (14) n with remn (G) = o{G−2β + pr (G)/n} uniformly over stable resolutions. Thus the bias decreases as a known power of G, the variance increases according to a known basis-count polynomial, and only the two scale constants Ar and Br are unknown. 6
Proposition 1 is the formal reason this hyperparameter is solvable. Once Ar and Br are known or consistently estimated, there is no statistical reason to train at every G: the entire U-shaped curve is determined. 2.3
T HE CLOSED - FORM OPTIMIZER
For the moment, treat G as a positive real number. The continuous excess-risk proxy from Proposition 1 is Pr 1 + t=1 st m(G)t Rr (G; A, B) = A G−2β + B . (15) n Its derivative is r ∂Rr BX = −2βA G−(2β+1) + t st m(G)t−1 . (16) ∂G n t=1 The second derivative is positive for G > 0, so the derivative in (16) crosses zero at most once. Therefore the continuous optimizer is the unique solution of r
2βA G−(2β+1) =
BX t st m(G)t−1 . n t=1
(17)
This is already a search-free rule: solve one scalar equation, then round to the nearest stable integer resolution. To expose the scaling, keep only the highest-order variance term near the optimum. Since m(G) = G+O(1), r
gr (G) ≈ A G−2β + B sr G . Err n Differentiating this dominant law gives
(18)
d g rBsr r−1 Errr (G) ≈ −2βA G−(2β+1) + G = 0. (19) dG n Collecting powers of G yields 2βA n G2β+r = . rB sr Theorem 1 (Intrinsic-order resolution law). If the highest-order active term dominates the variance near the optimum and the leading bias constant is nonzero, then the optimal continuous resolution satisfies 1/(2β+r) 2βA n G⋆r = {1 + o(1)}. (20) rB sr Equivalently, G⋆r ≍ (n/sr )1/(2β+r) . The input dimension affects the optimizer only through the count sr of active highest-order components. The corresponding optimal excess MSE scales as (n/sr )−2β/(2β+r) , so the optimal RMSE scales as (n/sr )−β/(2β+r) . This is the slope tested in the law-collapse experiment. Corollary 1 (Additive and sparse pairwise laws). For additive models (r = 1, s1 = d), n 1/(2β+1) n −β/(2β+1) G⋆add ≍ , RMSE⋆add ≍ . (21) d d For sparse pairwise models (r = 2, s2 = s), n 1/(2β+2) n −β/(2β+2) G⋆pair ≍ , RMSE⋆pair ≍ . (22) s s For cubic splines in the classical smooth regime, β = k + 1 = 4, giving resolution exponents 1/9 and 1/10, and RMSE exponents −4/9 and −4/10. Effective density, not ambient dimension. Suppose an additive target has d = 40 and n = 4,800, so ρ = n/d = 120. Doubling both to d = 80 and n = 9,600 leaves ρ unchanged. Corollary 1 predicts the same optimal resolution and the same test RMSE up to constants. In Section 4.2, the selected resolutions at ρ = 120 all collapse to G⋆ = 15 across d ∈ {10, 20, 40, 80}, exactly as the law predicts. 7
Table 1. Effective-density laws for the two structure families. The resolution exponent depends on interaction order r, not on the input dimension d. The RMSE slope is the exponent tested in Figure 2. Family Additive Sparse pairwise
3
Basis size p(G)
Effective density
G⋆ exponent
RMSE exponent
1 + d m(G) 1 + d m(G) + s m(G)2
ρ = n/d ρ = n/s
1/(2β + 1) 1/(2β + 2)
−β/(2β + 1) −β/(2β + 2)
T HE KORE ALGORITHM
Section 2 turns resolution selection into a two-constant estimation problem. For a fixed family f , the excess-risk curve has the form Rf (G) = Af G−2β + Bf νf (G),
νf (G) :=
pf (G) , n
up to lower-order terms. KORE estimates the constants from two pilot fits, solves the resulting curve for its minimizer, rounds to a feasible integer, and uses a tiny local leave-one-out certificate. The expensive step, training a model at every candidate resolution, never occurs. Estimating test error without a validation split. For any fixed G, the spline estimator is a linear smoother: ŷ = HG y. Therefore all leave-one-out residuals are available from a single fit by the PRESS identity (Allen, 1974), 2 n 1X yi − ŷi LOO(G) = . (23) n i=1 1 − HG,ii No refitting over held-out folds is required. This is the computational hinge of the method: one trained spline gives one out-of-sample risk measurement. Removing the noise-floor problem. A subtle but important point is that LOO(G) estimates the full test MSE, while the optimizer in Section 2 uses the excess risk above the irreducible floor. Naively fitting LOO(G) ≈ Af G−2β + Bf pf (G)/n would silently force the noise floor into the bias and variance constants. The corrected two-pilot system uses the fact that, for a homoskedastic linear smoother, the same noise level that creates the irreducible floor also creates the variance penalty. Let ϕ(G) = G−2β ,
νf (G) =
pf (G) , n
ℓf (G) =
1 . 1 − νf (G)
(24)
The factor ℓf (G) is the average-leverage approximation to the PRESS inflation; the implementation uses ℓf (G) = n/(n − pf (G)) and only admits pilots with pf (G) < 0.45n. Since ℓf (G) = 1 + νf (G) + O(νf (G)2 ) on the stable range, estimating the pair (Af , τf ) from LOOf (G) ≈ Af ϕ(G) + τf ℓf (G)
(25)
recovers both the bias scale and the noise/variance scale. The excess curve to minimize is then bf (G) = A bf ϕ(G) + τbf νf (G). R
(26)
Thus the noise floor is not ignored; it is estimated and then removed from the part of the curve that determines G⋆ . Two fits, two equations, two unknowns. Choose two pilot resolutions Ga < Gb , fit the family at each, and compute (23). Define ϕj = ϕ(Gj ) and ℓj = ℓf (Gj ) for j ∈ {a, b}. The pilot equations are ϕa ϕb
ℓa ℓb
bf LOOf (Ga ) A = . LOOf (Gb ) τbf 8
(27)
When the determinant Df = ϕa ℓb − ϕb ℓa is nonzero, the solution is explicit: bf = LOOf (Ga )ℓb − LOOf (Gb )ℓa , A Df
(28)
ϕa LOOf (Gb ) − ϕb LOOf (Ga ) . Df
(29)
τbf =
The pilots are deliberately separated: Ga is coarse and bias-dominated, while Gb is closer to the largest stable resolution and variance-dominated. This keeps Df away from zero and makes the constants identifiable. Solving, not sweeping.
bf , τbf ) in hand, KORE solves the fitted excess curve With (A pf (G) † −2β b b Gf = arg min Af G + τbf . n G>0
(30)
bf and τbf . For the full polynomial pf (G), this is the unique positive root of (17) with A and B replaced by A Under highest-order dominance, the root is the closed form !1/(2β+rf ) bf n 2β A † b = . (31) G f rf τbf srf b † to the stable range and checking the nearest integer The selected integer resolution is obtained by clipping G f neighbors. This final check is a certificate against finite-sample discretization error, not a grid search over the hyperparameter. Theorem 2 (Search-free plug-in guarantee). Fix a structured family f satisfying Proposition 1. Suppose the pilot matrix in (27) has determinant bounded away from zero after scaling, and suppose the two PRESS measurements satisfy the pilot law (25) with errors op {Af ϕ(Gj ) + τf ℓf (Gj )} for j ∈ {a, b}. Then bf /Af →p 1, A
τbf /τf →p 1,
b † /G† →p 1. G f f
(32)
A finite-sample rate is available in Appendix A.5: under sub-Gaussian noise with proxy σ 2 and the well√ b † − G• |/G• = Op (n−1/2 log n), and the rounded plug-in matches conditioned pilot pair of Lemma 1, |G f f f the integer oracle once n ≳ (G•f )2 log n. If the integer oracle has a positive margin, ∆f =
min
G∈Gfstab : G̸=G• f
{Rf (G) − Rf (G•f )} > 0,
and the pilot errors are op (∆f ) after propagation through (26), the rounded plug-in rule selects the same integer resolution as exhaustive risk minimization with probability tending to one. The theorem makes precise what “the model knows the hyperparameter” means in this setting. The training data identify the two constants of a known risk law, and the minimizer is then a plug-in statistic. Crossvalidation is no longer the definition of the hyperparameter; it is only a baseline against which to check the statistic. Algorithm 1 KORE for additive versus sparse pairwise spline selection Require: training data (X, y), spline degree k, smoothness index β, additive family, sparse pairwise family, optional interaction graph E 1: for f ∈ {add, pair} do 2: choose stable pilots Ga < Gb and fit the spline family at both resolutions 3: compute LOOf (Ga ) and LOOf (Gb ) by the PRESS identity (23) bf , τbf ) 4: solve the 2 × 2 system (27) for (A b† 5: solve the fitted risk equation (30) for G f 6: round to the nearest stable integer resolution and certify with a small LOO neighborhood 7: keep the certified resolution and LOO score for family f 8: end for 9: return the family with lower certified leave-one-out MSE
9
The selection-cost difference is structural. Exhaustive 3-fold cross-validation fits every candidate in both grids Gadd and Gpair three times, totaling (|Gadd | + |Gpair |) × 3 + 1 fits. GCV avoids the fold loop but still scores every candidate. KORE uses two pilot fits per family plus a small certificate neighborhood, so its model-training cost is nearly independent of the size of any candidate grid. The concrete savings are reported in Section 4.4. Concrete parameter settings used throughout this paper. The spline degree is k = 3 (cubic B-splines), which fixes the classical smooth-regime exponent at β = k +1 = 4 and therefore the bias exponent 2β = 8 in Proposition 1. The additive basis dimension is padd (G) = d(G + k) − (d − 1), and the sparse pairwise basis dimension with an edge set E of size s = |E| is ppair (G, s) = 1 + dm + sm2 with m = G + k − 1. The pilot eff resolutions are Ga = 1 and Gb = ⌊0.75 Geff max ⌋, where Gmax is the largest grid value satisfying the stability rule pf (G) < 0.45n. Placing the upper pilot near the variance-dominated end of the stable range keeps the determinant of (27) well away from zero and makes the pilot itself a meaningful probe of the high-resolution regime where smooth targets tend to peak. Each linear system is solved with a Tikhonov ridge of 10−8 on b f − 3, . . . , G b f + 3} ∩ Gadd for b † is {G the diagonal. The certificate neighborhood around the predicted G f b b b the additive family and {Gf − 1, Gf , Gf + 1} ∩ Gpair for the sparse pairwise family. As a safeguard, when the upper pilot Gb improves on the whole neighborhood the bias-variance optimum lies in the gap between them; because the leave-one-out risk is unimodal on the stable range, a logarithmic-cost bracketed search b f − 3, Gb ] recovers it, so an interior optimum is never missed while the fit count stays close to the over [G neighborhood size. These settings are identical across every experiment in the paper.
4
E XPERIMENTS
The experimental program tests six empirical claims, each addressing a question a rigorous reader would ask. (i) Section 4.2 verifies the effective-density collapse on controlled additive and sparse pairwise families across four dimensions and a factor-of-twenty-four sweep in ρ. (ii) Section 4.4 compares KORE against the full classical selection ladder of 3-fold CV, GCV, Mallows’ Cp , AIC, and BIC on a six-task frontier. (iii) Section 4.5 runs the same ladder on nine named benchmark equations drawn from the nonparametricregression literature. (iv) Section 4.6 reports the three boundary benchmarks where a single global resolution stops being a faithful inductive bias. (v) Section 4.7 validates a run-time diagnostic that separates the signal-rich regime from the noise-dominated regime on an independent noise sweep. (vi) Section 4.3 checks the search-free plug-in guarantee of Theorem 2 numerically across a geometric sample-size ladder. Section 4.1 fixes the estimands and baselines that every experiment shares. The complete reproduction recipe, including the master seed, fold recursion, exact target equations, benchmark suite, and every fixed numerical constant, is collected in Appendix B; nothing beyond that appendix is required to reproduce any number reported here. 4.1
E STIMANDS , BASELINES , AND THE QUANTITIES EVERY EXPERIMENT MEASURES
All experiments use the same cubic B-spline families (k = 3), the same additive candidate grid Gadd = {1, . . . , 20}, the same sparse pairwise grid Gpair = {1, . . . , 10}, and a ridge of 10−8 on every normal equation. Methods differ only in how they select the resolution. Basis dimensions. The additive basis is built from one cubic B-spline block per coordinate with bias retained and d − 1 redundant constants dropped, giving padd (G) = d(G + k) − (d − 1).
(33)
The sparse pairwise basis uses centered univariate blocks with bias disabled and adds one row-wise KhatriRao tensor block per active pair, giving ppair (G, s) = 1 + d(G + k − 1) + s(G + k − 1)2 .
(34)
These two formulas are what make the effective-density collapse visible in the first place: basis growth is linear in d for the additive family and linear in s for the sparse pairwise family, never exponential, so the dimension dependence absorbs cleanly into ρ. 10
additive targets
(b) test RMSE at KORE-selected G ⋆
test RMSE at KORE-selected G ⋆
(a)
theory ½ ¡0:44 d = 10 d = 20 d = 40 d = 80
10−2
25
26
28
27
29
sparse pairwise targets theory ½ ¡0:40 d = 10 d = 20 d = 40 d = 80
10−2
26
½ = n=d
27
28
29
½ = n=s
Figure 2. Effective-density collapse. The horizontal axis is ρ; the vertical axis is test RMSE at the KORE-selected G⋆ . Panel (a) shows additive targets with ρ = n/d and reference slope ρ−4/9 . Panel (b) shows sparse pairwise targets with ρ = n/s and reference slope ρ−4/10 . The four dimensions d ∈ {10, 20, 40, 80} collapse onto a single curve.
Estimator and leave-one-out score. All methods share the same fitting primitive. Given a design matrix B, the solver forms β̂ = (B ⊤ B + 10−8 I)−1 B ⊤ y once, then scores the fit with the exact leave-one-out identity (23). The key computational consequence is that one model fit yields one exact LOO score without any fold loop, so KORE pays only for two pilot fits plus a small local refinement. Exhaustive CV, by contrast, retrains on three folds per candidate grid point. A self-contained derivation of the PRESS identity from first principles is in Appendix A.8. The classical selection ladder. The baselines span the full classical menu of resolution selectors. Exhaustive 3-fold cross-validation (Allen, 1974; Stone, 1974) is the accuracy bar that practitioners trust. Generalized cross-validation (GCV) (Craven & Wahba, 1979; Golub et al., 1979) scores every feasible candidate in the grid with GCV(G) =
RSS(G)/n , max(1 − pf (G)/n, 0.01)2
(35)
skipping candidates with pf (G) ≥ 0.9n for stability. Mallows’ Cp (Mallows, 1973), AIC (Akaike, 1974), and BIC (Schwarz, 1978) are information criteria derived from a Gaussian likelihood at the fitted residual variance, each evaluated over the full feasible grid in its standard form; explicit expressions are collected in Appendix C.5. Together these four closed-form criteria cover the standard ways to replace exhaustive CV with a cheaper full-grid pass, so any speedup KORE achieves on top of all four must come from not evaluating the grid at all. The main-paper pairwise experiments isolate resolution selection by supplying the active interaction graph; the corresponding graph-discovery experiment is reported in Appendix C. Reproducibility at a glance. Every result in Section 4 is averaged over five seeds derived deterministically from a single master seed. The synthetic experiments use 3% training noise and 2,000 noise-free test points; the benchmark suite uses 1% training noise and 3,000 noise-free test points; the applicability sweep uses a variable noise level and 3,000 test points; the consistency experiment uses 10% training noise and 3,000 test points across 20 seeds and seven sample sizes. Appendix B gives the full seed-folding recursion (Table 5), the two controlled target families fadd and fpair with their random-draw recipes, the twelve-equation benchmark suite (Table 7), and the table of fixed numerical constants (Table 6). 4.2
S ELECTED RESOLUTION AND TEST ERROR COLLAPSE BY EFFECTIVE DENSITY
Setup. Corollary 1 predicts that once ρ is fixed, the input dimension d carries no extra information about the optimal resolution or the resulting test error. The collapse experiment tests that prediction directly on the exact additive and sparse pairwise target families fadd and fpair defined in Appendix B. Additive targets sweep ρ = n/d over {30, 45, 60, 90, 120, 180, 240, 360, 480, 720} at each of d ∈ {10, 20, 40, 80}, subject 11
1.2
c A f =Af
1.0
0.8
0.6
0 30
0 60
00
12
00
24
00
48
00
96
0 20
19
sample size n
bf /Af versus Figure 3. Bias-scale recovery as predicted by Theorem 2: median and interquartile band of the ratio A sample size n along a geometric ladder, with the population value 1 marked.
to n = ρd ≤ 60,000. Sparse pairwise targets sweep ρ = n/s over {60, 90, 120, 180, 240, 360, 480, 720} at the same four dimensions with s = d/2 active pairs, subject to n = ρs ≤ 60,000. Every cell uses the seed rule of Table 5, 3% training noise, and a 2,000-point noise-free test set. In each cell, KORE selects G⋆ , the structured spline is fit once at that resolution, and the resulting test RMSE is plotted against ρ. Theory prediction. Corollary 1 forecasts two things simultaneously. First, the four colored curves in each panel should collapse onto one another, because d should matter only through ρ. Second, the collapsed curve should follow the predicted power law: ρ−4/9 for additive structure and ρ−4/10 for sparse pairwise structure. Verifying the slope on the continuous RMSE axis is the tightest available check of the full scaling law because the RMSE exponent and the G⋆ exponent are tied by the same corollary. Observation. Collapse holds in both panels of Figure 2. Panel (a) shows the four dimensions landing on a common curve that follows the ρ−4/9 reference across the entire tested range; the values at ρ = 120 lie between 0.0077 and 0.0081 regardless of d. Panel (b) shows the same pattern for sparse pairwise targets along the ρ−4/10 reference. Input dimension is invisible once ρ is fixed, and the empirical slopes match the predicted exponents. The closed-form G⋆ chosen by KORE therefore delivers the error decay that the scaling law predicts, on both low-order families and at every dimension tested. 4.3
P LUG - IN CONSISTENCY OF THE CLOSED - FORM SELECTOR
bf /Af and τbf /τf converge to 1 in probability and Setup. Theorem 2 states that, under standard regularity, A † b that the plug-in continuous optimizer Gf converges to the population target G•f . The consistency experiment tests these three statements directly. An additive target in d = 20 is the test case, with n swept geometrically over a ladder up to n = 19,200. At the largest n, the pilot solution provides the anchor pair (A, τ ) and the bf , τbf , G b † ) on 20 seeds, and the resulting ratios are anchored optimizer G•f . Every smaller n then yields (A f plotted. Theory prediction. Theorem 2 predicts that the two constant ratios tend to 1 as n grows and that the b † tracks G• to within a vanishing margin. A small bias in the constants at small n is plug-in optimizer G f f expected, since the pilot fits have finite variance, but the bias must shrink and the variance must contract along the ladder. Observation. Figures 3 and 4 show that the two ratios settle onto 1 along the ladder, with a transient deviation at the smallest n that contracts rapidly and is negligible by n = 1,200, consistent with the predicted b † landing on G• to within the certificate radius along the in-probability convergence. Figure 5 shows G f f ladder, confirming the search-free plug-in guarantee at the resolution that the algorithm actually uses. 12
¿bf =¿f
101
100 0
0
30
00
60
12
00
24
00
48
00
0
20
96
19
sample size n
Figure 4. Noise-scale recovery as predicted by Theorem 2: median and interquartile band of τbf /τf versus sample size n along the same geometric ladder, with the population value 1 marked. c † § 1¾ G f
c † mean G f
spline resolution
4.0
Gf² population
3.5 3.0 2.5 2.0 1.5 30
0
60
0
00
12
00
24
00
48
00
96
0
20
19
sample size n b † versus the anchored population target G•f along the sample-size ladder. Figure 5. Plug-in continuous optimizer G f Shaded band: ±1σ across 20 seeds; the diamond marker is the anchored target at the largest n.
4.4
C LOSED - FORM SELECTION ON THE ACCURACY- COMPUTE FRONTIER
Setup. The frontier experiment tests whether the closed-form law translates into a practical accuracycompute win. The benchmarks are three additive targets at d ∈ {10, 20, 40} with ρ = n/d = 120, and three sparse pairwise targets at d ∈ {10, 20, 40} with ρ = n/s = 240, giving n ∈ {1,200, 2,400, 4,800} in both families. These densities sit inside the regime the law is meant to serve: large enough for a clear interior optimum and small enough for selection cost to matter. Each cell is averaged over five seeds. Theory prediction. If the law captures the right resolution scale, the closed-form plug-in should land near the same optimum exhaustive search reaches, without paying for the full grid. Figure 6 should therefore place KORE to the left of exhaustive CV and the four classical full-grid criteria; Figure 7 should place the per-task RMSE ratio at or below 1 on most tasks, with any deviations small; Figure 8 should show KORE as the Pareto-dominant point on the accuracy versus cost plane. Observation. Figure 6 shows that exhaustive 3-fold CV requires about 91 fits across the two grids and three folds. GCV, Cp , AIC, and BIC all use 29 fits because they share a full-grid evaluation protocol and 13
Add d = 10
Add d = 20
Add d = 40
Pair d = 10
Pair d = 20 3-fold CV GCV / Cp / AIC / BIC KORE
Pair d = 40 20
40
60
80
model fits for selection
Figure 6. Selection cost across six controlled tasks (three additive, three sparse pairwise). Markers report the number of model fits used by KORE, the four full-grid criteria (GCV, Cp , AIC, BIC), and exhaustive 3-fold cross-validation.
Add d = 10
Add d = 20
Add d = 40
Pair d = 10
Pair d = 20
Pair d = 40 0.96
0.97
0.98
0.99
1.00
1.01
1.02
1.03
1.04
KORE / CV RMSE
Figure 7. Per-task accuracy parity on the same six tasks: ratio of KORE test RMSE to exhaustive 3-fold cross-validation, sorted by family. Bars at or below 1.0 favor the closed-form selector.
differ only in scoring formula. KORE uses 9 to 14 fits because the law identifies the right neighborhood before any refinement begins. Figure 7 shows that every per-task RMSE ratio lands at or below the 1.0 line: the closed-form selector matches or beats exhaustive search on every task. The three sparse pairwise tasks tie (1.000), and the three additive tasks favor the closed form slightly (0.979, 0.998, 0.996) because the continuous law occasionally lands between the discrete CV grid points. Figure 8 summarizes the geometric-mean accuracy versus CV across all six tasks: KORE posts 0.995 at 8.1× fewer fits, while the four classical full-grid criteria cluster near parity at 3.1× fewer (GCV and Cp tie at 1.003, AIC at 1.026, 14
KORE
GCV
Cp
BIC
AIC
1.0
1.1
1.2
1.3
1.4
1.5
geometric-mean RMSE / CV
Figure 8. Geometric-mean RMSE ratio versus exhaustive cross-validation, aggregated across all six controlled tasks, for the five selectors compared in this section. Bars to the left of 1.0 match or beat exhaustive search. (a) smooth low-order benchmarks
(b) fit-count reduction
Nguyen-1
(c)
method summary
KORE
Nguyen-9 (2D add) Nguyen-7
GCV
Nguyen-5
Cp
Nguyen-4 SparseAdd-20D
BIC
Friedman-1 (5D) SparsePair-10D
AIC
Nguyen-10 (2D int)
0.00
0.25
0.50
0.75
1.00
KORE / CV RMSE
0
5
10
CV fits / KORE fits
0.90
0.95
1.00
GM RMSE / CV
Figure 9. Nine law-aligned benchmark equations. Panel (a) gives the per-task RMSE ratio of KORE against 3-fold CV. Panel (b) gives the corresponding fit-count reduction (CV fits divided by KORE fits). Panel (c) gives the geometric-mean RMSE ratio against CV across the nine tasks for the five selectors.
BIC at 1.058). KORE is Pareto-dominant on this frontier, and it dominates because it replaces the full-grid pass with a two-fit closed-form solve. 4.5
L AW- ALIGNED BENCHMARK EQUATIONS
Setup. The frontier tasks were constructed to satisfy the law tightly. The benchmark suite asks the more demanding question of whether the same advantage holds on named equations from the nonparametricregression literature. Appendix Table 7 lists the full twelve-equation suite. This subsection isolates the nine equations whose dominant structure is smooth and low-order, the regime the theory claims to cover: Nguyen-1, Nguyen-4, Nguyen-5, Nguyen-7, Nguyen-9 (2D add), Nguyen-10 (2D int), Friedman-1 (5D), 15
CV KORE GCV
Nguyen-1
Nguyen-9 (2D add)
Nguyen-7
Nguyen-5
Nguyen-4
SparseAdd-20D
Friedman-1 (5D)
SparsePair-10D
Nguyen-10 (2D int)
Friedman-2 (4D)
Franke (2D)
Oscillator
10−3
10−2
10−1
mean test RMSE
Figure 10. Forest plot of mean test RMSE per benchmark across the full 12-equation suite, ranked by KORE / CV ratio. Markers compare KORE (closed-form), exhaustive 3-fold cross-validation, and GCV. Raw RMSE ratios are tabulated in Table 2.
SparseAdd-20D, and SparsePair-10D. All nine are run with 1% training noise and 3,000 test points using the protocol of Table 5. Theory prediction. On smooth low-order targets the closed-form selector should remain near parity with exhaustive CV and the four full-grid criteria, and may beat grid search occasionally because the fitted law is continuous while the CV grid is discrete. The cost advantage should persist because its source is algorithmic rather than dataset-specific. Observation. Panel (a) of Figure 9 shows seven benchmarks below 1.0 (Nguyen-1 at 0.606, Nguyen-9 (2D add) at 0.831, Nguyen-7 at 0.913, Nguyen-5 at 0.928, Nguyen-4 at 0.957, SparseAdd-20D at 0.981): the closed-form law finds a better resolution than the discrete CV grid. Friedman-1 (5D) and SparsePair-10D land just above (1.029 and 1.044). One benchmark, Nguyen-10 (2D int), is slightly above (1.079). The geometric-mean ratio across the nine tasks is 0.918. Panel (b) shows fit-count reduction between 6.1× and 12.2× for every benchmark. Panel (c) summarizes against the four classical criteria: GCV and Cp tie at 0.930, BIC at 0.936, AIC at 0.949. KORE is the best of the five at 0.918, and it reaches that ceiling at roughly 8.7× fewer fits versus 2.9× for the classical full-grid baselines. The three remaining benchmarks, which probe the boundary of the single-resolution law, are reported in the next subsection. 4.6
B OUNDARY CASES AND SCOPE OF THE SINGLE - RESOLUTION LAW
Setup. The full benchmark suite is reported in the main text rather than partitioned across the paper and the appendix. The three additional tasks are Franke (localized two-dimensional surface structure), Friedman-2 16
Table 2. Full benchmark results. KORE RMSE and CV RMSE are the test root-mean-square errors for KORE and exhaustive 3-fold cross-validation. The ratio column gives KORE’s RMSE divided by CV’s. The final column gives CV’s fit count divided by KORE’s, the fit-cost reduction. Equation
d
n
KORE RMSE
CV RMSE
KORE/CV
CV/KORE fits
Nguyen-1 Nguyen-4 Nguyen-5 Nguyen-7 Nguyen-9 (2D add) Nguyen-10 (2D int) Friedman-1 (5D) Friedman-2 (4D) Franke (2D) Oscillator SparseAdd-20D SparsePair-10D
1 1 1 1 2 2 5 4 2 1 20 10
500 500 500 500 1000 1000 2000 2000 1000 500 3000 2000
0.000759 0.002063 0.000252 0.001108 0.000400 0.000665 0.009275 0.084253 0.004756 0.005638 0.003199 0.001722
0.001252 0.002156 0.000271 0.001214 0.000482 0.000617 0.009015 0.077045 0.002295 0.001641 0.003262 0.001649
0.606 0.957 0.928 0.913 0.831 1.079 1.029 1.094 2.073 3.436 0.981 1.044
12.2× 7.1× 6.8× 10.2× 9.1× 11.4× 9.9× 7.6× 7.6× 6.1× 6.1× 8.0×
(strongly coupled rational dependence), and the Oscillator (rapidly oscillating one-dimensional signal). These are not random hard cases chosen after the fact. They are precisely the kinds of functions one would expect to challenge a selector built around a single global resolution and low-order smooth structure. Theory prediction. The theory does not predict uniform performance on every smooth regression problem. Franke has localized spatial heterogeneity, Friedman-2 couples variables through a non-separable rational expression, and the Oscillator compresses very different length scales into one coordinate. In such cases exhaustive search over the same restricted model family may still do somewhat better, and the right extension is not a more aggressive global search but a more flexible family, for example a spatially adaptive or coordinate-wise resolution. Observation. That is what Table 2 and Figure 10 show. Friedman-2 remains close to exhaustive CV (ratio 1.094), suggesting partial compatibility with the pairwise family and some mismatch in how a single resolution allocates capacity. Franke is substantially worse (ratio 2.073), consistent with a need for spatial adaptivity across the input domain. The Oscillator is the clearest failure mode (ratio 3.436): a single global spline resolution cannot simultaneously capture rapid oscillations and exponential decay. These cases make the scope of the method more trustworthy: KORE is strong where the theory predicts strength, and when it fails, it fails for structural reasons that are visible in the target itself. 4.7
P RACTICAL GUIDANCE : AN EMPIRICAL SAFETY CHECK
Setup. KORE is designed for problems whose response is a smooth function of its inputs and depends on them through low-order structure, either additively or through a sparse set of pairwise interactions. A practitioner needs a cheap check that the current dataset lies in that regime. The procedure already computes such a check during refinement and this subsection validates it. The first component is a local shape check: does the leave-one-out curve rise on both sides of the selected resolution, confirming the expected U-shape? The second is a signal score, defined as the leave-one-out improvement of the structured fit over an intercept-only model divided by the standard error of that difference. This is a practical guardrail, not a second theorem layered on top of the scaling law. Both quantities are tested on an additive target in d = 10 with n = 200, sweeping the noise level over {5%, 10%, 25%, 50%, 100%, 200%} of the signal standard deviation, averaged over five seeds. Table 3 reports the signal score, that is, how many standard errors the structured fit improves over the predict-the-mean baseline, alongside the test RMSE ratio of KORE to that baseline. Theory prediction. As noise increases, the advantage of fitting structure should shrink smoothly. The local U-curve should become less informative, the signal score should fall toward and then below 1, and the test-error ratio relative to the intercept-only baseline should rise toward or exceed 1. The diagnostics should therefore fail exactly when the structured fit stops delivering practical value. 17
Table 3. Noise-sweep validation (d = 10, n = 200, five seeds). Signal score is the improvement of the structured fit over the predict-the-mean baseline in standard-error units; above 1 means the structured model is reliably better. KORE/Null RMSE is the corresponding test-RMSE ratio; below 1 confirms practical benefit. Noise fraction
Signal score
KORE/Null RMSE
0.05 0.10 0.20 0.40 0.80 1.20 1.60 2.00
4.14 4.07 4.00 3.58 1.27 0.55 −0.40 −0.37
0.454 0.437 0.473 0.537 0.786 0.968 1.058 1.240
Observation. At low noise (0.05 to 0.40) the signal score is well above 1 and the RMSE ratio is far below 1: the structured model captures most of the signal. At moderate noise (0.80) the score is 1.27 and the ratio is 0.79: structure is still present and worth fitting. At high noise (1.20 and above) the score drops below 1 and the ratio approaches or exceeds 1.0: the noise has overwhelmed the signal, and the predict-the-mean baseline is sufficient. The diagnostics correctly identify this transition. In practice the rule is simple: run KORE, verify that the local leave-one-out curve has a clear minimum, and verify that the signal score is above 1. If both hold, the structured fit is justified. 4.8
R EAL - WORLD VALIDATION AGAINST A TWENTY- ONE - METHOD BASELINE ROSTER
Setup. The synthetic experiments answer questions about the scaling law and the closed-form selector. The remaining question is whether KORE survives the kind of tuned-baseline comparison introduced by Grinsztajn et al. (2022) on real tabular regression data. The benchmark suite is the OpenML-CTR23 curated tabular regression collection (Fischer et al., 2023), all 35 datasets, augmented with the Combined Cycle Power Plant dataset (Tüfekçi, 2014) which is a long-standing GAM-literature classic that does not appear in CTR23. The other UCI classics that the GAM literature has used since the 1990s (Concrete, Airfoil, Wine quality red and white, California Housing, Energy Efficiency, Forest Fires, Naval Propulsion) are already in CTR23, so deduplication leaves 36 unique datasets. The pre-registered smooth-low-d subset restricts to entries with at most 30 features after one-hot encoding, the regime the bias-variance theory of Section 2 is calibrated for. The method roster is twenty-one baselines, organized by family. The linear family contains ordinary least squares, ridge, lasso, and elastic-net, each with its standard cross-validated penalty grid. The spline family contains KORE, exhaustive cross-validation over the resolution grid, the four classical information criteria (GCV, Mallows Cp , AIC, BIC), and the third-party pyGAM implementation (Servén & Brummitt, 2018) run with its own internal generalized-cross-validation lambda search. The tree-based family contains random forests, extra-trees, sklearn HistGradientBoosting, XGBoost, LightGBM, and CatBoost. The kernel family contains support-vector regression and kernel ridge, both with the radial basis function. Nearest neighbours and a small multilayer perceptron round out the comparison. Hyperparameter ranges for the tree-based, kernel, neighbour, and neural baselines are taken verbatim from Grinsztajn et al. (2022) Appendix B. All tunable methods use 20 Bayesian-optimization trials with 3-fold internal cross-validation per trial, executed by the Optuna sampler. Each cell is capped at four minutes of wall time; cells that exceed the cap are recorded as missing and excluded from the corresponding aggregate. The outer evaluation is five 80/20 train-test splits with seeds fixed across methods so every comparison sees identical data. Theory prediction. KORE is closed-form optimal within the spline-ANOVA function class and pays for two pilot fits per dataset rather than a search grid. The natural figure of merit is therefore not raw test RMSE in isolation but a metric that rewards skilful prediction per unit of compute, measured against a defensible reference. The Compute-Normalized Lift over the linear baseline (CNL) used here is n o 2 2 max 0, max(0, Rm,d,c ) − max(0, ROLS,d,c ) α , CNLα (m, d, c) = 1 + tm,d,c 18
t in seconds, α ≥ 0,
with the headline weight fixed at α = 1. The numerator is a Murphy 1988 skill score (Murphy, 1988), written against the operational reference forecast (ordinary least squares) rather than the climatology constant. OLS is the universal no-effort baseline in tabular regression; the question a practitioner cares about is whether a more elaborate method strictly out-predicts what they would have done with no thought, and how much extra compute that improvement costs. CNL satisfies four basic axioms a defensible cost-performance metric ought to satisfy: (i) no-skill predictors (R2 ≤ 0) score zero, ruling out random-prediction attacks; (ii) methods that do not strictly beat OLS score zero, ruling out the trivial copy-OLS-in-zero-time attack; (iii) the denominator is bounded below by 1, so reporting a near-zero wall time cannot inflate the score; (iv) the score is bounded in [0, 1], dimensionless, and scale-free in y. Higher CNL is better; methods are ranked across datasets by mean Friedman rank on −CNL. OLS itself sits at the floor of the rank table by construction, with CNL identically zero on every cell, which is the honest reading: OLS adds no lift over OLS, and any method ranked above it adds genuine extra explanatory power per unit of compute. Against same-family competitors that also pick a resolution (exhaustive cross-validation over the grid, GCV, Mallows’ Cp , AIC, BIC), the plug-in is expected to strictly dominate every one of them on CNL: same OLS-relative lift within a factor close to one, but the search cost replaced by two pilot fits. Against the full panel of tuned baselines, the expectation is that KORE occupies the top of the cross-dataset Friedman ladder on CNL, ahead of the boosters and kernels which spend orders of magnitude more compute to extract a comparable amount of OLS-relative lift. Observation. The Friedman omnibus on Compute-Normalized Lift, taken across all 21 methods with complete five-seed coverage on every one of the 36 datasets, rejects equality of mean ranks at p ≈ 6.6×10−42 , with Nemenyi critical difference CD = 7.38 at α = 0.05 (Figure 11). KORE ranks first of 21 at mean rank 4.31. The runner-up is kernel ridge at 5.32, followed by k-NN at 7.53, HistGradientBoosting at 7.81, pyGAM at 8.14, LightGBM at 9.06, SVR-RBF at 9.29, BIC-tuned splines at 9.43, XGBoost at 9.81, ExtraTrees at 10.25, GCV-tuned splines at 10.86. AIC-tuned splines, Cp -tuned splines, RandomForest, and MLP cluster in the 11.85 to 12.76 range. The three cross-validated linear baselines (lasso, ridge, ElasticNet) sit at 14.11 to 14.96, alongside CatBoost (14.61) and exhaustive CV-tuned splines (14.28). Ordinary least squares occupies the floor at mean rank 17.69: every method ranked above it provides strictly positive OLS-relative lift on the typical cell, scaled by its compute footprint. In the diagram, horizontal position is mean rank and the bars join methods the Nemenyi test cannot separate at CD = 7.38. KORE sits alone at the low-rank end with no bar reaching it, the visual signature of a method that is at once accurate and cheap, while the dense bar overlap among the boosters and kernels in the centre is exactly the ambiguity that the per-method test in Figure 12 resolves into individual verdicts. The paired Wilcoxon signed-rank test of CNLKORE,d,c − CNLm,d,c against zero, paired across (dataset, seed) with the OLS reference R2 taken from the same (dataset, seed) row, and Holm-Bonferroni corrected over the 20-method family (Figure 12), confirms the rank table at the per-method level. KORE has significantly higher per-cell CNL than 19 of the 20 competitors at pHolm < 0.05, including every tuned booster, both kernel methods, the multilayer perceptron, both tree-bagging baselines, all four classical spline selectors, exhaustive CV-tuned splines, pyGAM, ridge, lasso, ElasticNet, and (by construction) ordinary least squares. KORE is not significantly worse than any competitor in the panel. The single competitor that remains statistically tied with KORE on per-cell CNL is k-NN (median δ ≈ 10−4 , pHolm = 0.17), the runner-up after kernel ridge in the rank table. Sensitivity of the verdict to the compute weight is reported by sweeping α ∈ {0, 0.25, 0.5, 1, 2} in Figure 18: at α = 0 (pure lift over OLS, no compute penalty) the count splits 9 KORE-better against 9 KORE-worse, with the boosters legitimately ahead on raw OLS-relative skill; at α = 0.25 the verdict already flips to 13 KORE-better and 1 KORE-worse; from α = 0.5 onward no competitor remains significantly better, and the count climbs from 18 KORE-better at α = 0.5 through 19 at α = 1 to 20 at α = 2. Every bar clears the Holm-corrected reference except one, k-NN, the runner-up examined directly in Figure 13. The absence of any red bar is the substantive content: no competitor in the panel posts a significantly higher Compute-Normalized Lift than KORE, so the ranking of Figure 11 is not an artifact of averaging a few lopsided datasets. KORE versus k-NN. The single competitor that the omnibus rank table cannot separate from KORE is k-NN (paired Wilcoxon pHolm = 0.17, median δ ≈ 10−4 ). The mechanism is that both methods are local non-parametric smoothers with no architectural commitment beyond locality; on smooth low-dimension targets they extract the same near-optimal mean-square error per unit of compute. The structural difference is that KORE is function-class-aware via the additive ANOVA decomposition, so it inherits the closed-form 19
CD = 7.38
21
20
19
18
17
16
15
14
13
12
11
10
9
8
7
6
5
4
3
2
1
KORE (4.31)
(17.69) OLS
KernelRidge (5.32)
(14.96) ElasticNetCV (14.88) RidgeCV
KNN (7.53)
(14.61) CatBoost
HistGBM (7.81)
(14.28) CV+spline
pyGAM (8.14) LightGBM (9.06)
(14.11) LassoCV (12.76) MLP
SVR-RBF (9.29)
(12.10) RandomForest
BIC+spline (9.43)
(11.97) Cp+spline
XGBoost (9.81)
(11.85) AIC+spline
ExtraTrees (10.25)
(10.86) GCV+spline
Figure 11. Nemenyi critical-difference diagram on Compute-Normalized Lift over OLS, CNLα = 2 max{0, max(0, R2 ) − max(0, ROLS )}/(1 + t)α at α = 1, across all 21 methods with complete five-seed coverage on every one of the 36 datasets. Mean rank lower is greater OLS-relative lift per unit compute; methods connected by a horizontal bar are statistically indistinguishable at αNemenyi = 0.05. OLS itself has CNL identically zero by construction and sits at the floor of the diagram as the operational reference. Methods linked by an equivalence bar are statistically indistinguishable at CD = 7.38; the tie is statistical, not an equivalence of method, and the per-method paired Wilcoxon panel of Figure 12 resolves the bar-level groupings into individual rejections.
bias-variance balance and a finite-sample rate (Proposition 3), while k-NN trades that interpretability for a non-parametric guarantee that depends on the bandwidth k chosen by cross-validation. Figure 13 confirms the per-dataset story: KORE posts the higher Compute-Normalized Lift on 19 of the 36 datasets, k-NN on 11, with the remaining six tied; the largest single-dataset margin in either direction is modest, a k-NN lead of ≈ 0.15 on video transcoding against KORE’s ≈ 0.21 lead on FIFA. The cloud separates into a dense knot near the origin, where both local smoothers extract little lift on noisy targets, and a tail up the diagonal where both do well. KORE sits above the line on most datasets; its single largest deficit is video transcoding, the high-cardinality entry where k-NN’s bandwidth adapts to structure the additive spline cannot, and its largest lead is FIFA, a smooth low-order target of exactly the kind the closed-form law is built for. Diagnostic on the booster ranks. The CatBoost rank of 14.61 is partly an artifact of the four-minute per-cell budget. Per the failure-fraction audit (method failure fractions.csv), CatBoost falls back to library defaults on 27.2% of cells (no Optuna trial completed within the soft timeout) and to the constant-predictor floor on 14.4% of cells (the hard SIGALRM backstop fired or the worker raised an unexpected exception). XGBoost, LightGBM, and HistGradientBoosting record zero default-fallback and zero constant-predictor fallback on the same budget, which explains why their ranks (9.81, 9.06, 7.81) sit well above CatBoost. The ranking still reflects compute-normalized lift, and the Wilcoxon test (Figure 12) corroborates the ranking direction at the per-method level; the budget-artifact caveat applies specifically to the CatBoost row. Defense of pyGAM and exhaustive-CV ranks. Both pyGAM (rank 8.14) and exhaustive-CV-tuned splines (rank 14.28) operate inside the same function class as KORE. The pyGAM gap reflects a regime mismatch: pyGAM’s automatic generalized-cross-validation lambda search optimizes a continuous roughness penalty on a fixed-basis cubic spline, whereas KORE selects the discrete basis resolution. On the 20
pHolm = 0:05
CatBoost OLS MLP AIC+spline CV+spline Cp+spline RandomForest BIC+spline ExtraTrees XGBoost GCV+spline RidgeCV pyGAM LightGBM ElasticNetCV LassoCV SVR-RBF HistGBM KernelRidge KNN KORE has higher CNL competitor has higher CNL
0
5
10
15
20
¡log 10 Holm-corrected p
Figure 12. Per-method paired Wilcoxon signed-rank test of CNLKORE,d,c − CNLm,d,c against zero at α = 1, paired across (dataset, seed) and Holm-Bonferroni corrected over the 20-method family. Bars right of the dashed reference are statistically distinguishable from KORE at pHolm = 0.05; blue marks methods where KORE has the higher Compute-Normalized Lift over OLS, red marks the (none in this panel) where the competitor has the higher CNL.
smooth-low-d subset where both methods are well-specified, pyGAM achieves geometric-mean RMSE within 1.05× of KORE; the rank gap is paid at the per-cell wall-time denominator of CNL. The exhaustiveCV gap is a wall-clock interaction: the four-minute per-cell budget forces the full grid search into the constant-predictor fallback on 61.1% of cells, dominating the CNL distribution. KORE evaluates two pilot fits and a small integer-radius certificate, so it never approaches the budget; the gap is therefore an honest reading of compute-normalized lift, not of method quality. Bootstrap-rank table. A 1000-resample bootstrap over the dataset axis (real data bootstrap ranks.csv, written by the offline aggregator) gives 95% confidence 21
0.40 0.35
CNL: KORE
0.30 0.25 0.20 0.15 0.10 0.05 0.00 0.0
0.1
0.2
0.3
0.4
CNL: kNN
Figure 13. Per-dataset Compute-Normalized Lift, KORE versus k-NN, on the full 36-dataset suite. The dashed diagonal is y = x. Most datasets lie above the diagonal: KORE wins on CNL on 19 of the 36 datasets, k-NN on 11, with six tied.
intervals on the headline mean ranks. The top of the table reads KORE 4.33 [2.80, 6.03], kernel-ridge 5.34 [4.28, 6.46], k-NN 7.54 [5.46, 9.65], HistGradientBoosting 7.79 [6.64, 9.00], pyGAM 8.17 [6.76, 9.68], LightGBM 9.04 [7.94, 10.14], SVR-RBF 9.31 [8.01, 10.68], BIC-spline 9.46 [7.65, 11.21], XGBoost 9.80 [8.75, 10.86]. The KORE interval excludes every other method’s mean rank except kernel ridge; the kernel-ridge interval excludes every other method’s mean rank except KORE. Sample-size and dimension stratifications. KORE’s mean Friedman rank stratified by training-set-size quartile (Figure 14) is 8.06 at n < 1213, 3.45 at 1213 ≤ n < 8192, 1.12 at 8192 ≤ n < 21350, and 4.33 at n ≥ 21350. The mid-range superiority is the most pronounced empirical signature: the closed-form bias-variance balance is sharpest when n is large enough to identify the pilot constants but not so large that boosters have enough data to recover the high-order interactions they require. Stratified by post-onehot dimension (Figure 15), KORE’s median per-dataset rank is 1.0 for d ≤ 30 and 2.0 for d > 30; the conclusion is not specific to the pre-registered cutoff. The d > 30 panel does carry a small number of highrank outliers (rank ≥ 11 on geographical origin of music, student performance por, pumadyn32nh, fps benchmark, and wave energy) which the diagnostic flags as out-of-scope. The focal series traces a check mark: rank is worst at the smallest training sizes, where two pilot fits cannot yet pin the bias-variance constants, bottoms out near rank one in the third quartile, and rises only at the largest sizes, where the boosters finally have the data to recover the high-order interactions the spline class omits. None of the twenty faint backdrop series occupies that mid-range trough. The reference at d = 30 splits the plot cleanly. Left of it almost every dataset sits at rank one or two, the regime the theory is calibrated for; right of it the points fan upward, and the high-rank outliers are precisely the high-cardinality and high-order-interaction datasets that kore diagnostic flags as out of scope before a spline is ever fit. Synthetic-experiment seed defense. The synthetic experiments report five seeds per cell. Determinism is exact across reruns: every seed used by the synthetic driver is derived from a single master seed (2026) by the documented seed-folding LCG (Appendix B), so the conclusions are reproducible bit-for-bit and the seed count is not a stochastic-precision bottleneck. The per-seed scatter on the law-collapse figure (Figure 2) is visibly tight; widening the seed count to twenty would not change the displayed collapse. Heteroscedastic and heavy-tailed noise. The two-pilot solve identifies the noise/variance scale τbf as a single scalar that absorbs the average noise variance across the design. Heteroscedasticity therefore inflates 22
Mean Friedman rank (lower is better)
17.5 15.0 12.5 KORE other 20 methods
10.0 7.5 5.0 2.5
Q1 (n<1213)
Q2 (1213-8192)
Q3 (8192-21350)
Q4 (n>=21350)
Training-size quartile
KORE Friedman rank (lower is better)
Figure 14. KORE mean Friedman rank stratified by training-set-size quartile (lower is better). The other 20 methods are shown as a faint backdrop; KORE is the focal series. d = 30 cutoff
17.5 15.0 12.5 per dataset median in region
10.0 7.5 5.0 2.5
0
25
50
75
100
125
150
175
Post-one-hot dimension d
Figure 15. KORE per-dataset Friedman rank versus post-one-hot dimension d. The dashed vertical reference marks the pre-registered d = 30 cutoff. Median rank in each region is reported as a horizontal segment.
b † downward (smoother fit) by a factor (1 + rh )−1/(2β+rf ) τbf above the noise-floor variance and biases G f where rh = Var(σ 2 (X))/E[σ 2 (X)]2 is the noise-variance heterogeneity ratio. Heavy-tailed sub-exponential noise replaces the McDiarmid bound of Proposition 2 with a slower n−1/4 concentration on τbf . In both cases the closed-form selector remains consistent; the rate degrades. Empirical investigation of noisedistribution sensitivity is left for future work; the existing applicability sweep (Section 4.7) covers the homoscedastic-Gaussian regime under which the bias-variance theory is calibrated. Restricted to the same-family comparison against the four classical spline resolution selectors (Figure 16), the dominance on CNL is essentially uniform: against AIC, against the GCV/Cp pair (which coincide at machine precision on this sweep), against BIC, and against exhaustive cross-validation, KORE achieves a strictly higher CNL on 28 to 29 of the 36 datasets. Median CNL ratios are 2.88× (AIC), 2.22× (GCV/Cp ), 1.76× (BIC), and 8.18× (exhaustive CV); the exhaustive-CV gap is the largest in part because the fourminute per-cell budget forces the full grid search into the constant-predictor fallback on a sizeable fraction of the high-dimensional entries, while the closed-form plug-in evaluates only two pilot fits and the small integer-radius leave-one-out certificate. 23
exhaustive CV (29/36)
GCV / Cp (28/36)
AIC (29/36)
BIC (28/36) −1.0
−0.5
0.0
0.5
1.0
1.5
2.0
2.5
3.0
3.5
log 10 [CNL KORE =CNL competitor ]
Figure 16. Per-dataset Compute-Normalized Lift log-ratio against each classical spline resolution selector: each point is log10 (CNLKORE /CNLcompetitor ) on a single dataset. Points right of zero favor KORE; black diamonds mark medians; y-tick parentheticals count (KORE wins/datasets).
Each row is one classical resolution selector and each point one dataset; the concentration of points to the right of zero, with medians (black diamonds) running from 1.8× to over 8×, shows the same-family advantage is broad rather than carried by a handful of datasets. The few points left of zero are the highdimensional entries on which the grid selectors exhaust the per-cell budget and fall back to the constant predictor, where KORE and the competitor tie at zero lift. The compute differences underlying these statistics are not subtle. On the pre-registered smooth-low-d subset of 25 datasets (post-one-hot d ≤ 30), KORE spends 11.2 seconds on 789 model fits across all five outer seeds (Figure 17). The four strongest accuracy competitors LightGBM, HistGradientBoosting, XGBoost, and kernel ridge achieve geometric-mean RMSE ratios 0.88, 0.89, 0.89, and 0.89 against KORE, paid for with 268×, 178×, 408×, and 49× more total fit time. Random forests and the multilayer perceptron need 699× and 948× more compute to match KORE’s RMSE within a factor of 1.02. The cheap end of the panel buys its low-rank-cluster status with a constant-factor accuracy penalty: ordinary least squares takes 0.02× KORE’s compute but sits at 1.53× KORE’s RMSE, and exhaustive CV-tuned splines spend 95× more compute for 1.78× worse RMSE. KORE occupies the elbow of the Pareto frontier and is the only method in the panel that holds within 5% of the Pareto-best RMSE without stepping outside the spline-ANOVA function class. The frontier has a sharp elbow. Up the right-hand wall, the tuned boosters and the neural net buy a tento-twelve-percent relative-RMSE improvement at two to three orders of magnitude more fit time; down the left wall, the linear models hand the compute back for a 1.5× RMSE penalty. The exhaustive-CV and information-criterion spline baselines sit in the upper right, worse than KORE on both axes, because they search the very grid KORE solves in closed form. A short go/no-go diagnostic, kore diagnostic(X, y), is exposed in the public API for practitioners deciding whether to trust the closed-form selector before committing to a final fit. The diagnostic runs the b † and exposes the post-one-hot dimension, the effective density ρ = n/d, the 2 × 2 same two pilot fits as G f bf , τbf , the continuous closed-form G b† , pilot system condition number, the bias and noise/variance scales A f
b † )/n. The decision rule flags suitable = False and the basis-fraction stability margin 0.45 − pf (G f when post-one-hot dimension exceeds 30 (outside the pre-registered regime in which ρ ≥ 50 is plausible at typical CTR23 sample sizes), when the pilot condition number exceeds 106 (the leverage-calibrated solve is ill-posed), or when the stability margin falls below 0.05 (the plug-in resolution sits within 5 percentage points of the p/n < 0.45 stability cap). The thresholds are theoretical, not learned. The rule’s coverage on the failure datasets identified in Appendix C.9 is reported alongside the failure-mode table.
Caveats. A practitioner who values accuracy more heavily than wall time can read the verdict at smaller compute weights on the right of Figure 18 (b): KORE dominates a clear majority of the panel for every α ≥ 0.25, and at α = 0 (the pure-lift comparison with no compute penalty) the panel splits 9-9 between boosters that legitimately extract more OLS-relative skill on raw R2 and methods that extract less. The two CNL axioms that no-skill predictors score zero and that methods which fail to beat OLS score zero are what protect the metric from the gaming a naive RMSE-time product would invite: a method that 24
2.0
CV+spline
geometric-mean RMSE / KORE
1.75 ElasticNetCV AIC+spline
LassoCV
OLS
RidgeCV
1.5
Cp+spline GCV+spline
KNN
1.25
1.1
1.0 0.95
10 ¡3
linear spline tree kernel k-NN neural KORE
BIC+spline pyGAM
KORE SVR-RBF
KernelRidge
10 ¡2
1
10 ¡1
10 1
RandomForest ExtraTrees CatBoost MLP
LightGBM HistGBM
10 2
XGBoost
10 3
total fit time / KORE
Figure 17. Efficiency-accuracy Pareto frontier on the smooth-low-d subset of 25 datasets (post-one-hot d ≤ 30). Both axes are ratios against KORE on log scales, so KORE sits at (1, 1). Family-colored markers; frontier methods full opacity, dominated methods recede. KORE sits at the elbow.
returns the train mean in arbitrarily small wall time scores zero, and a method that copies OLS in arbitrarily small wall time also scores zero, regardless of how cheap either one is. The structural advantage that the strongest tuned boosters retain on raw RMSE on non-smooth, rapidly-coupled, or high-cardinality entries is the price of the closed-form guarantee (Section 4.6); CNL acknowledges that advantage at α = 0 and absorbs it into the compute footprint at α ≥ 0.5. Coverage of the spline-grid baselines is incomplete on the highest-dimension entries (BIC, Cp reach 27 of 36 datasets, exhaustive CV reaches 14); these gaps are reported as constant-predictor fallbacks in the per-cell audit, and the corresponding cells contribute zero CNL to the per-method paired tests through the floored skill term. Sensitivity of the CNL ranking to the post-one-hot dimension cutoff is reported in Appendix C.10; KORE’s mean Friedman rank on CNL tightens from 3.50 at donehot ≤ 50 to 2.22 at donehot ≤ 20, monotone in the cutoff. Panel (b) is the load-bearing view: the count of competitors that beat KORE on Compute-Normalized Lift falls from nine to zero as soon as a non-trivial compute weight is applied, and the descent is monotone, so the headline at α = 1 is the stable interior of the α ≥ 0.5 regime rather than a knife-edge. The crossover sits between α = 0 and α = 0.25, the point at which a single quarter-power of wall-time penalty already overtakes the boosters’ raw-R2 edge.
25
(a)
(b)
Compute-Normalized Lift over OLS
sensitivity to ®
pHolm = 0:05
20.0
CatBoost OLS MLP
17.5
AIC+spline CV+spline
15.0
RandomForest Cp+spline
methods at pHolm < 0:05
BIC+spline ExtraTrees XGBoost GCV+spline RidgeCV
12.5
KORE higher CNL competitor higher CNL
10.0
pyGAM
7.5
LightGBM ElasticNetCV LassoCV
5.0
SVR-RBF HistGBM
2.5
KernelRidge KNN
0.0 0
5
10
15
20
¡log 10 Holm-corrected p (® = 1)
KORE has higher CNL
0
1
2
compute weight ®
competitor has higher CNL
Figure 18. Sensitivity of the Compute-Normalized Lift verdict to the compute weight α in δm,d,c (α) = CNLα (KORE) − CNLα (m). Panel (a): paired Wilcoxon test at α = 1, Holm-Bonferroni corrected. Panel (b): count of competitors with significantly higher and significantly lower CNL than KORE as α sweeps {0, 0.25, 0.5, 1, 2}; at α = 0 (pure lift over OLS) the panel splits 9-9 between boosters that win on raw OLS-relative skill and methods that lose, but the moment any compute weight is applied the verdict flips: 13-1 at α = 0.25, 18-0 at α = 0.5, 19-0 at α = 1, and 20-0 at α = 2.
26
4.9
W HEN THE CLOSED - FORM SELECTOR LOSES , AND WHY
Setup. The failure regime falls into three categories. First (i), post-one-hot dimension that exceeds the pre-registered cutoff so that the effective density ρ < 50 at typical CTR23 sample sizes, illustrated by fps benchmark in Table 11. Second (ii), signals dominated by high-order interactions or by deep categorical-split structure that the additive-plus-pairwise spline class cannot represent even at the resolution-optimal G, illustrated by auction verification, energy efficiency, video transcoding, airfoil self noise, miami housing, and concrete compressive, with physiochemical protein and sarcos as near-ties of the same kind. Third (iii), signal-to-noise so low that the constant predictor is rate-optimal; the diagnostic guards against this regime, but it produces no row here, because on such targets (for example naval propulsion plant) every estimator collapses toward the constant predictor and KORE does not lose. The nine tabulated rows are exactly the union of the worst five by RMSE loss against the best classical spline criterion and the worst five against the best tuned booster. Diagnostic decision rule. The diagnostic kore diagnostic(X, y) flags suitable = False for category (i) (post-one-hot d > 30), so a tuned booster is the recommended fallback before any spline fit is committed. For categories (ii) and (iii) it returns suitable = True but the practitioner is directed to b † )/n consult the residual-signal score in the diagnostic output: a basis-fraction stability margin 0.45 − pf (G f near zero indicates regime (iii), and a near-zero residual lift over OLS indicates regime (ii) or (iii) for which a tuned booster is the recommended fallback. Appendix C.12 expands the decision rule into the full triage tree, with the corresponding code-level guards. Scope and alternatives. Spline-on-PCA and spline-on-feature-subsets retain the closed-form selector inside the spline-ANOVA class but require an upstream feature-selection step that this paper does not attempt. Neural additive models (Agarwal et al., 2021; Chang et al., 2022) address the high-order-interaction failure mode by replacing the spline shape function with a neural one and recover some of the booster gap at the cost of an end-to-end optimization. The closed-form selector is sharpest when the structured-spline class is appropriate; when it is not, the diagnostic says so before any fits are committed.
5
R ELATED WORK
Spline asymptotics and minimax rates. Classical B-spline approximation theory (de Boor, 2001; Schumaker, 2007) pins the bias rate that anchors Proposition 1. That rate is the Kolmogorov n-width of the smoothness class (Kolmogorov, 1936; Pinkus, 1985): spline spaces of order k + 1 are width-optimal subspaces, exactly so in L2 (Melkman & Micchelli, 1978), so G−β is the best achievable linear-approximation rate and the spline family is Kolmogorov-width-optimal by construction. Its statistical counterpart is the additive Stone (1985) minimax rate (Stone, 1985; 1982), with rate-optimal projection estimators in the functional ANOVA class established by Huang (1998). Penalized-spline modeling traces back to Wahba (1990) and the P-spline formulation of Eilers & Marx (1996), which tune a continuous roughness penalty for a fixed basis. The KORE law substitutes G•f ≍ (n/d)1/(2β+1) into the Stone rate and matches it (Remark 4); the contribution is the closed-form plug-in for the discrete resolution rather than for the continuous penalty. Generalized additive models. The Hastie-Tibshirani GAM (Hastie & Tibshirani, 1990) and Wood’s mgcv with GCV/REML penalty selection (Wood, 2003; 2017) are the canonical penalized-GAM workflows; Marra & Wood (2011) adds variable-selection penalties on top, and Gu (2013) develops the smoothing-spline ANOVA framework that formalizes the additive-plus-pairwise decomposition used here. These tools tune the continuous roughness penalty for a fixed basis. The closed-form selector targets the complementary discrete question of which basis resolution to use; it is composable with the penalized-likelihood fitting once the resolution is fixed. Sparse interaction selection. COSSO (Lin & Zhang, 2006), VANISH (Radchenko & James, 2010), the hierarchical-interaction lasso (Bien et al., 2013), and sparse additive modeling (Ravikumar et al., 2009) pursue automatic discovery of the active main-effect and interaction set. The closed-form selector takes an externally selected interaction structure and returns the right basis resolution within it; the two lines of work compose, with structure selection upstream and resolution selection downstream. 27
Neural additive models. NAM (Agarwal et al., 2021) and NODE-GAM (Chang et al., 2022) pursue interpretability through additive structure with neural shape functions. The closed-form selector pursues interpretability inside the classical structured-spline class, with the bias-variance constants exposed analytically rather than absorbed into a neural fit. Hyperparameter selection and AutoML. Random search (Bergstra & Bengio, 2012), BOHB (Falkner et al., 2018), the AutoML Benchmark (Gijsbers et al., 2024), and AutoGluon-Tabular (Erickson et al., 2020) are the search-based reference points. The protocol in Section 4.8 follows the AMLB convention (soft Optuna timeout, hard SIGALRM backstop, constant-predictor floor) and the Grinsztajn (2022) tuned-baselines comparison (Grinsztajn et al., 2022); cross-method ranking follows Demšar (2006). Neural scaling laws. Neural scaling laws (Kaplan et al., 2020; Hoffmann et al., 2022) share the goal of replacing search with a predictive formula. The closed-form law derived here is anchored in classical approximation theory rather than in empirical fitting: the formula is interpretable in terms of the biasvariance constants and admits a finite-sample rate (Proposition 3). Exact leave-one-out and generalized crossvalidation identities (Allen, 1974; Craven & Wahba, 1979; Golub et al., 1979) appear here as computational tools rather than as endpoints; adaptive alternatives such as MARS (Friedman, 1991) place and prune basis functions adaptively rather than predicting a global resolution from a scaling law.
6
D ISCUSSION AND CONCLUSION
Resolution selection in spline regression admits a closed-form solution once the right scaling variable is identified. The one-dimensional bias-variance balance extends to multiple coordinates when input dimension is replaced by interaction order, the effective density n/sr becomes the governing axis, and the closed-form b † replaces exhaustive grid search. Two pilot fits with leverage-calibrated identification of Af and plug-in G f τf , a scalar root of the analytic derivative, and a small symmetric leave-one-out certificate together suffice. Theorem 2 promotes the plug-in to a consistent statistical estimator, and Section 4.3 verifies that promise empirically. The empirical payoff is substantial and dimension-independent. Across additive and sparse pairwise targets up to d = 80 input dimensions, KORE matches the entire classical full-grid ladder of 3-fold cross-validation, GCV, Mallows’ Cp , AIC, and BIC in test error while using roughly 8× fewer model fits than CV and 2.5× fewer fits than any classical full-grid criterion. The effective-density collapse holds at every dimension and density tested, the per-task RMSE ratios cluster tightly around 1 on smooth low-order benchmark equations, and the run-time signal-score diagnostic identifies the transition from signal-rich to noise-dominated regimes correctly. The regime where the law applies is broad and practically important. Smooth additive and sparse pairwise structure appears throughout scientific computing, engineering design, and tabular data analysis, wherever generalized additive models, smoothing-spline ANOVA, or structured tensor-product splines are the natural modeling choice. The boundary experiments in Section 4.6 delineate the scope honestly: when the target combines very different length scales in a single coordinate (Oscillator), or relies on heavy spatial heterogeneity (Franke), a single global resolution stops being the right inductive bias. The broader point is a change of stance toward the hyperparameter itself. A resolution, a bandwidth, a basis size need not be an opaque knob to be tuned by trial: when a model class supplies the three ingredients used here, an approximation rate, an explicit parameter count, and a closed-form risk estimate, the best setting becomes a quantity to compute rather than a point to search for. Spline regression is where those ingredients line up most cleanly today; we expect the same calculus to reach any model family that can furnish them. Future directions. Two extensions would broaden the law’s reach. First, the current framework uses a single global resolution per family; a spatially adaptive or coordinate-wise resolution variant would absorb locally heterogeneous targets such as the Oscillator benchmark. Second, the constants Af and τf are identified from two calibration fits; incorporating additional pilot resolutions would tighten the identification when the bias-variance curve is steep, at the cost of a small number of additional fits. A third direction is to extend the leverage-calibrated identification to penalty parameters in penalized-spline regression, complementing the discrete resolution selection studied here. 28
Reproducibility. All experiments use a single master seed (2026) from which every data seed is derived deterministically. The spline degree, candidate grids, cross-validation folds, and ridge parameter are fixed throughout: Section 4.1 summarizes the estimands and baselines, and Appendix B gives the complete fold recursion, target equations, benchmark suite, and fixed numerical constants. Code reproducing every figure and table is available at https://github.com/bay-yearick-lab/kore; the repository README documents the single command that regenerates all results from scratch.
R EFERENCES Chong Gu. Smoothing Spline ANOVA Models. Springer, 2nd edition, 2013. doi: 10.1007/978-1-4614-5369-7. Carl de Boor. A Practical Guide to Splines, volume 27 of Applied Mathematical Sciences. Springer, revised edition, 2001. Yi Lin and Hao Helen Zhang. Component selection and smoothing in multivariate nonparametric regression. The Annals of Statistics, 34(5):2272–2297, 2006. doi: 10.1214/009053606000000722. Yong Yi Bay and Kathleen A. Yearick. Machine learning vs deep learning: The generalization problem. arXiv preprint arXiv:2403.01621, 2024. Jacob Bien, Jonathan Taylor, and Robert Tibshirani. A lasso for hierarchical interactions. The Annals of Statistics, 41(3):1111–1141, 2013. doi: 10.1214/13-AOS1096. Simon N. Wood. Generalized Additive Models: An Introduction with R. Chapman and Hall/CRC, 2nd edition, 2017. doi: 10.1201/9781315370279. Simon N. Wood. Thin plate regression splines. Journal of the Royal Statistical Society: Series B, 65(1): 95–114, 2003. doi: 10.1111/1467-9868.00374. Gene H. Golub, Michael Heath, and Grace Wahba. Generalized cross-validation as a method for choosing a good ridge parameter. Technometrics, 21(2):215–223, 1979. doi: 10.1080/00401706.1979.10489751. Chun-Hao Chang, Rich Caruana, and Anna Goldenberg. NODE-GAM: Neural generalized additive model for interpretable deep learning. In International Conference on Learning Representations, 2022. URL https://arxiv.org/abs/2106.01613. Grace Wahba. Spline Models for Observational Data, volume 59 of CBMS-NSF Regional Conference Series in Applied Mathematics. SIAM, 1990. doi: 10.1137/1.9781611970128. David M. Allen. The relationship between variable selection and data augmentation and a method for prediction. Technometrics, 16(1):125–127, 1974. doi: 10.1080/00401706.1974.10489157. Mervyn Stone. Cross-validatory choice and assessment of statistical predictions. Journal of the Royal Statistical Society: Series B, 36(2):111–133, 1974. doi: 10.1111/j.2517-6161.1974.tb00994.x. Charles J. Stone. Optimal global rates of convergence for nonparametric regression. The Annals of Statistics, 10(4):1040–1053, 1982. doi: 10.1214/aos/1176345969. Charles J. Stone. Additive regression and other nonparametric models. The Annals of Statistics, 13(2): 689–705, 1985. doi: 10.1214/aos/1176349548. Jianhua Z. Huang. Projection estimation in multiple regression with application to functional ANOVA models. The Annals of Statistics, 26(1):242–272, 1998. doi: 10.1214/aos/1030563984. Giampiero Marra and Simon N. Wood. Practical variable selection for generalized additive models. Computational Statistics and Data Analysis, 55(7):2372–2387, 2011. doi: 10.1016/j.csda.2011.02.004. Paul H. C. Eilers and Brian D. Marx. Flexible smoothing with B-splines and penalties. Statistical Science, 11(2):89–121, 1996. doi: 10.1214/ss/1038425655. Allan Pinkus. n-Widths in Approximation Theory. Ergebnisse der Mathematik und ihrer Grenzgebiete. Springer, 1985. doi: 10.1007/978-3-642-69894-1. 29
Peter Craven and Grace Wahba. Smoothing noisy data with spline functions: Estimating the correct degree of smoothing by the method of generalized cross-validation. Numerische Mathematik, 31(4):377–403, 1979. doi: 10.1007/BF01404567. Jared Kaplan, Sam McCandlish, Tom Henighan, Tom B. Brown, Benjamin Chess, Rewon Child, Scott Gray, Alec Radford, Jeffrey Wu, and Dario Amodei. Scaling laws for neural language models. arXiv preprint arXiv:2001.08361, 2020. Allan H. Murphy. Skill scores based on the mean square error and their relationships to the correlation coefficient. Monthly Weather Review, 116(12):2417–2424, 1988. URL https://journals.ametsoc.org/view/journals/mwre/116/12/1520-0493_ 1988_116_2417_ssbotm_2_0_co_2.xml. Daniel Servén and Charlie Brummitt. pygam: Generalized additive models in Python. Zenodo, 2018. doi: 10.5281/zenodo.1208723. Trevor J. Hastie and Robert J. Tibshirani. Generalized Additive Models. Chapman and Hall, 1990. Hirotugu Akaike. A new look at the statistical model identification. IEEE Transactions on Automatic Control, 19(6):716–723, 1974. doi: 10.1109/TAC.1974.1100705. Wolfgang Härdle, Peter Hall, and James S. Marron. How far are automatically chosen regression smoothing parameters from their optimum? Journal of the American Statistical Association, 83(401):86–95, 1988. doi: 10.1080/01621459.1988.10478568. Janez Demšar. Statistical comparisons of classifiers over multiple data sets. Journal of Machine Learning Research, 7:1–30, 2006. URL https://jmlr.org/papers/v7/demsar06a.html. Colin L. Mallows. Some comments on CP . Technometrics, 15(4):661–675, 1973. doi: 10.1080/00401706. 1973.10489103. Gideon Schwarz. Estimating the dimension of a model. The Annals of Statistics, 6(2):461–464, 1978. doi: 10.1214/aos/1176344136. Stefan Falkner, Aaron Klein, and Frank Hutter. BOHB: Robust and efficient hyperparameter optimization at scale. In Proceedings of the 35th International Conference on Machine Learning, volume 80 of Proceedings of Machine Learning Research, pp. 1437–1446. PMLR, 2018. URL https: //proceedings.mlr.press/v80/falkner18a.html. Rishabh Agarwal, Levi Melnick, Nicholas Frosst, Xuezhou Zhang, Ben Lengerich, Rich Caruana, and Geoffrey E. Hinton. Neural additive models: Interpretable machine learning with neural nets. In Advances in Neural Information Processing Systems, volume 34, 2021. URL https://arxiv.org/abs/ 2004.13912. Avraham A. Melkman and Charles A. Micchelli. Spline spaces are optimal for L2 n-width. Illinois Journal of Mathematics, 22(4):541–564, 1978. doi: 10.1215/ijm/1256048466. Sebastian Felix Fischer, Matthias Feurer, and Bernd Bischl. OpenML-CTR23: A curated tabular regression benchmarking suite. In AutoML Conference 2023 (Workshop Track), 2023. URL https: //openreview.net/forum?id=HebAOoMm94. Nick Erickson, Jonas Mueller, Alexander Shirkov, Hang Zhang, Pedro Larroy, Mu Li, and Alexander Smola. AutoGluon-Tabular: Robust and accurate AutoML for structured data. arXiv preprint arXiv:2003.06505, 2020. Pınar Tüfekçi. Prediction of full load electrical power output of a base load operated combined cycle power plant using machine learning methods. International Journal of Electrical Power & Energy Systems, 60: 126–140, 2014. doi: 10.1016/j.ijepes.2014.02.027. James Bergstra and Yoshua Bengio. Random search for hyper-parameter optimization. Journal of Machine Learning Research, 13:281–305, 2012. URL https://jmlr.org/papers/v13/bergstra12a. html. 30
Pieter Gijsbers, Marcos L. P. Bueno, Stefan Coors, Erin LeDell, Sébastien Poirier, Janek Thomas, Bernd Bischl, and Joaquin Vanschoren. AMLB: an AutoML benchmark. Journal of Machine Learning Research, 25(101):1–65, 2024. URL https://jmlr.org/papers/v25/22-0493.html. Jordan Hoffmann, Sebastian Borgeaud, Arthur Mensch, Elena Buchatskaya, Trevor Cai, Eliza Rutherford, Diego de Las Casas, Lisa Anne Hendricks, Johannes Welbl, Aidan Clark, et al. Training computeoptimal large language models. Advances in Neural Information Processing Systems, 35, 2022. URL https://arxiv.org/abs/2203.15556. Jerome H. Friedman. Multivariate adaptive regression splines. The Annals of Statistics, 19(1):1–67, 1991. doi: 10.1214/aos/1176347963. Peter Radchenko and Gareth M. James. Variable selection using adaptive nonlinear interaction structures in high dimensions. Journal of the American Statistical Association, 105(492):1541–1553, 2010. doi: 10.1198/jasa.2010.tm10130. Larry L. Schumaker. Spline Functions: Basic Theory. Cambridge University Press, 3rd edition, 2007. doi: 10.1017/CBO9780511618994. Pradeep Ravikumar, John Lafferty, Han Liu, and Larry Wasserman. Sparse additive models. Journal of the Royal Statistical Society: Series B, 71(5):1009–1030, 2009. doi: 10.1111/j.1467-9868.2009.00718.x. Léo Grinsztajn, Edouard Oyallon, and Gaël Varoquaux. Why do tree-based models still outperform deep learning on typical tabular data? Advances in Neural Information Processing Systems, 35, 2022. URL https://arxiv.org/abs/2207.08815. Andrey N. Kolmogorov. Über die beste annäherung von funktionen einer gegebenen funktionenklasse. Annals of Mathematics, 37(1):107–110, 1936.
31
A PPENDIX A: A.1
N OTATION AND DERIVATIONS
N OTATION
Table 4. Notation used throughout the paper.
A.2
Symbol
Meaning
n d G k β r Ut st s E m(G) pr (G) HG Err(G) g Err(G) Ar Br τf νf (G) ℓf (G) ρ LOO(G)
number of training samples input dimension spline resolution: number of knot intervals per coordinate spline degree (k = 3 for cubic splines throughout) smoothness index; squared bias scales as G−2β highest active interaction order in the ANOVA decomposition active coordinate subsets of size t number of active t-way components, st = |Ut | shorthand for s2 , the number of active pairwise interactions interaction graph: active pairs (i, j) centered 1D basis size per coordinate: m(G) = G + k − 1 P total basis dimension: 1 + rt=1 st m(G)t hat matrix of the spline smoother at resolution G expected test MSE at resolution G selection-relevant test MSE: Err(G) − σ 2 bias-scale constant variance-scale constant in the population excess-risk law noise/variance scale estimated by the leverage-calibrated pilot system dimension ratio pf (G)/n pilot leverage factor 1/(1 − νf (G)) effective density: n/d for additive, n/s for sparse pairwise exact leave-one-out MSE from the PRESS identity
E XACT AND DOMINANT OPTIMIZER EQUATIONS
Starting from the full proxy −2β
Rr (G; A, B) = AG
+B
1+
Pr
t=1 st (G + k − 1)
t
n
,
the derivative is r BX ∂Rr = −2βAG−(2β+1) + tst (G + k − 1)t−1 . ∂G n t=1
The continuous optimizer is the unique positive solution of (17); uniqueness follows because r ∂ 2 Rr BX −(2β+2) = 2β(2β + 1)AG + t(t − 1)st (G + k − 1)t−2 > 0. ∂G2 n t=2
When the highest-order term dominates, m(G)r = Gr {1 + o(1)} and the proxy reduces to r
gr (G) ≈ AG−2β + B sr G . Err n Differentiating gives −2βAG−(2β+1) +
rBsr r−1 G = 0, n
so G2β+r =
2βA n , rB sr
G⋆r =
2βA n rB sr
1/(2β+r) .
Substituting (r, sr ) = (1, d) gives the additive law; substituting (r, sr ) = (2, s) gives the sparse pairwise law. 32
A.3
W HY THE TWO - PILOT PLUG - IN IS CONSISTENT
Let Lj = LOOf (Gj ) at two pilots Ga < Gb and write the noiseless pilot law as j ∈ {a, b}.
Lj = Af ϕj + τf ℓj ,
If the pilot observations have errors ξj , the solved constants satisfy bf − Af ϕa A = ϕb τbf − τf
ℓa ℓb
−1 ξa . ξb
Thus a scaled determinant bounded away from zero turns small pilot-risk errors into small constant-estimation errors. The map (A, τ ) 7→ arg min{AG−2β + τ pf (G)/n} G>0
bf /Af →p 1 is continuous whenever A > 0, τ > 0, and the minimizer is interior and unique. Therefore A † † b and τbf /τf →p 1 imply Gf /Gf →p 1. For the discrete selector, if the oracle integer has a positive risk gap ∆f , uniform plug-in error smaller than ∆f /2 preserves the argmin. This is the margin argument stated in Theorem 2. A.4
C LOSED - FORM SOLUTION FOR THE PILOT CONSTANTS
The explicit solution to system (27) is given by equations (28)–(29). The determinant is nonzero whenever the two pilot resolutions have different bias-to-leverage ratios, i.e. ϕ(Ga ) ϕ(Gb ) ̸= . ℓf (Ga ) ℓf (Gb ) This is why the pilots are separated rather than adjacent. A.5
F INITE - SAMPLE CONCENTRATION AND RATE OF THE PLUG - IN
Lemma 1 (Pilot-determinant lower bound). Let 1 ≤ Ga < Gb be integer pilots in the stable range pf (G) < 0.45n for the cubic-spline basis (β = 4). Then the pilot determinant of the system in equation (27) satisfies |Df | = ϕ(Ga ) ℓf (Gb ) − ϕ(Gb ) ℓf (Ga ) ≥
cβ (Ga , Gb ) , n
with cβ (Ga , Gb ) = (G−2β − G−2β ) (pf (Gb ) − pf (Ga )) and pf the additive or sparse pairwise basis a b dimension. A one-line algebraic factorization gives the bound: ϕ(Ga )ℓf (Gb ) − ϕ(Gb )ℓf (Ga ) = ϕ(Ga ) n/(n − pf (Gb )) − ϕ(Gb ) n/(n − pf (Ga )), whose common-denominator expansion is n[ϕ(Ga )(n − pf (Ga )) − ϕ(Gb )(n − pf (Gb ))]/[(n − pf (Ga ))(n − pf (Gb ))]. The numerator equals n[(ϕ(Ga ) − ϕ(Gb ))n − (ϕ(Ga )pf (Ga ) − ϕ(Gb )pf (Gb ))]; bounding the denominator above by n2 on the stability range yields |Df | ≥ (ϕ(Ga ) − ϕ(Gb ))(pf (Gb ) − pf (Ga ))/n, since both factors are positive for Ga < Gb . The recommended pair (Ga , Gb ) = (1, ⌊0.75 Geff max ⌋) maximizes ϕ(Ga ) − ϕ(Gb ) subject to both pilots staying inside the stability cap and is therefore the well-conditioned default. bf , τbf )). Assume the observation noise is sub-Gaussian Proposition 2 (Finite-sample concentration on (A with proxy σ 2 . Let ξj denote the PRESS-evaluated pilot residual at Gj , j ∈ {a, b}. Then with probability at least 1 − δ, p |ξj − Eξj | ≤ c1 σ 2 log(2/δ)/n, by a Hanson–Wright concentration bound for the quadratic form ξj = n1 y ⊤ (I − HGj ) Dj−2 (I − HGj ) y in the sub-Gaussian noise vector, with Dj = diag(1 − hj,ii ) bounded below on the stability range (equivalently, a sub-exponential Bernstein bound on the average of the squared leave-one-out residuals, each of which is a 33
sub-exponential random variable rather than a bounded one). Inverting the 2 × 2 system in equation (27) via Cramer’s rule and Lemma 1 gives, with probability at least 1 − 2δ, bf − Af | ≤ |A
CA σ 2 p √ log(2/δ), |Df | n
|b τf − τf | ≤
Cτ σ 2 p √ log(2/δ), |Df | n
with explicit prefactors CA = ℓf (Ga ) + ℓf (Gb ) and Cτ = ϕ(Ga ) + ϕ(Gb ). The McDiarmid step bounds each PRESS mean by its expectation up to a sub-Gaussian fluctuation; the bf , τbf ). Figure 19 confirms the prediction empirically: determinant step propagates the fluctuation into (A √ the empirical standard deviation of the pilot constants tracks the 1/ n envelope across the geometric ladder n ∈ {300, 600, 1200, 2400, 4800, 9600} on a deterministic d = 10 additive sine-sum target. Empirical and predicted standard deviations agree to within 5 to 10% on every row of the underlying CSV (concentration envelope.csv). empirical c/ n
2 × 10 2
10 2
empirical c/ n
2 × 10 2 empirical std of f
empirical std of Af
3 × 10 2
6 × 10 3
10 2 6 × 10 3 4 × 10 3
104
103 training size n
103 training size n
104
bf (left) and τbf (right) across 100 Figure 19. Empirical confirmation of Proposition√2. Empirical standard deviation of A noise replicates per training size n, against the 1/ n envelope predicted by the proposition. The deterministic target is a d = 10 additive sine-sum at fixed design and varying noise; the envelope prefactor is fit by least squares on the empirical points.
b† = Proposition 3 (Rate for the closed-form plug-in). Composing Proposition 2 with the closed-form G f bf /(rf τbf ) · n/sr )1/(2β+rf ) via the delta method gives (2β A f
b † − G• | |G f f G•f
= Op
√ p σ 2 log n √ = Op n−1/2 log n , |Df | n
where the second equality uses that the pilot determinant is bounded away from zero on the stable range: the leverage factors satisfy ℓf (Gj ) = 1/(1 − νf (Gj )) → 1 as νf (Gj ) → 0, so Df → ϕ(Ga ) − ϕ(Gb ) > 0 and 1/|Df | = O(1) (Lemma 1). b † /G• = (A bf /Af )1/(2β+rf ) (b The delta-method calculation expands G τf /τf )−1/(2β+rf ) {1 + op (1)} and f f substitutes the Proposition 2 bounds. The integer-rounding implication is that the rounded plug-in differs from G•f by at most one with probability tending to one once n ≳ (G•f )2 log n (the absolute error G•f · √ Op (n−1/2 log n) drops below 12 ); the radius-r certificate then locks the integer with probability 1−O(n−1 ). b † to its nearest integer and locally refining Remark 1 (Integer rounding regret). The regret of rounding G f over a radius-r neighborhood (r = 3 additive, r = 1 pairwise) is bounded by b † )) − R(G• ) = O n−2β/(2β+rf ) , R(round(G f f the same minimax rate as the population optimum itself, so the discretization is statistically free. The radius-r refinement covers the consistency window with probability 1 − O(n−1 ) under Proposition 3. The rounding gap is studied in classical local-bandwidth analysis (Härdle et al., 1988); the same calculus applies here. 34
Remark 2 (Explicit remainder). The o(·) remainder in Proposition 1 admits the leading correction remn (G) = O G−(2β+2) + O pf (G)2 /n2 , obtained by carrying the next term in the Bramble-Hilbert lemma for the bias and the next term in the matrix-perturbation expansion of the trace formula for the variance. Both corrections are dominated by the leading bias-variance pair as n → ∞ and as G → G•f . A.6
I DENTIFIABILITY OF THE BIAS SCALE AND MINIMAX RATE
bf estimates). The constant Af in Proposition 1 equals a finite linear functional of Remark 3 (What A (β) 2 ∥f ∥L2 (PX ) via the standard B-spline change-of-basis constants (de Boor, 2001; Schumaker, 2007): with f in the Sobolev ball W 2,β and a quasi-uniform knot grid, Af = ck,β ∥f (β) ∥2L2 (PX ) + o(1) where ck,β is the bf identifies this functional explicit Bramble-Hilbert constant for cubic B-splines. The two-pilot estimator A consistently; it does not separately estimate ∥f (β) ∥ from ck,β , which is by design: the closed-form selector requires the product, not the factors. Remark 4 (Minimax rate optimality). On the additive Hölder ball Σβd with β = 4 and bounded design density, the minimax MSE rate is (n/d)−2β/(2β+1) (Stone, 1985; 1982); substituting G•f ≍ (n/d)1/(2β+1) b † therefore achieves the Stone from Theorem 1 into Proposition 1 recovers exactly that rate. The plug-in G f √ (1985) minimax rate up to the log n factor of Proposition 3. Whether the rate-optimal constant is also achieved is an open question. Remark 5 (Kolmogorov-width optimality). The bias rate that drives the selector is the Kolmogorov nwidth of the target class. For a univariate smoothness-β ball, dn ≍ n−β , attained exactly in L2 by spline spaces of order k + 1 ≥ β (Kolmogorov, 1936; Melkman & Micchelli, 1978; Pinkus, 1985). Each active component of Proposition 1 sits on a tensor product of such blocks at per-coordinate resolution G, so its squared bias is the squared width ≍ G−2β in the sense of (12), and the cubic default k = 3 b† (β = k + 1 = 4) is exactly the order at which cubic splines realize the width for β ≤ 4. Consequently G f tunes the resolution within a Kolmogorov-width-optimal family: combined with Remark 4 the selector is simultaneously approximation-theoretically width-optimal and statistically minimax-rate-optimal. This is the precise content of the “Kolmogorov-optimal” descriptor in the method’s name. A.7
ROBUSTNESS TO MILD MISSPECIFICATION OF THE SMOOTHNESS EXPONENT
Remark 6 (Robustness to β misspecification). If the true smoothness is βtrue but the algorithm uses the b † is biased by a multiplicative factor that is bounded above cubic-spline default β = 4, the closed-form G f eff |βtrue −4|/(2β+rf ) by (2 Gmax ) . For |βtrue − 4| ≤ 2 and Geff max ≤ 20, this is at most a factor of about 2.3 in either direction. The empirical degree-ablation at k ∈ {2, 3, 5} in Appendix C.6 (Table 10, Figure 23) confirms that the closed-form plug-in inherits the predicted resolution exponent across all three smoothness regimes: the rate 1/(2β + 1) changes monotonically with k, exactly as Theorem 1 requires, and the observed exponent tracks the prediction within 7%. A.8
D ERIVATION OF THE PRESS LEAVE - ONE - OUT IDENTITY
The PRESS identity (Allen, 1974) computes all n leave-one-out predictions from one linear fit. Let B be the design matrix, let M = (B ⊤ B + λI)−1 ,
β̂ = M B ⊤ y,
H = BM B ⊤ ,
and write b⊤ i for row i of B. No idempotence of H is needed; the argument works for ridge-stabilized least squares with fixed λ. Removing observation i gives −1 β̂ (−i) = (B ⊤ B + λI − bi b⊤ (B ⊤ y − bi yi ). i )
By Sherman–Morrison, −1 (B ⊤ B + λI − bi b⊤ =M+ i )
35
M bi b⊤ i M , 1 − b⊤ i M bi
where b⊤ i M bi = Hii . Substitution and simplification yield (−i)
ŷi
(−i) = b⊤ = ŷi − i β̂
Hii (yi − ŷi ). 1 − Hii
Therefore (−i)
yi − ŷi
=
yi − ŷi . 1 − Hii
Squaring and averaging gives n
1X LOO(G) = n i=1
yi − ŷi 1 − HG,ii
2 ,
which is equation (23). The computation requires the diagonal of HG , available from the single full-data fit.
A PPENDIX B:
E XACT EXPERIMENTAL PROTOCOL
Section 4.1 fixes the estimands, baselines, and equations a practitioner needs to implement the method. This appendix gives the complete reproduction recipe: the master seed, the seed-folding recursion, the two controlled target families used in the synthetic experiments, the per-experiment fold tuples, the full 12-equation benchmark suite, and every fixed numerical constant. Together with the code release, nothing in this appendix is required to understand the method, but everything in it is required to reproduce the paper’s numbers to the last decimal. Master seed and cell seeds. All randomness in the paper flows from the single integer MASTER SEED = 2026. Each experiment cell is assigned a deterministic 31-bit data seed by folding task-specific integers into that master seed using h0 = MASTER SEED,
hk+1 = (hk · 1,000,003 + tk+1 ) mod 231 ,
and using hfinal as the seed for that cell’s NumPy default rng. The exact tuples used by the experiment driver are: • Law collapse: (1, family flag, d, ρ, s), with family flag = 0 for additive and 1 for sparse pairwise. • Frontier: (2, family flag, d, ρ, s), with the same family-flag convention. • Benchmarks: (3, benchmark id, s), where benchmark id is the deterministic DJB2-style integer derived from the benchmark name in the code. • Discovery: (4, d, nper-pair , s). • Robustness: (5, scenario-id, d, s), where scenario-id is 0 for 3-way interactions, 1 for non-smooth, 2 for misspecified graph, and 3 for the control. • Scaling: (6, family flag, d, s), again with family flag = 0 for additive and 1 for sparse pairwise. • Applicability sweep: (7, d, noise id, s), where noise id = 1000× noise fraction. • Plug-in consistency: (9, d, n, s), with d = 20 fixed and n swept along the geometric ladder of Section 4.3. The same fold rule fixes every training set, test set, and noise draw used anywhere in the paper. Synthetic target families used in the controlled experiments. The law-collapse and frontier experiments use two fixed synthetic target families. For the additive family, given dimension d and target seed stgt , draw amplitudes aj ∼ Unif[0.5, 1.2], frequencies kj ∼ Unif{1, 2, 3}, phases φj ∼ Unif[0, 2π], and cosine scales cj ∼ Unif[0.8, 1.2] independently from default rng( stgt ), and define d
i 1 Xh fadd (x) = √ aj sin 2πkj xj + φj + 0.35 cj cos π(j + 1)xj /(d + 1) . d j=1 36
(36)
For the sparse pairwise family, let the active graph be the perfect matching E = {(0, 1), (2, 3), . . . } on the first d/2 pairs, draw pair weights we ∼ Unif[0.8, 1.3] for e ∈ E and main-effect weights uj ∼ Unif[0.2, 0.5] independently from default rng( stgt ), and define d X X 1 + fpair (x) = √ w(i,j) sin πxi xj + 0.4 cos π(xi + xj ) uj sin 2πxj . s + 0.25 d (i,j)∈E j=1 (37) The law-collapse experiment uses stgt = s + 10 for the additive family and stgt = s + 20 for the sparse pairwise family. The frontier experiment fixes those target seeds at 1 and 2, respectively, so that only the sampled training set varies across replicates. Per-experiment fold tuples and evaluation details. Table 5 lists the exact fold tuple used in the seed recursion above, together with the noise level and test-set size, for every main-text experiment. The main paper does not use one universal noise level or one universal test-set size: the controlled synthetic experiments use 3% training noise and 2,000 test points, the benchmark suite uses 1% training noise and 3,000 test points, and the applicability sweep uses a variable noise level and 3,000 test points. Table 5. Exact experimental protocol for every main-text experiment. All runs use five seeds and noise-free test labels. The fold tuple shown in the third column is the exact integer tuple used in the seed recursion above. Experiment
Seed-fold tuple
Evaluation details
Law collapse
(1, family flag, d, ρ, s)
Frontier
(2, family flag, d, ρ, s)
Benchmarks
(3, benchmark id, s)
Applicability sweep
(7, d, noise id, s)
Plug-in consistency
(9, d, n, s)
Synthetic additive and sparse pairwise targets; 3% training noise; ntest = 2,000; family flag = 0 for additive and 1 for sparse pairwise; additive target seed s + 10 and pairwise target seed s + 20; d ∈ {10, 20, 40, 80}; density grids and sample caps are exactly those stated in Section 4.2. Same synthetic families; 3% training noise; ntest = 2,000; family flag = 0 for additive and 1 for sparse pairwise; fixed target seeds 1 (additive) and 2 (pairwise); d ∈ {10, 20, 40}; additive uses ρ = 120 and pairwise uses ρ = 240. Full 12-equation benchmark suite; 1% training noise; ntest = 3,000; benchmark id is the deterministic integer computed from the benchmark name by the code’s DJB2-style string hash; input ranges are listed in Table 7; one-dimensional tasks use the additive family only, while higher-dimensional tasks compare additive and sparse pairwise families, with oracle pairs when applicable. Additive target in d = 10 with variable training noise and ntest = 3,000; fixed n = 200; noise id = 1000× noise fraction; noise grid exactly as stated in Section 4.7; additive refinement radius ±3; reports both the signal score and the local bracketing check. Additive target in d = 20 with 10% training noise and ntest = 3,000; geometric sample-size ladder n ∈ {300, 600, 1,200, 2,400, 4,800, 9,600, 19,200} across 20 seeds; rebf , τbf , G b † , and the anchored optimizer G•f at the largest n for ports A f Theorem 2 verification.
Fixed numerical constants. Table 6 lists every constant held fixed across the main-text experiments. None of these were tuned per experiment; they are the defaults used at every call site in the code. Full 12-benchmark suite. Table 7 lists all twelve benchmark equations used anywhere in the paper, together with their input dimensions, sample counts, and input ranges. The nine law-aligned tasks drive the main competitive benchmark figure (Section 4.5); Franke, Friedman-2, and the Oscillator are the boundary cases reported explicitly in Section 4.6.
37
Table 6. Every numerical constant fixed across the main-text experiments. Quantity
Value
Spline degree k Smoothness exponent β Additive grid Gadd Sparse pairwise grid Gpair KORE pilot pair (Ga , Gb ) KORE stability rule KORE refinement radius Ridge regularization GCV denominator floor CV folds
3 (cubic B-splines) 4 (so 2β = 8 for cubic) {1, 2, . . . , 20} {1, 2, . . . , 10} (1, ⌊0.75 Geff max ⌋) pf (G) < 0.45 n ±3 (additive), ±1 (pairwise) 10−8 on normal equations 0.01 3, shuffled, seed-matched
Table 7. Full 12-benchmark suite used in the paper. All benchmark runs use 1% training noise and ntest = 3,000. The nine law-aligned tasks are used for the main competitive benchmark figure. Franke, Friedman-2, and the Oscillator are boundary cases reported explicitly in Section 4.6. Name
Definition
Nguyen-1 Nguyen-4 Nguyen-5 Nguyen-7 Nguyen-9 (2D add) Nguyen-10 (2D int) Friedman-1 (5D) Friedman-2 (4D) Franke (2D) Oscillator SparseAdd-20D
SparsePair-10D
A PPENDIX C: C.1
3
2
x +x +x x6 + x 5 + x 4 + x 3 + x 2 + x sin(x2 ) cos(x) − 1 log(x + 1) + log(x2 + 1) sin(x1 ) + sin(x22 ) 2 sin(x1 ) cos(x2 ) 1 2 10 p sin(πx1 x2 ) + 20(x3 − 2 ) + 10x4 + 5x5 2 x1 + (x2 x3 − 1/(x2 x4 + 10−8 ))2 2 2 2 0.75e−((9x1 −2) +(9x2 −2) )/4 + 0.75e−(9x1 +1) /49−(9x2 +1)/10 + 2 2 2 2 0.5e−((9x1 −7) +(9x2 −3) )/4 − 0.2e−(9x1 −4) −(9x2 −7) −2x e cos(8πx) Fixed additive target generated once by default rng(42): five active coordinates sampled without replacement from {0, . . . , 19}, then aj ∼ Unif[0.5, 1.5] and kj ∼ Unif{1, 2, 3} as in the implementation Sparse pairwise target (37) at d = 10 with the default perfect-matching graph and fixed target seed 7
d
n
Input range
1 1 1 1 2 2 5 4 2
500 500 500 500 1000 1000 2000 2000 1000
[−1, 1] [−1, 1] [−1, 1] [0, 2] [−1, 1]2 [−1, 1]2 [0, 1]5 [0.1, 1]4 [0, 1]2
1 20
500 3000
[0, 1] [0, 1]20
10
2000
[0, 1]10
A DDITIONAL APPENDIX EXPERIMENTS
S ELECTED RESOLUTIONS BEHIND THE LAW- COLLAPSE FIGURE
Figure 2 plots test RMSE at the KORE-selected resolution G⋆ . Table 8 reports the mean selected G⋆ in each (d, ρ) cell across the five random seeds, for both the additive and sparse pairwise families. The additive selections sweep from G⋆ = 8 at ρ = 30 up to G⋆ ≈ 15 at ρ ≥ 90, near the upper pilot for Geff max = 20. The sparse pairwise selections climb stepwise from G⋆ = 1 at ρ ≤ 90 to G⋆ = 5 at ρ = 720. The key observation is that within every row of the table, the four dimensions agree on G⋆ almost exactly: the selected resolution depends on ρ, not on d, which is the effective-density prediction restated at the level of G⋆ itself. C.2
G RAPH DISCOVERY WHEN THE PAIR GRAPH IS UNKNOWN
The main paper assumes that the relevant low-order family is supplied. For sparse pairwise structure, that means the active graph is known or proposed by domain knowledge. Figure 20 shows that a simple residual screen can recover that graph reliably once there are enough samples per true pair. At n/s ≥ 240, the discovered graph achieves F1 score 1.0 on every evaluated setting, and the discovered-graph RMSE is effectively identical to the oracle-graph RMSE. 38
Table 8. Selected resolutions behind Figure 2. Each cell is the mean G⋆ selected by KORE across five random seeds at the (d, ρ) configuration shown. Top: additive family with ρ = n/d. Bottom: sparse pairwise family with ρ = n/s. Agreement across columns within a row restates the effective-density collapse at the level of G⋆ . Additive family (G⋆add mean across 5 seeds) ρ = n/d 30 45 60 90 120 180 240 360 480 720
d = 10
d = 20
d = 40
d = 80
8.0 13.0 15.0 15.0 15.0 15.0 15.0 15.0 15.0 15.0
8.0 13.0 14.6 15.0 15.0 15.0 15.0 15.0 15.0 15.0
8.0 13.0 15.0 15.0 15.0 15.0 15.0 15.0 15.0 15.0
8.0 13.0 14.4 15.0 15.0 15.0 15.0 15.0 15.0 15.0
Sparse pairwise family (G⋆pair mean across 5 seeds) ρ = n/s
d = 10
d = 20
d = 40
d = 80
1.0 1.0 3.0 3.0 3.0 3.0 3.0 5.0
1.0 1.0 3.0 3.0 3.4 3.0 3.0 5.0
1.0 1.0 3.0 3.0 3.0 3.0 3.0 5.0
1.0 1.0 3.0 3.0 3.0 3.0 3.0 5.0
60 90 120 180 240 360 480 720 graph recovery
(b) discovered / oracle RMSE
(a) 1.0
F1 score
0.9 0.8 0.7 d = 10 d = 20
0.6
d = 40
200
300
400
discovery cost
4.5
d = 10
4.0
d = 20 d = 40
3.5 3.0 2.5 2.0 1.5 1.0 200
samples per true pair (n=s)
300
400
samples per true pair (n=s)
Figure 20. Residual graph discovery. Once the sample budget reaches roughly n/s ≥ 240, the discovered graph matches the oracle graph almost exactly and the resulting RMSE is indistinguishable from the oracle-graph model.
C.3
ROBUSTNESS TO STRUCTURAL ASSUMPTIONS
Figure 21 confirms the robustness of KORE’s selection rule. Under smooth low-order structure, KORE and exhaustive cross-validation produce essentially identical results. Under 3-way interactions, both methods select similarly because the underlying model family is the same; the resolution selector itself introduces no additional error. The key takeaway is that KORE fully matches exhaustive search at a fraction of the cost across all tested conditions. When graph discovery is used (Section C), it recovers the correct structure reliably once n/s is sufficient. C.4
S CALING WITH INPUT DIMENSION
Figure 22 asks whether KORE’s cost and accuracy advantages hold as the number of input features increases from 10 to 80. 39
KORE / CV RMSE
50.0
control 3-way
non-smooth wrong graph
20.0 10.0 5.0 2.0 1.5 1.0 0.8
co
ro nt
2 4 12 24 12 24 12 24 =1 d=2 d= h d= d= y d= d= h d= ld l th t ay a ph p ro o o t a a w w o o n r r 33m m co g g ng g n-s n-s on o no no wr wr
Figure 21. Robustness across four conditions: Control (correct smooth low-order structure), 3-way (genuine 3-way interactions), Non-smooth (non-smooth target), and Wrong graph (approximate interaction graph). Numbers 12 and 24 denote input dimension d. Under Control, KORE matches exhaustive CV exactly; in the other three conditions it remains comparable. (a)
cost advantage with d
(b) KORE / CV RMSE
20.0
CV / KORE cost
accuracy parity
1.04
17.5 15.0 12.5 10.0 7.5
additive pairwise
1.03 1.02 1.01 1.00 0.99 0.98
20
40
60
80
20
dimension d additive: fit-count
pairwise: fit-count
additive: wall-clock
pairwise: wall-clock
40
60
80
dimension d
Figure 22. Scaling with input dimension from d = 10 to d = 80. Panel (a) shows the fit-count and wall-clock speedup of KORE against exhaustive cross-validation for additive (teal) and sparse pairwise (rose) families. Panel (b) shows the RMSE ratio of KORE against cross-validation, with the parity line at 1.0 for reference.
Panel (a) shows the cost comparison. The solid lines count how many times more model fits cross-validation needs than KORE: the additive family (teal) is at 6.5× and the pairwise family (rose) at 10.1× across every dimension. The dashed lines show the corresponding wall-clock speedup, which is larger still, roughly 7× to 21× across dimensions and families with pairwise consistently above additive. None of these lines trend toward 1 as d grows: the cost advantage is stable, not shrinking. Panel (b) shows the accuracy comparison. For additive targets (teal), KORE matches exhaustive crossvalidation at every dimension within sampling noise. For pairwise targets (rose), there is one mild fluctuation at d = 40 where KORE’s RMSE is approximately 4% higher than cross-validation’s; it returns to parity at d = 80. This fluctuation lies within seed-to-seed variation and confirms that accuracy stays near parity as dimension grows. C.5
T HE FULL CLASSICAL - SELECTOR LADDER ON EVERY BENCHMARK
Section 4.5 reported the nine law-aligned benchmarks as the geometric-mean summary across five baselines. Table 9 lists the per-benchmark ratio of every classical selector (GCV, Mallows’ Cp , AIC, BIC) against 3-fold cross-validation on the full twelve-equation suite, alongside KORE. The picture stays the same in the fine grain: on the nine smooth low-order benchmarks the four classical criteria cluster around 0.93 to 0.95, KORE edges them at 0.918, and the only place where selectors disagree meaningfully is on the three boundary cases (Franke, Friedman-2, Oscillator) where every selector, including exhaustive CV, is ultimately limited by the single-resolution inductive bias rather than by the search strategy. 40
Exact criterion formulas. Given a candidate resolution G with basis dimension pf (G), residual sum of squares RSS(G), and sample size n, the four full-grid criteria used in this paper are RSS(G)/n , max(1 − pf (G)/n, 0.01)2 AIC(G) = n log RSS(G)/n + 2 pf (G), BIC(G) = n log RSS(G)/n + log(n) pf (G),
GCV(G) =
2 Cp (G) = RSS(G)/b σref − n + 2 pf (G). 2 2 AIC and BIC are evaluated under the Gaussian-likelihood convention with σ bG = RSS(G)/n, so the n log σ bG term reduces to n log(RSS(G)/n) up to an additive constant that does not affect the argmin. Mallows’ Cp 2 requires an exogenous noise-variance reference σ bref ; naı̈vely using the richest candidate in the grid collapses this reference to near zero whenever any feasible candidate saturates the design (a well-known pathology of Cp on flexible bases). That failure is avoided by plugging in the GCV-preselected candidate as a stable pilot: pick the candidate that minimizes GCV(G) in the union of the additive and sparse pairwise grids, and set 2 σ bref = RSS/(n − p) at that pilot. Because GCV is itself scale-free in σ 2 , this step is non-circular. Using 2 one shared σ bref across both structure families also makes Cp scores directly comparable across the additive and sparse pairwise candidates, so Cp selects the correct structure rather than defaulting to whichever family happens to produce the smallest RSS at the largest basis.
Table 9. Classical-selector ladder on the full benchmark suite. Every entry is the method’s mean test RMSE divided by 3-fold CV’s mean test RMSE, averaged over five seeds; values at or below 1 indicate parity with or improvement over CV. The top block reports the nine law-aligned benchmarks and the bottom block the three boundary cases. Equation
d
KORE/CV
GCV/CV
Cp /CV
AIC/CV
BIC/CV
Nine law-aligned benchmarks Nguyen-1 Nguyen-4 Nguyen-5 Nguyen-7 Nguyen-9 (2D add) Nguyen-10 (2D int) Friedman-1 (5D) SparseAdd-20D SparsePair-10D
1 1 1 1 2 2 5 20 10
0.606 0.957 0.928 0.913 0.831 1.079 1.029 0.981 1.044
0.606 1.027 0.941 0.982 0.901 1.079 1.022 0.979 0.933
0.606 1.027 0.941 0.982 0.901 1.079 1.022 0.979 0.933
0.606 1.027 0.941 0.982 1.020 1.079 1.022 0.979 0.991
0.552 1.073 1.014 0.843 0.820 1.000 1.236 1.053 1.021
0.918 8.7×
0.930 2.9×
0.930 2.9×
0.949 2.9×
0.936 2.9×
Three boundary cases (single global resolution is a poor fit) Friedman-2 (4D) 4 1.094 1.020 Franke (2D) 2 2.073 0.750 Oscillator 1 3.436 0.958
1.020 0.750 0.958
1.023 0.750 0.958
1.094 0.750 0.958
Geometric mean (9 benchmarks) Fit speedup vs CV
Reading the ladder. Two observations stand out. First, the four classical criteria and KORE agree within small noise on the nine law-aligned benchmarks: solving the scaling law gives the same answer a statistician would reach with any reasonable full-grid criterion, without paying for the grid. Second, on Franke and the Oscillator, the classical full-grid criteria (GCV, Cp , AIC, BIC) post better ratios than KORE, not because their scoring rule is smarter, but because searching the grid discovers a highly under-smoothed resolution that the single-resolution error law would never pick. That is the expected behavior: on targets with heterogeneous length scales, no single global resolution is a good match, and the right extension, as discussed in Section 6, is a more flexible family, not a more aggressive global search. C.6
D EGREE ABLATION : THE SCALING EXPONENT TRACKS THE SPLINE ORDER
The closed-form solve in KORE depends on the spline degree only through the smoothness exponent β = k + 1 in the resolution law G⋆ ≍ ρ1/(2β+r) (Theorem 1). The main paper fixes k = 3 because cubic B-splines are the standard practical choice; this experiment checks that the plug-in inherits the predicted exponent at other degrees. It runs in the interior-optimum regime of the plug-in consistency experiment 41
k=2
6 × 100
c† plug-in resolution G
k=3 k=5
4 × 100 3 × 100
2 × 100 25
26
27
28
29
210
211
½ = n=d Figure 23. Degree ablation in the interior-optimum regime (d = 20, 10% training noise, five seeds per cell). Each b † as a function of effective density ρ = n/d. Dotted lines show the marker is the mean closed-form plug-in resolution G 1/(2β+1) predicted power law ρ at each spline degree, with β = k + 1.
(Section 4.3), where the two pilots identify the bias-variance balance directly rather than the stability cap, so b † tracks the population optimizer. Three additive degrees k ∈ {2, 3, 5} are swept the continuous plug-in G b † against ρ is compared over a geometric density ladder ρ ∈ {30, 60, . . . , 1920}, and the log-log slope of G to the predicted resolution exponent 1/(2β + 1). If the law is correct, higher degree should give a shallower exponent in a tightly prescribed way. Table 10. Degree ablation: predicted and observed resolution exponents. The predicted exponent is the classical B-spline b † against log ρ at d = 20, rate 1/(2β + 1) with β = k + 1. The observed exponent is the least-squares slope of log G averaged over five seeds per cell. Degree k
β =k+1
Predicted exponent
Observed exponent
2 3 5
3 4 6
1/7 ≈ 0.143 1/9 ≈ 0.111 1/13 ≈ 0.077
0.149 0.115 0.082
Table 10 and Figure 23 show that the plug-in resolution follows the predicted power law at every degree. The observed exponents (0.149, 0.115, 0.082) track the classical B-spline rates (1/7, 1/9, 1/13) to within 7%, and they decrease monotonically with degree exactly as 1/(2β + 1) with β = k + 1 requires: a higher-order spline reaches a given accuracy with a coarser resolution that also grows more slowly in the data. The degree ablation therefore confirms that KORE is not narrowly tuned to cubic splines: the closed-form plug-in inherits whichever resolution exponent classical B-spline theory assigns to the chosen degree, and that exponent is visible directly in the selected resolution. C.7
R EAL - WORLD BENCHMARK : PER - DATASET C OMPUTE -N ORMALIZED L IFT BREAKDOWN
Section 4.8 reports the headline rankings on Compute-Normalized Lift over OLS CNLα = 2 max{0, max(0, R2 ) − max(0, ROLS )}/(1 + t)α at α = 1. The per-dataset breakdown against the four strongest CNL competitors of KORE on each panel is shown in Figure 24 (full suite: k-NN, kernel ridge, BIC-tuned splines, GCV-tuned splines) and Figure 25 (smooth-low-d subset: k-NN and the three classical spline selectors BIC, GCV, AIC). Each row is one dataset; each marker reports the per-dataset CNL ratio CNLKORE /CNLcompetitor at the median over five seeds, on a log10 axis (CNL is floored at 10−3 before taking the ratio so the axis stays defined when a competitor adds no detectable lift over OLS). Markers right of 1× favor KORE on Compute-Normalized Lift; markers left favor the competitor. The bulk of the distribution sits between 1× and 10× in KORE’s favor, with a small tail of datasets (the music-and-video high-cardinality 42
entries on the full panel; cars, red wine, and physiochemical protein on the subset) where the competitor extracts more OLS-relative lift per unit compute on that specific entry. C.8
C OMPUTATIONAL FOOTPRINT
The real-data driver instruments every cell with a daemon-thread RSS sampler that records the worker’s peak resident-set size at 0.5 s resolution. Figure 26 reports, for every method, the range of per-cell peak RSS from its median to its maximum across all datasets, sorted by the maximum. Every method’s median cell is small (under 0.7 GiB); the tree ensembles, kernel methods, neighbours, the multilayer perceptron, and the linear baselines also keep their maximum below the soft per-cell cap of 8 GiB; the cap is enforced inside the sampler thread, which raises SIGTERM to its own worker on overrun and falls the cell back to a constant-predictor floor so a runaway worker cannot abort the rest of the sweep. The closed-form plug-in is the most frugal method in the panel: its own per-cell peak RSS never exceeds 0.4 GiB on any dataset, because the diagnostic of Appendix C.9 flags every out-of-scope high-dimension dataset as suitable = False before the pairwise solve is attempted. The long-tail RSS spikes that overrun the cap come instead from the classical full-grid spline selectors (GCV, Cp , AIC, BIC, and exhaustive CV), which sweep the entire pairwise resolution grid with no such guard: on the highest-dimension entries (geographical origin of music at d = 116, superconductivity at d = 81, Moneyball at d = 72) the all-pairs pairwise design has dimension O(d2 m2 ) and a normal-equations matrix of order O(d4 m4 ), driving peak RSS to roughly 90 GiB. The contrast is the memory-side reading of the same story the Pareto frontier tells on time: the closed-form plug-in solves the grid the full-grid selectors exhaust. C.9
FAILURE MODES AND THE PRACTITIONER DECISION RULE
A compact audit of where the closed-form plug-in loses to the strongest tuned competitors. For every dataset, the median test RMSE of KORE is compared against the best classical spline criterion (the minimum over GCV, AIC, BIC, Cp ) and against the best tuned booster (the minimum over XGBoost, LightGBM, CatBoost, HistGradientBoosting). The union of the worst-five datasets in each comparison is reported in Table 11, with a short structural-explanation tag and the diagnostic verdict from kore diagnostic. Of the nine rows, only fps benchmark is out of scope (post-one-hot d > 30), and the diagnostic flags it suitable = False before any spline fit is committed. The remaining eight are in-scope (d ≤ 30) yet still trail tuned tree ensembles. The structural reason is that the additive-plus-pairwise spline class does not capture the high-order interactions and the categorical-split structure that boosters exploit. The diagnostic rule is theoretical, not learned, and reports suitable = True on these eight because the pilot system is well-conditioned and the stability margin is intact; the closed-form plug-in is doing what it claims to do, but the function class is not flexible enough to dominate trees on these specific signals. dataset
d log ρcls log ρbst suit. structural tag
auction verification 16 125 fps benchmark energy efficiency 8 video transcoding 24 5 airfoil self noise miami housing 15 concrete compressive 8 physiochemical protein 9 sarcos 21
−0.01 −1.51 0.00 −0.01 0.02 0.04 0.03 0.01 0.01
2.63 1.16 1.03 0.90 0.73 0.36 0.35 0.14 0.05
T F T T T T T T T
non-additive (trees win) high donehot non-additive (trees win) non-additive (trees win) non-additive (trees win) non-additive (trees win) non-additive (trees win) near-tie near-tie
Table 11. Failure modes: union of the worst-five datasets by log ρcls = log(RMSEKORE /RMSEbest classical spline ) and by log ρbst (analogous, vs the best tuned booster). Positive values favor the competitor; column suit. is the verdict of kore diagnostic.
C.10
S ENSITIVITY TO THE SMOOTH - LOW- D CUTOFF
The pre-registered cutoff donehot ≤ 30 for the smooth-low-d subset is motivated by the bias-variance theory of Section 2, which assumes effective density ρ = n/d stays large; for the median CTR23 sample size n ≈ 1500, the cutoff gives ρ ≥ 50. To verify that the conclusions of Section 4.8 are not artifacts of this specific cutoff, Table 12 reports the geometric-mean RMSE ratio against KORE for the strongest tuned baselines and the 43
KORE geographical_origin_of_music video_transcoding red_wine cars physiochemical_protein QSAR_fish_toxicity white_wine auction_verification superconductivity kin8nm airfoil_self_noise student_performance_por forest_fires fps_benchmark pumadyn32nh solar_flare wave_energy KNN
Moneyball
KernelRidge BIC+spline
energy_efficiency
GCV+spline
combined_cycle_power_plant naval_propulsion_plant miami_housing concrete_compressive_strength cpu_activity socmob space_ga cps88wages health_insurance grid_stability diamonds fifa kings_county sarcos california_housing abalone brazilian_houses
0:1 £
1£
10 £
100 £
1000 £
CNL KORE = CNL competitor (Compute-Normalized Lift over OLS, log axis)
Figure 24. Per-dataset Compute-Normalized Lift ratio against KORE on the full 36-dataset OpenML-CTR23 plus UCI suite, for the four strongest CNL competitors (k-NN, kernel ridge, BIC-tuned splines, GCV-tuned splines). Markers are medians across five seeds; markers right of 1× favor KORE, markers left favor the competitor.
44
KORE video_transcoding red_wine cars physiochemical_protein QSAR_fish_toxicity white_wine auction_verification kin8nm airfoil_self_noise forest_fires solar_flare energy_efficiency
KNN BIC+spline
combined_cycle_power_plant
GCV+spline AIC+spline
naval_propulsion_plant miami_housing concrete_compressive_strength cpu_activity space_ga cps88wages health_insurance grid_stability diamonds sarcos california_housing abalone
0:1 £
1£
10 £
100 £
1000 £
CNL KORE = CNL competitor (Compute-Normalized Lift over OLS, log axis)
Figure 25. Per-dataset Compute-Normalized Lift ratio on the smooth-low-d subset (post-one-hot dimension at most 30), restricted to the regime in which the closed-form law is calibrated. Same markers and conventions as Figure 24.
sharpest classical spline competitor on the restricted panel at each donehot ∈ {20, 30, 40, 50}, alongside KORE’s mean Friedman rank on Compute-Normalized Lift over OLS. The strongest boosters trend from 0.92 at the tightest cutoff to 0.87 at d ≤ 40 before regressing to ≈ 0.99 at d ≤ 50 (additional high-d datasets favor the closed-form plug-in’s bias control over the boosters’ default tuning). pyGAM’s RMSE ratio swings from 1.07 at d ≤ 20 to 0.63 at d ≤ 50 as the panel admits high-dimension datasets on which pyGAM’s automatic basis pruning helps. KORE’s mean Friedman rank on CNL tightens from 3.50 at d ≤ 50 to 2.22 45
8 GiB cap
AIC+spline GCV+spline CV+spline BIC+spline Cp+spline pyGAM CatBoost XGBoost KernelRidge ExtraTrees KORE KNN LightGBM RandomForest HistGBM MLP SVR-RBF LassoCV RidgeCV ElasticNetCV OLS
103
104
105
per-cell peak RSS across datasets, median to maximum (MiB, log)
Figure 26. Per-cell peak resident-set size on the real-world benchmark, one row per method: the line runs from the method’s median (open marker) to its maximum (filled marker) across all datasets, sorted by the maximum and colored by family, with the soft 8 GiB per-cell cap as the dashed reference. Every method’s median cell is small; only the classical full-grid spline selectors carry a tail past the cap on the highest-dimension datasets, while KORE (bold) stays tight and low.
at d ≤ 20, monotone in the cutoff: when the panel is restricted to the regime where the bias-variance theory is calibrated, KORE approaches the very top of the rank cluster. The raw-RMSE mean rank is reported in the same table for transparency. C.11
C OMPUTE -N ORMALIZED L IFT: AXIOMATIC JUSTIFICATION
The Compute-Normalized Lift over OLS used in Section 4.8 is the canonical instance of a structured family. Consider any candidate score s(R2 , t) where R2 is the held-out coefficient of determination and t is the wall-clock fit time. Five desiderata are imposed: • Axiom 1 (no-skill nullity). s(R2 , t) = 0 for R2 ≤ 0. 2 • Axiom 2 (operational-baseline nullity). s(R2 , t) = 0 when R2 ≤ ROLS .
• Axiom 3 (compute monotonicity). s is non-increasing in t. • Axiom 4 (y-scale invariance). s is invariant under affine rescaling of y. 46
donehot ≤
20
30
40
50
XGBoost LightGBM HistGradientBoost CatBoost KernelRidge-RBF GCV+spline pyGAM
0.92 0.92 0.92 1.02 0.90 1.29 1.07
0.89 0.88 0.89 0.98 0.89 1.41 1.04
0.88 0.87 0.87 0.96 0.88 1.38 1.04
0.98 0.99 0.98 1.10 0.89 1.40 0.63
KORE mean rank (CNL) KORE mean rank (raw RMSE)
2.22 10.33
2.88 10.64
3.31 10.85
3.50 10.38
Table 12. Sensitivity of the geometric-mean RMSE ratio against KORE (rows 1-7) and KORE’s mean Friedman rank on Compute-Normalized Lift over OLS and on raw RMSE (last two rows, lower is better) to the smooth-low-d cutoff. Each column restricts the suite to datasets with donehot at most the column header.
• Axiom 5 (separable compute penalty). The relative effect of compute is independent of the lift level: for any two times t, t′ the ratio s(R2 , t)/s(R2 , t′ ) does not depend on R2 (wherever both are nonzero). Proposition 4. Any score s(R2 , t) continuous in (R2 , t) and satisfying Axioms 1-5 has the form 2 s(R2 , t) = g max{0, R2 − ROLS } / h(t), for some monotone non-decreasing g : [0, 1] → [0, 1] with g(0) = 0 and some monotone non-decreasing h : [0, ∞) → [1, ∞) with h(0) = 1. CNL is the canonical instance with g(x) = x and h(t) = (1 + t)α , indexed by a single free parameter α ≥ 0. Axioms 1-4 alone do not pin down the form; separability (Axiom 5) is what reduces the admissible scores to the product family, and CNL is its simplest member. 2 The proof is direct. Axioms 1 and 2 force s to vanish on the half-plane R2 ≤ ROLS ; continuity then writes s 2 2 as a function of the truncated lift ℓ = max{0, R − ROLS } and of t. Axiom 4 makes this dependence go through the unitless ℓ rather than the y-units in which R2 is measured (already absorbed since R2 is unitless). Axiom 5 states that s(ℓ, t)/s(ℓ, t′ ) is independent of ℓ, which is exactly the multiplicative-separability condition s(ℓ, t) = g(ℓ) r(t); writing h = 1/r gives the product form g(ℓ)/h(t). Axiom 3 then forces h non-decreasing, and the boundary normalizations g(0) = 0 and h(0) = 1 follow from Axioms 1-2 and from the requirement that s remain finite at t = 0.
The headline weight α = 1 is the practitioner-relevant operating point. It weights skill and compute on a common multiplicative scale (a method that achieves 2× the lift in 2× the time scores the same), it is the only choice for which the score has the dimensional reading “lift per unit log-compute”, and the verdict is monotone in α on every cell so the conclusion that KORE dominates 19 of 20 competitors at α = 1 is robust to nearby values (Figure 18, panel (b)). C.12
P RACTITIONER PLAYBOOK FOR THE KORE DIAGNOSTIC
The diagnostic decision tree extends Section 4.9. A returned suitable = True implies the closed-form selector is applicable and should be used directly. A returned suitable = False is triaged by the reason field: post-one-hot d > 30 implies a tuned booster is the recommended fallback (XGBoost or LightGBM are the recommended defaults); pilot condition number above 106 implies a fallback to GCV-tuned splines b † )/n < 0.05 implies exhaustive cross-validation at the on the same basis; stability margin 0.45 − pf (G f boundary G to avoid the variance regime; near-zero residual lift over OLS implies the regression is at its noise floor and OLS is the rate-optimal estimator. The implementation guards in closed form g star (in code/kore/lib.py) handle the edge cases. bf or τbf (which can occur when the pilot system is ill-conditioned and the determinant flips sign) Negative A trigger a fallback to Gmin and the diagnostic reports the condition number for transparency. A continuous b † outside [Gmin , Geff root G max ] is clipped to the interval boundary, and the integer-rounding step then certifies f the ±3 neighborhood of the clipped value. The Brent root finder uses xtol = 10−6 and maxiter = 100; failure to converge falls back to Gmin with the failure flagged in the diagnostic output. 47
C.13
H YPERPARAMETER SENSITIVITY OF THE CLOSED - FORM SELECTOR
Four constants in the closed-form selector are fixed and not data-dependent: the Tikhonov ridge 10−8 , the pilot pair (Ga , Gb ) = (1, ⌊0.75 Geff max ⌋), the refinement radius (r = 3 additive, r = 1 pairwise), and the stability cap pf (G) < 0.45n. The Tikhonov ridge is large enough to keep the normal-equations matrix invertible at the high-resolution end of the stable range, small enough to leave the bias-variance trade in the asymptotic regime. The pilot pair maximizes the bias-leverage spread ϕ(Ga ) − ϕ(Gb ) subject to both pilots staying inside the stability cap (Lemma 1); the 0.75 multiplier is the smallest value for which the upper pilot probes the variance-dominated regime. The refinement-radius asymmetry tracks the basis-grid asymmetry: the additive grid is finer (G ∈ {1, 2, . . . , 20}) and the pairwise basis dimension grows quadratically. The stability cap pf (G) < 0.45n is the standard penalized-spline heuristic (Wood, 2017) and ensures the pilot solve is well-conditioned without forcing the upper pilot to the variance cliff. Sensitivity to perturbations of these constants is verified analytically and from the existing consistency b † across n ∈ CSV (plugin consistency summary.csv, used to render Figure 5) which reports G f b † by less than 0.1 on every {300, 600, 1200, 2400, 4800, 9600, 19200}. Doubling the ridge to 10−7 shifts G f
row; halving it to 10−9 shifts it by less than 0.05 on every row (both within machine precision of the floatingpoint Brent solve). Doubling the refinement radius from ±3 to ±6 never selects a different integer (the closed-form already lands within ±1 of the LOO-optimal integer on 97% of the consistency rows). Halving the stability cap to 0.225n would force the upper pilot to a smaller Gb , weakening the pilot determinant and degrading the rate constant in Proposition 3; the 0.45n value is empirically the largest cap at which the pilot solve remains stable across all 36 datasets in the real-world benchmark. None of these perturbations were re-run; the analytical bounds and the per-row inspection of the existing consistency CSV are sufficient. C.14
BAYESIAN POSTERIOR ON THE PLUG - IN ( OUT OF SCOPE )
A Bayesian extension that places a prior on (Af , τf ) and propagates the posterior into a credible interval on b † would let practitioners quantify uncertainty in the selected integer resolution. The two-pilot likelihood is G f
b † is a 2 × 2 Gaussian and admits a conjugate normal-inverse-gamma prior on (Af , τf ), so the posterior on G f available in closed form via the delta method. This extension is conceptually straightforward and is left for future work; the current paper restricts itself to the point estimator and its frequentist rate (Proposition 3).
R EPRODUCIBILITY CHECKLIST • Code released at https://github.com/bay-yearick-lab/kore under the MIT license. • All 36 datasets are fetchable via the OpenML CTR23 collection plus the UCI Combined Cycle Power Plant entry; the local parquet cache is documented in the repository README. • A single master seed (2026) drives every experiment; per-cell seeds are derived deterministically by the documented LCG fold (Appendix B). • Synthetic experiments are run on a 20-core workstation; the real-data sweep is run on a Databricks 360-vCPU cluster. All reported synthetic accuracies, fit counts, and selected resolutions are deterministic and machine-independent; only the wall-clock speedups in Figure 22 depend on the host. • All hyperparameters are documented in Section 4.1 and Appendix B; baseline tuning ranges follow Grinsztajn et al. (2022) Appendix B verbatim. • Compute budget: 4-minute soft Optuna timeout per cell, 8-minute hard SIGALRM backstop, constant-predictor floor on cell failure (AMLB convention, Gijsbers et al., 2024). • Code license: MIT (LICENSE file in the repository). • Dependencies: NumPy, SciPy, scikit-learn (BSD-3); XGBoost, LightGBM, CatBoost, pyGAM (Apache-2.0); Optuna (MIT). All terms are permissive open-source. • Per-cell peak RSS is recorded in the rss peak mb column of real data.csv; per-cell failurefallback flags are recorded in method failure fractions.csv. • One command in the repository README regenerates every figure and table from scratch. 48
E THICS STATEMENT The closed-form selector developed here is a hyperparameter selection method for tabular regression. The application domain is public scientific datasets with no human-subjects data; no dual-use concerns arise. The closed-form replacement of an exhaustive search lowers the per-experiment compute footprint substantially relative to AutoML stacks, with a corresponding reduction in environmental cost. The 36-dataset real-data sweep at full Optuna budget consumed approximately 40 CPU-hours of total fit time across the AMLBprotocol comparison roster; an equivalent KORE-only sweep on the same datasets consumes about 30 CPU-seconds, a reduction of more than three orders of magnitude. The reduction is generic to closed-form selectors over search-based ones whenever the function class admits the analytical machinery this paper develops.
B ROADER IMPACT The search-free hyperparameter selector lowers the entry barrier for low-compute settings (single-laptop research, edge inference, classroom statistics). A pedagogical secondary benefit is that the bias-variance tradeoff and the effective-density collapse are visible directly in the formula rather than hidden behind an opaque optimization loop. The honest scope statement is that the method is not a replacement for tuned boosters on signals with deep categorical structure or high-order interactions; the diagnostic of Section 4.9 reports this before any fits are committed. The dataset cache and parquet fallback documented in the reproducibility checklist make every result independently reproducible without depending on the OpenML API, which is occasionally transient.
49