Fast prediction and evaluation of gravitational waveforms using surrogate models

[Abridged] We propose a solution to the problem of quickly and accurately predicting gravitational waveforms within any given physical model. The method is relevant for both real-time applications and in more traditional scenarios where the generation of waveforms using standard methods can be prohibitively expensive. Our approach is based on three offline steps resulting in an accurate reduced-order model that can be used as a surrogate for the true/fiducial waveform family. First, a set of m parameter values is determined using a greedy algorithm from which a reduced basis representation is constructed. Second, these m parameters induce the selection of m time values for interpolating a waveform time series using an empirical interpolant. Third, a fit in the parameter dimension is performed for the waveform's value at each of these m times. The cost of predicting L waveform time samples for a generic parameter choice is of order m L + m c_f online operations where c_f denotes the fitting function operation count and, typically, m<<L. We generate accurate surrogate models for Effective One Body (EOB) waveforms of non-spinning binary black hole coalescences with durations as long as 10^5 M, mass ratios from 1 to 10, and for multiple harmonic modes. We find that these surrogates are three orders of magnitude faster to evaluate as compared to the cost of generating EOB waveforms in standard ways. Surrogate model building for other waveform models follow the same steps and have the same low online scaling cost. For expensive numerical simulations of binary black hole coalescences we thus anticipate large speedups in generating new waveforms with a surrogate. As waveform generation is one of the dominant costs in parameter estimation algorithms and parameter space exploration, surrogate models offer a new and practical way to dramatically accelerate such studies without impacting accuracy.


I. INTRODUCTION
A direct detection of gravitational waves generated by the coalescence of a compact binary system is among the most anticipated discoveries to be made in gravitational wave physics. The signal from such an event will codify perhaps the only attainable information about the existence, dynamics, and underlying physics of the strongest gravitating objects in the universe. Currently, there are few, if any, direct observations pertaining to gravity in the strong field regime but there is enough data to show agreement with the predictions of general relativity when gravitational fields and speeds are not too large [1,2].
In the case of binary black holes, where the fields and speeds can be large, one must rely on numerical simulations of the Einstein equations to discover how these systems evolve. The resulting solution depends on the choice of initial data. The intrinsic parameter space of binary black holes in quasi-circular orbit is seven dimensional, consisting of the mass ratio and the three spin angular momentum components for each black hole. Different choices of parameters can lead to qualitatively different outcomes, such as the final speed of the merged black hole due to a "kick" from the asymmetric emission of gravitational waves [3][4][5][6][7][8][9][10][11][12][13][14]. In addition, potentially interesting effects due to strong precession from highly spinning black holes are waiting to be discovered and understood. Unfortunately, each numerical relativity (NR) simulation typically involves the use of large scale supercomputers, making an exploration of the parameter space a currently computationally intractable problem. For example, one might employ a uniform or random sampling strategy of the parameter space that, for a mere 4 points per dimension, requires 4 7 ≈ 16,000 expensive numerical solutions of binary black hole coalescences. This number, while still being a very coarse survey of the parameter space, is substantially greater (by more than an order of magnitude) than all the simulations performed by all of the numerical relativity groups to date [15][16][17][18].
To help alleviate this computational bottleneck, models of the inspiral, merger, and ringdown phases of a binary black hole (BBH) coalescence have been developed over the last decade [19][20][21][22][23][24][25][26][27][28][29]. The purpose of these phenomenological models is to provide a sufficiently accurate representation of a BBH waveform within some range of parameters by fitting certain coefficients and functions to a set of waveforms extracted from numerical simulations. In doing so, the models help to reduce the amount of information needed to represent NR waveforms. While these models are significantly faster than solving the Einstein field equations they remain computational bottlenecks for parameter estimation studies, arXiv:1308.3565v2 [gr-qc] 28 Feb 2014 which typically require generating millions of waveforms on the fly. Additionally, they still rely on waveforms computed from numerical simulations of binary black hole mergers and are thus unable, at least currently, to accurately model gravitational waveforms throughout the entire seven-dimensional parameter space, although efforts to attack this problem are underway [18,30].
Other important considerations come from precessing inspirals of compact binaries such as binary neutron stars. Generating the corresponding waveforms requires solving a set of ordinary differential equations (ODEs) and substituting the solutions into the post-Newtonian expressions for the phase and amplitude corrections. Given that around 520,000 to 860,000 waveforms are needed to build template banks for just non-precessing, slowly spinning binary neutron stars for advanced LIGO [31], which would already be a computational challenge, it follows that the large number of ODE solves would be prohibitively expensive in the general precessing case. Waveform generation for precessing compact binary inspirals constitutes the main computational bottleneck for both template bank construction and parameter estimation studies.
In this paper we offer a solution to the need for cheap and accurate generation of gravitational waveforms, that may otherwise be too expensive to compute for the application of interest. Alternative waveform prediction methods have recently been proposed [32][33][34][35] (see Appendix F for a brief discussion and Ref. [35] for comparison details, especially). These works have focused on gravitational waveform models known through closed-form expressions while the focus of this paper is on those described by differential equations. To achieve this, it is crucial to take advantage of the rich structure underlying the waveforms of interest. Importantly, our method builds accurate surrogate models that do not sacrifice the underlying physics but instead combines the efficiency and power of reduced order modeling techniques with high accuracy sparse representations and an offline-online decomposition of the problem.
Work over the last few years has shown that gravitational waveforms exhibit redundancy in the parameter space [36][37][38][39][40][41], suggesting that the amount of information necessary to represent a fiducial waveform model is smaller than might be anticipated. This reduction can be captured accurately using only a remarkably few number m of representative waveforms. These m representative waveforms can be found using a greedy algorithm and comprise a reduced basis [38] from which all other waveforms within the same physical model can be represented provided one can compute their projections onto the basis. In practice, this is neither feasible nor worthwhile because projecting onto the basis requires already knowing the waveform that one is seeking to represent in the first place. This is particularly the case for waveform families that are expensive to generate, such as those from numerical relativity (NR) simulations of the full Einstein equations. Instead, we aim to use only the information provided by the m representative waveforms of the reduced basis to predict waveforms accurately and cheaply for any desired parameter values.
To accomplish this goal we first build the reduced basis as mentioned above and described in more detail in Sec. III A and Appendix A. Second, we construct a temporal interpolant [42] whereby any fiducial waveform is fully specified through its evaluation at m appropriately chosen times. While this may seem surprising it is important to recall that there are only m independent pieces of information in the waveform family as indicated by the m waveforms that comprise the reduced basis. Indeed, we will show that the m reduced basis uniquely specify these m specially chosen times. The interpolation method outlined above, which is called empirical interpolation because it generates an interpolant specific to the given fiducial waveform family, takes advantage of this nearly optimal representation strategy in parameters to provide a corresponding representation strategy in time. See Sec. III B and Appendix B for more details. Finally, at each empirical interpolation time we perform a fit in the parameter dimension of the waveform's amplitude and phase. Evaluating these fits yield m time samples from which the waveform is accurately recovered through its empirical interpolant representation. Remarkably, the outlined method allows for a waveform within any physical model to be predicted for any parameter value of interest based solely on a knowledge of m fiducial waveforms.
Combining these pieces of information yields a surrogate model for the fiducial waveform family. The method to build the surrogate has several useful properties. First, the method is entirely hierarchical, i.e. the accuracy of the surrogate model can be improved, if necessary, by adding fiducial waveforms without discarding any of the previous ones. Second, the surrogate model can be evaluated using only O (mL + mc fit ) computational operations, where L is the number of time samples at which the model is evaluated and c fit is the typical fitting function operation count. This provides a significant speedup compared to the usual way that fiducial waveforms are generated, as we demonstrate below with a surrogate model for non-spinning Effective One Body (EOB) waveforms. The speedup compared to numerical simulations of the full Einstein equations is expected to be significantly larger.

