Numerical Calibration Without Overfitting

Sparse numerical data can calibrate a two-body model without replacing its analytic structure only when the fitted freedom is declared in advance, the loss compares shared invariants with numerical uncertainty accounted for, model choices are made on validation data, and a blind test set is opened once. Calibration estimates unresolved coefficients. It does not re-derive PN, PM, self-force, or EOB terms, and low training residual is not evidence of predictive accuracy. The two examples below contrast a successful fit to a known synthetic target with a failed fit to actual numerical-relativity scattering data.

Helpful background. The Einstein Equation supplies the governing dynamics; Weak Metric Perturbations and Gauge Freedom supplies analytic weak-field constraints; Presymplectic Potentials and Currents supplies canonical invariance; Computational Relativity supplies convergence and benchmark discipline; Gravitational Radiation and Waveforms supplies extraction conventions; and Relativistic Two-Body Dynamics supplies the model regimes and invariant variables.

Calibration as a constrained inverse problem

Write a dynamical prediction as

$$\mathbf y_{\rm model}(\lambda;\boldsymbol\theta) =\mathbf y_{\rm analytic}(\lambda) +\sum_{a=1}^{p}\theta_a\mathbf f_a(\lambda),$$

where $\lambda$ denotes invariant physical parameters such as $\nu$, dimensionless spins, frequency, energy, and angular momentum. The analytic part contains every coefficient already fixed by PN, PM, self-force, test-body, symmetry, and balance information. The basis functions $\mathbf f_a$ represent only unresolved higher-order or strong-field completion. Moving a known coefficient into $\boldsymbol\theta$ weakens the model and lets numerical data erase an analytic check.

For configuration $i$, compare a vector of invariants rather than coordinate trajectories:

$$\mathbf r_i(\boldsymbol\theta) =\mathbf y_{\rm model}(\lambda_i;\boldsymbol\theta) -\mathbf y_i^{\rm num}.$$

Examples are binding energy at fixed orbital frequency, scattering angle at fixed incoming $(E,J)$, periastron advance, radiated energy, or aligned waveform phase after a declared alignment. With numerical covariance $C_i$, a suitable objective is

$$\mathcal L(\boldsymbol\theta) =\sum_{i\in{\cal T}} \mathbf r_i^{\mathsf T}C_i^{-1}\mathbf r_i +\lambda_{\rm reg}\lVert L\boldsymbol\theta\rVert_2^2.$$

The training set ${\cal T}$ determines coefficients. The regularizer must encode a stated smoothness, naturalness, or asymptotic prior; it is not a device for hiding a poor fit. If extraction, truncation, or finite-radius errors are correlated across configurations, $C_i$ must be replaced by the appropriate block covariance. Treating repeated resolutions as independent physical systems artificially multiplies information.

Training, validation, test, and stress domains

A random row split is usually unsafe because adjacent times, resolutions, modes, or simulations from the same physical configuration are highly correlated. Split by configuration family before any feature construction:

Partition	Permitted use	Example blocked domain
Training	Estimate $\boldsymbol\theta$ for each predeclared candidate basis	Interior mass ratios and spins, with all resolutions of a configuration kept together
Validation	Choose basis, prior, regularization, optimizer tolerances, and stopping rule	Withheld mass-ratio bands, spin corners, and a distinct nonspinning family
Blind test	Report final predictive performance once	Unseen simulation IDs spanning both interpolation and edge configurations
Extrapolation stress set	Locate failure, not certify accuracy	Points beyond the convex hull in mass ratio, spin, eccentricity, or encounter energy

The blind test cannot be used to select a basis and remain a test. If its result motivates a change, the revised model needs a newly sequestered test set or must be labeled as development-stage. Likewise, an NR waveform used to tune an EOB coefficient cannot validate that coefficient. Validation must compare a physical invariant or waveform excluded from both coefficient estimation and model selection.

A published EOB example makes the distinction concrete. The SEOBNRv5 construction calibrated three declared dynamical parameters with 442 numerical-relativity waveforms, defined its mismatch and merger-time requirements, used posterior information rather than only best fits, and reserved 21 of 39 nonspinning simulations for validation Pompili et al. 2023, §IV.A–C, pp. 11–15, Eqs. (65)–(74), PDF. This is evidence for a documented workflow, not a claim that the model or split is universally optimal. The same paper reports difficult high-spin and high-mass-ratio regions and comparisons with independent approximants in its performance section Pompili et al. 2023, §V.A, pp. 15–18, PDF.

