Arby - Fast data–driven surrogates

Several problems arise in observational and theoretical contexts that demand the resolution of computationally intensive differential equations. From structural analysis in engineering to spacetime simulations in astrophysics, rapid and reliable evaluations of solutions to these equations are crucial due to the need for computing in real or near-real time quantities that depend on those solutions. Due to pervasive computational costs led by the inherent complexities present in many problems, achieving fast responses becomes an ubiquitous bottleneck.

As a case study, let us take the example of an ongoing problem in the field of gravitational wave research, namely the template bank problem (Field et al. (2012)). To be able to detect the very faint gravitational wave signals among the noisy background of ground-based interferometers, the LIGO-Virgo collaboration use match filtering against a bank of GW signal templates to maximize the signal-to-noise ratio (SNR) (Cutler & Flanagan (1994); Abbott et al. (2016, 2020)).

The observed time series is filtered to decide whether a binary coalescence took place, and the parameter estimation allows us to infer the properties of precursors of the remnant, such as masses and spins for binary black hole mergers.

It is convenient to count with large enough template banks of theoretical waveforms to fill the parameter space of GWs. However, this pose an exceedingly challenging task due to the complexity of solving the Einstein Equations of General Relativity. A single numerically generated GW waveform can take from days to weeks to become available for production (Lehner & Pretorius (2014)).

A palliative solution is the construction of approximate waveform models, which are more direct to build and deploy than numerical relativity ones. There are several methods to build these approximations. Analytical examples are the Post-Newtonian (Blanchet (2006)), Effective One Body (Damour & Nagar (2011)) and Phenomenological (Sturani et al. (2010); Hannam et al. (2014); Khan et al. (2019)) approximations. However, we focus here on a particular set of methods that have proven to be very fertile in GW research and produced some of the milestones in waveform modeling in recent years: Reduced Order methods.

Reduced Order Modeling (ROM) is an umbrella term that encompasses a variety of techniques for dimensional reduction, developed to address the problem of complexity in numerical simulations. In particular, Surrogate Models obtained through application of ROM to ground truth solutions are low resolution representations intended to be fast to build and evaluate without compromising accuracy. We take a data-driven approach, i.e. driven only by data, as opposed to more standard and intrusive ones in which reduction methods are coupled to differential solvers to build solutions (Quarteroni et al. (2015); Hesthaven et al. (2015)).

In waveform modeling, the combination of two ROM methods originally posed for intrusive problems and recreated later for data-driven ones led to significant success in constructing surrogate models for GWs. Those methods are the Reduced Basis (RB) (Boyaval et al. (2009); Field et al. (2011, 2012)) and the Empirical Interpolation (EI) (Barrault et al. (2004); Chaturantabut & Sorensen (2010)) methods, which we describe in the next section.

The primarily purpose of this work is to disseminate these tools in the astronomy/astrophysics community by introducing a single and user-friendly Python package for data-driven dimensional reduction: Arby.

Arby arises as a response to a lack of well documented, tested and actively maintained code for reduced basis and surrogate modeling in the scientific community while adhering to the data-driven and user-friendliness premise.

This section briefly describes the basics of a reduced-order pipeline for building surrogate models. As we stated above, it merges two main ingredients, the Reduced Basis and the Empirical Interpolation methods, for dimensional reduction of raw data. This pipeline was first introduced in (Field et al. (2014)) followed by the construction of several surrogate waveform models based on this method (Blackman et al. (2015, 2017b, 2017a); Varma et al. (2019a, b); Rifat et al. (2020)). See also (Tiglio & Villanueva (2021)) for a review.

In the first stage, we look for a low-resolution representation of the model $h_{\lambda}$. The RB method consist in representing a whole set of solutions, usually called the training set

by linear combinations of basis elements of the form

where