II. SURROGATE WAVEFORM MODELS
We denote the gravitational waveform produced from a fiducial model by h(t; λ). Here, t denotes time and λ is the waveform parameterization (e.g., mass ratio and spins). We denote the surrogate model of the fiducial waveform family by h S (t; λ) and describe its construction in this section.
When numerically generating waveforms, by solving partial or ordinary differential equations, one typically solves an initial (or initial-boundary) value problem for a fixed λ i thereby generating h(t; λ i ) on a densely sampled grid in time. In this paper we develop a procedure for building h S (t; λ) through judicious choices of λ i and the corresponding output h(t; λ i ) found by solving the relevant equations defining the fiducial problem. Crucially, given the complexity of existing numerical solvers, our approach to surrogate modeling is intentionally nonintrusive to legacy codes.
We seek a minimal number of λ i selections for a target accuracy such that the surrogate has a comparable or smaller error than that associated with the underlying waveform model. This is important both for the speed of evaluating the surrogate model and for overcoming computational challenges with building it in cases where one cannot generate h(t; λ i ) for arbitrarily many values of λ. Naturally, if more data is available it should be possible to include it and improve the surrogate's quality. This means that the surrogate model should be hierarchical by construction, improving as more simulations become available and without discarding previous ones.
The algorithm for building and evaluating a surrogate for a given fiducial family or model of gravitational waveforms is schematically depicted in Fig. 1 and outlined below: 1. (Offline) Described in Section III A. Select the most relevant m points in parameter space (shown as red dots in Fig. 1). The waveforms associated with these selections (shown as red lines) provide a nearly optimal reduced basis (RB) for this waveform family [38]. The resulting points and waveforms will be referred to as greedy data.

2.
(Offline) Described in Section III B. Identify m time samples of the full time series, which we call empirical nodes or times, to build an interpolant that accurately reconstructs any fiducial waveform. This step, called the Empirical Interpolation Method (EIM), only requires knowing the reduced basis. The number of empirical nodes m (shown as blue dots on the vertical axis in Fig. 1) exactly equals the number of basis elements m.
3. (Offline) Described in Section III C. At each empirical node perform a fit (e.g., least squares) in the parameter dimension for the amplitude and phase of the waveform using the greedy data from Step 1. The fits are indicated by blue lines in Fig. 1 We quantify the accuracy of the offline steps through the convergence rates in (9) and (20). The accuracy of the fast online step for the complete surrogate is estimated through the errors in (31) and (32). If each offline step is carried out with sufficiently good accuracy then the surrogate will satisfy for all t and λ in the given ranges and retain the physics of the original fiducial waveform family, whatever that might be. As discussed in Sec. IV, the waveform predictions by our surrogate model are indeed expected to have a small error with respect to the fiducial one.

III. SURROGATE MODEL BUILDING
The following four subsections expand on the steps outlined above. Each of these steps is illustrated with an application to non-spinning EOB waveforms. For simplicity we consider the (2, 2) mode of waveforms with mass ratios in the range q ∈ [1,2] and about 12,000M in duration. In Sec. VI we build surrogate models for astrophysical sources that include more cycles, cover larger mass ratio intervals, and contain higher spherical harmonics. Important technical details describing how these EOB waveforms were generated as well as our peak alignment scheme are discussed in Appendix E. Figure 2 shows the q = 1 EOB waveform. Despite its complicated structure, we shall demonstrate that waveforms such as this one can be represented accurately by relatively little information.
A gravitational waveform h(t; λ) is represented in terms of its two fundamental polarizations h + (t; λ) and h × (t; λ) by h(t; λ) = h + (t; λ) + ih × (t; λ). A natural inner product is given by the complex scalar product with an inherited norm given by h(·; λ) 2 = h(·; λ), h(·; λ) . Here, h * (t; λ) is the complex conjugate of h(t; λ). Other inner products might be more natural for different applications [43]. Throughout this paper we shall assume the waveforms are normalized such that h(·; λ) = 1. The overlap integral of two normalized waveforms, say, of a fiducial waveform and its surrogate model prediction, is given by Re h(·; λ), h S (·; λ) , This equality is useful to translate the error in approximating a fiducial waveform by its surrogate model prediction into an overlap integral that is used in some gravitational wave applications (c.f., Eq. (11)). A.
Step 1: Greedy selection of parameter samples and reduced basis We use a greedy algorithm (see Appendix A for more details) to select m parameter points {Λ i } m i=1 and corresponding waveforms h i (t) = h(t; Λ i ). The greedy algorithm provides a nearly optimal solution to the Kolmogorov n-width approximation problem [44,45], namely, given a set of waveforms where T denotes a compact parameter domain, find an m-dimensional function space that best approximates any h(t; λ) from this set.
and an associated set of waveforms that constitutes the reduced basis. The basis is hierarchical in the sense that if {h i } m i=1 is the basis for m < m then One of the key features of the greedy algorithm is its ability to select a small number of waveforms to serve as an accurate basis. For practical purposes of conditioning it is useful to use an orthonormal basis {e i } m i=1 , which spans the same approximation space as (6).
With the RB in hand, every waveform in the training set is well approximated by an expansion of the form whereas waveforms from T (even if not in the training set) continue to be well approximated by the RB if the training set is dense enough [38][39][40][41]. Since the waveform space is numerically finite dimensional [38], one can verify sufficiently dense training sets through convergence as M gets larger or by checking how well the basis represents randomly selected waveforms (see Appendix A). For an underlying model that requires prohibitively expensive numerical solves, one may use a simpler model to propose a training set building strategy. If there is sufficient similarity amongst the members of the original set then m M . This is found to be the case for gravitational waveforms [38,40,41].
Let be a user specified tolerance whose role is to guarantee that the approximation error for waveforms in the training set, which we will call the greedy error σ m , is bounded by , Then, the representation (8) is accurate to . The minimization over the coefficients {c i } in (9) is achieved by orthogonal projection P m h(t; λ) of h(t; λ) onto the span of the basis (see Appendix A for details) so that In Sec. III B we will find efficient approximations of the optimal projection representation in (10) that approximately retains its accuracy implied by (9). The error in (9) is directly related to the overlap between a waveform and its representation [46] which follows from (3). The quantity σ m quantifies the worst error of the best approximation by the basis. The greedy algorithm is nearly optimal in the sense that if the Kolmogorov n-width d m (defined as the smallest error (9) achieved by a best m-dimensional function space) decays exponentially then so does the greedy error [44,45], where D, a, b are positive constants andã = 2 −1−2b a. Recent work [38][39][40][41] has shown that for fixed but arbitrary physical and parameter ranges, a small number of basis functions is indeed sufficient to accurately represent any waveform of the same physical model and with an exponentially decaying greedy error (9). Such observations are expected for functions with smooth parameter dependence, as is the case with gravitational waveforms. To better understand these approximation properties one can make an analogy to the more familiar case of spectral methods. There, exponential decay with the number of basis elements is expected whenever there is smoothness with the physical dimension(s) (e.g., space or time).
Let us apply the greedy algorithm to build a reduced basis for our nominal EOB example introduced earlier. Figure 3 shows the exponential decay of the greedy error (9) over 501 waveforms in the training set, with only 19 RB waveforms needed to represent the EOB model to machine precision for the mass ratios considered. Errors of about 10 −3 are already achieved with as few as 5 RB waveforms. Later on (c.f., Fig. 9), we show that any waveform not present in the training set yields similarly small representation errors by the basis. This feature, due to a sufficiently well sampled training set, is essential for parameter estimation studies, which seek to explore the waveform continuum. The distribution of selected points is shown in Fig. 4. In Sec. III C we show how the greedy data from these parameter selections can be used to predict waveforms for any q in the range considered, including (and especially) values not in the original training set.   (9), over 501 EOB training set waveforms with mass ratios between 1 and 2. Labels at the dots indicate the selected mass ratios at each step in the greedy algorithm.

