Empirical Interpolation and Reduced Operators

A reduced basis is cheap only if its coefficients can be obtained without reconstructing every fine-grid field or evaluating every expensive nonlinear term. Empirical interpolation selects a small set of physical nodes and reconstructs the full reduced quantity from values there. The speedup is real when the offline construction, interpolation conditioning, parameter fits, and online cost are all measured separately. Exact agreement at selected nodes is a necessary algebraic check, not evidence of accuracy between them.

Required background. Probability and Inference for Gravity supplies sampling, held-out validation, and uncertainty assessment.

Helpful background. Formulations of General Relativity supplies the physical operators and constraints being approximated; Differential Equations for Relativity supplies residual and stability checks; Computational Relativity supplies convergence and reproducible-workflow controls.

Interpolation from selected nodes

Let $e_1(x),\ldots,e_m(x)$ span a reduced space for a field or nonlinear operator. Select distinct nodes $X_1,\ldots,X_m$ and define

$$V_{ji}=e_i(X_j).$$

For a new function $f$, solve

$$V c=f_X, \qquad (f_X)_j=f(X_j),$$

and reconstruct

$$(\mathcal I_m f)(x)=\sum_{i=1}^{m}c_i e_i(x).$$

If $V$ is nonsingular, $(\mathcal I_m f)(X_j)=f(X_j)$ exactly in exact arithmetic. A standard greedy construction chooses the first node at the largest component of $e_1$, subtracts the current interpolant from each new basis vector, and places the next node at the largest residual component. This creates a triangular construction and tries to keep each added direction visible.

The interpolation error obeys a bound of the form

$$\|f-\mathcal I_m f\| \le (1+\|\mathcal I_m\|) \inf_{v\in V_m}\|f-v\|,$$

Here the norm must make the sampling operator bounded: use, for example, the supremum norm of continuous functions on a compact domain, or a specified discrete vector norm. Point evaluation is not bounded on an arbitrary continuum $L^2$ space, where a function is defined only up to sets of measure zero. With that qualification, the operator norm is a Lebesgue-type amplification factor. The reduced-basis error can decay while $\|\mathcal I_m\|$ grows; then interpolation magnifies data, fit, or roundoff error. Never form $V^{-1}$ explicitly. Record the solve method, scaling, condition estimate, precision, and normwise backward residual

$$r_X= \frac{\|Vc-f_X\|} {\|V\|\,\|c\|+\|f_X\|}.$$

When the denominator vanishes, report the absolute residual instead and define the exactly zero case as $r_X=0$. This avoids a meaningless relative residual when the selected data vanish or nearly vanish.

The original empirical-interpolation waveform construction relates the interpolation error to greedy error and the amplification factor in Field et al. 2014, § III.B, Eq. (20) and Fig. 6. Discrete empirical interpolation applies the same logic to vector-valued nonlinear terms; its algorithm and error analysis are developed in Chaturantabut and Sorensen 2010, § 3, Algorithm 1 and Lemma 3.2.

Separating approximation from noisy node values

For any $v\in V_m$, reproduction gives $f-\mathcal I_mf=(f-v)-\mathcal I_m(f-v)$. The triangle inequality and minimization over $v$ yield the bound above. No assumption about a regression model enters this step.

Let $E$ map basis coefficients to the full discrete field, and suppose the predicted nodal values have error $\delta f_X$. Then

$$\widehat f-\mathcal I_mf=EV^{-1}\delta f_X, \qquad \|f-\widehat f\| \le\|f-\mathcal I_mf\|+\|EV^{-1}\|\,\|\delta f_X\|.$$

The inverse denotes the mathematical map; compute it by a solve. This identifies two separate tests: exact full-model node values test interpolation, while predicted node values test the regression and its amplification. A tiny nodal solve residual can coexist with a large nodal prediction error.

Offline–online reduction

Suppose a full discretization has state dimension $N$ and a reduced state $u\simeq U_n a$. A nonlinear term $F(U_na;\lambda)$ usually still costs $O(N)$. With a nonlinear basis $W_m$ and a selection matrix $P$ that extracts $m$ components, empirical interpolation gives

$$F(U_na;\lambda) \simeq W_m(P^TW_m)^{-1}P^TF(U_na;\lambda).$$

Let $M$ be the matrix of the discrete state-space inner product. Projecting onto the state basis yields

$$U_n^*MF \simeq \underbrace{U_n^*MW_m(P^TW_m)^{-1}}_{\text{offline}} \underbrace{P^TF(U_na;\lambda)}_{\text{online at }m\text{ entries}}.$$

