Reduced-Order and Surrogate Models

A reduced basis compresses a waveform family in time; empirical interpolation selects informative time nodes; a surrogate additionally learns how the waveform values at those nodes vary over source parameters. These are distinct approximation stages with distinct errors. A surrogate accelerates evaluation of a specified fiducial family—it does not improve the physical accuracy of that family or extend its calibration domain. This page develops the construction and a reproducible error hierarchy. Waveform production remains with Relativistic Two-Body Dynamics and Computational Relativity, while detector-specific inference verdicts remain downstream.

Helpful background. Einstein Equations supplies the underlying field theory; Waves, Green Functions, and Geometric Optics supplies waveform spaces and interpolation; Relativistic Two-Body Dynamics supplies analytic waveform families; Computational Relativity supplies numerical training waveforms; and Gravitational-Wave Data Analysis and Inference supplies application-specific accuracy criteria.

Offline compression and online prediction

Let $h(t;\lambda)$ be a fiducial mode family on a compact declared parameter domain $\mathcal D$. All training waveforms must share the same complex polarization $h=h_+-ih_\times$, unit-normalized spin-weight $-2$ harmonics, retarded-time grid or resampling rule, and named inertial or coprecessing frame. The usual four-stage construction is:

1. choose a training set in $\mathcal D$ and construct a reduced basis $\{e_i(t)\}_{i=1}^m$;
2. select $m$ empirical time nodes $T_i$ from that basis;
3. fit the amplitude/phase or real/imaginary waveform values $h(T_i;\lambda)$ over parameter space;
4. reconstruct the full time series from the fitted nodal values and precomputed basis functions.

Field and collaborators derive this reduced-basis/empirical-interpolation/surrogate pipeline in Field et al. 2014, §§II–III, Eqs. (4)–(33), PDF. The expensive waveform evaluations and basis construction are offline. The online cost is dominated by evaluating $m$ parameter fits and combining $m$ length-$L$ basis functions, scaling as $O(mL+mc_{\mathrm{fit}})$ for their representation Field et al. 2014, Abstract and §IV, PDF.

This scaling is conditional: evaluation of frame rotations, mode sums, start-frequency extension, or Fourier transforms can dominate in another interface. A speedup should name hardware, implementation, waveform length, modes, accuracy setting, and baseline; an operation count alone is not a timing benchmark.

Four errors that should not be merged too early

For a validation waveform, define the relative time-domain norm

$$\epsilon[a] =\frac{\lVert a-h\rVert_2}{\lVert h\rVert_2}, \qquad \lVert h\rVert_2^2 =\int_{t_{\min}}^{t_{\max}}\lvert h(t)\rvert^2dt.$$

The projection error compares $h$ with its orthogonal projection onto the reduced space. It asks only whether the basis is expressive. The empirical-interpolation error reconstructs the waveform from its exact nodal values and includes node conditioning. The regression-only waveform error passes fitted and exact nodal values through the same empirical interpolant, then measures their reconstructed difference in the displayed waveform $L^2$ norm. The total surrogate error compares the fitted-node reconstruction with truth. All four reported errors therefore use the same denominator and time norm. They do not generally add in quadrature, and a small projection error can coexist with a poor parameter fit.

The fiducial error lies outside this hierarchy: if the training waveform differs from general relativity or from an accurate NR product, a surrogate can reproduce that error extremely efficiently. For hybrid training data, alignment and blend errors are inherited too. Field et al. distinguish greedy, interpolation, fitting, and complete-surrogate errors in Field et al. 2014, §§III.A–III.D and Figures 3–9, PDF.

To see why interpolation is a distinct operation, let $V_{ji}=e_i(T_j)$. Exact nodal values give coefficients $c=V^{-1}h(T_j)$, so the interpolant agrees with $h$ at the chosen nodes whenever $V$ is invertible. It reproduces every waveform in the reduced span exactly, but need not be the orthogonal projection of a waveform outside that span. A poorly conditioned $V$ amplifies errors in both exact and fitted nodal values. The triangle inequality gives

$$\epsilon_{\rm total}\le\epsilon_{\rm EIM}+\epsilon_{\rm regression},$$

using the common denominator above. Adding projection error once more would double count part of the EIM error rather than provide an independent contribution.

Pinned synthetic family and split

The benchmark uses one inertial-frame synthetic $(2,2)$ mode on 512 uniform samples $t\in[-1,1]$ and parameter $\lambda\in[-1,1]$:

$$h(t;\lambda) =A(t;\lambda)e^{-i\phi(t;\lambda)},$$ $$A=e^{-0.15(t-0.2)^2} \left[1+0.08\lambda t+0.03\lambda^2(1-t^2)\right],$$ $$\phi=7t+1.4\lambda t^2+0.35\lambda^2t^3 +0.08\lambda^3\sin(\pi t).$$