B. Step 2: Greedy selection of time samples and empirical interpolation
Once a basis is built in Step 1 we can express any waveform evaluated at any time as a sum of m reduced basis elements. In Step 2, which is shown to significantly reduce the surrogate's evaluation cost in Appendix F, we now show how to leverage this knowledge to yield a temporal prediction scheme by recasting the problem as one of interpolation in time. Given a reduced basis {e i } m i=1 and m evaluations of a fiducial waveform at certain times {T i } m i=1 , we wish to recover the full fiducial waveform h(t; λ) with high accuracy for an arbitrary λ. A proper choice of these times {T i } m i=1 is crucial. Naively selected times, such as those randomly or equally spaced, do not guarantee that: i) the interpolation problem is well-conditioned or even has a solution, and ii) the interpolation error is minimized with a nearly optimal convergence rate.
A framework for finding a "good" set of times that achieve both criteria is provided by the Empirical Interpolation Method (EIM) [47,48]. These special times, which we call empirical times or nodes, are selected as a (sparse) subset of the waveform's given time series (or even the continuum). The empirical nodes are uniquely defined by the reduced basis waveforms and only these waveforms. Like the algorithm for building a reduced basis, the EIM is hierarchical and uses a greedy optimization strategy to select the most representative times. While the empirical times T i do not explicitly depend on parameters or their ranges, the parameter dependence is implicit, nevertheless, through the basis. For example, a reduced basis for spinning or precessing waveforms will exhibit different features, and the distribution of T i will reflect this. For the moment we shall assume that the empirical nodes are known; the precise algorithm for finding them is given in Appendix B. The empirical interpolant, which interpolates the waveform h(t; λ) in time for a given parameter λ, is denoted by I m [h](t; λ) and takes the form The coefficients {C i } m i=1 are defined by requiring the interpolant to equal the value of the waveform at the empirical nodes, which is equivalent to solving an m-by-m system is independent of the parameters λ.
The choice of empirical nodes given by the EIM algorithm together with the linear independence of the reduced basis ensure that V in (16) is as well-conditioned as possible and invertible [49] so that is the unique solution to (14). It then follows upon substituting (17) into (13) that the empirical interpolant is where and is independent of λ. Note that (18) is a linear combination of the fiducial waveform itself evaluated at the empirical times. The coefficients {B i } m i=1 are built directly from the reduced basis and provide a clean offline/online separation. Because of this the {B i } m i=1 can be pre-computed offline once the reduced basis is generated while the (fast) interpolation is computed during the online stage from (18) when the parameter λ is specified by the user. Evaluations of the fiducial waveform are still needed at the arbitrarily chosen parameter λ in order to construct the interpolant in (18). In the next subsection we explain how to estimate the fiducial waveform at any λ, thus approximating {h(T i ; λ)} m i=1 and completing the construction of the surrogate model.
The empirical interpolant satisfies [50] where σ m is the greedy error defined in (9) and Λ m is a computable Lebesgue-like quantity that changes slowly with m (see Appendix B). For problems with smooth dependence with respect to parameter variations we can expect an exponential decay of σ m with m and of the empirical interpolant's error.
Before describing how to estimate the values {h(T i ; λ)} m i=1 for arbitrary λ let us assume these values are known exactly and apply the EIM to build an empirical interpolant for our fiducial EOB example introduced earlier. Figure 5 shows all 19 empirical nodes set against a q = 1 waveform to compare with the structure of a typical waveform. Evaluating any q ∈ [1, 2] EOB waveform at these 19 nodes and computing (18) one can reconstruct the full time-series of the waveform with high accuracy. This is explicitly demonstrated in Fig. 6 where the solid black line denotes the largest empirical interpolation er- as a function of the number of reduced basis elements/empirical nodes for 1,000 randomly selected EOB waveforms drawn from q ∈ [1,2]. Notice that this error is remarkably close to the greedy error (dashed line) in (9) when using (10) for the coefficients. The bound in (20) (dashed-dotted line) guarantees an error better than 10 −8 , which is sufficient for many GW applications.  The solid line shows the maximum empirical interpolant error (21) taken over 1,000 randomly selected waveforms (i.e., not taken from the training set) for q ∈ [1,2]. The dash-dotted line shows the error bound provided by the right side of (20) and is based solely on the greedy error and Λm. All three errors display similar decay rates.