Here $*$ is the conjugate transpose for complex coordinates and $M=I$ recovers the Euclidean formula; $P^T$ remains the row-selection operator. The precomputed matrix is useful only when the selected components of $F$ can be evaluated without constructing the entire $N$-vector. Hidden global normalization, elliptic solves, or communication can destroy the promised online scaling. For nonpolynomial equations of state or fluxes, node evaluation must also retain physical-domain checks such as positive density and subluminal velocity.

Waveform surrogates use a closely related factorization. A reduced basis compresses the time dependence, empirical nodes choose a small set of times, parameter fits predict waveform values at those times, and interpolation reconstructs the series. Field and collaborators spell out the basis, node, fit, and online stages in Field et al. 2014, §§ III.A–III.D and Eqs. (20), (31)–(33). This exposes three distinct errors: basis projection, empirical interpolation, and parameter fitting.

Compact-binary operator benchmark

Freeze the same versioned compact-binary waveform family, parameter domain, mode convention, time interval, alignment, norm, and training/boundary/hole/out-of-domain split used for reduced bases. Store dataset and preprocessing checksums. Build a waveform basis $e_i(t)$, select empirical times $T_j$, and fit the complex waveform value at each $T_j$ over parameter space using a declared regression model.

For each $m$:

1. verify interpolation of every basis vector and training waveform at all selected nodes, with $r_X$ below $10^{-12}$ in binary64 unless conditioning predicts a larger floor;
2. report $\kappa(V)$, the basis greedy error, the actual empirical-interpolation error, and the computable amplification bound;
3. evaluate projection-only, interpolation-only, fit-only, and total errors on scattered, boundary, and contiguous-hole holdouts;
4. report time-domain relative error and detector-weighted mismatch against each held-out high-fidelity waveform;
5. compare direct full reconstruction time with online surrogate time after warm-up, while reporting offline basis, node, and fit costs separately.

A representative finite acceptance target is maximum in-domain held-out mismatch below $10^{-4}$, node residual below its conditioning-aware tolerance, and no unexplained increase in the worst boundary error as $m$ grows. The actual discrepancies and worst parameter locations must accompany the thresholds. Out-of-domain points test refusal or warning behavior; success there does not enlarge the validated domain.

Increase basis size, node count, training density, fit complexity, optimization tolerance, and arithmetic precision one at a time. When projection error continues to fall but total error plateaus, the parameter fit or interpolation matrix has become limiting. When the bound grows while actual error remains small, it is conservative but still diagnoses sensitivity. When both grow, adding nodes can make the model worse.

The environment record pins waveform revisions, split manifest, preprocessing, basis orientation, node ordering, fit architecture and priors, optimizer, seeds, linear-algebra and regression libraries, precision, compiler or interpreter, operating system or container, accelerator, thread counts, and timing method. Independently recompute selected overlaps with direct quadrature and reconstruct a subset from full-model evaluations that never enter fitting.

Conditioning, cost, and failure modes

The central failure modes are nearly dependent node rows, clustered nodes that miss a sharp feature, a parameter fit whose error is mistaken for interpolation error, data leakage through alignment or hyperparameter selection, and an online routine that secretly evaluates an $N$-dimensional object. Quoting speedup without hardware, batching, warm-up, sample rate, accuracy target, and amortized offline cost is also a failed comparison.

An empirical operator can preserve selected components yet violate a conservation law or constraint globally. Evaluate the original full residual and physical invariants on held-out reconstructed states. A small reduced residual computed with the same approximate operator is not independent evidence.

The complete rerun verifies inputs and environment; rebuilds the basis and ordered nodes; sweeps $m$, training coverage, fit tolerance, and precision; emits interpolation identities, conditioning, bounds, held-out errors, physical residuals, timing, and offline cost; and repeats a clean holdout comparison through an independent full-model path. Any node identity, conditioning, boundary, physical-residual, or total-error failure rejects the advertised reduced operator even if its online evaluation is fast.

Check your understanding

Take basis functions $e_1(x)=1$ and $e_2(x)=x$ on $[-1,1]$. The first empirical node is chosen at $x=-1$. After interpolating $e_2$ from $e_1$, where should the second node be placed, and what functions are then reproduced exactly?

Solution

Matching $e_2$ at $x=-1$ with a constant gives $-1$, so the residual is $x+1$, whose largest magnitude occurs at $x=1$. The interpolation matrix is

$$V=\begin{pmatrix}1&-1\\1&1\end{pmatrix},$$

which is nonsingular. Every affine function $a+bx$ is then reproduced exactly at all $x$, not merely at the two nodes. Choosing the same node twice would make the rows identical and the reconstruction undefined.

Continue with Gaussian-Process Emulators.

References

• Chaturantabut, Saifon, and Danny C. Sorensen. “Nonlinear Model Reduction via Discrete Empirical Interpolation.” SIAM Journal on Scientific Computing 32, no. 5 (2010): 2737–2764. DOI.
• 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. arXiv:1308.3565. DOI.