Paper deep dive

The Quadrilateral Loss: Additivity as a Measurable Behavior of Dense Neural Networks

Antonio Di Cecco

Year: 2026Venue: arXiv preprintArea: cs.LGType: PreprintEmbeddings: 92

Intelligence

Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 90%

Last extracted: 7/23/2026, 3:13:42 AM

Summary

This paper introduces the quadrilateral loss, a differentiable penalty that treats additivity as a measurable behavior in dense neural networks by penalizing second-order mixed differences. The authors demonstrate that this loss allows for the behavioral enforcement of additivity, which can be crystallized into exact structural additivity. They compare various routes to additivity (structural masks, backfitting, shared-section models, etc.) and introduce diagnostic tools like surrender curves to analyze interaction retention, finding that constraining behavior before structure is often superior to weight-space constraints.

Entities (10)

Quadrilateral Loss · method · 98%California Housing dataset · dataset · 96%Additive Model · concept · 95%TEAM · model · 93%Neural Additive Models (NAM) · model · 92%Backfitting · algorithm · 90%Shapley-GAM · concept · 90%Shared-Section Model · model · 88%WTA Decay · method · 87%Surrender Curves · metric · 85%

Relation Signals (7)

Quadrilateral Loss → measures → Additivity

confidence 95% · We introduce the quadrilateral loss, a differentiable penalty that treats additivity as a measurable behavior

TEAM → istypeof → Additive Model

confidence 92% · TEAM: A Tree Ensemble Additive Model... Additivity is exact and comes from the weak learner

Quadrilateral Loss → equalsexpectation → Shapley-GAM interaction mass

confidence 90% · equals in expectation the per-coordinate interaction mass of the interventional Shapley-GAM

WTA Decay → isproposedby → Di Cecco

confidence 90% · WTA decay (ours)... weight-space alternative of our design

Surrender Curves → reveals → Interaction Retention

confidence 88% · surrender curves... show... that pre-regularization interaction magnitude barely predicts what a regularized model retains

Quadrilateral Loss → enables → Crystallization

confidence 85% · crystallization, an optional post-processing that converts the penalized model into an exactly additive one

Backfitting → comparedwith → Quadrilateral Loss

confidence 80% · calibration against classical backfitting with neural smoothers, which quietly matches modern per-feature architectures

Cypher Suggestions (0)

No Cypher suggestions yet.

Abstract

Abstract:Additive models buy interpretability by forbidding feature interactions, a constraint that neural instantiations enforce architecturally. We introduce the quadrilateral loss, a differentiable penalty that treats additivity as a measurable behavior instead: a second-order mixed difference on pairs of training points swapping one coordinate, which vanishes if and only if the coordinate carries no interaction, remains informative for piecewise-linear networks, and equals in expectation the per-coordinate interaction mass of the interventional Shapley-GAM. The loss turns additivity into a dial - most learned interactions prove removable almost for free, and on small datasets a moderate penalty improves accuracy and additivity simultaneously - and into an online observable: its per-feature surrender curves show, across seeds and datasets, that pre-regularization interaction magnitude barely predicts what a regularized model retains, undermining post-hoc interaction rankings. Against this instrument we compare routes to exact additivity, spanning structural masks, behavioral penalties (optionally crystallized into exact structure), weight decay, backfitting, the shared-section model, and bagged boosted stumps: constraining behavior before structure dominates weight-space constraints, rankings reverse between data regimes, and converging routes agree on the shape functions themselves. Three silent failure modes we document share one anatomy: guarantees imported into settings that quietly void their preconditions.

PDF

Open source PDF →Open local PDF →

Trouble viewing inline? Open PDF directly →

Full Text

90,445 characters extracted from source content.

Expand or collapse full text