C. Step 3: Fitting at empirical nodes
The next step is to predict waveforms at the empirical nodes {T i } m i=1 for arbitrary parameter values λ based only on the knowledge of the fiducial waveforms at the greedy points {Λ i } m i=1 . To accomplish this, we fit h(T i ; λ) with respect to λ at each T i using only the following m values of the reduced basis waveforms: The accuracy of the fit using only this data relies, at least partially, on the fact that the reduced basis waveforms are chosen to be the most dissimilar from one another. Of equal importance is our choice of fitting function which, in principle, is arbitrary. We will focus on the choices most effective for our nominal EOB example while others could be more appropriate for different waveform families.
The behavior of most astrophysically relevant gravitational waveforms is highly oscillatory in time but the phase and amplitude themselves have a relatively simple structure. It is thus easier to perform high-accuracy fits of the phase and amplitude than of the complex waveform itself. The amplitude A and phase φ are defined through This third step then consists of finding 2m functions, , approximating the amplitude and phase of the waveform. Once these fitting functions have been found the approximation at each T i is Depending on the application some fitting functions might be more useful than others. Therefore, this third step in constructing a surrogate model is flexible in the way that the fitting is implemented and thus in how the surrogate is ultimately generated. This is quite a useful feature of the method that may be especially beneficial for building surrogate waveforms for highly precessing black hole binaries. Splines, rational polynomial, or weighted non-oscillatory fitting approaches could help limit the impact of numerical noise, for example.
We now return to our nominal EOB example and perform a least squares fit for both the amplitude and phase as a function of mass ratio at each empirical time using polynomials, where α i , β i < m are the degrees of the polynomials at the empirical time T i for i = 1, 2, . . . , m. Further details regarding how to select an optimal degree are provided in Appendix C.
Top: Amplitude (solid) and phase (dashed) of the fiducial EOB training space waveform at the fifteenth selected empirical time as a function of q along with the greedy data (circles). This empirical time is T15 = 28.5M after merger and corresponds to the largest pointwise relative error for the least squares fit to the amplitude as quantified by (26). Bottom: The pointwise least squares errors for the amplitude (red) and phase (blue) at T15 evaluated for 1,000 randomly selected waveforms. The dashed lines correspond to the maximum pointwise error for the second empirical node T2 = −2,367M , which has the smallest maximum error of all the nodes.
The top plot in Fig. 7 shows the amplitude and phase, along with the greedy data points, at the fifteenth empirical time node, T 15 , which is about 28.5M after merger. This node corresponds to the largest pointwise error for the relative amplitude for waveforms in the training set of our EOB test problem. T 15 also happens to correspond to the second largest difference for the phase, The bottom plot in Fig. 7 shows the pointwise errors (solid lines) of (26) and (27) as a function of mass ratio for 1,000 randomly selected waveforms. These errors are uniformly below 3 × 10 −3 . The horizontal dashed lines show the maximum errors for the empirical node for which (26) and (27) are smallest, which occurs for the second empirical time T 2 = −2,367M . These errors are of order 10 −5 . As we will discuss later on (see Fig. 9) all of this information translates into a mismatch of the surrogate model with respect to the underlying EOB family of < 10 −7 . The quality of a fit at each empirical node, using the greedy data, depends on the smoothness of those waveforms with respect to parameter variation. This is discussed in Appendix E. Here, it suffices to mention that the fitting errors depend sensitively on accurately aligning the waveforms at their peaks, which affects the fits most noticeably through merger and ringdown. This can be seen in the top panel of Fig. 8. Figure 8 shows the maximum of the pointwise differences from (26) and (27) for the relative amplitude (circles) and phases (crosses), respectively, evaluated at each empirical time. We see that the amplitudes are accurate to better than 10 −5 for the entire inspiral phase until the merger regime where the error increases to about 10 −3 after which it plateaus throughout the ringdown stage. The phase errors increase modestly during the inspiral and likewise plateau through ringdown with errors at the level of 10 −3 .  (26) and (27), maximized over the greedy mass ratios at each empirical time for our EOB example. The top panel shows these errors when using a polynomial least squares fit and the bottom panel when using a fitting function inspired by the post-Newtonian amplitude and phase. Both types of fits exhibit very low errors at all of the empirical times.
Instead of using polynomials for the fitting functions we next consider functions inspired by the expressions for the amplitude and phase through leading order and nextto-leading order, respectively, in the post-Newtonian expansion φ i (q) = a i,0 (q − 1) ai,1 q ai,2 1 + a i,4 (q + a i,5 ) ai,6 + a i,3 .
The bottom panel of Fig. 8 shows the maximum of the pointwise differences from (26) and (27) using these post-Newtonian-inspired fitting functions. These fitting functions have a least squares fitting error comparable to the polynomial errors shown in the top panel. In both cases, the fit quality decreases rapidly at the merger but still exhibit very low errors at all of the empirical times. We thus see in this example that the third offline step for building the surrogate is flexible in the choice of fitting functions. This insight could be useful for other fiducial models such as waveforms with precession. D.
Step 4: Completing the surrogate model Finally, our complete surrogate model h S (t; λ) for the fiducial waveform family is given by substituting the fitting approximation (24) into the empirical interpolant (18), which yields This is the culmination of the offline steps. Only the m reduced basis waveforms evaluated at the m empirical times are needed to build the surrogate model and to predict an approximation for a fiducial waveform at any time and parameter value. In addition, the are computed once and for all offline; only the fitting functions for the amplitude and phase need to be evaluated during the online stage once λ is specified.

IV. ASSESSING THE SURROGATE MODEL
One of the errors of interest for the complete surrogate model is a discrete version of the normed difference between a fiducial waveform and its surrogate, which is, for L equally spaced time samples, where ∆t = (t max − t min )/(L − 1). We will sometimes refer to this as the surrogate error. Recall, from (3) and (11) that the square of the normed difference between two waveforms is directly related to their overlap. Other errors of interest are the pointwise ones for the phase and amplitude, Figure 9 shows a variety of comparisons between the surrogate and fiducial model for our EOB test case, using L = 16,384 time samples [51]. The top plot shows that the surrogate error (31) is uniformly below 10 −7 , where the mass ratio q = 1.068 corresponds to the largest error. The middle panel of Fig. 9 shows the fiducial EOB and surrogate waveforms for q = 1.068. Both waveforms are visually indistinguishable and, from the bottom panel of the same figure, we see that both amplitude and phase pointwise errors (32) are indeed very small. The largest errors are 10 −3 and are smaller than: i) the differences, for the same quantities, between the EOB model and the NR simulations used to calibrate the former [26], and ii) the numerical error of those NR simulations (see, e.g., [52]) and of more recent state-of-the-art simulations [18], as quantified through self-convergence tests. As discussed in Sec. III C and App. E, these maximum errors for the surrogate take place shortly after merger and are directly related to the accuracy with which one can determine the peak amplitude of the fiducial waveforms used to build the surrogate.  (31), which is related to the overlap error through (3), for 1,000 randomly selected mass ratios. The mass ratio yielding the largest surrogate model error is q = 1.068. Middle: The fiducial EOB waveform and its surrogate prediction for q = 1.068. There is visual agreement throughout the entire duration of ≈ 12,000M . Bottom: The fractional errors (32) in the amplitude and the phase difference between the fiducial EOB waveform and its surrogate model prediction for q = 1.068. The differences are smaller than the errors intrinsic to the EOB model itself as well as those of state-of-the-art numerical relativity simulations.
In Appendix D we derive the following error bound for the discrete norm (31), This bound identifies contributions from two sources. The first term in (33) describes how well the empirical interpolant (i.e., the basis and empirical nodes) represents h(t; λ). The expected exponential decay of the greedy error σ m with m along with a slowly growing Lebesgue constant Λ m results in this term being very small. The term Λ m σ m corresponds exactly to the curve labeled "EIM Bound" in Fig. 6. The second term in (33) is related to the quality of the fit. Incidentally, the fitting step has the dominant source of error in the surrogate model compared to the first two steps of generating the reduced basis and build the empirical interpolant (see also the discussion in Section III C).