Numerical uncertainty also sets a floor. The SXS catalog compares the two highest resolutions and checks whether that difference decreases relative to the next pair; it warns that resolution changes can slightly change physical parameters and that extrapolation and truncation errors are not automatically separable Boyle et al. 2019, §4, pp. 24–28, Eqs. (22)–(24), PDF. A calibration residual below this floor cannot be interpreted as increased physical accuracy.

Reproducible calculations

A reproducible fit keeps the input data and their versions, the physical definitions, every fixed and adjustable coefficient, the split and loss, and the actual numerical environment. It also retains predictions and residuals, so another reader can distinguish a fitting error from a poor model. A nonlinear optimizer needs checks from different starts; a full-rank linear least-squares problem has a unique minimum that can instead be checked with independent QR and SVD solvers. The physical example below includes downloadable source, a fixed design, results and a complete rerun command.

Benchmark and acceptance criteria

Before fitting physical data, use a synthetic target with known coefficients. One minimal benchmark is

$$Q_\star(z,\nu) =Q_{\rm analytic}(z,\nu) +\nu^2\left(0.80\,z^3-0.35\,\nu z^2\right),$$

where $z$ is an abstract dimensionless input, not the circular redshift. The $\nu^2$ factor leaves prescribed test-mass and first-order self-force terms untouched. This target checks a fitting procedure; it is not NR data or a calibrated gravitational observable.

Here is a complete example using only the Python 3 standard library. It fits $y=Q_\star-Q_{\rm analytic}$ with independent Gaussian errors of known standard deviation $\sigma=2\times10^{-5}$. The basis and all three partitions are fixed before fitting; validation and test points are never used to change the model. Save the code as calibration-example.py and run python3 calibration-example.py.

from math import sqrt
from random import Random

rng = Random(1729)
sigma = 2e-5
truth = (0.80, -0.35)
domains = {
    "train": [(z, nu) for z in (.4, .6, .8, 1.)
              for nu in (.10, .17, .24)],
    "validation": [(.5, .12), (.5, .21), (.9, .12), (.9, .21)],
    "test": [(.45, .14), (.55, .19), (.75, .14), (.95, .19)],
}

def features(z, nu):
    return (nu**2*z**3, nu**3*z**2)

data = {}
for name, points in domains.items():
    data[name] = []
    for z, nu in points:
        f = features(z, nu)
        y = sum(a*b for a, b in zip(truth, f)) + rng.gauss(0, sigma)
        data[name].append((f, y))

train = data["train"]
a = sum(f[0]**2 for f, y in train)
b = sum(f[0]*f[1] for f, y in train)
d = sum(f[1]**2 for f, y in train)
c = [sum(f[k]*y for f, y in train) for k in (0, 1)]
det = a*d-b*b
theta = ((d*c[0]-b*c[1])/det, (a*c[1]-b*c[0])/det)
sd = (sigma*sqrt(d/det), sigma*sqrt(a/det))
for k in (0, 1):
    print(f"theta[{k}] = {theta[k]:.6f} +/- {sd[k]:.6f}")
for name, rows in data.items():
    residuals = [(sum(t*v for t, v in zip(theta, f))-y)/sigma
                 for f, y in rows]
    rms = sqrt(sum(r*r for r in residuals)/len(rows))
    print(f"{name}: n={len(rows)}, standardized RMS={rms:.6f}")

The checked result is $\theta_0=0.800512\pm0.001245$ and $\theta_1=-0.352285\pm0.005058$, with one-standard-deviation errors from $\sigma^2(F^{\mathsf T}F)^{-1}$. The standardized RMS values are

Partition	Points	Standardized RMS
Training	12	0.787011
Validation	4	0.735561
Test	4	1.081949

Both injected coefficients lie within one reported standard deviation in this realization. With only four test points, a score near one is a small demonstration, not a coverage study. A physical pipeline should repeat known synthetic targets across noise realizations and include deliberately misspecified targets. The explicit two-by-two inverse is transparent here; larger or poorly conditioned fits should use QR or SVD instead of forming normal equations.

A production calibration then needs never-fit invariants or waveforms. A defensible acceptance rule is:

• all qualifying optimizer starts give validation loss within 1% and predictions differing by less than one numerical-data standard deviation;
• enlarging training through nested, prespecified subsets changes validation predictions by less than the numerical floor; the blind set is assessed only after this choice is fixed;
• deleting any one configuration family does not move a prediction by more than its declared combined uncertainty;
• the selected basis beats the analytic-only baseline on validation and blind test, not merely on training;
• extrapolation failures are reported as a domain boundary rather than averaged into a favorable global score.

These numbers are protocol thresholds, not universal constants. They must be fixed before the blind set is opened and justified against the numerical covariance and intended observable.