The baseline training grid has 33 equally spaced parameters. A fixed NumPy random seed 1729 selects 200 uniform held-out validation parameters; these points are not used to fit. A singular-value basis is truncated at rank $m$, empirical nodes are selected greedily from successive basis residuals, and degree-eight polynomials fit real and imaginary nodal data. The norm is the unweighted complex $L^2$ norm on the declared time interval, not a detector match.

In the executable benchmark, this norm is discretized with equal sample weights, $\lVert h\rVert_{\rm grid}^2=\Delta t\sum_j\lvert h_j\rvert^2$; the common $\Delta t$ cancels in each ratio. The SVD projection is orthogonal in that discrete norm. Table entries are maxima over the same 200 validation points, not a continuum-domain supremum, and no time-grid convergence is inferred merely from increasing basis rank.

The exact code was run with Python 3.9.6 and NumPy==2.0.2. Increasing basis complexity gives the required error hierarchy:

Basis rank	Projection error	EIM error	Regression-only waveform error	Total surrogate error
2	1.024259e-1	1.528075e-1	7.695581e-6	1.528064e-1
4	4.899394e-3	1.169924e-2	6.015496e-6	1.169971e-2
6	4.747071e-5	7.773693e-5	4.363972e-6	7.776104e-5
8	1.130405e-6	2.355175e-6	3.550956e-6	3.880782e-6

The rank-eight total no longer tracks projection error alone because regression is becoming relevant. That is a feature of the diagnostic, not a reason to report only the smallest column.

Sampling, fit, and extrapolation tests

Held-out results for controlled variations are:

Variation	Projection error	EIM error	Regression-only waveform error	Total error
17 training points	5.376819e-5	8.801412e-5	5.749182e-6	8.864736e-5
33-point baseline	4.747071e-5	7.773693e-5	4.363972e-6	7.776104e-5
65 training points	6.387302e-5	9.808489e-5	3.885611e-6	9.980491e-5
Degree-four fit	4.747071e-5	7.773693e-5	3.093455e-3	3.083251e-3
Rank 8, degree 10	1.130405e-6	2.355175e-6	1.073218e-7	2.363555e-6

The 65-point SVD is not monotonically better at fixed rank under this discrete training measure; the common held-out set catches that fact. The degree-four row isolates regression failure. At $\lambda=1.20$, outside the declared domain, the relative error is $1.769152\times10^{-3}$—about 22.75 times the baseline in-domain maximum. No basis tolerance certifies extrapolation.

For production families, vary parameter density and placement, basis rank or greedy tolerance, empirical-node count and conditioning, fit family and regularization, waveform resolution, and validation split. Varma et al. use training-point cross-validation and mismatches across mass, orientation, and parameter space for a hybrid NR surrogate Varma et al. 2019, §§VI–VII, PDF. Their reported performance is conditional on their hybrid family, 104 NR waveforms, parameter range, modes, alignments, mass scaling, one-sided detector PSDs, and validation design; it is not a generic guarantee for “surrogates.”

Executable rerun

import numpy as np

SEED = 1729
TIME = np.linspace(-1.0, 1.0, 512)
VALIDATION = np.random.default_rng(SEED).uniform(-1.0, 1.0, 200)

def waveform(parameter):
    p = np.asarray(parameter)[..., None]
    t = TIME[None, :]
    amplitude = (1 + 0.08*p*t + 0.03*p**2*(1-t**2))*np.exp(-0.15*(t-0.2)**2)
    phase = 7*t + 1.4*p*t**2 + 0.35*p**2*t**3 + 0.08*p**3*np.sin(np.pi*t)
    result = amplitude*np.exp(-1j*phase)
    return result[0] if np.ndim(parameter) == 0 else result

def eim_nodes(basis):
    nodes = [int(np.argmax(abs(basis[0])))]
    for k in range(1, len(basis)):
        matrix = basis[:k, nodes].T
        coefficients = np.linalg.solve(matrix, basis[k, nodes])
        residual = basis[k] - coefficients @ basis[:k]
        nodes.append(int(np.argmax(abs(residual))))
    return np.array(nodes)