V. COST AND SPEEDUP FOR SURROGATE MODEL PREDICTIONS
Next we discuss the cost (in terms of operation counts) to evaluate a surrogate model. We also present the large speedups that can be achieved when evaluating a surrogate model for our nominal EOB example compared to generating a fiducial waveform using the EOB solver as implemented in the LAL software package, which we refer to as the EOB-LAL code.
The complete surrogate model is given in (30) where the m coefficients B i (t) in (19) and the 2m fitting functions {A i (λ)} m i=1 and {φ i (λ)} m i=1 are assembled offline as described in Sections III A, III B, and III C. In order to evaluate the surrogate model for some parameter λ 0 we only need to evaluate each of those 2m fitting functions at λ 0 , recover the m complex values {A i (λ 0 )e −iφi(λ0) } m i=1 , and finally perform the summation in (30). Each B i (t) is a complex-valued time series with L samples. Therefore, the overall operation count to evaluate the surrogate model at each λ 0 is (2m − 1)L plus the cost to evaluate the fitting functions. Figure 10 shows timing results for the nominal EOB test case with m = 10 and a surrogate error (31) uniformly below 10 −7 for all mass ratios between 1 and 2. The top panel confirms that the cost of evaluating the surrogate model is linear in the number of samples L, as discussed above.
Depending on the sampling rate, the speedup in evaluating the surrogate model compared to generating an EOB waveform with the EOB-LAL code is between two and almost four orders of magnitude. For a sampling rate of 2 11 = 2,048 Hz, which is the rate used in the S5 and S6 searches for gravitational waves from binary black holes by the LIGO-VIRGO-GEO600 collaboration Average time to generate a single fiducial EOB waveform from a standard EOB code (circles) and through evaluation of its surrogate (crosses). Here we show results for the nominal example when using polynomial least squares fits for the amplitudes and phases. Bottom: The speedup, defined as the ratio of waveform generation times for EOB-LAL code to the surrogate model. [27,53], the speedup is ≈ 2,300 as shown in the bottom panel of Fig. 10. This is about three orders of magnitude faster than the EOB-LAL code.
The speedups indicated here are not an artifact of studying waveforms from binaries with nearly equal masses. Repeating these experiments for waveforms with mass ratios from 9 to 10 (chosen so that the typical duration ≈ 11,000M and number of waveform cycles ≈ 80 are comparable to our nominal EOB example), we find that only m = 15 reduced basis waveforms are needed to span the space with σ m = 10 −11 . The resulting surrogate model has an error from (31) of 8 × 10 −9 with a corresponding speedup in the online stage of about 5,000 at a sampling rate of 2,048 Hz. Again, the speedup is about three orders of magnitude.
As already mentioned in Sec. IV, the fitting step for building the surrogate potentially introduces the largest errors in the surrogate model. For the EOB example, these largest errors are still small (see Fig. 9) and suggest that one does not need to include all 19 basis waveforms/empirical times in order to yield a sufficiently accurate approximation. The top panel of Fig. 11 shows the surrogate error in (31), maximized over 1,000 randomly selected waveforms, as a function of the number of selected RB waveforms m. After m = 7 there is little to be gained by including more basis waveforms because the surrogate error is roughly constant until m = 19 while, from the bottom panel of Fig. 11, its evaluation time continues to grow with m. The dash-dotted line in the top panel shows the expected error computed by averaging the surrogate's error bound (33) over q. Taking the average (maximum) of (33) over q we are guaranteed surrogate errors of better than 10 −5 (5 × 10 −5 ), which is sufficient for many GW applications. The actual errors, which might be inaccessible for some fiducial waveform models, are better than 10 −7 (c.f. Fig. 9 and the solid curve in the top panel of Fig. 11).

VI. ASTROPHYSICAL SURROGATES
For pedagogical considerations we have primarily focused on the (2, 2) mode of non-spinning EOB waveforms in the range q ∈ [1, 2] and about 12,000M in duration. In this section we build surrogate models for a variety of astrophysical sources relevant for detection templates and parameter estimation with gravitational wave detectors. Typically six or fewer digits of accuracy suffice for these applications. We therefore build surrogates here with this criteria in mind by considering less ambitious error requirements of ≈ 10 −6 instead of ≈ 10 −9 . The surrogate models presented here also have more cycles, cover larger mass ratio intervals, and include higher spherical harmonics. Surrogates built for these more challenging scenarios continue to be both accurate and fast to evaluate. Most importantly, we can apply exactly the same method described earlier in Sec. III. Table I provides a summary of the surrogate models presented here, which we discuss in more detail below. Surrogates 1 and 2 were discussed earlier in Sec. III. In Sec. III we show that the time-domain overlap error (i.e. the mismatch) is one-half the L 2 error measure we use in Table I. Since our goal is to directly match the output of the EOB-LAL code we do not minimize over intrinsic or extrinsic parameters to compute the error. Hence, both the faithfulness and effectualness diagnostics will be even smaller than those implied by Table I. A surrogate model needs to be computed "once and for all time" for a given set of specifications and so we always strive to make surrogates with the highest possible accuracy unless otherwise indicated. By working to high accuracies one has guaranteed results equivalent to using the full underlying model (in this case EOB waveforms) without the need for special case-by-case studies of systematic biases. For particular applications reduced accuracy could be acceptable, especially at the benefit of faster model evaluations. We shall not pursue such application specific optimizations here.  Top: The greedy error in (9) computed for waveforms of fixed duration and mass ratios in [1, qmax] with qmax = 2, 4, 6, 8 and 10. The curves correspond to cases 1, 3, 4, 5 and 6 from Table I. Bottom: Greedy error for mass ratios in [1,2] with different durations in time, thus numbers of cycles. The curves correspond to cases 1, 7, and 8 from Table I.

A. Larger mass ratios (Cases 1-6)
The top panel in Fig. 12 shows that the number of reduced basis waveforms needed to approximate larger mass ratios accurately increases only mildly. In particular, the number of basis functions to achieve 6 × 10 −8 accuracies grows from 19 to 40 when q min = 1 and q max is raised from 2 to 10. Furthermore, surrogates built for the intervals [9,10] and [1,2] use a total of 33 basis functions, which is nearly the amount (40) needed for the entire [1,10] range. This feature is typical of global approximation methods (such as reduced basis and empirical interpolation) since they tend to promote sparseness whenever the underlying model is sufficiently smooth.  The bottom panel in Fig. 12 shows the number of RB functions needed to accurately cover waveforms with longer durations (i.e., more cycles) also increases mildly. For the longest, most accurate surrogate, 8B, the evaluation time is as large as 10 −3 seconds, which is an order of magnitude larger than the shorter, but otherwise equivalent, case in surrogate 1. However, the EOB-LAL code also runs slower. Thus, the overall speedup is found to be about 750. The lower accuracies required for gravitational wave detection templates (as opposed to the higher accuracy standards for parameter estimation) imply that the speedup can be improved to about 1,100 (see case 8A).