The Quadrilateral Loss: Additivity as a Measurable Behavior of Dense Neural Networks Antonio Di Cecco Università “G. d’Annunzio” Chieti–Pescara, and School of AI. Author of Explainable AI with Python[20]. (July 2026) Abstract Additive models buy interpretability by forbidding feature interactions, a constraint that neural instantiations enforce architecturally. We introduce the quadrilateral loss, a differentiable penalty that treats additivity as a measurable behavior instead: a second-order mixed difference on pairs of training points swapping one coordinate, which vanishes if and only if the coordinate carries no interaction, remains informative for piecewise-linear networks, and equals in expectation the per-coordinate interaction mass of the interventional Shapley-GAM. The loss turns additivity into a dial—most learned interactions prove removable almost for free, and on small datasets a moderate penalty improves accuracy and additivity simultaneously—and into an online observable: its per-feature surrender curves show, across seeds and datasets, that pre-regularization interaction magnitude barely predicts what a regularized model retains, undermining post-hoc interaction rankings. Against this instrument we compare routes to exact additivity, spanning structural masks, behavioral penalties (optionally crystallized into exact structure), weight decay, backfitting, the shared-section model, and bagged boosted stumps: constraining behavior before structure dominates weight-space constraints, rankings reverse between data regimes, and converging routes agree on the shape functions themselves. Three silent failure modes we document share one anatomy: guarantees imported into settings that quietly void their preconditions. 1 Introduction Additive models of the form f(x)=β0+∑i=1dfi(xi)f(x)= _0+ _i=1^df_i(x_i) occupy a privileged point on the accuracy–interpretability spectrum: each coordinate’s contribution is a one-dimensional function that can be plotted, audited, and constrained, while the component functions themselves remain arbitrarily nonlinear [1]. Neural Additive Models [2] realize this class with one subnetwork per feature, and their extensions with selected pairwise terms [4, 5] recover much of the accuracy lost to the additivity constraint. In all of these approaches the constraint is structural: the architecture makes interactions impossible by construction. Structural enforcement is clean but rigid: it commits to the additive hypothesis before training, discards any pretrained dense network, and offers no middle ground between “fully additive” and “fully connected.” Yet it is far from the only road. One can penalize a dense network’s interactions behaviorally and crystallize the structure afterwards; decay its cross-feature weights during training; run the classical backfitting algorithm with neural smoothers; read an additive predictor directly out of a dense network’s mean-conditioned sections; or leave neural networks aside and bag depth-one boosted trees. All of these produce (exactly or asymptotically) the same functional class. The question this paper asks is the one a practitioner actually faces: among the many routes to an additive model, which ones work, when, and do they arrive at the same functions? Answering that question requires an instrument before it requires a method: a way to measure how interacting a network is. The natural functional measure faces an obstacle that is easy to underestimate. The mixed partial derivative is identically zero almost everywhere for piecewise-linear networks—their interactions live in the gradient’s discontinuities, invisible to any pointwise test—so a useful measure must be nonlocal. We use the elementary second-order finite difference on a quadrilateral of inputs, which detects interactions across activation boundaries and costs four forward passes; the same statistic serves throughout as diagnostic, as training penalty, and as the criterion by which one-dimensional response plots earn the right to be called shape functions (Section 3). Our contributions, in the order the paper develops them: (i) the quadrilateral loss, a new differentiable penalty, with a proof that it characterizes additivity on product domains (Proposition 1), survives piecewise linearity where mixed-derivative penalties fail, and—via a two-line Möbius argument—equals, in expectation, the total interaction mass of the interventional Shapley-GAM containing each coordinate, so that driving it down also certifies the faithfulness of ordinary SHAP explanations (Section 3); (i) the λ-sweep, showing a strongly asymmetric accuracy–interaction trade-off in which most interactions are removable almost for free (Figure 2), together with a regularization sweet spot on small data—seed-universal on Wine—where moderate λ improves accuracy and additivity simultaneously, a dimensionless-λ argument and balance heuristic λbal _bal for choosing it, and a tuning experiment showing why λ should not be selected on accuracy alone (Sections 5.1 and 7); (i) surrender curves—per-feature interaction tracked online during training—whose five-seed, four-dataset replication shows that pre-regularization magnitude is only a weak, unstable predictor of which interactions a regularized model retains (Section 5.2); (iv) the observation that one-dimensional response sections only become well-defined shape functions once measured interaction is driven down (Figure 6); (v) crystallization, an optional post-processing that converts the penalized model into an exactly additive one at no accuracy cost relative to a from-scratch NAM, and whose comparison against a weight-space alternative of our design (winner-take-all decay) shows that constraining behavior before structure dominates (Section 6.1); (vi) calibration against classical backfitting with neural smoothers, which quietly matches modern per-feature architectures (Section 6.2); (vii) a shared-parameter section model on which sequential backfitting diverges for identifiable structural reasons but joint training of the same object—the shared-section model—is stable and yields our best exactly additive neural model (Section 6.3); and (viii) a six-dataset benchmark including TEAM and the soft penalized model itself, showing that route rankings are regime-dependent and that the near-additive model, at the heuristic default λ=1λ=1, leads the additive family on three datasets—plus three open clinical cohorts on which every additive route beats the dense network, the regime where additivity pays outright (Sections 7–7.1). Figure 1 organizes the routes by where the constraint acts—structure, behavior, fitting algorithm, or read-out—and marks the paper’s contributions within that ontology. Additive predictor f(x)=β0+∑ifi(xi)f(x)= _0+ _if_i(x_i)where is the constraint imposed?StructureBehaviorFitting algorithmRead-outNAMmasksTEAMstumpsQuad loss⋆→ crystallizeWTA decay⋆→ pruneBackfittingpartial residualsShared sections⋆mean-conditionedsoft/hard pair, shared target⋆introduced in this paper bold border: ours TEAM: Di Cecco [19] Figure 1: Ontology of the routes compared, organized by where the additivity constraint acts: on the structure (the architecture cannot express interactions), on the behavior (a dense model is penalized, then optionally crystallized), on the fitting algorithm (components are estimated additively), or on the read-out (an unconstrained network is queried through its one-feature sections). The quadrilateral loss additionally serves as the measuring instrument for all four families, and the two routes we introduce are the soft and hard treatments of one constraint, sharing a target but not a parametrization (dashed; Section 6.3). We emphasize at the outset what this paper does not claim. All experiments use a single tabular dataset and small multilayer perceptrons; we make no claim of generality across domains or scales, and we position this work as a carefully analyzed case study of a simple idea rather than a benchmark sweep. Related penalties on Hessian off-diagonals and functional decompositions exist in the literature (Section 2); training-time penalization of explanations also predates us [28, 30]; the quadrilateral loss itself—a finite-difference penalty targeting global additivity with an exact characterization—is, to our knowledge, new, as are the λ-interpolation analysis and the surrender-curve decomposition built on it. 2 The Additive Modeling Landscape An additive model predicts f(x)=β0+∑ifi(xi)f(x)= _0+ _if_i(x_i) with arbitrary univariate components: each fif_i is a plottable, auditable object, and the price is the exclusion of interactions. The classical instantiation fits the components with smoothers via backfitting [1, 16]. The tree lineage replaces smoothers with depth-one ensembles—boosted stumps are additive by construction [17]—running through the intelligible models of Lou et al. [3], GA2M’s screened pairwise terms [4], and the Explainable Boosting Machine [22]; TEAM [19, 21], described in Section 4, is a compact member, and on tabular data this lineage remains hard to beat [18]. Its canonical deployment is clinical risk: Caruana et al. [36] famously used an intelligible additive model on pneumonia data to surface the rule asthma implies lower mortality risk—a triage artifact that a dense model would have silently exploited—and that episode is the template for what our surrender curves are built to do online: expose which interactions a model is retaining, so that a domain expert can decide whether they are signal or artifact. The neural lineage gives each feature its own network [2], with successors sharing bases across features [6], modeling interactions via low-rank polynomial decompositions [23], building components from differentiable trees [5], adding structured pairwise terms [24], or treating the sub-networks in a Bayesian fashion, with credible intervals and interaction ranking [37]; our shared-section model (Section 6.3) is the limiting case of basis sharing, with the whole network as basis and no per-feature parameters. Adjacent threads include Kolmogorov–Arnold Networks [25], whose one-layer form with identity outer functions is exactly a GAM, and concurvity regularization [26], which penalizes correlation between fitted components to make decompositions identifiable under dependent features—complementary to our penalty on interactions within a dense model. On the measurement side, our statistic descends from the functional-ANOVA tradition [12, 7], the partial-dependence-based H-statistic of Friedman and Popescu [38], interaction detection by restricted retraining [39], and post-hoc interaction attribution [8, 9, 10]; Tsang et al. [11] constrain interactions through architecture discovery, and training-time penalization of explanations predates us in the form of gradient penalties [28], attribution priors [29], and contextual decomposition penalization of practitioner-specified attributions [30]. Pointwise second-derivative penalties are ill-posed for piecewise-linear networks, which motivates our finite-difference formulation; the quadrilateral loss itself—a nonlocal penalty targeting global additivity with an exact characterization—is, to our knowledge, new, as are the λ-interpolation analysis and the surrender-curve decomposition built on it. Our crystallization step borrows saliency from magnitude–gradient pruning [13] and group sparsity [14]. 3 The Quadrilateral Loss We introduce the quadrilateral loss, a differentiable training penalty whose population version exactly characterizes additivity and which remains informative for piecewise-linear networks. Adjacent constructions exist—mixed differences as post-hoc diagnostics in the functional-ANOVA tradition [12], and training-time penalization of local, practitioner-specified attributions [28, 30, 29]—but to our knowledge the statistic below has not previously been used as a training objective, nor tied to a global additivity guarantee. Let fθ:ℝd→ℝf_θ:R^d be a network and let x,x~x, x be two inputs. For a coordinate i, define the single-coordinate swaps x(i)=(x~i,x−i)x^(i)=( x_i,x_-i) and x~(i)=(xi,x~−i) x^(i)=(x_i, x_-i), where x−ix_-i denotes all coordinates except i. The quadrilateral residual is Qi(x,x~)=fθ(x)−fθ(x(i))−fθ(x~(i))+fθ(x~).Q_i(x, x)\;=\;f_θ(x)\;-\;f_θ(x^(i))\;-\;f_θ( x^(i))\;+\;f_θ( x). (1) This is the second-order mixed difference of fθf_θ along the coordinate-i direction versus all remaining directions, evaluated at the four corners of an axis-aligned quadrilateral. If fθf_θ is additive, every term fj(⋅)f_j(·) for j≠ij≠ i appears once with each sign and cancels, as does fi(⋅)f_i(·), so Qi≡0Q_i≡ 0. Conversely: Proposition 1. Let =1×⋯×dX=X_1×·s×X_d be a product domain and f:→ℝf:X . Then Qi(x,x~)=0Q_i(x, x)=0 for all x,x~∈x, x and all i if and only if f(x)=∑ifi(xi)f(x)= _if_i(x_i) for some functions fi:i→ℝf_i:X_i . Proof. Sufficiency is the cancellation above. For necessity, fix a reference point c∈c and define fi(xi)=f(xi,c−i)−d−1df(c)f_i(x_i)=f(x_i,c_-i)- d-1df(c). Setting x~=c x=c in Qi(x,c)=0Q_i(x,c)=0 gives f(x)=f(x(i)∣x~=c)+f(xi,c−i)−f(c)f(x)=f(x^(i)\!\! _ x=c)+f(x_i,c_-i)-f(c), i.e., replacing coordinate i of any point by cic_i changes f by an amount depending on xix_i alone. Applying this coordinate by coordinate telescopes f(x)f(x) down to f(c)f(c) plus a sum of single-coordinate increments, which is the claimed additive form. ∎ Two properties matter in practice. First, (1) is kink-robust: it integrates interaction across activation boundaries rather than sampling the gradient at a point, so it detects interactions in ReLU networks for which ∂2f/∂xi∂xj=0∂^2f/∂ x_i∂ x_j=0 almost everywhere. Second, it is cheap: four forward passes, no double backpropagation, and the four corner evaluations batch trivially. Training objective. We penalize the empirical second moment of the residual over data pairs and coordinates, ℒ(θ)=[(fθ(x)−y)2]⏟task+λ1|S|∑i∈Sx,x~[Qi(x,x~)2],L(θ)\;=\; E [(f_θ(x)-y)^2 ]_task\;+\;λ\, 1|S| _i∈ S\;E_x, x [Q_i(x, x)^2 ], (2) where x~ x is obtained by permuting the minibatch and S is a random subset of coordinates resampled each step (|S|=4|S|=4 throughout). Sampling x~ x from the data distribution rather than a product measure means the penalty is enforced on-manifold: the network may retain interactions in regions unsupported by data. For interpretability applications this is usually the desired semantics; for worst-case guarantees it is not, and we return to this limitation in Section 10. Connection to the Shapley-GAM. The statistic has an exact game-theoretic meaning. Bordt and von Luxburg [33] show that every subset-compliant value function induces a unique functional decomposition f=∑S⊆[d]fS(xS)f= _S [d]f_S(x_S) via the Möbius transform—the Shapley-GAM—and that for the interventional SHAP value function v(x,S)=z[f(z)|do(xS)]v(x,S)=E_z[f(z)\,|\,do(x_S)] the pairwise component is precisely a quadrilateral: fij=[f|do(xi,xj)]−[f|do(xi)]−[f|do(xj)]+[f]f_ij=E[f|do(x_i,x_j)]-E[f|do(x_i)]-E[f|do(x_j)]+E[f]. Our statistic is the two-point Monte Carlo version of this construction, and the correspondence is exact: writing v(x,S)=∑L⊆SfLv(x,S)= _L Sf_L and taking the expectation of (1) over the partner point x~∼ x with x fixed gives x~[Qi(x,x~)]=f(x)−v(x,[d]∖i)−v(x,i)+v(x,∅)=∑L∋i,|L|≥2fL(xL)E_ x[Q_i(x, x)]=f(x)-v(x,[d]\! \!i)-v(x,\i\)+v(x, )= _L i,\,|L|≥ 2f_L(x_L), the sum of all interaction components of the interventional Shapley-GAM that contain coordinate i. The quadrilateral loss therefore penalizes the ℓ2 ^2 mass of the Shapley-GAM’s interaction components, coordinate by coordinate, and driving it to zero is equivalent to the Shapley-GAM being of order one. Two consequences follow (a third, structural one—that additive functions are the fixed points of a mean-anchored section read-out, making our penalty the soft counterpart of a hard projection—is developed with the section model in Section 6.3). First, by the recovery theorem of Bordt and von Luxburg [33], once a model is (near-)additive its ordinary interventional SHAP values coincide with its shape functions even under arbitrarily dependent features—so λ is also a dial for the faithfulness of standard SHAP explanations, quantified by the same statistic. Second, the surrender curves of Section 5.2 become a tractable window onto an otherwise exponential object: the per-coordinate interaction mass of the Shapley-GAM, which Bordt and von Luxburg [33] estimate by enumerating all 2d2^d subsets, is probed here at O(d)O(d) forward passes per checkpoint, online, during training. Higher-order probes and uncertainty. The construction generalizes beyond single coordinates: swapping a subset S of coordinates between x and x~ x and taking the alternating sum over the 2|S|2^|S| corners yields a one-sample Monte Carlo estimate of ∑L⊇SfL(x) _L Sf_L(x), the joint interaction mass of all Shapley-GAM components containing S—for pairs, an 88-point statistic that screens candidate interaction terms for a GA2M at a handful of forward passes per pair, in place of exponential subset enumeration. Bayesian treatments of NAMs also rank candidate interaction pairs [37]; our methods acquire uncertainty by composition rather than redesign—the penalty is a loss term through which any Laplace or ensemble treatment of the dense network trains unchanged, the section read-out inherits credible bands by evaluating the posterior on the cross, and TEAM’s bagged replicates already form a bootstrap ensemble whose between-bag spread is a confidence band at no extra cost. We use only the single-coordinate probe in this paper. The knob, honestly stated. As λ ranges over [0,∞)[0,∞), (2) interpolates between an unconstrained dense network and an approximately additive one. We stress that intermediate λ does not correspond to a structured model class such as a GA2M of fractional order: the intermediate models are dense networks with uniformly small but nonzero interactions of all orders. The interpolation is a trade-off curve in behavior space, not a nested family of hypotheses. Where a structured intermediate class is required, the penalty extends naturally by restricting the sum in (2) to coordinate swaps across prescribed feature blocks, leaving intra-block interactions free; we leave a systematic study of this block variant to future work. 4 Experimental Setup All experiments use the California Housing dataset [40] (20,433 block groups after dropping rows with missing bedroom counts; eight standardized features; target median house value in units of $100,000), with a fixed random 80/20 train/test split. Unless stated otherwise the architecture is a two-hidden-layer MLP (64+6464+64 units, GELU) trained with Adam [15] at learning rate 10−310^-3, batch size 512, 120 epochs, with λ annealed linearly from zero after a 15-epoch warmup; the λ-sweep of Section 5.1 uses a single-hidden-layer network (H=64H=64, 80 epochs, 3 seeds per point) to keep the sweep cheap, and all comparisons in that section are against the same architecture at λ=0λ=0. We measure interactions with the test-set statistic Q2¯=1d∑i[Qi2] Q^2= 1d _iE[Q_i^2] estimated by batch permutation (3 replicates); for reference, the test target variance is 1.311.31. Routes compared. For self-containedness we describe every method at the level needed to reimplement it. NAM from scratch. The same MLP with block-diagonal binary masks fixed at initialization: first-layer units are pre-assigned to features in balanced groups (unit j may connect only to feature jmoddj d), second-layer units connect only to first-layer units of their own group, and the linear output layer sums all groups. The masks are re-applied after every optimizer step, so each hidden path processes exactly one feature and the model is additive by construction [2]. We use standard GELU units rather than the exp-centered (ExU) activations of the original NAM; we return to this choice, which is not innocent, in Section 5.3. WTA decay (ours). Dense training in which, after a 15-epoch warmup, every non-maximal input weight of each first-layer unit is multiplied by (1−ε)(1- ) at each step (ε annealed to 0.020.02), so each unit’s receptive field contracts toward a single self-selected feature; a grouped analogue acts at the second layer with a 2.5×2.5× stronger rate. Training is followed by hard pruning to the winner topology (with the capacity-balancing repair of Section 6.1) and fine-tuning under the masks. Quadrilateral penalty (ours). Dense training with (2), λ annealed to 20, followed by the same pruning-and-fine-tuning protocol; first-layer feature assignments are extracted from an exponential moving average of the weight–gradient saliency |w⊙∇wℒ||w _wL| [13]. Backfitting. Backfitting [1] uses single-feature networks (1→64→11→64→1, GELU) as smoothers, fitting them cyclically on partial residuals y−α−∑j≠ifj(xj)y-α- _j≠ if_j(x_j), eight inner epochs per component per sweep, warm-started across 40 sweeps, with each component mean-centered after its update and the intercept α absorbing the shifts. Shared-section model (ours). A single ordinary dense network gθg_θ read out through Eq. 3: the additive predictor is the sum of gθg_θ’s mean-conditioned one-feature sections, trained end-to-end by backpropagating the task loss through all d+1d+1 masked forward passes jointly. TEAM. A Tree Ensemble Additive Model [19, 21]: a bagged ensemble (B bootstrap replicates) of gradient-boosted regression trees of depth one. Additivity is exact and comes from the weak learner, not from any constraint machinery: a depth-one tree (stump) is a step function of a single feature, so a boosted sum of stumps is a sum of per-feature step functions, ∑mνhm(xim)=∑ifi(xi) _mν\,h_m(x_i_m)= _if_i(x_i) with piecewise-constant fif_i; bagging averages additive models and is therefore additive. Boosting chooses at every stage which feature to split, so capacity is allocated across features adaptively by the functional gradient rather than fixed per branch; bagging reduces the variance of that greedy allocation. This is also where TEAM departs from the EBM lineage of Lou et al. [3], Nori et al. [22]: EBMs fit shape functions by cyclic gradient boosting, visiting every feature in round-robin with a small learning rate to wash out order effects, which costs a full pass over the feature set per boosting round; TEAM lets each stage select its feature greedily by standard split search and delegates the de-biasing of that greed to bagging, whose replicates are embarrassingly parallel. The result is a simpler and faster training loop at the price of coarser per-feature budget control. Shape functions are read out by the same mean-conditioning used in Eq. 3. We report TEAM at B=5B=5 bags of 1,000 stumps (learning rate 0.10.1). Figure 2: Accuracy–interaction trade-off under the quadrilateral penalty (single hidden layer, H=64H=64; mean ± s.d. over 3 seeds). Left: test MSE versus λ; the dashed line is the identical architecture at λ=0λ=0 and the dotted line a NAM of matched width trained from scratch. Right: residual interaction Q2¯ Q^2 versus λ (log–log). Interaction falls by nearly four orders of magnitude, from 0.860.86 to 1.1×10−41.1× 10^-4, while MSE rises from 0.3700.370 to 0.4700.470; at λ=0.1λ=0.1, 91%91\% of the interaction mass is removed for +0.008+0.008 MSE. The right curve saturates near 10−410^-4, the noise floor of the stochastic estimator rather than exact additivity. Reproducibility. Code reproducing every experiment, table, and figure—the loss, all training routes, the balance heuristic, the Hyperband protocol, and the surrender-curve instrumentation—is available at github.com/AntonioDiCecco/quadrilateral-loss. 5 What the Loss Reveals: a Dial and an Observable The single coefficient λ turns additivity into a training-time dial; the coordinate-wise decomposition of the statistic turns interaction into an online observable; and driving the statistic down determines when one-dimensional plots earn an interpretation. This section develops the three in turn. 5.1 The Trade-off is Strongly Asymmetric Figure 2 reports the sweep λ∈0,0.1,0.3,1,3,10,30,100,300λ∈\0,0.1,0.3,1,3,10,30,100,300\. Three observations. Most interactions are cheap to remove. At λ=0.1λ=0.1 the interaction statistic drops from 0.8550.855 to 0.0730.073—a 91%91\% reduction—while test MSE moves from 0.37030.3703 to 0.37800.3780. By λ=10λ=10 interactions are down three orders of magnitude (1.0×10−31.0× 10^-3) at MSE 0.4230.423, still short of the full additive cost. The dense optimum evidently encodes its function with substantial interaction redundancy: alternative, nearly-additive parameterizations of almost the same function exist nearby in function space, and a weak penalty suffices to select them. The economically meaningful interactions—those the network genuinely needs—are concentrated in the last order of magnitude of suppression, which is where the accuracy price is paid. The useful dynamic range of λ is narrow. Both curves are steep only for λ∈[0.1,10]λ∈[0.1,10]; below, the penalty is inert, and above, the interaction estimate saturates at ∼10−4 10^-4. This saturation is a property of the estimator, not the function: with finite batches and |S|=4|S|=4 sampled coordinates per step, the stochastic penalty cannot certify exact additivity, only drive its estimate to the noise floor. Exactness requires the crystallization step of Section 6.1. Near-additive dense networks can match structural NAMs. At λ=1λ=1 the single-layer penalized network attains MSE 0.3990.399 with residual interaction 1.2×10−21.2× 10^-2—statistically indistinguishable from the from-scratch NAM (≈0.39≈0.39, Table 2). The “almost additive” regime is thus a legitimate operating point rather than a no-man’s-land, and there is a principled reason to prefer it that goes beyond accuracy. Exact additivity pushes complexity into the components: the best additive projection of a genuinely interacting target can demand shape functions of extreme roughness—the observation that led the NAM authors to engineer ExU activations for “jumpy” univariate fits [2]—and a shape function jagged enough to absorb what interactions used to express is formally one-dimensional but no longer humanly readable, which defeats the purpose of the constraint. This is the expressivity paradox of additive models: perfect additivity and interpretable components can be in tension. The paradox has a classical pedigree. The Kolmogorov–Arnold superposition theorem [34] shows that every continuous multivariate function is a finite superposition of univariate functions—so one-dimensionality of the building blocks, by itself, carries no interpretability guarantee whatsoever—and the univariate functions that make the theorem true are necessarily wild: they cannot in general be taken smooth [35], and the standard constructions produce continuous but fractal-like inner functions. The theorem is thus the extreme statement of the trade our sweep exhibits in miniature: multivariate complexity does not disappear when arity is constrained, it is displaced into the regularity of the components. A GAM is far weaker than a full superposition—one layer, no outer function—so it cannot absorb everything; but its best additive projection of an interacting target already sits on the same slope, and the roughness that ExU activations were engineered to reach is that displaced complexity surfacing. Soft enforcement dissolves the dilemma from the side: at moderate λ the model retains a small interaction budget precisely where forcing additivity would deform the components most, and the resulting sections stay smooth while the quadrilateral statistic certifies, quantitatively, how far from additive the model still is. Near-additivity with a measured residual can be more interpretable than exact additivity with wild components. A regularization sweet spot on small data. The California Housing sweep is monotone: every unit of additivity costs accuracy. But Section 7 shows that on smaller datasets the additive class outperforms the dense network—which predicts that there the trade-off curve should dip below the unconstrained baseline before rising. It does, with a strength that itself depends on the dataset (Figure 3; ten seeds, paired per seed since penalty runs share split and initialization with their baseline). On Wine the effect is unambiguous: test MSE falls from 0.4800.480 at λ=0λ=0 to 0.4540.454 at λ=3λ=3, the paired improvement holds on ten of ten seeds for every λ∈[0.3,3]λ∈[0.3,3] (−0.023±0.009-0.023± 0.009 at λ=1λ=1), the interaction statistic drops an order of magnitude, and even λ=30λ=30 remains at or below the baseline on seven of ten seeds. On Concrete a ten-seed replication demotes our initial three-seed reading: the dip survives only at the smallest penalty (λ=0.3λ=0.3, −0.8±1.4-0.8± 1.4 paired, eight of ten seeds) and λ≥3λ≥ 3 already costs accuracy—the earlier, larger dip was partly seed luck. The honest summary is therefore conditional rather than universal: in the small-n regime there can exist λ at which the model is simultaneously more accurate and more additive than the unconstrained one, decisively so on Wine, marginally on Concrete, and where it happens the quadrilateral penalty is not a constraint purchased with accuracy but a regularizer paid for by nothing—the interactions were the overfitting. Figure 3: Win–win regime of the quadrilateral penalty on the two datasets where the additive class beats the dense network (2-hidden-layer net, 10 seeds, bands ± 1 s.e.m.). On Wine the dip below the unconstrained baseline (dashed) is seed-universal for λ∈[0.3,3]λ∈[0.3,3]; on Concrete it is marginal and confined to the smallest penalty. Where the dip exists, the penalty regularizes rather than trades. 5.2 Surrender Curves: Interaction as a Training-Time Observable Because the quadrilateral statistic decomposes by coordinate, the penalty doubles as a free online diagnostic: tracking [Qi2]E[Q_i^2] per feature during training yields surrender curves showing which features’ interactions resist the regularizer and which capitulate (CalHousing panel of Figure 5). This is information no structurally constrained model can emit—a NAM never had interactions to surrender—and it turns post-hoc interaction detection [7, 8] into a live training signal. Feature [Qi2]E[Q_i^2] at warmup end [Qi2]E[Q_i^2] at convergence Reduction Rank at convergence MedInc 0.0900.090 8.0×10−48.0× 10^-4 113×113× 2 HouseAge 0.0730.073 4.8×10−44.8× 10^-4 153×153× 7 AveRooms 0.5920.592 8.9×10−48.9× 10^-4 667×667× 1∗ AveBedrms 0.6620.662 6.6×10−46.6× 10^-4 1006×1006× 4 Pop 0.0220.022 3.2×10−43.2× 10^-4 69×69× 8 AveOccup 0.0110.011 1.9×10−31.9× 10^-3 6×6× 1 Lat 0.8990.899 6.7×10−46.7× 10^-4 1342×1342× 3 Lon 0.8700.870 5.7×10−45.7× 10^-4 1518×1518× 5 Total 3.223.22 6.3×10−36.3× 10^-3 511×511× Table 1: Per-feature interaction [Qi2]E[Q_i^2] on the test set at the end of the unconstrained warmup (epoch 15, the peak of freely grown interaction) and at convergence (mean of the last five checkpoints), with the reduction factor achieved by the penalty. The largest interactions (Lat, Lon, the ratio features) are suppressed by three orders of magnitude, while AveOccup—smallest at warmup—is reduced only 6×6× and ends as the largest residual (∗AveRooms ranks first only through its transient spikes; on the smoothed tail AveOccup leads). On this run (seed 0) warmup magnitude and final residual are uncorrelated (ρs=−0.05 _s=-0.05); Figure 4 shows the seed-replicated picture. The curves corrected our own prior. We expected latitude and longitude—the carriers of the dataset’s dominant interaction—to be the last holdouts; instead they surrender to ∼7×10−4 7× 10^-4, while the stubborn features are average occupancy (∼2×10−3 2× 10^-3, drifting upward late) and median income. Two readings, not mutually exclusive. First, [Qi2]E[Q_i^2] is a variance-weighted measure evaluated on random data pairs: the heavy-tailed ratio features (rooms, bedrooms, occupancy) produce extreme swap pairs whose quadrilateral residuals dominate the average, which also explains their outsized pre-penalty peaks. Second, the geographic interaction, while large, is apparently easy to trade away—the network finds nearby near-additive substitutes for it—whereas the occupancy–income interactions are entangled with the task gradient in a way that keeps regenerating them, visible as the late upward drift and the transient re-emergence spikes around epochs 59 and 103. Table 1 quantifies the effect: reduction factors range from 1518×1518× (longitude) down to 6×6× (occupancy). We caution against reading the strong positive correlation between warmup magnitude and reduction factor as a finding: since the factor is the ratio peak/final, it correlates with the peak mechanically whenever the final residual is independent of the peak—and independence is precisely what we observe (below). Either way, the operational point stands: which interactions are expensive to remove is an empirical property of the optimization, not something readable off post-hoc importance scores, and the surrender curves measure it directly during training. A dynamical-systems analogy. We note, with the caution due to any cross-field metaphor, that the structure of the finding echoes KAM theory [27]. There, an integrable Hamiltonian—the mechanical analogue of a separable, interaction-free system—is perturbed, and the question is which invariant tori survive; the celebrated answer is that survival is governed not by any notion of a torus’s size but by an arithmetic property of its frequencies (a Diophantine non-resonance condition), with resonant tori destroyed at arbitrarily small perturbation. Our setting runs the perturbation in the opposite direction—the penalty pushes an interacting system toward separability—but exhibits the same signature: which interactions survive the perturbation is decided not by their magnitude but by a structural property of how they are entangled with the task gradient, a property the magnitude ranking barely predicts. We do not claim a formal correspondence; we record the analogy because it correctly reframes the question the surrender curves answer—from “how big is this interaction?” to “how resonant is it with what the loss needs?”—and because it suggests that a quantitative stubbornness criterion, playing the role of the Diophantine condition, may be identifiable. Replication across datasets and seeds: magnitude is a weak, unstable predictor of stubbornness. Repeating the experiment on Wine, Abalone, and Boston (same protocol, per-dataset batch sizes as in Section 7) and—crucially—across five seeds per dataset revised our own first reading, twice. On the initial seed, the rank correlation between a feature’s warmup interaction magnitude and its final residual appeared to vanish on all four datasets (|ρs|≤0.12| _s|≤ 0.12), suggesting clean decoupling. Five seeds show that this was a fortunate draw: per-seed correlations range from −0.26-0.26 to +0.81+0.81, with per-dataset means of +0.37±0.31+0.37± 0.31 (CalHousing), +0.20±0.26+0.20± 0.26 (Wine), +0.39±0.27+0.39± 0.27 (Abalone), and +0.05±0.14+0.05± 0.14 (Boston); pooling all twenty runs, the grand mean is ρs=+0.25±0.29 _s=+0.25± 0.29, significantly positive (Wilcoxon p=0.002p=0.002) but weak (Figure 4). The corrected claim is therefore not decorrelation but unreliability: initial interaction magnitude carries some signal about final stubbornness, yet so little—and so seed-dependent—that on any single training run the magnitude ranking routinely inverts (seed-level ρs _s crosses zero on three of four datasets). The identities of the stubborn features remain dataset-specific and plausible (shucked and whole weight on Abalone, whose near-collinear relation to total weight is genuinely interactive; rooms and %lower-status on Boston), consistent with stubbornness being a real feature-level property whose relation to raw magnitude is loose. The practical consequence for anyone using post-hoc interaction scores to decide which pairwise terms a GA2M should keep survives in weakened but usable form: the magnitude ranking explains at best a small fraction of the retention ranking’s variance, and only an online measurement like the surrender curves, run on the actual training in question, reveals which interactions that run will keep. Figure 4: Seed replication of the magnitude–stubbornness relation. Each dot is one training run’s Spearman correlation between per-feature warmup interaction magnitude and final residual; red diamonds: per-dataset mean ± s.d. over five seeds. The relation is weakly positive on average (+0.25±0.29+0.25± 0.29 pooled, Wilcoxon p=0.002p=0.002) and crosses zero within seeds on three of four datasets: a single run’s ranking is unreliable, which is itself the operational finding. Figure 5: Surrender curves on four datasets (three most stubborn features highlighted per panel; shaded band = penalty annealing). Panel titles report ρs _s between warmup magnitude and reduction factor; panel curves are shown for seed 0; the seed-replicated statistics are in Figure 4. 5.3 When Do Sections Become Shape Functions? Figure 6 shows, for seven values of λ, the one-dimensional sections xi↦fθ(xi,x¯−i)x_i f_θ(x_i, x_-i) with the remaining coordinates at their means. These are the objects practitioners routinely plot as “the effect of feature i.” The figure makes the epistemic status of such plots visible. For MedInc and HouseAge the sections are essentially invariant across the entire λ range, including λ=0λ=0: their effects are additive already in the unconstrained model, and the penalty has nothing to remove. For latitude and longitude—the interaction carriers—the λ=0λ=0 section is the most irregular curve in each panel and its shape drifts as λ grows, converging only around λ≳10λ 10 to a stable profile shared across penalization strengths. The average-ratio features (AveRooms, AveBedrms, AveOccup) show a related effect in their sparse tails, where apparent structure at λ=0λ=0 flattens under the penalty; caution is warranted there regardless, as the standardized tails beyond ±2σ±2σ contain few data points. The interpretation is not that the penalty “smooths” the sections. It is that at λ=0λ=0 a section is a one-dimensional slice through a strongly interacting function, and its shape depends on the arbitrary choice of conditioning point x¯−i x_-i—slicing at a different context would give a different curve, which is precisely what interaction means. Only as Q2¯→0 Q^2→ 0 does the section become invariant to the conditioning point and thus interpretable as a shape function fif_i. The λ-convergence visible in Figure 6 is therefore a direct visualization of sections acquiring well-definedness. It is the training-time dual of a phenomenon Bordt and von Luxburg [33] exhibit post-hoc: their Shapley-value partial-dependence scatter collapses onto the component functions as the explanation order n grows toward d, ours as the interaction budget shrinks toward zero—two paths to the same identified object. A sharper version of this diagnostic—plotting, at each λ, a bundle of sections conditioned at multiple context points and watching the bundle collapse—would quantify section validity per feature; we propose it as a cheap sanity check for any practitioner reading partial-dependence-style plots off a dense model. 6 Routes to Exact Additivity Near-additivity suffices for many uses, but exactness is sometimes the point. This section covers how to get it: a brief account of crystallization as post-processing for the penalized model, then the two full modeling routes—classical backfitting with neural smoothers, and reading an additive predictor out of a single shared network, which will turn out to be not a separate design but the hard-projection counterpart of the quadrilateral penalty itself. 6.1 From Near to Exact: Crystallization as Optional Post-Processing The penalized model is additive in behavior, not in structure: the stochastic estimator bottoms out near 10−410^-4, so exactness, where required, needs one further step. Crystallization converts the penalized dense network into a masked NAM topology—first-layer units assigned to features by the saliency |w⊙∇wℒ||w _wL|, second-layer units grouped under per-feature capacity caps, cross-group weights zeroed, then fine-tuned under the masks—after which the additive decomposition holds to machine precision. On California Housing this recovers test MSE 0.395±0.0060.395± 0.006 over three seeds, matching a NAM trained from scratch (0.388±0.0030.388± 0.003, Table 2): pruning is violent in parameter space but benign in function space, because the function being re-encoded is already additive. The contrast that makes this a finding rather than a recipe is our weight-space alternative, WTA decay, which constrains parameters during training (non-winning weights decay even while carrying signal) and plateaus at 0.4390.439 (seed-matched single run) after the identical pruning protocol: when the target is a functional property, constraining behavior first and structure second dominates constraining the weights directly. Two silent failure modes surfaced here—Adam’s stale moments moving pruned weights off their masks (gap inflating from 10−610^-6 to 3×10−23× 10^-2; fix: re-project after every step), and unconstrained winner-take-all grouping collapsing onto a single feature (fix: capacity caps with a minimum quota)—and are dissected in Section 9. We emphasize the demotion implied by this subsection’s title: crystallization is an optional exactness post-processing, not a modeling route, and we exclude it—together with its weight-space ablation, WTA decay, which exists only as the comparator for the behavior-first claim—from the cross-dataset benchmark accordingly; on small datasets the multi-stage pipeline is fragile, and part of the deficit we initially attributed to it traced, on re-examination, to a batch-size protocol artifact of the kind Section 7 documents for Boston. 6.2 Classical Backfitting with Neural Smoothers A classical baseline predating neural additive modeling by decades calibrates what joint fitting buys. Backfitting. The naive alternative to joint fitting—train d single-feature networks independently on y and stack them linearly—fails badly (0.4990.499, barely above linear regression) despite living in the same additive class, because each branch learns the marginal [y|xi]E[y|x_i], which under correlated features double-counts effects that the additive class needs as partials; the fitted stacking weights betray it with large negative coefficients subtracting the overlap. The classical repair is the backfitting algorithm of Hastie and Tibshirani [1]: cyclically refit each component on the partial residual y−α−∑j≠ifj(xj)y-α- _j≠ if_j(x_j), with mean-centering of each component for identifiability. Instantiated with the same single-feature networks as smoothers (8 inner epochs per component per sweep, warm-started across sweeps), backfitting descends monotonically from the stacking level, crosses the from-scratch NAM at roughly sweep 20, and plateaus at 0.380±0.0060.380± 0.006 over three seeds (Table 2). The trajectory is exactly the theory: partial residuals progressively reassign shared signal among correlated features, converging to the joint additive fit—backfitting is coordinate descent in the space of additive functions, and each sweep provably decreases the training objective because updating fif_i leaves every other component exactly fixed. We do not read the final 0.3800.380 vs. 0.3880.388 margin as a systematic advantage over the NAM: the backfitting ensemble has roughly 8×8× the first-layer parameters per feature, and the difference is within seed and tuning noise. The robust finding is the shared plateau: masked NAM, penalty-then-crystallize, and backfitting—three entirely different procedures—all land at ≈0.39≈0.39, identifying the additive limit of the dataset itself. 6.3 A Single Network as Its Own Additive Model The routes above either dedicate parameters per feature (NAM, backfitting) or train a dense network and discard its cross-feature weights (crystallization). A third possibility uses one dense network gθ:ℝd→ℝg_θ:R^d of ordinary architecture and defines the additive predictor through its mean-conditioned sections: F(x)=gθ(x¯)+∑i=1d[gθ(xi,x¯−i)−gθ(x¯)],F(x)\;=\;g_θ( x)\;+\; _i=1^d [\,g_θ(x_i, x_-i)-g_θ( x)\, ], (3) where x¯ x is the training mean. Each bracketed term depends on a single coordinate, so F is exactly additive by construction, for any θ; the question is only how to fit θ. Sequential backfitting on shared weights diverges. The direct transplant of backfitting—cyclically fit section i on partial residuals by training gθg_θ on mean-masked inputs (xi,x¯−i)(x_i, x_-i)—reaches 0.540.54 within ∼15 15 sweeps and then diverges catastrophically (test MSE >103>10^3 by sweep 40). The failure is structural rather than a matter of step sizes. Classical backfitting converges because its components are parameter-disjoint: refitting fif_i leaves every fjf_j, j≠ij≠ i, untouched, which is what makes the procedure a genuine coordinate descent with stationary partial residuals within each update. With shared weights, training section i moves all sections simultaneously; the partial residuals computed at the start of a sweep are stale before the sweep ends, the error is baked into the next component’s regression target, and the feedback loop amplifies. The convergence guarantee of Buja et al. [16] dies precisely where the coordinates cease to be orthogonal in parameter space. Penalize or project: the section model and the quadrilateral loss are two faces of one idea. The two routes we introduce are more closely related than their descriptions suggest, and the relation is an identity. For d=2d=2, the gap between the dense network and its section read-out is exactly a quadrilateral anchored at the mean: g(x)−F(x)=g(x1,x2)−g(x1,x¯2)−g(x¯1,x2)+g(x¯1,x¯2)=Q(x,x¯)g(x)-F(x)=g(x_1,x_2)-g(x_1, x_2)-g( x_1,x_2)+g( x_1, x_2)=Q(x, x); for general d it is the sum of the anchored decomposition’s terms of order two and higher. Additive functions are precisely the fixed points of the read-out (g additive ⇒F[g]=g F[g]=g, as noted above), and the quadrilateral statistic characterizes that fixed-point set. The two methods are then the classical soft and hard treatments of the same constraint: the quadrilateral loss penalizes the distance from the fixed-point set, with the partner point sampled from the data so that the penalty weights interactions by the interventional measure; the section model composes with the projection, discarding the interaction component by fiat rather than suppressing it. In the λ→∞λ→∞ limit they share the same target—minimum MSE over the additive class, whose population optimum is the additive L2L^2 projection—but they do not coincide as procedures, and the distinction is visible in our numbers. The penalty forces the network itself to become additive, a severe constraint on a dense parametrization; the projection leaves the network free and additivizes only the read-out. At finite width the two parametrizations of the additive class have different effective capacity, and the section model’s empirical edge over the crystallized penalty route (0.3760.376 vs. 0.3950.395 on California Housing) is precisely that capacity difference showing. The remaining asymmetries are consequences of soft versus hard, not separate designs: the projected model earns the exact guarantee everywhere but its gradients only ever flow on the cross of one-dimensional manifolds through the anchor, while the penalized model enforces Q=0Q=0 only on the data’s swap distribution, trains on the full manifold, and buys accuracy with its small measured budget (Table 2). Joint training of the shared-section model is stable and wins. The repair is to abandon the sequential structure altogether and minimize the joint objective [(F(x)−y)2]E[(F(x)-y)^2] directly, backpropagating through all d+1d+1 masked forward passes of (3) in every step. There are no stale residuals by construction. With the same 64+6464+64 architecture and budget as the dense baseline, this attains test MSE 0.376±0.0040.376± 0.004 over three seeds—the best exactly additive neural model in our study, ahead of both the from-scratch NAM (0.388±0.0030.388± 0.003) and separate-network backfitting (0.380±0.0060.380± 0.006), in an architecture-matched comparison with the dense reference (Table 2). Why does an exactly additive predictor built from shared weights beat one built from disjoint branches? Our working hypothesis is implicit basis sharing: the shape functions of this dataset share qualitative structure (monotone segments, saturations), and a single network can learn that structure once and reuse it across sections, whereas disjoint branches must each rediscover it from their share of the data. On this reading, (3) is the limiting case of basis sharing anticipated in Section 2: the whole network is the basis, the sections are the read-outs. The price is computational: d+1d+1 forward passes per example, linear in d, which is acceptable at d=8d=8 and untested at scale. A further property we note but do not yet exploit: the trained θ supports two read-out modes—the additive F(x)F(x) of (3) and the ordinary dense forward gθ(x)g_θ(x)—and the relationship between them (in particular, whether training through sections implicitly additivizes the dense forward pass, as the quadrilateral penalty does explicitly) is a natural follow-up experiment. 7 Across Datasets: Rankings Are Regime-Dependent The single-dataset picture invites overgeneralization in both directions. We therefore repeated the converging routes on five further tabular benchmarks—Wine Quality (red), Abalone, Boston Housing, Concrete compressive strength, and daily Bike rentals—under a uniform protocol (same architectures, Adam 10−310^-3, three seeds, 80/20 splits; batch size reduced on the smallest datasets—to 64 on Boston and Bike, 128 on Concrete—where on Boston ntrain=404n_train=404 would otherwise yield one gradient step per epoch and spuriously catastrophic neural results, a protocol artifact worth flagging because it initially masqueraded as method failure). We additionally include TEAM [19, 21], a Tree Ensemble Additive Model: a bagged ensemble of depth-one gradient-boosted trees. Because every stump splits on a single feature, each boosting stage—and hence the bagged ensemble—is exactly additive by construction, in the boosted-stumps-as-GAM lineage of Friedman [17] and the EBM family [3]; its mean-conditioned marginal-effect read-out coincides with the section read-out of Eq. 3, making TEAM the tree-based twin of our shared-network section model. We report TEAM at 5 bags of 1,000 stumps; the repository default of 100 stumps underfits noticeably (e.g. 0.480.48 on California Housing), so capacity matters for fair comparison. CalHousing WineRed Abalone Boston Concrete Bike n=20433,d=8n=20433,d=8 n=1599,d=11n=1599,d=11 n=4177,d=8n=4177,d=8 n=506,d=13n=506,d=13 n=1030,d=8n=1030,d=8 n=728,d=13n=728,d=13 Dense (non-additive) 0.282±.0020.282±.002 0.534±.0080.534±.008 4.74±1.144.74± 1.14 9.90±2.319.90± 2.31 64.2±2.764.2± 2.7 0.373±.0070.373±.007 Quad soft, λ=1λ=1 (ours, near-additive) 0.328±.0040.328±.004 0.448±.0290.448±.029 4.40±.484.40±.48 10.65±3.0810.65± 3.08 70.4±3.770.4± 3.7 0.377±.0160.377±.016 Quad soft, λ tuned (ours; median λ∗λ^*) 0.301±.0180.301±.018 (0.06) 0.444±.0340.444±.034 (2.6) 4.38±.494.38±.49 (0.75) 10.22±2.6110.22± 2.61 (0.81) 71.4±4.071.4± 4.0 (0.75) 0.382±.0160.382±.016 (0.74) NAM from scratch 0.388±.0030.388±.003 0.500±.0450.500±.045 5.05±.565.05±.56 15.91±3.9415.91± 3.94 80.1±12.280.1± 12.2 0.478±.0100.478±.010 Backfitting 0.380±.0060.380±.006 0.418±.0550.418±.055 4.55±.444.55±.44 12.13±2.6812.13± 2.68 48.0±4.448.0± 4.4 0.399±.0200.399±.020 Shared sections (ours) 0.376±.0040.376±.004 0.441±.0580.441±.058 4.41±.514.41±.51 14.45±3.2014.45± 3.20 51.9±2.751.9± 2.7 0.419±.0130.419±.013 TEAM (bagged stumps) 0.338±.0060.338±.006 0.412±.0560.412±.056 4.71±.384.71±.38 13.27±3.6113.27± 3.61 32.3±1.332.3± 1.3 0.425±.0090.425±.009 Table 2: Test MSE (mean ± s.d., 3 seeds) across six tabular datasets. Bold: best model of the additive family per dataset (ties within one s.d. both bolded); italics: dense reference where it wins. The first quad-soft row is near-additive—a fixed λ=1λ=1 for all datasets, no per-dataset tuning, with a small measured interaction residual (Table 3)—so its comparison against the exactly additive rows trades a weaker guarantee for accuracy. The second tunes λ per dataset by Hyperband-style successive halving on a held-out validation split (nine log-uniform candidates, 9→3→19→ 3→ 1, retrained on the full training set at the selected value; parenthetical numbers are the median selected λ∗λ^*): it matches the fixed default within noise everywhere except California Housing, where it wins by selecting λ∗≈0.06λ^*\!≈\!0.06—that is, by largely abandoning additivity. Target variances: 1.331.33, 0.650.65, 10.410.4, 8484, 279279, 3.73.7. On Concrete the best additive model beats the dense network by a factor of two. Method Constraint acts on Additivity guarantee ⟨Q2¯/Q2¯dense⟩6 datasets Q^2/ Q^2_dense _6 datasets Dense network — none 1.001.00 (reference) NAM structure (masks) exact, by architecture 0 TEAM structure (stumps) exact, by architecture 0 Backfitting fitting algorithm exact, by construction 0 Shared sections (ours) read-out exact, by read-out 0 Quad → crystallize (ours) behavior, then masks exact after pruning ∼10−6 10^-6 Quad soft, λ=1λ=1 (ours) behavior near-additive, measured 0.150.15 Table 3: Additivity contracts across methods, with the interaction residual in a single comparable column: the quadrilateral statistic normalized by the dense model’s (Q2¯ Q^2 is in units of y2y^2, so the raw values are not comparable across datasets; the ratio is dimensionless), averaged over the six benchmark datasets. Exactly additive methods are 0 identically and need no measurement; the soft model is the only method whose additivity is a measured quantity, and at the default λ=1λ=1 it retains on average 15%15\% of the dense interaction mass (per dataset: CalHousing 0.0160.016, Abalone 0.0270.027, Bike 0.170.17, Concrete 0.210.21, Boston 0.220.22, Wine 0.260.26), i.e. removes 7474–98%98\%, while paying the accuracy shown in Table 2. Table 2 also includes the penalized model itself, taken soft: the dense network trained with the quadrilateral loss at a fixed λ=1λ=1 for every dataset—the central recommendation of the balance heuristic developed alongside the table (median λbal≈1.5 _bal≈ 1.5 across the six datasets); we deliberately avoid per-dataset λ selection, which without a validation protocol would leak test information. It is the strongest model of the additive family on California Housing (where it also edges TEAM), Boston, and Bike (where it matches the dense network to within noise), and tied-strongest on Abalone, at the price of a weaker contract: its additivity is near-exact and measured rather than structural. This is also where our approach parts ways with the GA2M lineage [4, 24]. GA2M-style models prescribe the interaction order from outside—main effects plus an explicitly screened list of pairwise terms—so the analyst decides in advance which interactions may exist and of what order. The quad-soft model inverts this: no order is prescribed; the optimization retains whatever interaction mass the task loss defends against a global, continuously tunable budget, and the retained mass is then quantified (per coordinate, online) by the same statistic that enforced the budget. One buys structural semantics and enumerable terms; the other buys adaptivity and a certificate. The higher-order probes of Section 3 suggest the two can meet—using the quadrilateral statistics to discover which pairs a GA2M should include—but we do not pursue that here. Choosing λ: a dimensional argument and a balance heuristic. Two observations remove most of the guesswork. First, the task loss and the penalty have the same units—QiQ_i is a signed combination of function values, so [Qi2]E[Q_i^2] is measured in y2y^2 exactly like the MSE—hence λ is dimensionless, and a fixed O(1)O(1) value can transfer across datasets whose target variances differ by orders of magnitude (here 0.650.65 to 279279), as the sweet-spot figure confirms. Second, the natural scale within that O(1)O(1) range is the balance point λbal=MSE(0)/Q2¯(0) _bal=MSE(0)/ Q^2(0), the coefficient at which the two loss terms weigh equally for the unconstrained model: it is cheap to compute (one λ=0λ=0 run plus one pass of the statistic) and encodes the right diagnostic—a model whose interaction mass is large relative to its error is using interactions load-bearingly and warns against pushing, while a small ratio flags them as prunable. Across our six datasets λbal _bal spans 0.180.18 (CalHousing, where indeed any λ>0.1λ>0.1 already costs accuracy) through 1.051.05 (Bike), 1.151.15 (Boston), 1.91.9 (Abalone), ≈2.2≈2.2 (Concrete), to 4.44.4 (Wine), and on every dataset with a measured sweep the optimum falls at or below the balance point, in [0.1,1]⋅λbal[0.1,1]· _bal (Wine 0.7λbal0.7 _bal, CalHousing ∼0.5λbal 0.5 _bal, Concrete at the low end, ∼0.15λbal 0.15 _bal). The benchmark’s fixed λ=1λ=1 is the heuristic’s central recommendation (median λbal≈1.5 _bal≈ 1.5 across the six datasets); our earliest experiments had used λ=3λ=3, the aggressive upper edge of the Wine sweet-spot plateau, and switching to the heuristic default improved or matched every dataset (California Housing most visibly, from 0.3640.364 to 0.3280.328). As an out-of-sample test of the heuristic rather than a fit to it: Bike’s λbal=1.05 _bal=1.05 predicts λ∗≈0.5λ^*\!≈\!0.5–11, and λ=0.5λ=0.5 indeed yields 0.375±0.0160.375± 0.016, matching the dense network (0.3730.373). Per-dataset heuristic selection would improve the benchmark row further still; we keep the single fixed value to preserve its no-tuning reading. λ prices a constraint; do not tune it on accuracy. Table 2’s tuned row makes the comparison explicit and teaches two things. Practically, accuracy-tuned λ (Hyperband on a validation split) buys nothing over the heuristic default on the five small datasets—it matches within noise and, on two, validation noise makes it slightly worse—while the selected λ∗λ^* track the ordering of λbal _bal. Conceptually, the one dataset where tuning wins is the cautionary tale: on California Housing the tuner selects λ∗≈0.06λ^*\!≈\!0.06 and recovers most of the dense model’s accuracy by dissolving the very property the penalty exists to buy. λ is not a nuisance parameter: it prices a constraint, and tuning it on accuracy alone amounts to asking the market how much it will pay for a constraint it is free to ignore. The answer—nothing, where interactions pay—is correct and useless. The right protocol fixes the interpretability budget first (via the heuristic, or a target residual Q2¯ Q^2) and lets accuracy adjust, not the reverse. For calibration, the gap between the additive plateau (≈0.44≈0.44 for this width; ≈0.39≈0.39 for the deeper two-layer additive models of Table 2) and the unconstrained dense network (0.2840.284 for the two-layer model) reflects genuine interaction structure in the data, dominated by the joint dependence of house value on latitude and longitude, which no additive model can capture; this is consistent with the known geography of the dataset [40]. Figure 6: One-dimensional response sections as a function of λ, without pruning: fθf_θ evaluated along each feature with all others held at their (standardized) mean, curves centered for comparability. Red dashed: λ=0λ=0; viridis dark→ : increasing λ. Features with genuinely additive effects (MedInc, HouseAge) have λ-invariant sections; interaction-carrying features (Lat, Lon) have sections that change shape with λ and stabilize only for λ≳10λ 10. At λ=0λ=0 these sections are conditioning-point artifacts, not feature effects. Table 2 revises several single-dataset conclusions, which is precisely why it exists. The additivity tax is not universal. On California Housing additivity costs ∼35% 35\% MSE, but on three of the six datasets the best additive model beats the dense network—most dramatically on Concrete, where TEAM halves the dense error: with n≲4,000n 4,000 noisy samples the additive constraint acts as a regularizer, and where the true shape functions carry sharp thresholds (curing age, water–cement ratio) a well-matched additive learner can dominate outright. The framing “additivity trades accuracy for interpretability,” which our own Section 5.1 adopted, is a large-n statement. Trees remain formidable on tabular data. TEAM is the best exactly additive model on three of six datasets—and on Concrete the best model overall by a wide margin—and is never far from the front; only the near-additive quad-soft row, with its weaker contract, edges it on the largest dataset. This is consistent with the broader pattern that gradient-boosted trees dominate tabular benchmarks [18]; the interesting refinement is that the dominance survives the restriction to depth one, i.e. to the exactly additive subclass. The comparison between TEAM and the shared-section network is the cleanest in the table—identical functional class, identical read-out, differing only in the function approximator—and it splits: trees win where shape functions carry sharp axis-aligned structure (California’s geography, Concrete’s thresholds), the network wins on Abalone’s smooth monotone relationships. Why TEAM is hard to beat. Four ingredients, none exotic, compound. Stumps fit piecewise-constant shape functions, which capture the thresholds, plateaus, and saturations that dominate tabular relationships—and that a smooth GELU branch must approximate with finite curvature; California’s geographic effects, full of sharp administrative and coastal boundaries, are the clearest beneficiary. Split-based fitting is invariant to monotone feature transforms and therefore indifferent to the heavy tails that visibly destabilize the neural routes on the ratio features. Stagewise boosting is a functional-gradient coordinate descent that allocates its next unit of capacity to whichever feature currently has the largest residual signal—an adaptive budget the fixed per-branch widths of a NAM cannot imitate, and a one-at-a-time discipline that handles correlated features in the same spirit as backfitting. And the whole pipeline has essentially no optimization hyperparameters to mistune: no learning-rate schedule, no batch size, no seed-sensitive initialization, which is worth a great deal precisely in the small-n regime where Table 2 shows the elaborate neural pipelines faltering. The flip side, visible on Abalone, is that when the true shape functions are smooth and monotone, step-function bias wastes capacity on staircase artifacts and a smooth approximator wins. A distillation probe confirms that the advantage is representational rather than informational: training our penalized network on TEAM’s predictions instead of the labels performs worse than training on the labels directly (0.4110.411 vs. 0.4010.401 on California Housing), the smooth student’s inability to cheaply reproduce the teacher’s sharp axis-aligned steps swallowing the entire teacher margin. What TEAM knows, a network can learn from the data; how TEAM represents it, a smooth network cannot afford. Variance dwarfs method gaps on small data. On Boston, seed-to-seed standard deviations (2.72.7–3.93.9) exceed most pairwise method differences; the honest reading is that backfitting leads and the two crystallization routes trail, with everything else statistically entangled. Small-n tabular comparisons without seed bars are not informative, ours included wherever bars overlap. 7.1 Clinical Tabular: the Regime Where Additivity Pays The additive lineage earned its reputation in clinical risk, so we close the benchmark with three open clinical cohorts under the identical protocol: the scikit-learn diabetes progression task (regression, n=442n=442, d=10d=10), Pima Indians diabetes (n=768n=768, d=8d=8), and UCI Heart Disease (Cleveland, n=303n=303, d=13d=13); the two binary cohorts are trained as proper classification: every method fits an additive—or, for the dense reference and the penalized model, unconstrained—logit under cross-entropy, and we report test cross-entropy and AUROC. Dense NAM Backfitting Sections (ours) Quad soft λ=1λ=1 (ours) TEAM DiabProg MSE 0.320±.0470.320±.047 0.284±.0190.284±.019 0.294±.0230.294±.023 0.296±.0250.296±.025 0.292±.0350.292±.035 0.332±.0350.332±.035 Pima CE 0.565±.0570.565±.057 0.485±.0410.485±.041 — 0.477±.0500.477±.050 0.482±.0380.482±.038 0.476±.0340.476±.034 AUROC 0.812±.0300.812±.030 0.837±.0270.837±.027 — 0.845±.0280.845±.028 0.838±.0270.838±.027 0.839±.0240.839±.024 Heart CE 0.642±.2620.642±.262 0.355±.0610.355±.061 — 0.352±.0630.352±.063 0.354±.0700.354±.070 0.411±.0970.411±.097 AUROC 0.895±.0400.895±.040 0.923±.0220.923±.022 — 0.926±.0230.926±.023 0.924±.0240.924±.024 0.905±.0310.905±.031 Table 4: Clinical cohorts. DiabProg is regression (test MSE, 3 seeds; bold: best). Pima and Heart are trained as proper classification—additive logit, cross-entropy—and report test cross-entropy and AUROC over ten seeds. No entry is bolded on the binary rows deliberately: the additive routes are statistically indistinguishable there (paired differences within one seed s.d.; even the nominal winner flips between seed subsets), while every one of them dominates the dense reference, whose Heart fits are additionally unstable (±.26±.26). Classical backfitting under a logistic link is the local-scoring algorithm, which we do not implement (dashes). The result is the paper’s smallest-data regime taken to its conclusion: every additive route beats the dense network on every cohort and every metric (Table 4), and the win–win regime survives the change of loss with room to spare. At λ=1λ=1 the penalized model improves on its unconstrained counterpart on regression (0.2920.292 vs. 0.3100.310 MSE, λbal=6.6 _bal=6.6, the largest in the study) and, on the link scale, in classification—cross-entropy 0.4820.482 vs. 0.5650.565 on Pima and 0.3540.354 vs. 0.6420.642 on Heart over ten seeds, AUROC rising in both cases—while retaining a fraction of a percent of the logit-scale interaction mass. Within the additive family the binary rows carry no bold deliberately: at ten seeds the additive routes are statistically tied—paired cross-entropy differences of 0.0030.003–0.0080.008 against seed spreads several times larger, win counts of 66–77 out of 1010—and the tie is robust in an instructive way: TEAM, last on Pima in an earlier three-seed run, is nominally first there at ten, the kind of winner-flip that is exactly what a tie looks like through a small-sample keyhole. What survives every replication is the regime claim, not a ranking: at a few hundred samples every additive route dominates the dense reference, and there are too few samples to separate the additive routes from one another. The penalty lives on the link scale. Nothing in the construction is tied to squared error: for classification the network outputs a logit zθ(x)z_θ(x), the task loss is cross-entropy, and the quadrilateral penalty is applied to zθz_θ directly—driving it to zero makes the logit additive, i.e. the model a logistic GAM, exactly as classical additive modeling prescribes for generalized links; the same recipe extends to any GLM link. Table 4’s binary rows are produced this way, and the penalized logit removes 99.7%99.7\% of the logit-scale interaction mass on both cohorts while improving both metrics. These are open UCI-class cohorts rather than credentialed ICU data; the MIMIC-scale study that the asthma episode of Caruana et al. [36] calls for—and that Bayesian NAMs have begun [37], with credible intervals on MIMIC-I mortality shape functions—is the dedicated follow-up flagged in Section 10, where surrender curves would audit which clinical interactions a model refuses to give up. 8 Do the Routes Agree on the Functions? Matching test MSE does not imply matching models: additive decompositions are identifiable only up to constants under a product measure, and different procedures could in principle attain similar accuracy with different shape functions when features are correlated. Figure 7 overlays the centered shape functions of the benchmark’s four additive routes together with the penalized model’s mean-conditioned sections at four values of λ. Figure 7: Agreement and convergence in one picture, centered per feature. Four procedures with disjoint mechanisms—fixed masks (NAM, dotted grey), cyclic partial-residual refitting (backfitting, dashed green), the shared-section model (red), and bagged boosted stumps (TEAM, purple)—recover mutually consistent functions on the high-density regions of every feature, diverging only in sparse tails where the data cannot constrain them. The blue gradient overlays the penalized dense model’s sections at λ=0.3,1,3,10λ=0.3,1,3,10: at small λ residual interactions leak into the sections and pull them off the common curves, and as λ grows the sections collapse onto the consensus—the λ-convergence of Section 5.3 rendered directly against the independent routes that define the target. The agreement is striking given how different the mechanisms are—TEAM’s piecewise-constant staircases trace the same profiles the smooth neural fits draw, despite sharing no machinery with them. On median income all five curves coincide over the full range; on house age, latitude, and longitude they coincide over the data-dense regions, including the nontrivial multi-modal profiles of the two geographic coordinates—which, being the one-dimensional surrogates of a fundamentally two-dimensional effect, had no a priori reason to be stable across procedures. Disagreement concentrates exactly where it should: in the extreme tails of the heavy-tailed ratio features, where standardized values beyond a few σ correspond to a handful of block groups and the components are unconstrained by data. The WTA route (the weakest by MSE among the six) is also the visible outlier where curves separate, consistent with its soft phase distorting components that fine-tuning cannot fully repair. We read this as evidence that on this dataset the additive projection of the regression function is a well-identified object: the plateau at MSE ≈0.38≈0.38–0.400.40 is not a family of accuracy-equivalent but functionally distinct solutions, but a single solution reached from five directions. A practical corollary: when two independent additive-fitting procedures disagree on a shape function in a data-dense region, at least one has not converged. Identifiability under correlated features. The tail disagreements deserve a sharper diagnosis than data scarcity alone, because a structural issue lurks beneath them. When features are strongly dependent—average rooms and average bedrooms here, with correlation above 0.80.8—the additive sum fi(xi)+fj(xj)f_i(x_i)+f_j(x_j) is well determined on the data manifold, but the individual components are not: adding g(xi)g(x_i) to one and subtracting its best predictor g^(xj) g(x_j) from the other changes each shape function while leaving predictions nearly unchanged wherever the two features co-vary. This is concurvity, the nonlinear analogue of collinearity, and it means that in strongly correlated directions our agreement figure shows five estimates of an under-determined quantity—their consistency in dense regions reflects shared inductive biases toward smoothness as much as identification by the data. Two families of remedies exist. Decompositions can be defined to be identifiable under dependence: the weighted functional ANOVA of Hooker [12] imposes orthogonality relative to the joint distribution, accumulated local effects [31] build components from local derivatives so that correlated marginals do not leak into one another, and purification [32] post-hoc moves mass between components into canonical form. Alternatively, regularization can select among the equivalent decompositions during training: the concurvity penalty of Siems et al. [26] shrinks the correlation between fitted components, pinning down the representative with least mutual cancellation. The two regularizers compose naturally—the quadrilateral penalty controls interactions within the model while the concurvity penalty controls redundancy between its additive components—and their combination, quadrilateral for additivity plus concurvity for identifiability, is in our view the most promising follow-up experiment this paper does not run. 9 Anatomy of the Failure Modes Three qualitatively different failures surfaced across the routes, and juxtaposing them reveals a common structure. Adam’s stale moments moved pruned weights off their zeros, silently inflating the additivity gap by four orders of magnitude (Section 6.1); unconstrained winner-take-all group assignment collapsed onto a single feature, killing the branches for median income and both geographic coordinates (Section 6.1); and sequential backfitting on shared weights diverged because updating one section moved all the others, invalidating the partial residuals mid-sweep (Section 6.3). In each case, a guarantee was imported from a setting whose preconditions the new setting quietly violated. Mask-based sparsity assumes the entire update rule respects the mask, but Adam’s state is part of the update rule and predates the mask. Winner-take-all assignment inherits its intuition from balanced clustering, but nothing in arg-max dynamics enforces balance, and rich-get-richer feedback does the rest. Backfitting’s convergence proof [16] assumes parameter-disjoint components, and weight sharing removes exactly that hypothesis. None of the three failures announces itself as a violated precondition: the first manifests as a functional property silently degrading while accuracy looks fine, the second as a model inexplicably worse than linear regression, the third as training that improves for fifteen sweeps and then explodes. The practical lesson we take is that constraints on learned functions must be audited against the whole training system—optimizer state, assignment dynamics, staleness of cached quantities—not just against the loss and the architecture. Each fix was cheap once the failure was understood (re-project after every step; cap group capacities; replace sequential residual fitting with a joint objective), and each failure is generic enough that we expect it to be present, and silent, in existing pipelines that combine masking, learned grouping, or alternating optimization with adaptive optimizers. 10 Limitations The λ-sweep, section analysis, and crystallization comparison are conducted on California Housing only; Section 7 extends the converging routes to three further datasets and indeed overturns two single-dataset readings (the universality of the additivity tax, and the competitiveness of the crystallization pipeline), but six small tabular benchmarks remain a narrow base; the λ-sweep itself is single-dataset, while the surrender-curve analysis has now been replicated over five seeds on all four datasets (Section 5.2), which is what exposed the fragility of its initial single-seed reading. The penalty is enforced on the data manifold, so no claim is made about interactions off-distribution, and the stochastic estimator has a noise floor (∼10−4 10^-4 here) below which it cannot certify additivity—exactness always requires crystallization. Intermediate λ yields behaviorally near-additive but structurally unconstrained models, which is a feature for the trade-off analysis and a bug for anyone hoping the knob traverses a nested family of interaction orders; the block-restricted variant of the penalty, and its comparison against GA2M-style baselines with a learned (Lat, Lon) block, is the natural next experiment and is not performed here. The shared-net section model leads the neural additive routes on two of four datasets, but the implicit-basis-sharing hypothesis we offer for its advantage has not been isolated by ablation, and its d+1d+1 forward passes per example are untested beyond d=13d=13. The KAM analogy of Section 5.2 is offered as a reframing, not a result. Finally, the per-step cost of the penalty is four forward passes times |S||S| sampled coordinates, roughly tripling training cost at our settings; scaling behavior in d is untested. The natural next deployment is clinical tabular risk—MIMIC-style ICU mortality and readmission cohorts, where the additive lineage earned its reputation and where a measured interaction budget has audit value—which requires credentialed data access we deliberately leave to a dedicated study. 11 Conclusion We compared routes from dense networks to additive predictors under a fixed budget, on six tabular datasets. A four-point finite-difference statistic serves throughout as instrument, penalty, and certificate: it turns additivity into a continuous training-time dial, reveals that most learned interactions are redundant and removable almost for free, and—because it decomposes by coordinate—exposes interaction as a quantity that evolves under optimization pressure. Watching that evolution taught us something we did not expect and then, on replication, corrected our first reading of it: which interactions survive an additivity regularizer is only weakly and unstably predicted by how large they were, so post-hoc interaction rankings are unreliable guides to what a regularized model will keep, and only online measurement on the actual run reveals it. The route comparison delivered its own reversals. Behavior-first-then-crystallize dominates weight-space constraints, classical backfitting with neural smoothers quietly matches modern per-feature architectures, one ordinary dense network read out through its own mean-conditioned sections—the shared-section model, the hard-projection counterpart of the quadrilateral penalty, sharing its λ→∞λ→∞ target though not its parametrization—is the strongest exactly additive neural model we trained, and the same network trained soft with the penalty is stronger still where near-additivity suffices, and bagged boosted stumps—adaptive capacity allocation, transform invariance, nothing to mistune—remain the method to beat on tabular data. None of these rankings is universal: in low data, additivity itself becomes a regularizer—the quadrilateral penalty then improves accuracy and interpretability simultaneously—while multi-stage pipelines grow fragile; and the coefficient that prices the constraint should be set by the interpretability budget, not tuned on accuracy, which where interactions pay will simply buy the constraint away. Conceptually, the soft regime earns a defense that exact additivity cannot give: since perfect additive projections may demand shape functions too wild to read—the statistical shadow of Kolmogorov’s superposition theorem, whose univariate representations are necessarily non-smooth, and the tension that forced ExU activations into the original NAM—a dense model with a small, measured interaction residual and smooth sections can be the more interpretable object. And the three silent failures we hit along the way (Section 9) taught a single lesson we suspect generalizes: each had imported a guarantee whose preconditions the new setting had quietly voided—and none of them announced it. References Hastie and Tibshirani [1986] T. Hastie and R. Tibshirani. Generalized additive models. Statistical Science, 1(3):297–310, 1986. [2] R. Agarwal, L. Melnick, N. Frosst, X. Zhang, B. Lengerich, R. Caruana, and G. Hinton. Neural additive models: Interpretable machine learning with neural nets. In Advances in Neural Information Processing Systems, 2021. Lou et al. [2012] Y. Lou, R. Caruana, and J. Gehrke. Intelligible models for classification and regression. In KDD, 2012. [4] Y. Lou, R. Caruana, J. Gehrke, and G. Hooker. Accurate intelligible models with pairwise interactions. In KDD, 2013. [5] C.-H. Chang, R. Caruana, and A. Goldenberg. NODE-GAM: Neural generalized additive model for interpretable deep learning. In International Conference on Learning Representations, 2022. [6] F. Radenovic, A. Dubey, and D. Mahajan. Neural basis models for interpretability. In Advances in Neural Information Processing Systems, 2022. [7] I. M. Sobol’. Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates. Mathematics and Computers in Simulation, 55(1–3):271–280, 2001. [8] S. M. Lundberg and S.-I. Lee. A unified approach to interpreting model predictions. In Advances in Neural Information Processing Systems, 2017. [9] J. D. Janizek, P. Sturmfels, and S.-I. Lee. Explaining explanations: Axiomatic feature interactions for deep networks. Journal of Machine Learning Research, 22(104):1–54, 2021. [10] M. Tsang, D. Cheng, and Y. Liu. Detecting statistical interactions from neural network weights. In International Conference on Learning Representations, 2018. Tsang et al. [2018] M. Tsang, H. Liu, S. Purushotham, P. Murali, and Y. Liu. Neural interaction transparency (NIT): Disentangling learned interactions for improved interpretability. In Advances in Neural Information Processing Systems, 2018. Hooker [2007] G. Hooker. Generalized functional ANOVA diagnostics for high-dimensional functions of dependent variables. Journal of Computational and Graphical Statistics, 16(3):709–732, 2007. [13] P. Molchanov, S. Tyree, T. Karras, T. Aila, and J. Kautz. Pruning convolutional neural networks for resource efficient inference. In International Conference on Learning Representations, 2017. [14] W. Wen, C. Wu, Y. Wang, Y. Chen, and H. Li. Learning structured sparsity in deep neural networks. In Advances in Neural Information Processing Systems, 2016. [15] D. P. Kingma and J. Ba. Adam: A method for stochastic optimization. In International Conference on Learning Representations, 2015. Buja et al. [1989] A. Buja, T. Hastie, and R. Tibshirani. Linear smoothers and additive models. The Annals of Statistics, 17(2):453–510, 1989. Friedman [2001] J. H. Friedman. Greedy function approximation: A gradient boosting machine. The Annals of Statistics, 29(5):1189–1232, 2001. Grinsztajn et al. [2022] L. Grinsztajn, E. Oyallon, and G. Varoquaux. Why do tree-based models still outperform deep learning on typical tabular data? In Advances in Neural Information Processing Systems, 2022. Di Cecco [2024] A. Di Cecco. TEAM: Tree ensemble additive models. GitHub repository, https://github.com/AntonioDiCecco/TEAM, 2024. Di Cecco and Gianfagna [2025a] A. Di Cecco and L. Gianfagna. Explainable AI with Python. Springer, Cham, 2025. ISBN 978-3-031-92229-9. Di Cecco and Gianfagna [2025b] A. Di Cecco and L. Gianfagna. Additive models for interpretability. In Explainable AI with Python, chapter 6. Springer, Cham, 2025. Nori et al. [2019] H. Nori, S. Jenkins, P. Koch, and R. Caruana. InterpretML: A unified framework for machine learning interpretability. arXiv preprint arXiv:1909.09223, 2019. Dubey et al. [2022] A. Dubey, F. Radenovic, and D. Mahajan. Scalable interpretability via polynomials. In Advances in Neural Information Processing Systems, 2022. Yang et al. [2021] Z. Yang, A. Zhang, and A. Sudjianto. GAMI-Net: An explainable neural network based on generalized additive models with structured interactions. Pattern Recognition, 120:108192, 2021. Liu et al. [2024] Z. Liu, Y. Wang, S. Vaidya, F. Ruehle, J. Halverson, M. Soljačić, T. Y. Hou, and M. Tegmark. KAN: Kolmogorov–Arnold networks. arXiv preprint arXiv:2404.19756, 2024. Siems et al. [2023] J. Siems, K. Ditschuneit, W. Ripken, A. Lindborg, M. Schambach, J. S. Otterbach, and M. Genzel. Curve your enthusiasm: Concurvity regularization in differentiable generalized additive models. In Advances in Neural Information Processing Systems, 2023. Arnold [1963] V. I. Arnold. Proof of a theorem of A. N. Kolmogorov on the preservation of conditionally periodic motions under a small perturbation of the Hamiltonian. Russian Mathematical Surveys, 18(5):9–36, 1963. Ross et al. [2017] A. S. Ross, M. C. Hughes, and F. Doshi-Velez. Right for the right reasons: Training differentiable models by constraining their explanations. In International Joint Conference on Artificial Intelligence, 2017. Erion et al. [2021] G. Erion, J. D. Janizek, P. Sturmfels, S. M. Lundberg, and S.-I. Lee. Improving performance of deep learning models with axiomatic attribution priors and expected gradients. Nature Machine Intelligence, 3(7):620–631, 2021. Rieger et al. [2020] L. Rieger, C. Singh, W. J. Murdoch, and B. Yu. Interpretations are useful: Penalizing explanations to align neural networks with prior knowledge. In International Conference on Machine Learning, 2020. Apley and Zhu [2020] D. W. Apley and J. Zhu. Visualizing the effects of predictor variables in black box supervised learning models. Journal of the Royal Statistical Society: Series B, 82(4):1059–1086, 2020. Lengerich et al. [2020] B. Lengerich, S. Tan, C.-H. Chang, G. Hooker, and R. Caruana. Purifying interaction effects with the functional ANOVA: An efficient algorithm for recovering identifiable additive models. In International Conference on Artificial Intelligence and Statistics, 2020. Bordt and von Luxburg [2023] S. Bordt and U. von Luxburg. From Shapley values to generalized additive models and back. In International Conference on Artificial Intelligence and Statistics, pages 709–745. PMLR, 2023. Kolmogorov [1957] A. N. Kolmogorov. On the representation of continuous functions of many variables by superposition of continuous functions of one variable and addition. Doklady Akademii Nauk SSSR, 114:953–956, 1957. Vitushkin [1954] A. G. Vitushkin. On Hilbert’s thirteenth problem. Doklady Akademii Nauk SSSR, 95:701–704, 1954. Caruana et al. [2015] R. Caruana, Y. Lou, J. Gehrke, P. Koch, M. Sturm, and N. Elhadad. Intelligible models for healthcare: Predicting pneumonia risk and hospital 30-day readmission. In ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, 2015. Bouchiat et al. [2024] K. Bouchiat, A. Immer, H. Yèche, G. Rätsch, and V. Fortuin. Improving neural additive models with Bayesian principles. In International Conference on Machine Learning, pages 4416–4443. PMLR, 2024. Friedman and Popescu [2008] J. H. Friedman and B. E. Popescu. Predictive learning via rule ensembles. The Annals of Applied Statistics, 2(3):916–954, 2008. Sorokina et al. [2008] D. Sorokina, R. Caruana, M. Riedewald, and D. Fink. Detecting statistical interactions with additive groves of trees. In International Conference on Machine Learning, 2008. Pace and Barry [1997] R. K. Pace and R. Barry. Sparse spatial autoregressions. Statistics & Probability Letters, 33(3):291–297, 1997.

Instant research discovery