defines the inner product between training functions, $\bar{e_{i}}$ is the complex conjugate of $e_{i}$ –if we deal with complex functions– and $\Omega$ is the physical domain. The set $\{e_{i}\}_{i=1}^{n}$ is called the reduced basis and is composed by a subset of optimally chosen solutions from the training set. The construction of the reduced basis is iterative: at each step of the algorithm, the most dissimilar (orthogonal) element from the training set ${\cal K}$ joins to the current basis, and the process stops when an user-specified tolerance is met. This tolerance is related to the maximum error of the difference between solutions and approximations. In consequence, the different spaces $\mathcal{X}_{i}(i=1,2,\ldots)$ spanned by each reduced basis built at each step are nested, i.e. $\mathcal{X}_{1}\subset\mathcal{X}_{2},\ldots$ The addition of a new element to the basis implies a fine tuning of the previous approximation space.

The RB method supposes a compression in the parameter space of solutions. The next step is to achieve compression in time. To this we turn to interpolation replacing the projection-based approach described so far by an interpolation scheme. We pose the problem by looking for an efficient linear interpolation operator ${\cal I}_{n}$ such that

subject to

for strategically selected nodes $X_{i}\,(i=1,\ldots,n)$ out from the physical domain. The EI method (Maday et al. (2009); Barrault et al. (2004); Chaturantabut & Sorensen (2010)) gives us an algorithm for building such interpolant. The Empirical Interpolation (EI) algorithm, as described in (Field et al. (2014)), selects iteratively the nodes $\{X_{i}\}$ from the physical domain following a local optimization criteria (Tiglio & Villanueva (2020)).

The EI algorithm receives as unique input the reduced basis and selects the interpolation nodes for building the interpolant. Note that there is no need of the whole training set ${\cal K}$ since the RB algorithm already did the introspection of it, and we assume that the relevant information about ${\cal K}$ is hard-coded in the reduced basis.

It is possible to show that, under some conditions, the interpolation error is similar to the projection one, which in most applications has exponential decay in the number of basis elements (see (Tiglio & Villanueva (2021)) and citations therein). This leads to an efficient and (in most cases of interest) accurate representation of the training set by means of an empirical interpolation.

To summarizing: first, a reduced basis is built from a training set using the RB method, which leads to the linear representation in (2). This step completes a compression in the parameter domain. Next, an empirical interpolant is built solely from the reduced basis. This step completes a compression in the physical domain. Finally, we end up with an empirical interpolation (4) which provides efficient and high-accuracy representations of all functions in the training set.

Let us rewrite the interpolation in (4) as

where

and

Recall the approach is data-driven, so we do not fill the training set with more solutions to approximate newer ones. Instead, we predict them by performing fits along data that we already know, that is, along $h_{\lambda}(X_{i})$ $(i=1,\ldots,n)$. For $1$-D problems ($\lambda\in\mathbb{R}$), deciding for problem-agnostic fitting methods that are well-suited for most cases can be a challenging task, not to mention the high dimensional case, which remains an open question (Tiglio & Villanueva (2021)). For a rough classification, we refer to regression and interpolation methods. The former deals with calibrating free parameters of some model by optimizing an error function; the latter, with solving an interpolation problem consisting essentially in solving algebraic systems which are possibly subject to constraints (e.g. Eqs. (4,5)).

Here we address the second approach for parametric fits. The procedure consists in interpolating the values

along parameter samples for each empirical node $X_{i}\,,\,i=1,\ldots,n$. As we describe in Section 4, as of now, Arby uses splines for this step, though we expect to generalize this. Once fits are done, we end up with a surrogate model $h_{\lambda}^{surr}(x)$ which can represent and predict $h_{\lambda}$ at any $\lambda$ in the parameter domain with high accuracy and low computational cost.