A failed fit to physical scattering data

Now consider the public scattering simulations of Long et al. The target is the mean of the two individual deflections in the initial nominal center-of-mass frame, $\theta_{\rm avg}=(\theta_1+\theta_2)/2$. For unequal masses this is not generally the relative angle supplied by an EOB model: radiative recoil makes the distinction physical Long et al. 2025, §II.D and Eqs. (26)–(27), PDF.

We deliberately choose an incomplete comparison model, so this is a test of failure diagnosis rather than a proposed EOB calibration. Let $\Gamma=E/M$, $\nu=m_1m_2/M^2$, $\gamma=1+(\Gamma^2-1)/(2\nu)$, and $j=J_{\rm orb}/(GM\mu)$, with $\mu=\nu M$ and $c=1$. Starting with the exact Schwarzschild scattering angle, define

$$\begin{aligned} \theta_{\rm base}(\gamma,j,\nu) &=\theta_{\rm Schw}(\gamma,j) +\frac{3\pi(5\gamma^2-1)}{4j^2} \left(\frac1\Gamma-1\right),\\ \theta_{\rm model} &=\theta_{\rm base}+\nu x^3 P(x,t,\nu), \qquad x=\frac4j,\quad t=\frac{\gamma-1}{0.1}. \end{aligned}$$

The correction restores the finite-mass 2PM coefficient, while the polynomial starts at 3PM. Thus the model preserves exact test-mass scattering and the full 1PM and 2PM angles; it does not impose all known higher-order physics Damour 2020, Eqs. (1.4), (1.9)–(1.11), (2.45)–(2.46), PDF. The three predeclared bases for $P$ are $(1)$, $(1,x,t,\nu)$, and $(1,x,t,\nu,x^2,xt,t^2)$.

The fixed split groups whole nominal mass ratios $q=m_1/m_2$: training uses $q=1,4,7$, validation uses $q=3,6$, and test uses $q=2,5,8,9,10$. Eligibility requires nominally nonspinning data, $1.0199\leq\gamma\leq1.16$, and a scattering Schwarzschild baseline. This leaves 18, two and five cases, respectively. The fits use the authors’ asymmetric error widths as positive scales; no Gaussian likelihood or independence of the two angle errors is assumed. The design and code were fixed before the test residuals were evaluated.

Polynomial	Training RMS, degrees	Validation RMS, degrees
Constant	73.791934	61.888296
Linear	46.376724	22.618076
Quadratic	25.522674	0.013501

The smallest validation RMS of residual divided by its positive source scale selects the quadratic. Its five-case test RMS is also small: $0.055791^\circ$, or 1.0471 source error scales. Nevertheless the model is rejected. Its largest training error is 50.2%, far above the predeclared illustrative requirement of at most 1% error on every training, validation and test case. Independent QR and SVD calculations agree on this result. The training errors already warranted rejection before opening the test; its frozen scores are shown here to demonstrate how a narrow test domain can also give false reassurance.

Why do the withheld scores look so good? They sample only $\gamma\simeq1.02$, $j\simeq4.8$, whereas training reaches $\gamma\simeq1.09$ and $j\simeq7.33$. Holding out mass ratios does not create coverage in energy and angular momentum. The scaled design condition number is about 3717; changing the training subset or adding weak regularization substantially moves validation predictions. The one otherwise eligible lower-energy stress case misses by $42.98^\circ$. Fourteen spinning cases are unsupported, and twenty other rows have no scattering Schwarzschild baseline; that last limitation is not a claim of binary capture. No model choice is changed after inspecting these results.

This is model misspecification combined with inadequate validation coverage, not a rule that seven coefficients always overfit. A good calibration must first represent its training domain, then predict genuinely new configurations throughout the domain for which accuracy is claimed.

To reproduce the calculation, save the Python script and fixed design in one folder, then run:

python3.12 -m venv .venv
.venv/bin/python -m pip install numpy==2.2.6 scipy==1.18.1 mpmath==1.3.0
.venv/bin/python calibrate-scattering.py --fetch --input data/ScatteringAngles.dat --output-dir results

The instructions and provenance give the source checksum, extraction definitions, checks and comparison tolerances; the recorded results include coefficients, conditioning, residuals and sensitivity tests. The script downloads the original table and checks its bytes. Repeating this published exercise reproduces an already opened test, not a new blind test.

Convergence, identifiability, and stability