def build(training_count=33, rank=6, fit_order=8):
    parameters = np.linspace(-1.0, 1.0, training_count)
    training = waveform(parameters)
    _, _, right_vectors = np.linalg.svd(training, full_matrices=False)
    basis = right_vectors[:rank]
    nodes = eim_nodes(basis)
    matrix = basis[:, nodes].T
    nodal_training = training[:, nodes]
    real_fits = [np.polyfit(parameters, nodal_training[:, j].real, fit_order)
                 for j in range(rank)]
    imag_fits = [np.polyfit(parameters, nodal_training[:, j].imag, fit_order)
                 for j in range(rank)]

    truth = waveform(VALIDATION)
    projection = (truth @ basis.conj().T) @ basis
    exact_nodes = truth[:, nodes]
    eim_coefficients = np.linalg.solve(matrix, exact_nodes.T).T
    interpolant = eim_coefficients @ basis
    predicted_nodes = np.column_stack([
        np.polyval(real_fits[j], VALIDATION) + 1j*np.polyval(imag_fits[j], VALIDATION)
        for j in range(rank)
    ])
    surrogate_coefficients = np.linalg.solve(matrix, predicted_nodes.T).T
    surrogate = surrogate_coefficients @ basis

    def maximum_error(approximation):
        numerator = np.sum(abs(approximation-truth)**2, axis=1)
        denominator = np.sum(abs(truth)**2, axis=1)
        return np.sqrt(np.max(numerator/denominator))

    regression = np.sqrt(np.max(
        np.sum(abs(surrogate-interpolant)**2, axis=1)
        / np.sum(abs(truth)**2, axis=1)
    ))
    return (maximum_error(projection), maximum_error(interpolant),
            regression, maximum_error(surrogate), nodes,
            basis, matrix, real_fits, imag_fits)

print("basis_rank,projection_error,eim_error,regression_waveform_error,total_error")
for rank in (2, 4, 6, 8):
    values = build(rank=rank)
    print(f"{rank},{values[0]:.6e},{values[1]:.6e},{values[2]:.6e},{values[3]:.6e}")

print("variation,projection_error,eim_error,regression_waveform_error,total_error")
for name, count, rank, order in (
    ("sparse_training", 17, 6, 8),
    ("baseline", 33, 6, 8),
    ("dense_training", 65, 6, 8),
    ("low_fit_order", 33, 6, 4),
    ("high_complexity", 33, 8, 10),
):
    values = build(count, rank, order)
    print(f"{name},{values[0]:.6e},{values[1]:.6e},{values[2]:.6e},{values[3]:.6e}")

projection, eim, regression, total, nodes, basis, matrix, real_fits, imag_fits = build()
outside = 1.20
truth = waveform(outside)
nodal = np.array([np.polyval(real_fits[j], outside) + 1j*np.polyval(imag_fits[j], outside)
                  for j in range(len(nodes))])
surrogate = np.linalg.solve(matrix, nodal) @ basis
outside_error = np.linalg.norm(surrogate-truth)/np.linalg.norm(truth)
print(f"outside_parameter,{outside:.2f},relative_error,{outside_error:.6e}")
print("empirical_nodes", ",".join(str(index) for index in nodes))

Save the Python program above as surrogate_benchmark.py. The program emits the rank, variation, extrapolation, and empirical-node results. The exact rerun is:

python3 -m venv .venv
. .venv/bin/activate
python3 -m pip install numpy==2.0.2
python3 surrogate_benchmark.py

Failure modes and open questions

Stop if validation points overlap training through leakage, basis or fit errors are reported only on the training set, node matrices become ill-conditioned, error concentrates near merger or parameter boundaries without targeted sampling, frame alignment changes with the split, or extrapolated values are returned without a domain flag. Also stop if the surrogate error is smaller than an undocumented training-waveform uncertainty: that comparison is uninterpretable, not impressive.

An online model should return its version, domain, modes, frame, and warnings with the waveform. A detector match, if later used, must state its one-sided PSD, frequency band, total-mass scaling, modes, orientation, and time/phase maximization. This chapter establishes construction and validation evidence only.

Exercise

Using the rank-six and rank-eight rows, identify the dominant listed error at each rank. Then compute the factor by which the out-of-domain error at $\lambda=1.20$ exceeds the rank-six baseline total.

Solution

At rank six, the empirical-interpolation error $7.773693\times10^{-5}$ is much larger than the regression-only waveform error $4.363972\times10^{-6}$ and nearly equals the total $7.776104\times10^{-5}$. Basis/interpolation is therefore the limiting listed stage. At rank eight, the EIM error is $2.355175\times10^{-6}$ while the comparable regression-only waveform error is $3.550956\times10^{-6}$; regression has become the larger listed stage and the total is $3.880782\times10^{-6}$. The total is not the quadrature sum because the reconstructed differences are correlated.

The extrapolation factor relative to the baseline is

$$\frac{1.769152\times10^{-3}}{7.776104\times10^{-5}} =22.75.$$

The in-domain validation maximum therefore provides no safe bound at $\lambda=1.20$, even though the extrapolation is only $0.20$ beyond the training interval.

References

• Field, Scott E., Chad R. Galley, Jan S. Hesthaven, Jason Kaye, and Manuel Tiglio. “Fast Prediction and Evaluation of Gravitational Waveforms Using Surrogate Models.” Physical Review X 4 (2014): 031006. DOI. Open PDF.
• Varma, Vijay, Scott E. Field, Mark A. Scheel, Jonathan Blackman, Lawrence E. Kidder, and Harald P. Pfeiffer. “Surrogate Model of Hybridized Numerical Relativity Binary Black Hole Waveforms.” Physical Review D 99 (2019): 064045. DOI. Open PDF.