C. Higher harmonics (Cases 1 and 9-12)
Gravitational waveforms have multiple spin-weighted spherical harmonic modes. Surrogate models must accommodate these multi-mode functions in order to maximize their usefulness. One direct approach is to build a surrogate for each mode separately using exactly the same steps described earlier in Sec. III. The resulting multi-mode surrogate model is then defined by the set of single-mode surrogates. Some q = 1 modes are identically zero, and while the reduced basis can exactly approximate zero modes they slightly complicate the treatment of parametric fits (e.g., the phase is undefined): we construct fits on an open interval q ∈ (1, 2]. Our surrogates are thus defined on this open interval with q = 1 modes given by zero. To asses the error of each multimode surrogate model we continue to draw 1,000 random waveform samples from the closed interval [1,2]. As a demonstration, we have built surrogate models for the (2, 1), (2,2), (3,3), (4,4) and (5,5) modes in the same physical and parametric ranges used for the nominal EOB example problem (surrogate 1) considered throughout this paper. These five modes exhaust the currently known ones provided by the EOB model. Compared to the (2, 2) mode, surrogates built for these higher harmonics are more sensitive to peak alignment (see App. E) which translates into larger surrogate errors. These errors are still small and, furthermore, higher harmonics typically have less contribution to the overall gravitational wave strain measurement. Another way to build multi-mode surrogate models starts by integrating the complex scalar product in (2) over the angles (θ, φ) on the 2-sphere. The orthonormality of the spin-weighted spherical harmonics implies a sum over the scalar product in (2) for each mode, which is used for constructing the norm in the greedy error in (9). When the RB-greedy algorithm is performed using the scalar product in (34), all modes contribute to selecting the relevant parameters in the space. Likewise, the residual used for finding the empirical nodes in Step 6 of Algorithm 2, and other quantities required by the EIM, can also be integrated over the angles of the 2sphere so that all modes can contribute to building the global empirical interpolant. We will not cover here all possible variations for multimode waveforms since the choices made are largely dependent upon the specific application and waveform model studied. This is typical in problems that require learning from data, such as this one.

D. Building costs
Ignoring training set generation, the surrogates listed in Table I typically took between 5 and 20 minutes to build. However, we have found the main cost, both in terms of computational and memory requirements, to be in creating the training set. These costs are significantly greater than what might be expected from any particular surrogate's properties (e.g. sampling rate and duration). As discussed further in App. E, the error of an EOB surrogate is dominated by the error in resolving and localizing the waveform's peak. Consequently, we must generate training data having well-resolved waveform peaks. For example, the EOB-LAL code was called using a sampling rate of 2 20 Hz to generate the training set data for cases 1 through 6 in Table I. Surrogate 6 was trained on 2001 EOB waveforms, which took nearly 8 hours to generate. Waveform generation times quoted throughout this paper, for both the EOB-LAL code and surrogate evaluations, do not depend on these settings in anyway whatsoever.
For longer waveforms, such as cases 7 and 8, we are unable to maintain these high sampling rates. A single waveform cannot be produced (on a personal computer) due to the larger memory overhead. In lieu of using higher memory nodes we instead decreased the sampling rate to 2 18 Hz to make the problem more manageable. As a result, the pointwise (maximum) waveform errors increase but remain acceptably small in many cases.

VII. CONCLUDING REMARKS AND OUTLOOK
We introduced a solution to the problem of quickly and accurately generating predictions for a given family of gravitational waveforms. The solution constructs a surrogate for this fiducial set of waveforms in three offline steps. In the first step, a reduced basis is generated that spans the space of waveforms in the given range of parameters. In the second step, an application-specific (i.e., empirical) interpolant is constructed using only these m reduced basis waveforms. The empirical interpolation method selects a corresponding set of m times that are used to build the interpolant but requires knowing the fiducial waveform at any parameter value at those times in order to evaluate the interpolant. In the third step, we complete the offline part by implementing a fit for the parametric dependence of the waveform's phase and amplitude at each empirical time. In this way, the value of the fiducial waveform at each empirical time can be estimated and then fed into the empirical interpolant. The result of these three offline stages is an accurate surrogate model (30) for the underlying family of waveforms that is cheap to evaluate for any parameter value in the considered range.
Surrogate models offer a new and complementary approach to other modeling endeavors. Indeed, our goal is to clone the input-output functionality of an existing waveform generation code thereby permitting fast evaluations for any task of interest. Consequently, we have intentionally worked in a detector-independent context, to very high accuracies, and without regard for the systematic errors of the underlying waveform family. To ensure the surrogate can be used in place of an underlying model (without introducing bias) it is best to be as accurate as possible. However, for particular applications one may wish to sacrifice accuracy at the benefit of even faster surrogate model evaluations.
The standard paradigm for fast online evaluation of new solutions within reduced order modeling frameworks (see, e.g., [54] for a review) is to numerically solve a small problem that is essentially a projection of the original problem onto the basis built in the offline stage. Nonlinear terms or non-affinely parametrized problems can be dealt with using the EIM [55]. This approach has some advantages. For example, for many problems of interest rigorous error bounds can be guaranteed for the resulting output, which is often referred to as a certified approach.
In this paper we deviated from this standard course and sought a different and more heuristic one for two major reasons specific to gravitational waveforms. First, the complexity of projecting the full nonlinear Einstein equations onto a basis to obtain a certified approach is highly nontrivial. Second, our goal has been to develop a non-intrusive approach that does not resort to manipulating, in any way, the original equations and codes that generate the fiducial waveform model. Of course, such equations have to be used to generate the fiducial waveforms in the offline stage in order to build the reduced basis to start the construction of the surrogate model. However, the approach introduced in this paper does not intrude upon or require editing those codes.
In order to demonstrate the basic ideas and methods in this paper, we have focused on surrogate models for single-mode, non-spinning black hole binary EOB waveforms [56]. For mass ratios in [1,2], we find that evaluating the surrogate is three orders of magnitude faster than generating EOB waveforms in the standard way. However, the construction of the surrogate model is not limited to such a short range of mass ratios, to nonspinning binaries, nor to single-mode waveforms. We demonstrated this in Sec. VI by building surrogates for astrophysically motivated problems relevant for template bank generation and parameter estimation studies with gravitational wave detectors. Regarding the range of mass ratios (or other parameters), depending on the application and the target accuracy, a partitioning of the parameter space might provide faster online queries. This issue is familiar when solving differential equations where one may choose to use a single domain or utilize a domain decomposition, as with a spectral or hp-element approach (see, e.g., [57,58]). Similar tools for parame-ter space subdomain decomposition, known as hp-greedy algorithms, have been employed as an adaptive sampling strategy for large problems (see [59][60][61][62] for further details).
A preliminary cost-benefit analysis of domain decomposition is provided by Table I, which summarizes all surrogates considered in this paper. Taking the 19 basis functions for the [1,2] range (surrogate 1 in Table I) as indicative of the reduced basis size needed for each successive integer range of mass ratios up to q = 10, a naive scaling with a domain decomposition approach would suggest ≈ 19 × 9 = 171 basis elements. Compare this number to the 40 elements needed for the whole range [1,10]. While the latter gives fewer basis functions for representing EOB waveforms across the whole [1,10] range, the cost to evaluate surrogate models increases with m. For example, if one were interested in only waveforms with q from 9 to 10 then surrogate 2 in Table I would be preferable in terms of speedup since m = 15, not 40. Such optimizations are application specific and, as such, were not pursued in this paper.
Finally, the method presented in this paper for building a surrogate model can be applied to other waveform families, including precessing inspiral waveforms and multimode inspiral-merger-ringdown waveforms such as those from NR simulations of binary black hole coalescences. We anticipate extremely large speedup factors for predicting a NR waveform with a surrogate model compared to solving the Einstein equations for the same parameters because the cost of evaluating the surrogate is independent of the offline costs required to build it. Given that a single production-quality simulation for a non-spinning equal mass binary takes around ∼ 10 4 − 10 5 hours and predicting a single-mode waveform with a NR-based surrogate model takes about 10 −4 seconds (as implied by Fig. 11), it follows that one may expect speedup factors of ∼ 10 11 or more.