The pipeline for building surrogates described in the previous section is valid in any number of parameter and physical dimensions if we consider arbitrary fittings through parameter space. For surrogate modeling, Arby supports in its present version 1-D parameter and physical domains (real number intervals of the form $[a,b]$ for $a,b\in\mathbb{R}$) and real-valued functions. Again, this restriction is only for building surrogate models. On the other hand, Arby supports multidimensional parameter domains (although still restricted to 1-D domains in the physical dimension) and complex-valued functions for building reduced bases and empirical interpolants.

Below we summarize the algorithm for surrogate modeling. We refer the reader to the Appendices for technical details about the RB and EI algorithms used in intermediate stages.

The inputs are the training set ${\cal K}=\{h_{\lambda_{i}}\}_{i=1}^{N}$, the parameter set ${\cal T}:=\{\lambda_{i}\}_{i=1}^{N}$, and the greedy tolerance $\epsilon\in\mathbb{R}$.

Let’s make some remarks on Alg. 1.

• •

  The training set ${\cal K}$ is built from a discretization ${\cal T}$ of the parameter domain. In the current version of Arby, ${\cal T}$ is a discretized real interval ${\cal T}:=\{\lambda_{1},\ldots,\lambda_{N}\}$.

• •

  The RB algorithm used for building the reduced basis is fully described in Alg. 2, see the Appendices. It selects from ${\cal T}$ $n$ points $\lambda_{i}=\Lambda_{i}\,(i=1,\ldots,n)$, called the greedy points, labeling those functions in the training set that conform the reduced basis. For conditioning purposes, the basis is orthonormalized at each step, so the algorithm’s final output is a set of orthonormal basis elements along with the set of greedy points. So we use the term reduced basis interchangeably for both, the basis conformed by greedy solutions and its orthonormal version, due to its equivalence (they span the same space).

  The number $n$ depends on the greedy tolerance $\epsilon$. In Arby we must specify a discretization of the physical domain $\Omega$ so as to be able to do integrals (see (3)). In this context, $\Omega$ is a real interval $[x_{a},x_{b}]$ and Arby must receive as input an equispaced discretization of it.

• •

  Step 3 implements the EI algorithm described in Alg. 3, see the Appendices. It receives the reduced basis as unique input and finds $n$ empirical nodes $X_{i}\,(i=1,\ldots,n)$ to build the interpolant. In practice, the interpolant is specified by assembling the $n$ functions $B_{i}(x)$ defined in Eqs. (6,7).

• •

  To achieve predictability Steps 4-6 perform parametric fits along training values for each empirical node $X_{i}$. Let’s illustrate this by looking at the first iteration of the loop in Steps 4-6. For the first node $X_{1}$ collect all values

  $$\{h(\lambda_{1};X_{1}),h(\lambda_{2};X_{1}),\ldots,h(\lambda_{N};X_{1})\}$$

  and perform a fit along them. This results in a function $h^{fit}_{1}(\lambda)$ that is continuous in the interval $[\lambda_{1},\lambda_{N}]$. This is repeated $n$ times for each empirical node. The resulted functions $h(\lambda)_{i}^{fit}$ constitute along the reduced basis the building blocks for the final surrogate assembly. The current Arby version implements splines (J.H. Ahlberg & Walsh (1967)) for parametric fits, i.e., piecewise polynomial interpolation of some degree arbitrarily set by the user.

• •

  From eqs. 6 and 7 the empirical interpolant is defined through the $n$ functions $B_{i}(x)$. They comprise all the RB-EI information. Combining the functions $B_{i}(x)$ with the fits $h^{fit}_{i}(\lambda)$ built on previous steps, Step 7 leads to the desired surrogate $h^{surr}(\lambda;x)$ which is continuous in $\lambda$ inside the real interval $[\lambda_{1},\lambda_{N}]$.

Arby is a Python package for data-driven surrogate modeling satisfying standard software compliance along quality assurance (see Section 4.4). It allows the user for building reduced basis and empirical interpolants at any number of parameter dimensions. At the current release, Arby builds surrogate models for 1-D domains.