Optimization convergence is not model convergence. Repeat the fit from widely separated starts, inspect the Hessian or posterior covariance, and plot residuals against every invariant input. A small Hessian eigenvalue signals a nearly unconstrained coefficient combination. Individual parameters may then vary greatly while predictions remain stable; report both parameter and predictive uncertainty. If predictions also vary, reduce the basis, add genuinely informative data, or impose a physical prior whose influence is shown by an ablation.

Increase the training set in nested blocks and vary the basis and regularization. Stable training loss with degrading validation loss is overfitting. Oscillatory coefficients under modest basis changes reveal an ill-conditioned representation. Correlated residuals along mass ratio or spin reveal missing structure even when their mean is zero. A posterior over calibration parameters can propagate coefficient uncertainty, but it covers only the assumed model family and data likelihood. Recent work models SEOBNR calibration posteriors and checks them against waveforms not used to train the uncertainty model; it also finds that calibration uncertainty alone does not cover all higher-mode or out-of-domain modeling error Pompili, Buonanno, and Pürrer 2026, §§II–IV and Appendix B, PDF.

Failure modes

Data leakage. Modes, time samples, resolutions, or near-duplicate simulations from one configuration cross partitions. Group before splitting and publish the group manifest.

Parameter nonidentifiability. Several coefficients produce the same invariant over the sampled domain. Use Hessian/posterior eigenvectors, not optimizer success, to diagnose it.

Optimizer trapping. Different starts reach distinct basins with materially different validation predictions. Report the ensemble and do not select a basin using the blind test.

Misweighted data. Ignoring numerical covariance lets densely sampled or noisy configurations dominate. Whiten residuals with the documented covariance and test covariance sensitivity.

Basis instability. Added high-order terms cancel each other inside the training hull and explode outside it. Compare analytic-only, minimal, and extended bases through ablations and stress tests.

Extrapolation. A fit is evaluated beyond the convex hull in mass ratio, spin, eccentricity, or frequency. Give the distance or domain violation explicitly; a smooth curve is not evidence.

Underestimated numerical error. Extraction, resolution, and alignment effects are treated as exact. Repeat the calibration across plausible numerical variants and propagate their shift.

Solved model-selection diagnostic

Two candidate completions have standardized root-mean-square residuals

Model	Training	Validation	Blind test at high mass ratio
Analytic plus two fitted terms	0.20	0.90	1.10
Analytic plus six fitted terms	0.05	0.85	4.10

The six-term fit also has Hessian eigenvalues whose ratio is $10^9$. Before opening either blind result, a predeclared simplicity and conditioning rule selected the two-term model because the validation improvement was too small to justify the extra freedom. What does the blind test establish, and may it be used to tune a seventh term?

Solution

The selected two-term model has a blind residual compatible with the quoted data scale at the tested edge. The alternative’s failure illustrates the danger already suggested by its ill-conditioned fit; it did not select the two-term model. The Hessian ratio diagnoses a poorly constrained parameter combination only when its absolute curvature is small enough relative to the noise and prior scales. The ratio alone flags numerical conditioning in the chosen parameter scaling, not a proof that the underlying physics is wrong. If the choice were instead made after seeing both blind scores, those scores would become model-selection evidence and fresh test data would be required.

The blind result may diagnose a high-mass-ratio boundary, but it cannot tune a seventh term while remaining a blind validation. If the model is revised, that point becomes development data and a new test family must be sequestered. Until then, report the simpler model’s tested domain and the alternative model’s blind-test failure separately.

Calibration supplies inputs to the uncertainty budget, but its posterior must not be mistaken for PN truncation, shared analytic discrepancy, or numerical extraction error.

References

• Boyle, Michael, et al. “The SXS Collaboration Catalog of Binary Black Hole Simulations.” Classical and Quantum Gravity 36, 195006 (2019). DOI. Open PDF.
• Damour, Thibault. “Classical and Quantum Scattering in Post-Minkowskian Gravity.” Physical Review D 102, 024060 (2020). DOI. Open PDF.
• Long, Oliver, et al. “Highly Accurate Simulations of Asymmetric Black-Hole Scattering and Cross Validation of Effective-One-Body Models.” Physical Review D 112, 124039 (2025). DOI. Open PDF.
• Pompili, Lorenzo, Alessandra Buonanno, and Michael Pürrer. “Accounting for Numerical-Relativity Calibration Uncertainty in Gravitational-Wave Modeling and Inference.” arXiv:2410.16859v3 (2026). Open PDF.
• Pompili, Lorenzo, et al. “Laying the Foundation of the Effective-One-Body Waveform Models SEOBNRv5: Improved Accuracy and Efficiency for Spinning Nonprecessing Binary Black Holes.” Physical Review D 108, 124035 (2023). DOI. Open PDF.