VIII. ACKNOWLEDGMENTS
We thank Frank Herrmann and Evan Ochsner for help during this project, including some software tools, as well as Yi Pan, Alessandra Buonnano, and Collin Capano for helpful discussions about the EOB model and its generation using the LAL code. This work was supported in part by NSF grants PHY-1208861 and PHY-1005632 to the University of Maryland and by NSF grant PHY-1068881 and CAREER grant PHY-0956189 to the California Institute of Technology.

Appendix A: The reduced basis method
We use a greedy algorithm to build a reduced basis (RB), which accurately approximates any fiducial waveform within the given parameter ranges (see, e.g., Ref. [38]). The greedy algorithm, outlined in Algorithm 1, takes as inputs a discretization of the parameter space T ≡ {λ i } M i=1 (or the training space) and the associated waveforms, an arbitrary parameter Λ 1 ∈ T (or seed), and a threshold error for a target representation accuracy (or greedy error). The output consists of the m RB waveforms and m greedy points. 10: The naive implementation of the classical Gram-Schmidt procedure can lead to a numerically illconditioned algorithm. This is related to the fact that the Gramian matrix, which would have to be inverted, can become nearly singular [63]. To overcome this we use an iterated Gram-Schmidt algorithm or a QR decomposition in step 9. See [64,65] for discussions about the conditioning and numerical stability of different orthonormalization procedures.
As mentioned in Sec. III A, minimization over the coefficients {c i } in (8) is satisfied by orthogonal projection P m h(t; λ) of h(t; λ) onto the span of the basis. For example, for an orthonormal basis which takes its global minimum when After applying the greedy algorithm to build a reduced basis and find the greedy points, we check that the basis accurately approximates the continuum space of waveforms for the given parameter range by verifying at a randomly chosen set of test points.
The Empirical Interpolation Method (EIM) provides a sparse subset of empirical time (or frequency) nodes from which it is possible to reconstruct the waveform at any other time with very high accuracy using an applicationspecific interpolant. The selection of the empirical time nodes and the construction of the empirical interpolant proceeds using a greedy algorithm, which is hierarchical and is applicable to unstructured meshes in several dimensions.
Consider a basis {e i } m i=1 (e.g., a RB) whose span approximates the functions of interest. Let {t i } L i=1 denote a set of L time samples and define the L-vector t = (t 1 , t 2 , . . . , t L ) † . For compactness of notation, denote other functions evaluated at these time samples as vectors so that, for example, h(λ) := h( t; λ) and e i := e i ( t).

Given an input of m evaluated basis functions
the output of the EIM algorithm is a set of m empirical nodes The empirical interpolant is constructed in step 5 of Algorithm 2. At the j th iteration the empirical interpolant is built from the first j basis functions and nodes, where the C i coefficients are solutions to the j-point interpolation problem for all λ and where k = 1, . . . , j.

Algorithm 2 The Empirical Interpolation Method
2: i = argmax| e 1 | (argmax returns the largest entry of its argument).
for L equally spaced time samples. The empirical interpolant's error is then directly related to the greedy error (9) through [50] h where the first equality follows from I m [P m h] = P m h, I is the identity matrix, (I−I m ) 2 d = I m 2 d holds whenever the operator norm is induced by the vector norm (as is the case here, see Refs. [66,67]) and is a computable Lebesgue-like quantity that generally changes slowly with m. For problems with smooth dependence with respect to parameter variations we can expect an exponential decay of σ m with m and, from the left side of (20), of the EIM's error.
In practice, Λ m is computed from the matrix representation of B from (19), where each column of E = [ e 1 , . . . , e m ] is an evaluated reduced basis function and V is the interpolation matrix defined in (16). The matrix operator B, as written above, acts on an m-vector h( T ; λ) whose components are evaluations of h at the empirical nodes.

Appendix C: Details of polynomial least squares
When performing a least squares fit we must select the degree n LS of each of the 2m least squares polynomials, balancing accuracy and stability of the resulting fit. For each fit there are m greedy data points so n LS < m. A small value of n LS would result in a low accuracy fit while too large of a value can exhibit Runge's phenomenon [68]. Furthermore, a large value of n LS can fit (numerical) noise thereby leading to low quality fits (this is sometimes called overfitting [69]). Reference [70] provides a computable expression for the largest n LS that avoids this phenomenon and gives an error estimate for the resulting fit.
For our nominal EOB example we proceed in a straightforward way. We construct m separate fits (using only the greedy data) for all degrees 0 ≤ n LS < m and we select the one that minimizes the sum of the squared residuals relative to the training set data. This additional offline work guarantees, in a simple way, that each polynomial fit has the optimal degree. Figure 13 shows the results for our EOB test problem. We see that for empirical times in the early inspiral the optimal polynomial degrees are relatively large and decrease until merger and ringdown. This is a consequence of noisy data stemming from discrete uncertainties in locating the amplitude peak (see App. E for further details).

Appendix D: Surrogate error estimates
In this appendix we derive the error bounds shown in (33) for the surrogate model. We differentiate between the surrogate waveform model h S (t; λ), whose computation requires an estimate for the waveform at each empirical node T i , from the empirical interpolant I m [h](t; λ), whose computation assumes the exact (fiducial) values h. For any λ we have, with Λ m being the same constant defined in (B6). The first equality follows from I m [h S ( T )] = h S ( t). The second line follows from the empirical interpolant's matrix representation (B7). The error in approximating an underlying model h(t; λ) by the surrogate h S (t; λ) is, for any λ, which follows from the error bounds (D1) and (B5) (or (20)) as well as the triangle inequality. Notice that Λ m and σ m are computable quantities as are the differences h(T i ) − h S (T i ), which are only due to least square fitting errors.
Appendix E: On generating the fiducial EOB waveform family In this paper we demonstrated how to build a surrogate using an EOB model of non-spinning binary black hole coalescence waveforms. Here, we discuss some of the technical details regarding how these EOB waveforms were generated.
The specific version of the model that we used is from Ref. [26] and implemented in the routine EOBNRv2 as part of the publicly available LIGO Analysis Library (LAL) Suite [71]. Other versions and models are equally applicable (e.g. [28]). In its simplest description, the code takes as input a starting frequency f min and the mass components m 1 and m 2 . From initial conditions, determined through post-Newtonian expressions, the EOB differential equations are solved to give the system's orbital evolution until merger, which is defined to be the time at which the orbital frequency begins to decrease. From the compact binary system's orbit a gravitational wave is generated up to the time of merger, after which quasinormal modes are attached.
Our nominal EOB example uses a training space of mass ratios q ∈ [1,2]. We sampled this parameter range with 501 equally spaced points, solving the original model at each q using the aforementioned code. We checked that this number of training set samples was dense enough to reach the convergent regime for building a faithful reduced basis representation.
We generated the EOB waveforms with f min = 9Hz and m 1 + m 2 = 80M , which corresponds to roughly 65 − 70 waveform cycles before merger in the (2, 2) mode. We avoided generating short waveforms (where the initial radial separation is less than 20M ) because the ODE initial data could become less accurate. The waveform's coalescence phase was determined implicitly through initial data instead of specifying a particular value [72]. The relevant (2, 2) modes h 22 + (t) and h 22 × (t), as opposed to their spin weighted values, comprised the training set.
Waveforms generated by EOBNRv2 are automatically aligned at f min and thus have lengths that depend on their mass ratio. A typical example is shown in the top panel of Fig. 14. We have found that, when applied to a set of waveforms with varying lengths, the greedy error (9) has a very slow decay rate as indicated by the bottom panel of Fig. 14.
To overcome this, we shift each waveform in time so that their peak amplitudes are aligned. We first align all waveforms in the training set in this way and then "chop off" the beginning portions so that all waveforms have a length (from start to peak amplitude) equal to that of the shortest waveform (here, q = 1). Next, we adjust each waveform's phase (23) to be initially zero. The benefits of waveform alignment are evident from the curve in Fig. 3, which should be compared with the pre-alignment case shown in the bottom panel of Fig. 14. For example, to achieve a greedy error of 10 −7 one needs ≈ 7 (400) with (without) peak alignment. Top: EOB waveforms for q = 1, 2 starting at the same initial frequency but not aligned at the peak amplitudes. Bottom: Not aligning the waveforms results in more reduced basis elements needed to accurately span the space of waveforms. Here we see that nearly all 501 points in the training space are selected whereas only 19 points are required if the waveforms are aligned at the peak amplitude (compare with Fig. 3).
Aligning the waveforms in the manner discussed above is expected to depend smoothly on the mass ratio q since the time of maximum amplitude, measured from the start of an orbital evolution with a fixed f min , is expected to depend smoothly on q. In practice, waveforms are only known at time intervals ∆t so that each waveform's peak time is determined within ∆t. Consequently, aligning discrete waveforms introduces some degree of "nonsmoothness." We initially found the surrogate's error to be dominated by this effect. To overcome this difficulty, we generated each waveform on a temporal grid of spacing ∆t fine , which allowed the time of maximum amplitude to be resolved within ∆t fine . Next, we downsampled each waveform to a sampling rate of interest, say 2,048Hz, such that the peak was located on the downsampled grid. Once a downsampled waveform is generated, neither building nor evaluating the surrogate carries a cost that depends on ∆t fine in any way. Such observations are not unique to surrogate modeling. Indeed, other applications that align waveforms, especially those that need (or expect) some degree of smoothness with parametric dependence, will encounter similar issues. Figure  15 shows how the surrogate error in (31) for our EOB example changes with the grid spacing ∆t fine . In this paper, the smooth parameter dependence and the aligning of the waveforms at the peak amplitude combine to give fast convergence of the surrogate model to the fiducial one.  Fig. 9) as a function of the resolution (i.e., time steps) of the peak amplitude. The trend is linear in ∆t/M . The resolution leads to an uncertainty in estimating the peak amplitudes and thus into aligning the waveforms. This is the dominant source of error in the surrogate model that translates directly into errors in the fits of the last offline step for building the surrogate.

Appendix F: Other approaches for waveform prediction
In this paper we provided a three-step solution for quickly and accurately predicting gravitational waveforms within any given physical model. Here, we discuss a few other approaches that could have been taken instead, which include: i) interpolating the projection coefficients {c i (λ)} m i=1 , defined from (10), in λ, ii) interpolating the (complex) waveforms {h(T i ; λ)} m i=1 at each empirical time, and iii) fitting the amplitudes and phases at all times. The first approach is an alternative to the empirical interpolation in Step 2 and the fitting in Step 3, the second approach is an alternative to the fitting in Step 3 (Sec. III C), and the third approach is an alternative to empirical interpolation in Step 2 (Sec. III B). We consider these in turn.
The first alternative is to build an interpolating (e.g., Chebyshev) grid in λ for each c i (λ). This approach was carried out in Refs. [32,34,35] for inspirals (in the stationary phase approximation in the frequency domain) and phenomenological waveforms for large chirp masses [73]. Problems with such an approach include: i) waveforms from binaries with many GW cycles require increasingly dense interpolation grids [35], ii) the number of grid points scales exponentially with the number of parameter dimensions, iii) standard grid-based interpolation is not hierarchical and dictates sampling locations at predetermined points that are not tailored to the waveform family of interest, and iv) grids that are essential to resolving one projection coefficient may not be useful for resolving another projection coefficient. Finally, the projection coefficients can be poorly behaved as functions of λ because they stem from nontrivial overlaps between waveforms and the basis. This is confirmed using our nominal EOB example, as shown in Fig. 16 where some of the coefficients become noisy, and also in Ref. [35]. Furthermore, using only the m greedy points, which comprise an unstructured and under-sampled grid for this problem, exacerbates many of the aforementioned problems. The second alternative is to interpolate in λ the complex waveforms at each empirical time. This approach has the same problems as interpolating the projection coefficients discussed above. Figure (17) shows the structure of the waveforms as a function of mass ratio at several empirical times in our nominal EOB example.
The third alternative is to perform fits for the waveform amplitude and phase at all time samples instead of the ones dictated by the EIM. It is instructive to compare the operation counts for the online evaluation between this all-times fitting alternative and our EIMbased method. If c fit is the operation count of the fitting functions at each time, taken to be constant for simplicity, then the dominant operation count is 2c fit L for the all-times fitting and 2m(L + c fit ) using EIM and fitting at each empirical time. Therefore, making the reasonable assumption that m L, the EIM-based approach is more efficient whenever c fit m .
In one parameter dimension, the standard way of eval- FIG. 17. The curves depict the values of the real parts of the waveforms along with the greedy data as a function of the mass ratio for our EOB example case. Only curves at a few representative empirical times are shown. While there is less structure here than appears in the coefficients shown in Fig. 16, the majority of functions still require additional sampling to be accurately resolved by global polynomial fits.
uating a polynomial fit of degree n is through Horner's algorithm [74], which has an optimal operation count of 2n. It would then seem to follow from (F1) that the online evaluation cost of the EIM-based approach is comparable to fitting at all L times. However, operation counts can be misleading as they do not take into consideration other aspects of an algorithm's implementation that are also relevant for the total execution time. We conducted numerical experiments with our nominal EOB example and found that, for our particular implementation, fitting at all L ≈ 10,000 samples is between 20 and 1,000 times slower. These timing experiments depend sensitively on both the number of surrogate basis/nodes as well as using "vectorized" for-loops. Therefore, the actual online evaluation cost in the examples considered in this paper are consistently an order of magnitude or more faster than what a naive operation count would suggest. The operation count for evaluating polynomials grows with the dimensionality. While the most efficient scheme for evaluating multivariate polynomials is not presently known [75], it is an active area of research. In general, (F1) is easily met by surrogate models in higher parameter dimensions and we expect the EIM-based surrogate approach to be more efficient than one based on fitting at all output times. In addition, the cost to construct L separate fits in higher parameter dimensions could make this offline step prohibitively expensive.