# Reduced Order and Surrogate Models for Gravitational Waves

Manuel Tiglio111email<EMAIL_ADDRESS>Aarón Villanueva222email:
<EMAIL_ADDRESS>
(Facultad de Matemática, Astronomía, Física y Computación,
Universidad Nacional de Córdoba, Córdoba (5000), Argentina
)

###### Abstract

We present an introduction to some of the state of the art in reduced order
and surrogate modeling in gravitational wave (GW) science. Approaches that we
cover include Principal Component Analysis, Proper Orthogonal Decomposition,
the Reduced Basis approach, the Empirical Interpolation Method, Reduced Order
Quadratures, and Compressed Likelihood evaluations. We divide the review into
three parts: representation/compression of known data, predictive models, and
data analysis. The targeted audience is that one of practitioners in GW
science, a field in which building predictive models and data analysis tools
that are both accurate and fast to evaluate, especially when dealing with
large amounts of data and intensive computations, are necessary yet can be
challenging. As such, practical presentations and, sometimes, heuristic
approaches are here preferred over rigor when the latter is not available.
This review aims to be self-contained, within reasonable page limits, with
little previous knowledge (at the undergraduate level) requirements in
mathematics, scientific computing, and other disciplines. Emphasis is placed
on optimality, as well as the curse of dimensionality and approaches that
might have the promise of beating it. We also review most of the state of the
art of GW surrogates. Some numerical algorithms, conditioning details,
scalability, parallelization and other practical points are discussed. The
approaches presented are to large extent non-intrusive and data-driven and can
therefore be applicable to other disciplines. We close with open challenges in
high dimension surrogates, which are not unique to GW science.

###### Contents

  1. 1 Introduction
  2. 2 Reduced Order Modeling and Gravitational Waves
  3. 3 Mathematical Preliminaries
  4. I Representation and Compression
    1. 4 Principal Component Analysis
    2. 5 Proper Orthogonal Decomposition
      1. 5.1 PCA and POD
    3. 6 Spectral expansions
      1. 6.1 Spectral methods
        1. 6.1.1 Fourier Expansions
        2. 6.1.2 Jacobi Polynomials
    4. 7 Parametrized problems and optimal approximations
    5. 8 Reduced Basis
      1. 8.1 Introduction
      2. 8.2 The Training set
      3. 8.3 Greedy algorithms
        1. 8.3.1 Convergence rates and near-optimality
        2. 8.3.2 Complexity, scaling, computational aspects
  5. II Predictive Models
    1. 9 Polynomial interpolation
      1. 9.1 Representation versus prediction
      2. 9.2 Prediction through polynomial interpolation
      3. 9.3 Convergence rates, collocation points and Runge’s phenomenon
      4. 9.4 Discrete expansions and interpolation
      5. 9.5 Multiple dimensions
    2. 10 The Empirical Interpolation Method
      1. 10.1 From projection to interpolation
      2. 10.2 EIM algorithm
      3. 10.3 Accuracy and conditioning of the EIM
    3. 11 Surrogate models
      1. 11.1 Surrogate models for components
      2. 11.2 Empirical interpolant based surrogates
    4. 12 Surrogates of compact binaries
      1. 12.1 Numerical relativity binary black holes
      2. 12.2 Numerical relativity hybrid binary black holes
      3. 12.3 Extreme mass ratio inspirals
        1. 12.3.1 Eccentric inspirals
      4. 12.4 Effective One Body
      5. 12.5 Post-Newtonian (PN)
      6. 12.6 Ringdown
  6. III Data Analysis
    1. 13 Reduced Order Quadratures
      1. 13.1 ROQ, other quadrature methods, dependence on dimensionality
    2. 14 Accelerating Parameter Estimation with ROQ
      1. 14.1 Constructing the Reduced Order Quadrature
      2. 14.2 Implementation in LALInference
    3. 15 Challenges, open issues and comments
    4. 16 Acknowledgments

#### Notation

  * •

$h_{\lambda}$: a function (for the purposes of this article, a GW waveform)
associated with parameter value $\lambda$. It can be, for example, a time
series $h_{\lambda}=h_{\lambda}(t)$, or a frequency domain function
$h_{\lambda}=h_{\lambda}(f)$.

  * •

$\lambda$: a (usually multi-dimensional) parameter for $h_{\lambda}$.

  * •

$\tt dim$: number of parameters of $\lambda$.

  * •

$L$: number of time/frequency samples.

  * •

${\cal F}$: an abstract space of functions of interest. Typically, ${\cal
F}:=\\{h_{\lambda}\,|\,\lambda\in\Phi\\}$, with $\lambda$ in some compact
region (continuous or discrete) $\Phi$. We also refer to ${\cal F}$ as the
fiducial or underlying model.

  * •

$N$: size of the training set. Also, the number of points in standard
quadratures, polynomial interpolation, etc.

  * •

${\cal T}$: training set of parameter points ${\cal
T}:=\\{\lambda_{i}\\}_{i=1}^{N}$. It is assumed to be compact.

  * •

${\cal K}$: training set of functions in ${\cal F}$, ${\cal
K}:=\\{f_{i}\\}_{i=1}^{N}=\\{f(\lambda_{i},\cdot)\\}_{i=1}^{N}\subset{\cal
F}$.

  * •

RB: Reduced Basis as a framework.

  * •

$\Lambda_{i}$: selected (typically, by a greedy algorithm) parameter values.

  * •

rb: a specific reduced basis.

  * •

$t,f$: time, frequency domains.

  * •

EIM: The Empirical Interpolation Method.

  * •

$T,F$: selected (typically, by the EIM) time/frequency interpolation nodes.

  * •

$n$: number of basis elements and of EIM nodes (they are equal by
construction). The goal of ROM, if the problem is amenable to dimensional
reduction, is to find an accurate basis such that $n\ll N$.

  * •

$\Lambda_{n}$: the Lebesgue constant for an approximation with $n$ basis
elements/interpolation points. This involves some ambiguity: we use the same
symbol for greedy selected points, but it should be clear from context which
one we refer to.

  * •

$h_{\tt s}$: a surrogate function for $h$.

  * •

Boldfaces are used for matrices, for example $\bf A$ has elements $A_{ij}$.

  * •

$\langle a,b\rangle$: scalar product between two vectors or functions $a,b$.

###  1 Introduction

Gravitational wave (GW) science has reached a level of maturity in which many
tools from areas such as modern approximation theory, data science, machine
learning, and artificial intelligence are being incorporated into the field.
These attempt to address challenges such as dealing with complex modeling,
analysis, and handling of big data. A common feature of these challenges is
the computational cost involved, which in many cases can be prohibitive, to
the point that it cannot be solely overcome with larger or faster
(super)computers, specialized hardware such as GPUs, or software optimization.
This is particularly the case for parametrized problems, where each query
depends on multiple input parameters that might only be known at run time.
This is exacerbated as the number of parameters grow, usually resulting in the
curse of dimensionality. This refers to the complexity of the problem (here
leaving the term complexity ambiguous on purpose) growing fast, sometimes
exponentially, with the number of parameters. In the case of gravitational
waves from binary systems, parameters can be intrinsic or extrinsic. The
former relate to parameters such as the mass and spin of the binary
components, the initial separation and eccentricity of the system, and
equations of state if matter is present. Extrinsic parameters include distance
of the source to Earth, sky position, orientation, time and phase of arrival.

One of the first challenges of a parametrized problem is sampling it. With
standard methods (for example, equally spaced, using the metric approach in
GWs, or stochastic sampling) the accuracy of such catalogs for GW detection
increases in many cases at best linearly with the number of samples, and their
sizes typically increase exponentially with the number of parameter dimensions
dim as ${\cal O}(N^{\tt dim})$. In fact, for the metric approach, the number
of templates grows as $(1-\text{MM})^{-\text{dim}/2}$, with MM the minimal
match. Producing such catalogs, as well as their storage and analysis, can
become challenging and, again, even not tractable through raw computational
power and storage.

One approach to this problem is decreased fidelity: the modeling of the
problem is simplified by approximations to Einstein’s General Relativity
equations which are cheaper to solve for and analyze. Decreased fidelity is a
delicate approach, though, since the accuracy of the approximations might not
be known without access to the high fidelity models for an arbitrary query,
which if not available could become an issue when attempting to assign error
bars and a given precision needs to be guaranteed for statistical purposes.
Decreased fidelity models can also lead to missed signals or biases in
parameter estimation of detected GWs.

Thus, another challenge is to simulate the GWs emitted by, for example, the
coalescence of compact binary objects in real time, without any physical
approximation. Here the gold standard is high accuracy numerical relativity
(NR) solutions to the full Einstein field equations. It might appear
unattainable to replace online (that is, not evaluated in advance), for an
arbitrary query, a supercomputer simulation which might take per query
hundreds to tens of thousands of hours of computing time with a substitute or
surrogate model of equal accuracy but which can be evaluated in real time on a
commodity laptop. As we show throughout this review, this can be – and is
being – done, though there are still outstanding challenges left and the
problem is not completely solved.

Finally, another challenge of equal importance is to perform parameter
estimation on the sources of any detected GW signal in quasi real time.
Meaning, fast enough so as to process the large number of detections by modern
GW laser interferometers. And, most important, in the case of sources with
electromagnetic counterparts, fast enough to allow for rapid telescope
followups. This requires both online fast evaluation of the GW waveforms, as
well as rapid likelihood evaluations.

In this review we discuss some state of the art approaches which attempt to
(and in some, but not all, cases manage to) obtain accurate and fast-to-
evaluate and analyze surrogate models of gravitational waves emitted by binary
systems, accomplishing some of the aforementioned challenges. The common
aspect underlying all these methods is reduced order modeling (ROM). The
approaches here reviewed are Principal Component Analysis (PCA), Singular
Value Decompositions (SVD), Reduced Basis (RB), the Empirical Interpolation
Method (EIM), and derivatives, such as Reduced Order Quadratures (ROQ). The
goal is intendedly not to be a definite survey, since it is a very active
field. Instead, we attempt to provide an introduction to some of the
approaches that are being used in practice and do deliver on some or several
of the above challenges, along with some basic theory underlying each of them.
The review is divided into three parts, each of which builds upon the previous
one:

  1. 1.

Representation

One example is the generation of compact catalogs or banks of gravitational
waves for searches. Another one is to analyze a system and look for
redundancies.

  2. 2.

Predictive Models

The most ambitious goal here is to build surrogate models that can be
evaluated in real time and are indistinguishable from numerical relativity
supercomputer simulations of the Einstein equations, without any physical
approximation. The target is the evaluation of a waveform in the order of
milliseconds per angular mode; that is, a speed up of at least $\sim 10^{8}$
with respect to NR and without loss of accuracy.

  3. 3.

Data Analysis

One of the main goals of ROM and other efforts in GW science is to achieve
very fast parameter estimations, in particular so that real time alerts can be
sent for searches of electromagnetic counterparts. From an astrophysical point
of view, the target is from months using standard methods to around $10$
minutes using Focused Reduced Order Quadratures (discussed in Section 14),
including millions to tens of millions waveform evaluations and likelihood
evaluations.

We encourage the reader to provide us with feedback, including topics to cover
in future versions of this Living Review article. For briefness we have to
skip many references, we apologize in advance for any and all omissions in
such an active field.

###  2 Reduced Order Modeling and Gravitational Waves

This article reviews approaches which attempt to solve many of the
aforementioned problems in GW science through Reduced Order Modeling (ROM),
also known as Dimensional or Complexity Reduction, and the related field of
Surrogate Modeling. ROM as a field has been around for a long time, but over
the last decade and a half there have been major theoretical advances followed
by a rapid raise in the number of applications and pace at which they both
take place, in many areas of science, engineering and technology.

This is in part, again, due to powerful new approaches and results, from
approximation theory to numerical analysis and scientific computing, but also
due to the recognition of the power of dimensional reduction in many important
problems, some of them long-standing. That is, problems which are only now
being seen in the light of the latest developments of ROM. On top of that,
fields such as data science (DS), Machine Learning (ML) and Artificial
Intelligence (AI) are considerably benefiting from ROM to either eliminate
redundancies in big data, making them more amenable to analysis, and/or
identifying relevant features for further studies using techniques from these
disciplines. In fact, depending on the definitions of these fields, some
include ROM as a subdiscipline.

There is a big difference, though. There are many books on established DS, ML,
and AI practices and theory, while the literature on modern approaches to ROM
for parametrized problems, with some important exceptions that we mention in
this review, is largely composed of technical papers, which could be difficult
to grasp for practitioners as introductory material. Also, these notable
exceptions usually focus on time-independent partial differential equations
(PDEs), and even as introductory material they might be hard to absorb by non-
mathematicians. In contrast, this review covers approaches which are purely
data-driven and attempts to be amenable as an introduction for GW science
practitioners, keeping theory to the minimum necessary to build intuition on
why a given approach works as it does, what are the challenges left, and
possible approaches to them.

One of the reasons why ROM for parametrized system has been largely devoted to
time-independent PDEs is because in that case many rigorous results, such as a
priori error bounds on the reduced model, can be proved; these are called
certified approaches. Here we move beyond that particular arena to ROM for
general parametrized systems. This includes problems which might not involve
any differential equation, be time-dependent, to purely data-driven problems
such as analysis and handling of large amounts of data. This broadening is at
the expense, in many cases, of more heuristics and less rigor, such as a
posteriori validation as opposed to rigorous a priori error bounds. The
rationale for this is simply that many problems of interest in GW science are
far too complex for existing detailed rigorous theorems. On the other hand, GW
problems share many similarities with others for which approaches with proven
properties have been developed and therefore certain algorithms can be adapted
to cases of interest in GWs. In other instances it is quite the opposite:
available rigorous results are quite generic, and can be abstracted from any
previous application in which they were introduced.

Half-rigorous, half-heuristic approaches are not as pessimistic as it might
seem. Many widely used techniques in DS, ML and AI are somewhat heuristic in
nature and require a posteriori validation. Similarly with half-rigorous,
half-heuristic cutting-edge ROM approaches as those discussed here: they can
(and should) always be a posteriori validated.

Our presentation keeps technicalities and analyses to the minimum necessary
for building intuition of why a technique works as it does, under which
conditions it does so, and what can be expected given the properties of the
problem of interest. At the same time, it is not a survey of all the work done
in the field.

There are two important books about Reduced Basis; though focused on partial
differential equations there is a lot of valuable information (well beyond the
scope of this short introduction): [70] and [111]. A website for ROM with
several resources is [9], while [11] has the speakers slides from a workshop
on ROM in General Relativity; even though the latter is several years old and
there has been much progress since then, the topics covered are still highly
relevant. A very recent and lengthier, year-long program at ICERM (Brown
University) focused on ROM in GR can be found on the website [6], including
many of the slides and talks of each workshop.

###  3 Mathematical Preliminaries

This section serves as a brief recap of some basic mathematics and to
introduce some notation used hereon. Since in this review we deal with model
reduction for parametrized problems, we next consider an example amenable to
it.

Example 1.

Consider the model

$h(t;\lambda)=\exp(\lambda t)\,.$ (1)

where $\lambda$ is a complex parameter and $\operatorname{Im}({\lambda})>0$.
This represents oscillatory, exponentially damped functions; it can be seen as
a toy model for the ringdown of a black hole (discussed in Section 12.6).

In the notation that we use throughout this review, $h$ is a function of the
physical variable $t$ parametrized by $\omega$ and for compactness we often
use the following type of shortcuts:

$h=h_{\lambda}=h(\lambda)=h(t;\lambda)\,.$

Suppose now that we have a set of $N$ samples of $\lambda$, leading to a
training set

${\cal
K}=\\{h_{i}\\}=\\{h_{\lambda_{i}}\\}=\\{h(\lambda_{i})\\}=\\{h(t;\lambda_{i})\\}\,,\quad
i=1\ldots N.$

Since the functions of interest (1) are known in closed-form, building such a
training set is straightforward, though in general this is not the case (as
when solving the Einstein equations is required).

One of the goals of ROM is to find a reduced basis, that is, a subset of
${\cal K}$ with number of elements $n\leq N$ – with, hopefully, $n\ll N$ –,
such that its span (the span of a set of vectors is its set of linear
combinations) represents ${\cal K}$ with arbitrarily high accuracy. The
compression rate is then

$C_{r}:=N/n\,.$

Furthermore, in practice one also needs to discretize time (or frequency). So
let’s sample these training set elements at an arbitrary set of times
$\\{t_{i}\\}_{i=1}^{L}$, where $L$ stands for Length. Using, for example, the
Empirical Interpolation Method (EIM, discussed in Section 10), one can
subsample these time samples to a subset $l\leq L$. In fact, in EIM, if the
number of reduced basis is $n$, then by construction $l=n$.

That is, the initial set being of size $N\times L$ can now be reduced to
$n^{2}$, with a double compression rate (both in parameter and time domains),

$C_{r}:=\frac{N\times L}{n^{2}}\,.$

#### Inner products and norms for functions or vectors

Let $\Omega$ denote the physical domain. Throughout this review it represents
time or frequency interval, $[t_{\rm min},t_{\rm max}]$ or $[f_{\rm
min},f_{\rm max}]$ respectively, though in general it could be space, space-
time, or some more abstract arena. For the sake of discussion we consider time
intervals. For any two complex-valued functions, $h(t)$ and $g(t)$, we
consider inner/scalar/dot products and their corresponding norms of the form

$\langle
h_{1},h_{2}\rangle:=\int_{\Omega}\bar{h}_{1}(t)h_{2}(t)\omega(x)dx,\quad\|f\|^{2}:=\langle
f,f\rangle\,,$ (2)

where the bar over $h_{1}$ denotes complex conjugation and $\omega(x)$ is a
generic weight function. These are referred to as weighted $L_{2}$ scalar
products. In data analysis/signal processing and usually in the frequency
domain, $S:=\omega^{-1}$ characterizes the sensitivity of the detector and is
referred to as the power spectral density (PSD).

We consider also discrete inner products and norms of the form

$\langle
h_{1},h_{2}\rangle:=\sum_{i=1}^{L}\bar{h}_{1}(t_{i})h_{2}(t_{i})\omega_{i}\,\quad\|h\|^{2}:=\langle
h,h\rangle\,.$ (3)

This discrete version of an $L_{2}$ scalar product is usually denoted as
$\ell_{2}$, with a lowercase to distinguish it from the continuum case.

Whenever (3) is the discrete approximation of an integral (2) the values
$t_{i}\in\Omega$ and $\omega_{i}$ are respectively referred to as quadrature
nodes and weights. Together, $\\{t_{i},\omega_{i}\\}_{i=1}^{L}$ is called a
quadrature rule. When $\omega_{i}=1$, (3) is referred to as the Euclidean
scalar product. Throughout this review we also use the infinity, or max, norm

$\|h\|_{\infty}:=\max_{1\leq i\leq L}|h(t_{i})|\,.$

of a vector $f\in\mathbb{C}^{L}$. This norm is not induced by an inner
product. That is, it can be shown that there is no scalar product
$\langle\cdot,\cdot\rangle$ for which $\|h\|_{\infty}=\langle h,h\rangle$.

Formally speaking, (2) is defined for functions while (3) is for vectors. We
use the same notation $\langle\cdot,\cdot\rangle$ for both cases; with the
hope that the distinction is clear from the context. Some proofs are sometimes
more convenient in a continuous setting while numerical computations are
restricted to discrete values. Depending on the context, we sometimes switch
between these two settings using the same notations for both of them.

Example 2.  Polynomials.

Consider the space of degree $n$ polynomials defined on $\Omega=[-1,1]$. Any
element

$f(t)=\sum_{i=0}^{n}c_{i}t^{i}\,,$ (4)

can be written as a sum of $(n+1)$ terms and such space is a linear one of
dimension $(n+1)$. The first $(n+1)$ normalized Legendre polynomials
$\\{P_{i}(t)\\}_{i=0}^{n}$ form an orthonormal basis with respect to the
scalar product (2) with $\omega(t)\equiv 1$.

#### Matrices

Consider a complex $L\times N$ matrix $\textbf{H}\in\mathbb{C}^{L\times N}$,
where $L$ is the number of rows (“length of the time series”) of H and $N$ is
the number of columns (“number of parameter samples”),

${\bf H}:=\begin{pmatrix}h_{11}&h_{12}&\cdots&h_{1N}\\\
h_{21}&h_{22}&\cdots&h_{2N}\\\ \vdots&\vdots&\ddots&\vdots\\\
h_{L1}&h_{L2}&\cdots&h_{LN}\end{pmatrix}\,.$ (5)

As one might imagine, for gravitational waves $h_{ij}=h_{j}(t_{i})$. Each
column is referred to as a snapshot. A matrix ${\bf H}$ of shape $L\times N$
appears often throughout this review, where $N$ is the number of samples in
the training set (or $n$, the dimensionality of the reduced basis) and $L$ is
the length of each time series, so it is a very concrete example.

A square matrix ${\bf H}$ is said to be non-singular if its inverse – denoted
by ${\bf H}^{-1}$ – exists, i.e. ${\bf H}{\bf H}^{-1}={\bf H}^{-1}{\bf H}={\bf
I}$, where

${\bf I}=\begin{pmatrix}1&0&\cdots&0\\\ 0&1&\cdots&0\\\
\vdots&\vdots&\ddots&\vdots\\\
0&0&\cdots&1\end{pmatrix}=:diag(1,1,\ldots,1)\,.$

The (Hermitian) transpose of ${\bf H}$, denoted by $H^{\dagger}$ is defined as
the matrix with elements

$\left({\bf H}\right)^{\dagger}_{ij}=\bar{{\bf H}}_{ji}\,,$

the bar again indicates complex conjugation, and we say that ${\bf H}$ is
symmetric (hermitian or self-adjoint, in the current context), if ${\bf
H}={\bf H}^{\dagger}$.

Eigenvectors and eigenvalues:

A vector $x\neq 0$ is an eigenvector of ${\bf H}$ with eigenvalue $\alpha$ if

${\bf H}x=\alpha x\,.$

Eigenvectors are defined up to a normalization constant; that is, if $x$ is an
eigenvector with eigenvalue $\alpha$ so is $ax$, for any non-zero $a$.

##### Matrix Norms

The Frobenius norm. Imagine “unpacking” ${\bf H}$ into a long vector of size
$L\times N$. Measuring this vector with the Euclidean norm defines the
Frobenius norm

$\|{{\bf H}}\|^{2}_{F}:=\sum_{i}^{L}\sum_{j}^{N}\,\left|h_{ij}\right|^{2}$ (6)

of a matrix ${\bf H}$.

Later we will see that the Frobenius norm plays a key role when discussing
approximation by proper orthogonal decompositions.

Induced norm.

The matrix can be viewed as a linear operator mapping vectors from
$\mathbb{C}^{N}$ to $\mathbb{C}^{L}$. Given norms for both $\mathbb{C}^{L}$
and $\mathbb{C}^{N}$, which need not be the same, the induced norm of a matrix
is defined as

$\|{{\bf H}}\|=\max_{\|x\|=1}\|{{\bf H}}x\|\,.$ (7)

This norm characterizes the maximal possible “amplification” from the
application of ${\bf H}$ to $x$. More precisely: the value of $\|{{\bf H}}\|$
is the smallest positive number $c$ such that

$\|{{\bf H}}x\|\leq c\|x\|\,.$ (8)

Notice that, following standard practice, we use the same notation – namely,
$||\cdot||$ – for both the norm of a vector ($\|x\|$ and $\|{\bf H}x\|$ in the
above definition) and a matrix ($\|{\bf H}\|$). It should be clear from the
context, though, which one we are referring to.

##### Rank and Kernel

The range of a complex $L\times N$ matrix ${\bf H}$ is

$\mbox{range}({\bf H})=\\{y\in\mathbb{C}^{L}|\,y={\bf H}x\;\;\mbox{for some
}x\in\mathbb{C}^{N}\\}\,,$

its rank

$\mbox{rank}({\bf H})=\dim\left(\mbox{range}({\bf H})\right)\,,$

its kernel

$\mbox{ker}({\bf H})=\\{x\in\mathbb{C}^{L}|\,{\bf H}x=0\\}\,,$

and its nullity

$\text{nul}({\bf H})=\text{dim}(\text{ker}({\bf H}))\,.$

As a linear algebra exercise it can be shown that

$\text{rank}({\bf H})+\text{nul}({\bf H})=N\,.$

For square $N\times N$ matrices, the following properties are equivalent:

  1. 1.

${\bf A}$ is non-singular

  2. 2.

$\mbox{det}({\bf A})\neq 0$

  3. 3.

$\mbox{ker}({\bf A})=\\{{\bf 0}\in\mathbb{C}^{N}\\}$

  4. 4.

$\mbox{rank}({\bf A})=N$

  5. 5.

${\bf A}$ has linearly independent columns and vectors.

#### Approximation by projection

##### The least squares problem

It is rather easy to approximate one function by other functions, such as
Fourier, wavelets, polynomial expansions, or somewhat physically inspired
bases. In this article we focus on representing gravitational waves by
themselves. Intuitively, this should be (and is) more efficient. In this
subsection we briefly introduce projection-based approximations, with more
details given in Section 6.

To motivate the problem consider an $n$-dimensional vector space $W_{n}$ which
is itself a subspace of a Hilbert’s one ${\cal H}$ (for the purposes of this
article, this means that a scalar product is assumed to exist). A common
approximation criteria is a least squares (LS) one. That is, one seeks to
approximate $h\in{\cal H}$ by $h^{(n)}\in W_{n}$ which is the solution to

$h^{(n)}={\tt argmin}_{\tilde{h}}\|h-\tilde{h}\|^{2}\,,$ (9)

where $\tilde{h}\in W_{n}$. Equation (9) means finding the element of $W_{n}$
that minimizes the squared norm of the representation error.
If $\\{e_{i}\\}_{i=1}^{n}$ is an orthonormal basis of $W_{n}$ and
$\langle.,.\rangle$ the scalar product associated with ${\cal H}$, the
following is the (unique) solution to the LS problem,

$h^{(n)}=\sum_{i=1}^{n}c_{i}e_{i}\,,\quad\text{with}\;c_{i}=\langle
e_{i},h\rangle.$ (10)

The solution $h^{(n)}$ is the orthogonal projection of $h$ onto $W_{n}$,
therefore denoted as

${\cal P}_{n}h:=h^{(n)}\,.$

What it means is that if $h\in W_{n}$ then

${\cal P}_{n}h=h\,,$ (11)

and the residual $\delta h={\cal P}_{n}h-h$ satisfies

$\langle{\cal P}_{n}h,\delta h\rangle=0\,.$ (12)

The solution in 10 is basis-independent and $h^{(n)}$ is uniquely defined. The
orthogonal projection onto any linear space is a geometric construction, so it
is independent of the basis used to represent it: If $\\{e_{i}\\}_{i=1}^{n}$
and $\\{\tilde{e}_{i}\\}_{i=1}^{n}$ are any two orthonormal bases, and the
projection coefficients $c_{i}$ and $\tilde{c}_{i}$ are computed according to
(10) (replacing in the second case $e_{j}$ by $\tilde{e}_{j}$), then

${\cal
P}_{n}h=\sum_{i=1}^{n}c_{i}e_{i}=\sum_{i=1}^{n}\tilde{c}_{i}\tilde{e}_{i}\,.$
$\sum_{i=1}^{n}\langle e_{i},v\rangle
e_{i}=\sum_{i=1}^{n}\langle\tilde{e}_{i},v\rangle\tilde{e}_{i}\,,$ (13)

leading to the same solution of the LS problem. One might ask why the emphasis
on orthogonal or orthonormal bases; it is due to a conditioning issue which we
discuss this in Section 6. For the time being the summary is that one should
always use orthogonal or orthonormal bases. An exception is when building
predictive models, as discussed in Section 11.

##### Representing a set of elements: collective error

Suppose we seek to approximate not one function but a set of them:

${\cal F}=\\{h_{1},h_{2},\dots,h_{N}|h_{i}\in{\cal H}\\}$ (14)

From the above discussion in Section 3 we know that ${\cal P}_{n}h_{i}$ is the
the best representation of any particular $h_{i}$. The question of how one
measures the approximation error of the collective set naturally arises. Two
reasonable and standard criteria are

$\frac{1}{N}\sum_{i=1}^{N}\|h_{i}-{\cal P}_{n}h_{i}\|^{2}\,,$ (15)

and

$\max_{i}\|h_{i}-{\cal P}_{n}h_{i}\|^{2}\,,$ (16)

which have the interpretations of the mean and maximum errors, respectively.
Whether the notion of error ${\rm(1)}$ or ${\rm(2)}$ is chosen for an error
minimization criteria will lead to the main two ROM approaches discussed in
this review. Namely, Proper Orthogonal/Singular Value Decompositions, or the
Reduced Basis-greedy approach.

#### Further reading

The material discussed can be found in standard linear algebra/numerical
analysis books. There are many introductory good ones, we suggest [129] and
[130].

## Part I Representation and Compression

###  4 Principal Component Analysis

Principal Component Analysis (PCA) is perhaps one of the most used tools when
seeking for redundancy or a hierarchy of importance of variables in
statistical analysis. It has a close relationship with Proper Orthogonal
Decomposition (POD) as we discuss in Section 5.1. With this approach one seeks
to determine the most statistically relevant variables and to potentially
dimensionally discard from the problem the least relevant ones. In order to do
this, we recall the definition of the covariance between two stochastic
variables $X,Y$ as given by

$\displaystyle{\rm Cov}(X,Y)$ $\displaystyle=\big{\langle}(X-\langle
X\rangle)(Y-\langle Y\rangle)\big{\rangle}$ $\displaystyle=\langle
XY\rangle-\langle X\rangle\langle Y\rangle\,,\,$ (17)

where the brackets denote expectation values. The covariance between two
variables provides a measure of the degree to which their fluctuations are
correlated. A smaller (larger) covariance implies lower (higher) correlation.
In particular, the covariance of a variable with itself is its variance (i.e.,
the standard deviation squared) and measures deviations from the mean value.

When there are multiple stochastic variables $X_{i}$ ($i=1,\ldots,n$) one can
construct their associated covariance matrix $\bf{C}$ with components
$C_{ij}={\rm Cov}(X_{i},X_{j})$. This matrix is symmetric, non-negative
definite and can therefore be diagonalized with an orthogonal transformation.
Consider the $i^{\rm th}$ normalized eigenvector $\hat{{\bf V}}_{i}$. If we
set ${\vec{\bf X}}:=(X_{1},\ldots,X_{n})$ then the principal components (PCs)
are the associated eigenmodes,

${\cal E}_{i}={\vec{\bf X}}\cdot\hat{\bf V}_{i}\,.$ (18)

The PCs are new variables representing directions in which the $X_{i}$ have
different levels of variance. They are uncorrelated, a consequence of the
orthogonality of the eigenvectors, and their associated eigenvalues $\lambda$
are their variances,

${\rm Cov}({\cal E}_{i},{\cal E}_{j})=\lambda_{i}\delta_{ij}\,.$ (19)

The fact that, by construction, principal components are uncorrelated with
each other is important since they provide independent pieces of statistical
information.

The smaller an eigenvalue $\lambda_{i}$, the more likely that the
corresponding linear combination ${\cal E}_{i}$ will not deviate from its
average value for a randomly set of variables. Therefore, if there exist small
eigenvalues then the associated PCs are largely conserved in a statistical
sense. Conversely, the larger an eigenvalue is then the more relevant the
associated PC is in describing the dynamics and variations in the problem.

There are two related but different senses in which for a parametrized time
series a principal component with small variance can be semi-conserved. The
first is being constant as a function of time for an arbitrary but fixed set
of parameters. The second one is in the statistical sense that deviations of a
principal component from the mean value are small for arbitrary but fixed
initial and final times over a set of runs with the initial configurations. A
small variance automatically implies approximate conservation in the second
sense but not necessarily in the first one. The interest is not only in those
PCs which have the smallest variances (and thus identify semi-conserved
quantities in the second sense) but also in those with the largest variances,
which encode the most information about the system dynamics. This should be
clarified through the following example.

Example 3.  PCA for spin dynamics.

One may wonder if given a uniform initial spin orientation distribution, after
a while those orientations turn into some preferred orientation, such as being
aligned into some preferred direction. This was studied in [57] through
massive numerical simulations. We only briefly review some of the results of
that reference. The interest at the time was motivated by the unexpectedly
large kicks found in numerical relativity simulations and whether they were of
a generic nature in a statistical sense; for a review on recent comprehensive
studies on this topic see Section 12.1. In order to perform a large enough
analysis the Post-Newtonian equations were used, up to 3.5PN order in the
angular frequency and 2PN with the covariant spin supplementary condition.

Next, in any PCA study one has to define which variables to analyze. In order
for the quantities to be invariant under rotations of the system of reference,
at least in a Newtonian sense, the first quantities chosen were scalar
products between the normalized spins of each binary component and the orbital
angular momentum of the system. Moreover, their differences between initial
and final values:

$\displaystyle\Delta({\bf\hat{S}}_{1}\cdot{\bf\hat{L}})$ $\displaystyle=$
$\displaystyle{\bf\hat{S}}_{1}\cdot{\bf\hat{L}}|_{f}-{\bf\hat{S}}_{1}\cdot{\bf\hat{L}}|_{i}=:{{\Delta}_{1\text{L}}}$
(20) $\displaystyle\Delta({\bf\hat{S}}_{2}\cdot{\bf\hat{L}})$ $\displaystyle=$
$\displaystyle{\bf\hat{S}}_{2}\cdot{\bf\hat{L}}|_{f}-{\bf\hat{S}}_{2}\cdot{\bf\hat{L}}|_{i}=:{{\Delta}_{2\text{L}}}$
(21) $\displaystyle\Delta({\bf\hat{S}}_{1}\cdot{\bf\hat{S}}_{2})$
$\displaystyle=$
$\displaystyle{\bf\hat{S}}_{1}\cdot{\bf\hat{S}}_{2}|_{f}-{\bf\hat{S}}_{1}\cdot{\bf\hat{S}}_{2}|_{i}=:{{\Delta}_{12}}\,,$
(22)

where hats stand for unitary vectors. The choice of initial and final values
is something to define and is discussed in [57], the summary is that the main
results do not depend on these choices, which gives insight into the fact that
the binary problem in GR is highly redundant.

The orbital angular momentum and spin orientations naturally become correlated
due to spin-orbit and spin-spin interactions as each of these binary black
hole configurations evolve in time. However, at least within the PN
approximation here considered, the orbital angular momentum and spin vectors
remain perfectly uniformly distributed [68]. For example, a Kolmogorov-Smirnov
test for a representative configuration returns a p-value of $\sim 10^{-5}$
when testing for lack of uniformness [68]. Higher PN expansions might
introduce small biases [89] but if so they appear to be at a level in which
approximating the mean of the above scalar products at any instant of time by
zero is a very good approximation.

Ref. [57] starts with a simple case and makes contact with previous
conservation results. The authors start building towards the more general case
by first doing a PCA using only the two spin-orbit (SO) variables in (20) and
(21),

${{\Delta}_{1\text{L}}}=\Delta({\bf\hat{S}}_{1}\cdot{\bf\hat{L}})\,,\quad{{\Delta}_{2\text{L}}}\mathrel{\mathop{:}}=\Delta({\bf\hat{S}}_{2}\cdot{\bf\hat{L}})\,.$
(23)

However, spin-spin interactions in both the numerical simulations and in the
analytical calculations are included in the PN equations of motion used when
solving for the evolution of each configuration.

For mass and spin magnitudes ($m_{j},\chi_{j}$) of each black hole the
covariance matrix for the variables (23) is

${\bf C}=\left(\begin{array}[]{cc}{\rm
Cov}({{\Delta}_{1\text{L}}},{{\Delta}_{1\text{L}}})&{\rm
Cov}({{\Delta}_{1\text{L}}},{{\Delta}_{2\text{L}}})\\\ {\rm
Cov}({{\Delta}_{2\text{L}}},{{\Delta}_{1\text{L}}})&{\rm
Cov}({{\Delta}_{2\text{L}}},{{\Delta}_{2\text{L}}})\end{array}\right)\,,$ (24)

where the entries can come either from numerical simulations or from what the
authors call the instantaneous approximation. The matrix $\bf C$ is then
diagonalized to find the principal components.

From numerical simulations the authors find that, sampling across many random
initial spin orientations, each of the principal components has zero mean over
time (to numerical accuracy), $\langle\Delta{\cal
E}_{j}^{\text{SO}}\rangle=0$, a consequence of the spin orientation
distributions remaining highly uniform during the inspiral.

Furthermore, they find $\lambda_{2}$ to be in the range $\sim 10^{-9}-10^{-4}$
for the parameters sampled and that it grows with both spin magnitudes, which
is expected from physical intuition, but also that it increases as the equal
mass case is approached, which was unexpected.

As an example with $(m_{1},m_{2},\chi_{1},\chi_{2})=(0.4,0.6,1.0,1.0)$, Figure
1 shows a graphical representation of the principal components overlaid on a
scatter plot of the ${{\Delta}_{1\text{L}}}$ and ${{\Delta}_{2\text{L}}}$ data
from $1,\\!000$ out of $100,\\!000$ numerical simulations using random initial
spin orientations. Notice that the first PC, which points along the direction
of the eigenvector ${\bf\hat{V}}_{1}$ with the largest eigenvalue
$\lambda_{1}$, captures the largest variation in the data while the second PC,
pointing along ${\bf\hat{V}}_{2}$, indicates that there is very little spread
in the data in that direction, which is also implied by the smallness of
$\lambda_{2}$ relative to $\lambda_{1}$. Therefore, for the time interval
considered, the second PC is largely irrelevant. This figure is almost an
ideal example of dimensionality reduction through PCA. For a detailed analysis
see [57].

Figure 1: A graphical representation of the principal components for the spin-
orbit variables and the numerical data of ${{\Delta}_{1\text{L}}}$ and
${{\Delta}_{2\text{L}}}$ for a binary black hole system with, as illustration,
$m_{1}=0.4$ and maximal spin magnitudes. One can see that PC2 is largely
irrelevant. For details see [57].

#### Further Reading

Principal Component Analysis is a very well known tool to look for
redundancies or categorize the relevance of different variables. One of its
weaknesses is that it depends on which variables to look for. The case study
that we chose is due to the fact that if looking in more detail at the
original reference reviewed, there are strong indications that there are at
least three (out of eight) redundant quantities in the problem, and therefore
the problem is amenable to an unsupervised dimensional reduction approach, as
presented throughout this review. Reference [76] is a very pedagogical book
solely on PCA, with many case studies and modern developments.

###  5 Proper Orthogonal Decomposition

Given the snapshots $H_{1},\ldots,H_{N}\in\mathbb{C}^{L}$ (the training
samples), the snapshot matrix is

${\bf H}=[H_{1},\ldots H_{N}]\in\mathbb{C}^{L\times N}$

with, in general, $L\geq N$. The goal here is to find $n\leq N$ orthonormal
vectors $\\{\phi_{i}\\}_{i=1}^{n}$ in $\mathbb{C}^{L}$ minimizing the average
approximation error (15)

$\frac{1}{N}\sum_{i=1}^{N}\|H_{i}-{\cal P}_{n}H_{i}\|^{2}\,,$ (25)

where ${\cal P}_{n}$ is the orthogonal projector onto
$W_{n}:=\text{span}\\{\phi_{i}\\}_{i=1}^{n}$ and $\|\cdot\|$ is the Euclidean
norm. The minimizer subspace of (25) is called a Proper Orthogonal
Decomposition (POD) of rank $n$, and its generating basis the POD one.

Notice that in general the POD basis elements $\phi_{i}$ are not members of
$\\{H_{i}\\}_{i=1}^{N}$. In the gravitational wave case, this means that the
reduced basis is not a subset of waveforms, neither does the method provide a
set of the “most representative” points in the space of waveforms for building
a posteriori bases through numerical relativity simulations.

Another disadvantage is that the problem, involving a potentially very large
linear algebra system, is not easily parallelizable in distributed memory
architectures. Still, the method is straightforward to implement and can, for
example, quickly provide insight into whether a problem is amenable to
dimensionality reduction: after performing a POD decomposition of the snapshot
matrix, one can look at the decay rate of singular values as a function of $n$
in order to know if the problem can be codified in a space of lower
dimensionality than the original one. This can be turned into an actual
strategy: a POD for a small subset of the problem can provide a quick insight,
and if the singular values decay fast enough with $n$, the larger problem of
interest can be tackled with a greedy approach, as discussed in Section 8.3.

The optimal solution to the minimization problem defined by Eq. (25) can be
accomplished by means of a Singular Value Decomposition (SVD) procedure. In
this framework the POD basis is given by the first $n$ left singular vectors
of the snapshot matrix $\bf{H}$ (this can be proved by considering the first-
order optimality conditions for the minimization problem).

Consider then the following singular value decomposition of H:

$\bf{H}=\bf{U\Sigma V}^{\dagger}\,,$ (26)

where ${\bf U}=[u_{1},\ldots,u_{L}]\in\mathbb{C}^{L\times L}$, ${\bf
V}=[v_{1},\ldots,v_{N}]\in\mathbb{C}^{N\times N}$ are orthogonal matrices and

${\bf\Sigma}=\begin{bmatrix}{\bf D}\\\ {\bf
0}\end{bmatrix}\in\mathbb{C}^{L\times N},{\bf
D}=diag(\sigma_{1},\ldots,\sigma_{N})\in\mathbb{C}^{N\times
N},\sigma_{1}\geq\ldots\geq\sigma_{N}\geq 0.$ (27)

Then, the rank-$n$ POD basis is given by $\phi_{i}=u_{i},\,i=1,\ldots,n$. This
basis provides a low-rank approximation to H and represents the minimizer
basis for the minimization problem stated in (25). To see this, consider from
(26) the relation

$u_{i}^{\dagger}{\bf H}=\sigma_{i}v_{i}^{\dagger}\,,$

valid for all column vectors of ${\bf U}$. Next, multiply by $u_{i}$ on the
left and sum

$\sum_{i=1}^{n}u_{i}u_{i}^{\dagger}{\bf
H}=\sum_{i=1}^{n}\sigma_{i}u_{i}v_{i}^{\dagger}\,.$ (28)

The l.h.s. of (28) is exactly the orthogonal projector associated to the basis
$\\{u_{i}\\}_{i=1}^{n}$ acting on $\bf H$ and the r.h.s. represents its
rank-$n$ approximation ${\bf H}_{n}$. Therefore we can rewrite (28) as

${\cal P}_{n}{\bf H}={\bf H}_{n}\,.$ (29)

Remarks on the SVD decomposition

  * •

The numbers $\sigma_{i}$ are known as the singular values of the matrix H and
correspond to the positive square roots of the eigenvalues of the associated
matrix $\bf{K}=\bf{H}^{\dagger}\bf{H}$. They are usually chosen in descending
order in practice to facilitate the recognition of the most relevant principal
components.

  * •

For the matrices U and V, only the first $\mathrm{rank}(\mathbf{H})$ columns
are unique, whereas the remaining ones are arbitrarily extended such that
orthogonality is maintained. Since $\sigma_{i}=0$ whenever
$i>\mathrm{rank}(\mathbf{H})$ the factorization (26) does not depend on the
choice of this extension, since such extension is annihilated by null entries
of $\bf\Sigma$.

  * •

The first $N$ columns of U and V are known as the left- and right-singular
vectors of H respectively. The right-singular vectors, $v_{i}$, are the
normalized eigenvectors of ${\bf H}^{\dagger}{\bf H}$ and the left-singular
vectors, $u_{i}$, are the normalized eigenvectors of ${\bf H}{\bf
H}^{\dagger}$.

  * •

The number $K$ of non-zero elements in $\mathbf{\Sigma}$ is exactly
$\text{rank}(\mathbf{H})$. Consequently, Eq. ((26)) is sometimes called a
rank-revealing factorization.

Error of a POD approximation

The accuracy of the low rank-$n$ approximation ${\bf H}_{n}={\cal P}_{n}{\bf
H}=[{\cal P}_{n}H_{1},\ldots,{\cal P}_{n}H_{N}]$ is given by the following
lemma.

Lemma 1. Define the approximation error by

$\epsilon=\frac{1}{N}\sum_{i=1}^{N}\|H_{i}-{\cal P}_{n}H_{i}\|^{2}.$ (30)

It can be shown that $\epsilon$ satisfies

$\epsilon=\frac{1}{N}\sum_{i=n+1}^{N}\sigma_{i}^{2}.$ (31)

In terms of the Frobenius norm,

$\epsilon=\frac{1}{N}\|{\bf H}-{\bf H}_{n}\|^{2}_{F}.$ (32)

Intuitively, the square of the Frobenius norm of the difference between a
matrix and its low-rank approximation represents the total squared difference
of the rows of ${\bf H}$ due to omitting the last $(N-n)$ singular values when
forming ${\bf H}$. As discussed above, the projection ${\bf H}_{n}$ represents
the best rank-$n$ approximation of H in the Frobenius norm.

In summary, going back to section 3, one can see that the minimization problem
(1) related to the average error of a set of functions ${\cal F}$ is optimally
solved by the POD/SVD decomposition of the matrix associated to those
functions.

In the context of a gravitational-wave template bank the rows of the snapshot
matrix are the waveforms evaluated at different time/frequency values, and the
different rows correspond to different intrinsic parameters. The inner product
of rows with themselves, in turn, correspond to the total power in the bank.
Thus, the squared Frobenius norm of a template bank is the total power
contained within all the templates in it. For normalized templates, this is
simply the size of the bank. It can be seen how Eq. (32) directly corresponds
to the total power in a template bank lost from a low-rank approximation to
it. More precisely, the error measure is directly related to the average
fractional SNR loss which, up to a constant, is the squared Frobenius norm of
the difference between the full template bank and its low-rank approximation.
Fixing the average fractional SNR loss thus determines the total number of
non-zero singular values which must be retained to guarantee that the rank-
reduced bank remains effective. This is discussed in practical terms through
the following case study.

Example 4.  SVD for gravitational waves.

In Ref. [35] the authors studied the application of SVD to gravitational wave
templates to reduce the redundancies in the bank and build an orthogonal basis
to represent the whole set. As a proof of concept, the authors applied a SVD
approach to a set of CBC waveforms corresponding to a sliver of the BNS
parameter space. Next, we summarize some of the results of this reference,
closely following its notation, which might be different from the rest of this
review.

In order to detect a GW signal, the common choice of the minimal match between
an arbitrary point in parameter space and its nearest point of the template
bank is $97\%$.

In order to compare and filter the data against the entire template bank, an
approximation to the matched filtering
$\rho_{\alpha}=\langle\bar{h}_{\alpha},s\rangle$ is sought for, where
$h_{\alpha}$ is a complex waveform vector and $s$ is the data vector (the
presumed signal). In this way the number of evaluation of inner products can
be reduced as well as its computational cost. Let’s define the $N\times L$,
where $N=2M$ and $M$ stands for the number of complex waveform (there is no
obvious reason why to perform an SVD on real and imaginary parts of a
waveform, given that the POD approach can handle complex snapshots) template
matrix

${\bf H}=\\{h_{1}^{R},h_{1}^{I},\ldots,h_{M}^{R},h_{M}^{I}\\}\,,$ (33)

where $h_{i}^{R,I}$ are the real and imaginary part of the $i$-template
waveform, each one corresponding to the rows of $\bf H$. Applying an SVD
decomposition to ${\bf H}$ and writing it in component form (here we follow
the notation of [35])

${\bf
H}_{\mu\nu}=\sum_{\kappa=1}^{N}v_{\mu\kappa}\sigma_{\kappa}u_{\kappa\nu}\,,$
(34)

one can define the truncated sum

${\bf
H}^{\prime}_{\mu\nu}=\sum_{\kappa=1}^{N^{\prime}}v_{\nu\kappa}\sigma_{\kappa}u_{\kappa\nu}\,,$
(35)

and approximate $\rho_{\alpha}$ by

$\rho^{\prime}_{\alpha}=\langle
H^{\prime}_{2\alpha-1}-iH^{\prime}_{2\alpha},s\rangle=\sum_{\nu=1}^{N^{\prime}}(v_{{2\alpha-1}\,\nu}\sigma_{\nu}-iv_{{2\alpha}\,\nu}\sigma_{\nu})\langle
u_{\nu},s\rangle\,,$ (36)

where the $H^{\prime}_{j}$ are the vector rows of ${\bf H}^{\prime}$.

Figure 2 shows a representation of the matrix of waveforms H and its
associated SVD-basis. One disadvantage of this kind of dimensional reduction
is that the remaining basis barely resembles the structure of the original
template. This can be fixed with the RB-greedy approach that is presented in
Section 8.

Figure 2: Left: Representation of template waveforms corresponding to $\bf H$
against time. Right: Representation of the first four SVD-basis. Notice that
since the SVD basis functions do not correspond to actual waveforms, they
display non-physical behavior which might be difficult to fit for when
building a surrogate model as discussed in Section 11. Figures taken from
[35].

Accuracy. An approximation to the data snapshot matrix is sought such that
$\|{\bf H}_{\mu}-{\bf H}^{\prime}_{\mu}\|\sim 1\%$. Following this
requirement, the fractional SNR

$\Big{\langle}\frac{\delta\rho}{\rho}\Big{\rangle}:=\frac{1}{M}\sum_{\alpha=1}^{M}\frac{\delta\rho_{\alpha}}{\rho_{\alpha}}$
(37)

can be approximated by:

$\Big{\langle}\frac{\delta\rho}{\rho}\Big{\rangle}=\frac{1}{2N}\sum_{\mu=N^{\prime}+1}^{N}\sigma_{\mu}^{2}\,,$
(38)

resulting from a Taylor expansion valid in the range
$\langle\delta\rho/\rho\rangle<10\%$. Note that this approximation is
proportional to the squared Frobenius norm of the truncation error of ${\bf
H}$,

$\Big{\langle}\frac{\delta\rho}{\rho}\Big{\rangle}=\frac{1}{2N}\|{\bf H}-{\bf
H}^{\prime}\|_{F}^{2}\,.$ (39)

As a case of application, [35] shows an SVD analysis to gravitational waves
emitted by a CBC-BNS, with chirp masses $1.125M_{\odot}\leq
M_{c}<1.240M_{\odot}$ (that is, a rather small sliver in parameter space) and
component masses $1M_{\odot}\leq m_{1},m_{2}<3M_{\odot}$. In order to satisfy
a minimal match of $96.8\%$, a number of templates $M=456$ ($N=912$) to cover
the parameter space was found. In Fig. 3 (plot of
$\langle\delta\rho/\rho\rangle$ vs. # of SVD-basis elements) it can be seen
that, to obtain $\langle\delta\rho/\rho\rangle=10^{-3}$, the number of basis
elements needed to reconstruct the whole template bank to that accuracy can be
reduced from $N=912$ to $N^{\prime}=118$. Though POD/SVD is a good starting
point for ROM, in following sections we will develop a more modern framework
for modeling reduction.

Figure 3: SNR (Eq. (37)) for a CBC-BNS gravitational waves as a function of
the number of SVD-basis elements – valid for $SNR<10\%$ –. Figure taken from
[35].

####  5.1 PCA and POD

Principal Component Analysis and Proper Orthogonal Decomposition are closely
related from a mathematical point of view, though usually applied in different
contexts. PCA is in general used in statistical analysis. It represents a
solution to an optimization problem; namely, that one of finding the
uncorrelated directions of maximum statistical variability in a stochastic
data set. POD is usually applied in dimension reduction modeling. As in PCA,
POD also represents a solution to an optimization problem: finding the best
low-rank approximation to a data set matrix. This is achieved by minimizing
the Frobenius norm of the error matrix (Eq. (32)).

Suppose we have a data set matrix of size $L\times N$, $N$ being the number of
variables and $L$, the number of samples corresponding to each variable. In
PCA one is not interested in the “length” of each variable and their
individual components, but rather on their correlations. Put differently, if
there are $N$ stochastic variables, the covariance matrix, regardless of the
number of components in each stochastic variable, is in all cases of size
$N\times N$. Compression can be achieved by ignoring the directions in which
the stochastic data have little variance.

Mathematical relationship between PCA and POD

Although conceptually different, both methods are mathematically related. To
see this, we make a simple observation concerning the snapshot matrix and its
covariance matrix. Let ${\bf X}:=[X_{1},\ldots,X_{N}]\in\mathbb{C}^{L\times
N}$ be the snapshot matrix with each column $i$ storing $L$ independent
observations of the random variable $X_{i}$. We assume that each $X_{i}$ is
identically distributed. Subtracting to each column $X_{i}$ the corresponding
mean value $\langle X_{i}\rangle$, $X_{s,i}\rightarrow X_{s,i}-\langle
X_{i}\rangle\,,\,s=1,\ldots,L$, and replacing $\bf{X}$ by this “centered”
matrix version, the associated covariance matrix $\bf{C}$ can be written as

${\bf C}=\frac{1}{L-1}{\bf X}^{\dagger}{\bf X}\,.$ (40)

Without loss of generality, we will suppose $\bf X$ has full rank. Therefore,
being an hermitian, non-negative definite matrix, $\bf C$ can be diagonalized
as

${\bf C}={\bf V}{\bf\Lambda}{\bf V}^{\dagger}\,,$

where ${\bf\Lambda}=\text{diag}(\lambda_{1},...,\lambda_{N})$ is the matrix of
non-zero descending-ordered real eigenvalues and $\bf V$ is the unitary matrix
storing eigenvectors of $\bf C$. The columns of $\bf V$ are called principal
axes or principal directions in the PCA framework. PCs correspond to the
projection of the data onto these principal directions, the columns of the
matrix product ${\bf X}{\bf V}$.

Now, perform a POD/SVD decomposition of the snapshot matrix X

${\bf X}=\mathbf{U}\mathbf{\Sigma}\mathbf{\tilde{V}}^{\dagger}$

and rewrite Eq. (40) as

${\bf
C}=\frac{1}{L-1}\,\,(\mathbf{U}\mathbf{\Sigma}\mathbf{\tilde{V}}^{\dagger})^{\dagger}\mathbf{U}\mathbf{\Sigma}\mathbf{\tilde{V}}^{\dagger}=\mathbf{\tilde{V}}\,\frac{\mathbf{\Sigma}^{\top}\mathbf{\Sigma}}{L-1}\,\mathbf{\tilde{V}}^{\dagger}\,.$
(41)

The similarity relation between $\bf C$ and
$\frac{\mathbf{\Sigma}^{\top}\mathbf{\Sigma}}{L-1}$ implies that both have the
same eigenvalues:

$\lambda_{i}=\frac{\sigma^{2}_{i}}{L-1}\,,\quad i=1,\ldots,N\,.$

Another benefit of this mathematical correspondence is that PC’s are
straigthforward to compute as ${\bf X}{\bf V}={\bf U}\mathbf{\Sigma}$.

#### Further reading

The reader may consult [76, 140] for clear expositions about the classical
matrix factorization methods presented in this section. QR decompositions form
another class of matrix factorizations which can be used for low-rank
approximations. A theoretical framework for rank revealing QR decompositions
of a matrix was developed in the the late 1980’s and 1990’s [38, 73],
motivated partly by the SVD’s high computational cost. QR decompositions form
the basis for many modern fast algorithms [66, 64, 47, 42]. Actually,
motivated by GW research an important resuIt has been shown: the RB-greedy
algorithm is completely equivalent to a certain type of QR decomposition [16].

###  6 Spectral expansions

In this section we discuss fast converging classical linear approximations,
since they naturally lead to Reduced Basis when considering parametrized
problems.

Consider a set of complex functions,

${\cal F}:=\\{h:\mathbb{C}\rightarrow\mathbb{C}\\}\,.$

A standard linear approximation consists of the following: a basis with $n$
(non necessary orthogonal) functions is somehow chosen,

$\\{e_{i}\\}_{i=1}^{n}\,.$

Next, a function $h$ in ${\cal F}$ is approximated by a linear combination of
the chosen basis elements,

$h(z)\approx h^{(n)}(z):=\sum_{i=1}^{n}c_{i}e_{i}(z)\,,\quad
c_{i}\in\mathbb{C}.$ (42)

There are several criteria to choose the expansion coefficients $\\{c_{i}\\}$.
We recap least-squares (LS), in Section 10 we will discuss interpolation. In
the LS approach the coefficients $\\{c_{i}\\}$ are chosen such that the
representation error is minimized with respect to some chosen norm,

$\left\|h-\sum_{i=1}^{n}c_{i}e_{i}\right\|^{2}\,.$ (43)

As was stated in Section 3, the solution to the LS problem (Eq. (9), Section
3) is the orthogonal projection onto the span of the basis,

$h^{(n)}={\cal P}_{n}h\,.$ (44)

That is, the approximation (42) minimizes the error (43) when the coefficients
$\\{c_{i}\\}$ are chosen such that the residual

$\delta^{(n)}h=h-h^{(n)}$

satisfies $\langle\delta h,h^{(n)}\rangle=0$, which implies

$\left\langle
e_{j}(\cdot),h(\cdot)-\sum_{i=1}^{n}c_{i}e_{i}(\cdot)\right\rangle=0\quad\quad\text{
for }j=1,\ldots,n\,.$ (45)

The solution to (45) is

$c_{i}=\sum_{j=1}^{n}({\bf G}^{-1})_{ij}\langle e_{j},h\rangle\,,$ (46)

where ${\bf G}^{-1}$ is the inverse of the Gram matrix or Gramian ${\bf G}$,
with entries

${\bf G}_{ij}:=\langle e_{i},e_{j}\rangle\,.$

If the basis is orthonormal, this matrix is the identity and one recovers the
familiar expression

${\cal P}_{n}h=\sum_{i=1}^{n}\langle e_{i},h\rangle e_{i}\,.$ (47)

The Gramian matrix can be very ill-conditioned in general, meaning that the
calculation of its inverse, needed in (46), can have large numerical errors
(see, for example [136]). Therefore, from a numerical conditioning point of
view, it is convenient to work with an orthonormal (or orthogonal) basis since
the Gram matrix is the identity.

In iterative approximations one defines the convergence rate as the rate at
which the representation error

$\|\delta^{(n)}h\|=\left\|h-{\cal P}_{n}h\right\|$ (48)

decreases as $n$ increases for any given $h$. In the context of RB, when one
is dealing with parametrized systems, this error will depend also on the
parameters of the system.

####  6.1 Spectral methods

In terms of accuracy and optimal convergence rates for approximation of a
space ${\cal F}$, we have not yet discussed two related aspects:

  1. 1.

The choice of a “good” basis or, more precisely, the approximation space
$W_{n}$.

  2. 2.

The optimal choice of how many basis elements to use. One would think that the
larger, the better. This is related to the regularity of the functions in
${\cal F}$. That is, how smooth or differentiable they are.

The first comment might be puzzling, after all we have emphasized that given
an approximation space $W_{n}$, the LS approximation is uniquely defined in a
geometric way. The question really is what the approximation space $W_{n}$
should be to minimize the representation error.

#####  6.1.1 Fourier Expansions

We start with the simplest and most familiar case: that one of periodic
functions in $[0,2\pi]$, unit weight, $\omega(x)=1$ 333The choice of norm is
crucial, see the next Section about Jacobi polynomials for more insights about
this., and Fourier modes as (orthonormal) basis,

$e_{j}(x)=\frac{1}{\sqrt{2\pi}}e^{ijx}\,,\quad j\in\mathbb{Z}$

Assuming for simplicity that $n$ is even, we have

${\cal P}_{n+1}h(x)=\sum_{-n/2}^{n/2}\hat{h}_{j}e_{j}(x)\,,$

where $\hat{h}_{j}$ are the Fourier coefficients of $h$.

One can show that if $h$ has $s$ derivatives, then there exists a constant
$C>0$ independent of $n$ such that

$\|h-{\cal P}_{n+1}h\|\leq C(n+1)^{-s}\left\|\frac{d^{s}h}{dx^{s}}\right\|$
(49)

for all $n\geq 0$. In particular, if $h$ is smooth ($h\in C^{\infty}$), then
the representation error decays faster than any power law with
$n\rightarrow\infty$, which is referred to as spectral convergence. Under
further conditions, in particular the case in which $h$ is analytic, the error
decay is actually exponential. The Fourier case gives a very intuitive way of
how this happens just by using integration by parts, for more details see
Chapter 9 of [120].

The summary here is that Fourier modes are not a good basis just because they
are periodic functions, but they provide fast convergence as a representation
space, as fast as the regularity of the function(s) being represented has.
This is the core idea of spectral methods and, for parametrized systems,
reduced basis. A brief summary of the non-periodic case follows.

#####  6.1.2 Jacobi Polynomials

In the case of non periodic functions a similar result to that one discussed
for periodic functions holds, whether the domain is bounded or infinite.
Again, the choice of weight is crucial. Any interval can be mapped into
$[-1,1]$ or $(-1,1)$, that is why we usually refer to those intervals. The
open interval case is because in some cases (such as Chebyshev) the weights
actually diverge at the $\pm 1$ end points – there is a reason for this but we
shall skip it.

Polynomials are a natural basis to use. Why? Just because after centuries we
understand them well; going beyond that is one way of looking at Reduced
Basis. Given a maximum degree, the span of polynomials is the same, so the
question is only what kind of weights guarantee fast convergence upon
regularity of the function to be represented. The following family of weights
is a sufficient class, thought not a necessary one,

$\omega(x)=(1-x)^{\alpha}(1+x)^{\beta}\,\quad x\in(-1,1)$ (50)

with $\alpha,\beta>-1$.

Under these assumptions, the following holds in a Jacobi Polynomial
approximation:

$\|h-{\cal P}_{n+1}h\|\leq
C(n+1)^{-s}\left\|(1-x^{2})^{s/2}\frac{d^{s}h}{dx^{s}}\right\|$ (51)

for all $n>(s-1)$ and $C$ independent of $n$. If $h$ is smooth there is
asymptotic spectral convergence. Under additional conditions on the smoothness
(or, but not necessarily, analyticity), the convergence is in fact
exponential. In this review we will loosely associate smoothness with
asymptotic exponential convergence of application-specific spectral
expansions, without discussing these additional assumptions.

Most of spectral methods literature is based on solving numerical problems,
prominently differential equations. In this review we shift this focus to that
one of a more fundamental representation problem, after which all the calculus
for approximating quadratures, taking derivatives, and solving differential
equations follows in a rather straightforward way. This is not just a matter
of taste, but will make the introduction of Reduced Basis and all its
associated calculus and applications very natural and almost trivial.

#### Further reading

The weights (50) lead to Jacobi polynomials, which are solutions to a singular
Sturm-Liouville problem, the properties of which guarantee the above mentioned
spectral convergence. Standard examples are Legendre and Chebyshev
polynomials. Their span is the same, it is with respect to which scalar
product they guarantee fast convergence and their discrete version and
relation to interpolation, discussed in Section 10. They are also
orthonormalized with respect to their own scalar products, which helps avoid
the typical conditioning problem of inverting the associated Gramian matrix.
For more details on and a quick glimpse at spectral methods, from a physicist
or practitioner perspective, see Chapter 9 of Ref. [120]. Reference [139] is a
surprisingly compact and efficient book to get started with spectral methods
(despite its title, MATLAB is not necessary at all to digest the content). A
long classic reference, especially for practitioners, is [29]. The book is
legally available for free from the author’s webpage, though with some typos.
A great book, in some sense targeted at ordinary differential equations is
[56]. For a modern presentation and the latest results, from theory to current
approaches to beat Gibb’s phenomena, see [70] and references therein.

###  7 Parametrized problems and optimal approximations

We are interested in approximating some abstract space of parametrized
functions ${\cal F}$, which in general is not linear (the sum of two waveforms
does not need to be a waveform) but we assume that it can be embedded in a
Hilbert one ${\cal H}$ (for example, the one of integrable functions in the
$L_{2}$ sense). We denote the underlying space of parameters as $\Phi$, which
we assume to be compact. For example, ${\cal F}$ can be the space of
$\lambda$-parametrized solutions $u_{\lambda}(x)$ of a partial differential
equation representing the dynamics of a physical system. In this review we
place emphasis in gravitational waves, parametrized for example by the mass
and spin of each black hole in a binary collision.

The question discussed next is how well one can theoretically approximate all
of ${\cal F}$ by a set of $n$ basis elements of ${\cal H}$ in a linear and
most compact way. This leads to the Kolmogorov n-width 444There are other
widths, see [107]. $d_{n}$ of ${\cal F}$ with respect to ${\cal H}$,

$\displaystyle d_{n}=d_{n}({\cal F},{\cal
H}):=\min_{\\{e_{i}\\}_{i=1}^{n}\in{\cal
H}}\max_{\lambda\in\Phi}\min_{c_{i}\in\mathbb{C}}\bigg{\|}h(\cdot;\lambda)-\sum_{i=1}^{n}c_{i}(\lambda)e_{i}(\cdot)\bigg{\|}^{2}\,.$
(52)

Comment 1.  We explain the meaning of (52), from right to left, with respect
to the $\min,\max,\min$ properties.

  * •

The first $\min$ implies that the optimal representation with respect to the
(so far arbitrary) basis $\\{e_{i}\\}_{i=1}^{n}$ is used. That is, given a
basis, the best representation (in the induced norm $\|.\|$) is considered. We
already discussed that this is the orthogonal projection ${\cal P}_{n}$ onto
the span of the basis $\\{e_{i}\\}_{i=1}^{n}$, so we can replace (52) by

$\displaystyle d_{n}:=\min_{\\{e_{i}\\}_{i=1}^{n}\in{\cal
H}}\max_{\lambda\in\Phi}\bigg{\|}h(\cdot;\lambda)-{\cal
P}_{n}h(\cdot;\lambda)\bigg{\|}^{2}\,.$ (53)

  * •

Next, the largest error in parameter space for such a best approximation is
picked. That is the “worst best” approximation given a choice of basis.

  * •

Finally, a choice of basis (more precisely, the subspace $W_{n}$) which
minimizes this worst best error is chosen and the associated error is by
definition the n-width. The compactness of $\Phi$ is important to guarantee
that these minimum and maximum exist in the search through the parameter
space.

In other words, the $n$-width is an upper bound of “what is the best that one
can do” if one could optimally choose an optimal basis.

This is mostly a theoretical problem, since solving for such a basis is
impractical (or, actually, in most cases, intractable) from a computational
point of view because it carries combinatorial complexity (all elements of the
basis have to be simultaneously chosen). It is more of a theoretical upper
bound against which to benchmark the quality of any computable approximation
which seeks for an approximate solution to the $n$-width problem.

Parameter regularity and fast convergence

Exact expressions for the $n$-width in general can only be achieved in a few
cases (see, for example, [91]). In general, the best one can do is to set up
upper bounds under specific assumptions and infer the dependence of the
$n$-width with respect to $n$. Indeed, Kolmogorov calculated this distance for
a special class of functions [91] with its first $(r-1)$ derivatives with
respect to parameter variation absolutely continuous and obtained a power law
dependency for the $n$-width,

$d_{n}\sim n^{-r}\,.$

If the set of functions considered were $C^{\infty}$ with respect to parameter
variation – the functions themselves can be discontinuous – the rate of
convergence of its optimal basis representation would be better than any power
law. This is called spectral convergence and will be relevant in next chapters
at the moment of quantifying the fast convergence of a Reduced Basis approach.
This is exactly the case in many scenarios of GWs, since they depend smoothly
on the parameters of their sources and one therefore expects an optimal
approach to have very fast (in fact, exponential) convergence, as opposed to
random sampling, for example, for which the convergence is sublinear. This is
the main reason underlying the extreme high accuracy of very compact surrogate
models based on reduced bases. That is, one expects in most GW scenarios

$d_{n}\sim e^{-an^{b}}\,,$

and an optimal basis should be extremely accurate and compact; in fact in some
cases super-exponential convergence ($b>1$) has been found in the GW context
[53]. The question then turns into how to find an approximate basis which is
not only computable but is also nearly optimal with respect to the n-width.

It is “common knowledge” that parametrized problems with regularity with
respect to parameter variation show spectral convergence of the n-width in
practice. One might ask how this observation can be possible at all when we
have mentioned that the n-width is not computable in practice: we will return
to this point when we discuss the greedy algorithm. But in fact, up to our
knowledge there is no rigorous proof of this expectation for general
parametrized systems. But there is a very compelling argument. Asymptotic
exponential convergence is also observed in practice in GWs and, being
physicists, that is good enough for us. The argument is as follow (courtesy of
Albert Cohen):

The convergence rate of the n-width depends on the smoothness/regularity of
the functions with respect to the parameters of the problem. We now argue that
if there is regularity with respect to parameter variation, as in many cases
of interest, then one can expect fast (in fact, up to exponential or even
super-exponential) decay of the approximation error. The argument is a
spectral standard one (spectral methods were discussed in the previous
section) applied to parameter variation and is as follows:

If the parametric map

$\lambda\rightarrow h_{\lambda}(\cdot)$

is smooth enough with respect to $\lambda$, then it can be very well
approximated in some appropriate basis in the $\lambda$ parameter variable,
for instance by Fourier or Jacobi polynomials. This means that there is an
expansion of the form

$c_{1}(\lambda)h_{1}(\cdot)+\ldots+c_{n}(\lambda)h_{n}(\cdot)+\ldots\,,$ (54)

with $h_{i}\in{\cal F}$, that converges fastly towards $h_{\lambda}(\cdot)$ in
${\cal F}$ uniformly in $\lambda$.

Now, a partial sum of the form (54) is a member of the span of
$\\{h_{1},\ldots,h_{n}\\}$, which means that this linear space approximates
well all functions in ${\cal F}$.

What “well” means depends on the approximation result used. If the number of
parameters is finite and the dependence is analytic, then exponential rates of
the form $exp(-n^{1/d})$ can easily be proved. If the dependence is only
$C^{s}$, then a rate $n^{-s/d}$ can also easily be proved. This argument
establishes, then, that there are bases in which expansions of the form (54)
converge very fast, and as a consequence the n-width can only decay faster.

To summarize, for problems with smooth parametric dependence, fast convergence
in terms of greedy reduced bases can be expected and it is not surprising that
such global methods outperform local ones.

We mentioned that the functions themselves can be discontinuous, which might
be confusing so it is worthwhile explaining it in more detail, even if
qualitatively. Consider a problem of fluid dynamics, where the solutions might
develop shocks at a finite time, and imagine that one is solving the partial
differential equations for a parametrized family of initial data. This is
completely fine in terms of the fast convergence of the n-width, since the
location, shape, etc., of the shock depend smoothly on the initial data; the
time series (in this case) itself does not need to have regularity with
respect to the physical variable (time in this example) but with respect to
parameter variation.

###  8 Reduced Basis

In this section we present the RB-greedy framework for parametrized systems in
order to address the resolution of the Kolmogorov problem in a quasi-optimal
way. First mathematical conventions related to the RB-scheme are presented,
followed by a discussion of its near-optimality with respect to different
$n$-width behaviors.

####  8.1 Introduction

In previous sections we have discussed linear approximation of functions by
orthogonal projection onto the span of a basis. The latter was taken to be a
generic, problem-independent one, such as Fourier modes or polynomials. The
Reduced Basis (RB) approach is a framework for efficiently solving
parametrized problems, representing the solutions in a compact way, and
predicting new ones based on an offline-online decomposition.

In parametrized problems one is interested in functions of the form

$h_{\lambda}=h_{\lambda}(\cdot)=h(\cdot;\lambda)\,.$

where is (in general a multidimensional) parameter $\lambda$, the
discretization of which will define the training space.

In fields related to scientific computing and data science one is usually
interested in multiple evaluations and analyses of functions in real time. The
approach of RB is especially tailored to that one in which some numerical
problem has to be solved in order to obtain each function, and obtaining such
numerical solutions is very expensive. The approach is also very powerful when
large existing data sets are known and a sparse representation is needed
and/or multiple, fast operations on them.

Then a minimal, nearly optimal set of representative such solutions
(waveforms, for the purpose of this review) is sought for in order to
construct a reduced basis (rb) for the whole solution set. These rb solutions
constitute an application-specific basis. That is, a set

$\\{h_{\Lambda_{i}}\\}_{i=1}^{n}$ (55)

of functions in the space of interest for carefully chosen parameter values
$\Lambda_{i}$ is used as a basis itself, as opposed to generic basis such as
polynomials or Fourier modes.

The scalar product in this parametrized space is taken at fixed parameter
values; that is, with respect to the physical dimension(s),

$\langle h_{i},h_{j}\rangle=\langle
h_{\lambda_{i}},h_{\lambda_{j}}\rangle=\langle
h({\cdot;\lambda_{i}}),h({\cdot;\lambda_{j}})\rangle\;.$

The same results of Section 6 follow through for the parametrized and
application-specific case. Namely, the optimal solution to the least-squares
approximation is

$h(\cdot;\lambda)\approx{\cal
P}_{n}h(\cdot;\lambda)=\sum_{i=1}^{n}c_{i}(\lambda)e_{i}(\cdot)\,,$ (56)

with the coefficients $c_{i}$ given by the equivalent of Eq. (46), in the
physical dimension(s) and the reduced basis 55 relabeled as
$\\{e_{1},\ldots,e_{n}\\}$. Namely,

$c_{i}(\lambda)=\sum_{j=1}^{n}({\bf G}^{-1})_{ij}\langle
e_{j}(\cdot),h(\cdot;\lambda)\rangle\,,$ (57)

where the entries of the Gramian matrix are

${\bf G}_{ij}:=\langle e_{i},e_{j}\rangle\,.$

Clearly, the Gramian coefficients depend on the special parameters
$\Lambda_{i}$ associated to each basis element.

So far, the span of the reduced basis uniquely determines the rb
representation. This reduced space depends primarily on the choice of the
selected parameter points in order to define a starting basis. As of a
convenient basis itself, for the same reasons discussed in Section 6, from a
conditioning numerical perspective, it is convenient to work with an
orthonormal basis built out of the rb solutions through a simple
orthonormalization procedure. This does not change the span of the rb and is
numerically convenient.

####  8.2 The Training set

A commonly used approach to construct a basis, described below in Section 8.3,
is through a greedy algorithm. In its simplest version, the algorithm
identifies a set of $n$ points in parameter space out of a representative
enough set of functions of interest that are actually known. We call this set
of functions the training set:

${{\cal K}}:=\\{h_{\lambda_{i}}\\}_{i=1}^{N}\,.$ (58)

The subset $\\{h_{\Lambda_{i}}\\}_{i=1}^{n}$ constitutes a nearly optimal
basis for application-specific spectral expansions of any function in ${\cal
K}$ in a precise sense discussed in Section 8.3.1. If there is partial
redundancy/similarity in the latter, then $n<N$ or even $n\ll N$. Note here
that we have used the symbol ${\cal K}$ to represent a discretization of the
space of interest ${\cal F}$ in order to perform actual computer calculations.

The training set can be constructed by any means, including simple random or
uniform sampling, more sophisticated stochastic methods [95], the metric
approach [102], or those of Ref. [92], for example. For large problems,
resampling it while constructing the basis can be critical; for an application
in the case of GWs see for example [28]. Regardless of the method used to
populate the training set, the RB-greedy formalism produces a compact and
highly accurate representation of the training space catalog. The reduced
basis is used to approximate other functions in the space of interest, whether
they were in the training set catalog or not, through linear combinations that
represent an orthogonal projection onto its span,

$h_{\lambda}={\cal P}_{n}h_{\lambda}+\delta h_{\lambda}\,,$ (59)

where the RB approximation is ${\cal P}_{n}h_{\lambda}$ and satisfies, by
construction, $\langle{\cal P}_{n}h_{\lambda},\delta h_{\lambda}\rangle=0$.

####  8.3 Greedy algorithms

We mentioned that a commonly used way to generate a reduced basis is through a
greedy algorithm. In its simplest form, such as when the waveforms are
inexpensive to compute or the data is already somehow known, the greedy
algorithm, outlined in Algorithm 1, has as input a discretization of the
parameter and solution space,

${\cal T}:=\\{\lambda_{i}\\}_{i=1}^{N}\quad{{\cal
K}}=\\{h_{\lambda_{i}}\\}_{i=1}^{N}$ (60)

with the elements of ${\cal T}$ usually called training points. To put
emphasis on structure instead of form, all functions in the training set in
this review are normalized. How to recover the norm, particularly for
parameter estimation of any detected signal, is discussed in Section 14.1.

The scheme needs an arbitrary seed $\Lambda_{1}\in{\cal T}$ to initialize it,
and a threshold error $\epsilon$ for a target representation accuracy (or
greedy error). Part of the output of the algorithm is a sequential selection
of $n$ parameter greedy points

$\\{\Lambda_{1},\Lambda_{2},\ldots,\Lambda_{n}\\}\subset{\cal T}$

and their associated waveforms

$\\{h_{\Lambda_{1}},h_{\Lambda_{2}},\ldots,h_{\Lambda_{n}}\\}\subset\\{h_{\lambda_{i}}\\}_{i=1}^{N}\,.$

The set of waveforms $\\{h_{\Lambda_{i}}\\}_{i=1}^{n}$ constitutes the reduced
basis. As was already discussed, for numerical conditioning purpose it is
sometimes convenient (see Eq. (10)) to work with an orthonormalized set
$\\{e_{i}\\}_{i=1}^{n}$ instead of directly the
$\\{h_{\Lambda_{i}}\\}_{i=1}^{n}$.

Another output of the algorithm is precisely the set of projection
coefficients for functions in the training space catalog, whereas coefficients
for any other known function can be computed through projection onto the
basis. In Section 10 we discuss how to approximate these coefficients through
interpolation in the physical dimension (time in the case of gravitational
waveforms), and in Section 11 a predictive approach for accurate and fast
evaluation of new (unknown) functions in ${\cal F}$ through surrogate models.
These are predictive models: as opposed to a known function being projected
into a compact basis, they predict new solutions (waveforms).

Algorithm 1 Greedy algorithm for reduced basis

1:Input: $\\{\lambda_{i}\,,h(\cdot;\lambda_{i})\\}_{i=1}^{N}$, $\epsilon$

2:Initialize $i=0$ and define $\sigma_{0}=1$

3:Seed choice (arbitrary): $\Lambda_{1}\in{\cal T}$,
$e_{1}=h(\cdot;\Lambda_{1})$

4:rb = $\\{e_{1}\\}$

5:while $\sigma_{i}\geq\epsilon$ do

6: $i=i+1$

7: $\sigma_{i}=\max_{\lambda\in{\cal T}}\|h(\cdot;\lambda)-{\cal
P}_{i}h(\cdot;\lambda)\|^{2}$

8: $\Lambda_{i+1}=\text{argmax}_{\lambda\in{\cal T}}\|h(\cdot;\lambda)-{\cal
P}_{i}h(\cdot;\lambda)\|^{2}$

9: $e_{i}=h(\cdot;\Lambda_{i+1})-{\cal P}_{i}h(\cdot;\Lambda_{i+1})$ (Gram-
Schmidt)

10: $e_{i+1}=e_{i+1}/\|e_{i+1}\|$ (normalization)

11: rb = rb $\cup\,e_{i+1}$

12:end while

13:Output: rb $\\{e_{i}\\}_{i=1}^{n}$ and greedy points
$\\{\Lambda_{i}\\}_{i=1}^{n}$

Comment 2.

  * •

Greedy-type algorithms are global optimization procedures used in contexts
outside reduced basis or even dimensional reduction.

  * •

Being a global optimization algorithm, the choice of the seed $\Lambda_{1}$ is
largely irrelevant. This was explicitly discussed as a sidenote in Ref. [36],
see Fig. 4 below. We also refer to Section 12.6 for a discussion of Ref. [36]
in the context of RB for multi-mode black hole ringdown.

  * •

What the greedy algorithm does in Step 8 is to select the waveform for which
its representation error onto the existing basis with $i-1$ elements is worst,
and in Step 11 adding it to the enrichment of the rb representation.

  * •

In steps 9 and 10 the rb waveforms are orthonormalized to avoid ill-
conditioning of the computation of the projection (see the discussion in
Section 6, before 6.1).

  * •

Given an arbitrary user-defined tolerance error $\epsilon$, the algorithm
stops when the approximation (56) meets the error tolerance – introduced as an
input in the greedy algorithm –,

$\|h_{\lambda}-{\cal
P}_{n}{h_{\lambda}}\|^{2}\leq\epsilon\,\,\,\,\forall\,\,\lambda\in{\cal T}.$

  * •

The expected exponential convergence of the method for the problems of
interest implies that $\epsilon$ can be made arbitrarily small with a
relatively small number $n$ of basis element, with $n<N$ and, in many cases,
$n\ll N$.

Figure 4: Taken from [36]. This figure shows the representation error as a
function of the number of reduced basis waveforms for a single QNM catalog.
The authors iterate over all possible seed waveforms in the training set. The
dark line represents the average, and the shaded area the maximum dispersion
around it. This numerical experiment confirms that the seed choice becomes
nearly irrelevant due to the global nature of the the greedy algorithm at each
step, as intuitively expected being a global optimization.

.

In practice, one uses a better conditioned orthogonalization algorithm in step
$9$ than the standard (i.e., “classical”) Gram-Schmidt one. The naive
implementation of the classical Gram-Schmidt procedure is actually ill-
conditioned. This is related to the fact that the Gramian matrix, which would
have to be inverted, can become nearly singular [136]. To overcome this one
can use an iterated Gram-Schmidt algorithm or a QR decomposition. A popular
alternative to the classical Gram-Schmidt, largely used, and called the
modified Gram-Schmidt, is also ill-conditioned, contrary to common perception
since its ill conditioning only becomes apparent for large large data sets.
See [117, 63, 131] for discussions about the conditioning and numerical
stability of different orthonormalization procedures.

Error definitions. We define the greedy error of the rb
$\\{{e_{i}}\\}_{i=1}^{n}$ as

$\displaystyle\sigma_{n}:=\max_{\lambda}\bigg{\|}h_{\lambda}-{\cal
P}_{n}h_{\lambda}\bigg{\|}^{2}\leq\epsilon\,,$ (61)

where $\epsilon$ is the user-defined tolerance error, which depends on the
number of basis elements.

The quantity $\sigma_{n}$ represents the largest error in the parameter space
of the best approximation by the greedy-reduced basis.

As discussed in [53], in the limit of sufficiently dense training spaces the
greedy error is comparable to the minimal match (${\mathrm{MM}}$) through

$\sigma_{n}\sim 1-{\mathrm{MM}}{\rm~{}~{}as~{}~{}}N\rightarrow\infty.$ (62)

Since the RB framework allows us to compress the information presented in the
training set, it is useful to introduce the quantity

$C_{r}:=N/n\,,$ (63)

called the compression ratio [119], to measure it.

#####  8.3.1 Convergence rates and near-optimality

The greedy algorithm chooses, in a precise sense, a nearly-optimal basis for
the function spaces ${\cal K}$ or ${\cal F}$, depending on whether one is
discussing the discrete or continuum cases. In order to quantify this near-
optimality, recall the definition of the $n$-width from Section 7, Eq. (53),
but write it in a more geometric way:

$\displaystyle d_{n}=d_{n}({\cal F},{\cal H})=\min_{W_{n}\in{\cal
H}}\max_{\lambda}\bigg{\|}h_{\lambda}-{\cal P}_{n}h_{\lambda}\bigg{\|}^{2}\,.$
(64)

The search space here is the whole Hilbert space and one looks for an optimal
$n$-dimensional subspace $W_{n}$ for approximating the function manifold
${\cal F}$. As was stated in Section 7, $d_{n}$ is a theoretical upper bound
to any practical algorithm to perform a linear approximation. In practice,
solving such optimization problem becomes unfeasible due to its intrinsic
combinatorial complexity. The nested nature of the greedy algorithm becomes
crucial for reducing the complexity of the Kolmogorov problem. This means that
each $H_{m}={\tt Span}\\{e_{i}\\}_{i=1}^{m}$ satisfies $H_{1}\subset
H_{2}\subset\ldots\subset H_{n}$ and this feature dramatically reduces the
search space in (64). One expects that this reduction of the search space is
at the expense of losing completely – if not a modicum – the $n$-width
optimality but, as we discuss in the next paragraph, this expectation is the
opposite from being true.

In which sense does the reduced basis-greedy procedure degrade the optimality
in the Kolmogorov sense? To answer this, lets summarize two important results
in relation with the greedy algorithm: if $d_{n}$ decays exponentially with
$n$ then so does $\sigma_{n}$ [48],

$d_{n}\leq
D\mathrm{e}^{-an^{\alpha}}\implies\sigma_{n}\leq\sqrt{2D}\gamma^{-1}\mathrm{e}^{-a^{\prime}_{\alpha}n^{\alpha}}\,,$
(65)

where $D$, $a$, $\alpha$ are positive constants.

Similarly, if the n-width has polynomial decay then so does the greedy error,

$d_{n}\leq Dn^{-\alpha}\implies\sigma_{n}\leq
D^{\prime}_{\alpha}n^{-\alpha}\,,$ (66)

where $D,\alpha>0$. More generally, for any decay rate of the n-width,

$\sigma_{n}\leq 2\gamma^{-1}d_{n/2}\,.$ (67)

The factor $\gamma$ in Eqs. (65) and (67) is a constant in $(0,1]$ [48]. In
recent years there have been efforts to improve these bounds. For details see
[148, 100].

In light of these results, we see that the reduced basis-greedy procedure
inherites the optimality of the $n$-width. If the latter has exponential
convergence, so it does the greedy error. In this precise sense the reduced
basis-greedy approach is nearly-optimal: we cut down the complexity of the
Kolmogorov problem at the very low expense of losing quality which, in most
practical applications, becomes insignificant.

The convergence rate of the $n$-width (and, in consequence, of the greedy
error) depends on the parametric smoothness/regularity of the
functions/waveforms. Indeed, if the functions are analytic and can be extended
to a complex region, exponential decay for the greedy error can be proven (see
[111], section 5.5, and citations therein). Therefore, the rb expansion is
expected to have very fast convergence to the original waveform as the number
of basis elements is increased. This is similar in spirit to the standard fast
convergence of spectral methods, but here smoothness in the parameters of the
problem is exploited, and the basis are elements themselves of the space of
functions of interest. For this reason RB is sometimes referred to as an
application-specific spectral expansion.

To summarize, for problems with smooth parametric dependence, fast convergence
in terms of greedy bases can be expected and it is not surprising that such
global methods outperform local ones.

Comment 3.  Not only the choice of basis is important, but a also global
expansion. To illustrate this consider again the example of the space of
smooth periodic functions from section 6.1.1. If one is interested in a family
of them, and is able to choose $n$ Fourier modes to represent them as best as
possible, one could:

  * •

build a discretization of the space of periodic functions and choose $n$
Fourier modes which when compared against the whole training set gives the
best global approximation error; this approach is a global variation of the
well known best m-term approximation problem, in which one is interested in
approximating some function and wants the best set of $m$ elements chosen from
a dictionary of functions to do it;

  * •

follow the standard approach, in which the functions of interest are projected
onto the span of the first $n$ Fourier modes. If the functions of interest are
smooth, this approach gives exponentially decaying representation errors with
the number of Fourier modes. This is, of course, not a predictive model, since
the function to be projected has to be known.

The first path would return a global and problem-dependent basis by
construction. The second one, on the contrary, will return a standarized basis
with no further feeling for the problem than the minimization of the maximum
error due to the projection.

#####  8.3.2 Complexity, scaling, computational aspects

Besides its near-optimality the greedy algorithm has several implementation
advantages:

  1. 1.

It is simple to implement (see Algorithm 1).

  2. 2.

It follows a stricter aproximation criteria (see section 3) than POD/SVD,
namely, the minimization of the worst error, as opposed to the average one, in
parameter space due to the projection onto the rb basis. It provides a
stricter error control.

  3. 3.

It is embarrassingly parallel. Each greedy sweep, defined by Step 7 of Alg. 1,
can be carried out simultaneously and independently for different values of
$\lambda\in{\cal T}$.

  4. 4.

As we remarked in Sec. 8.3.1, the greedy points and associated basis are
hierarchical (nested). This is particularly important if the algorithm is to
be used to guide which numerical relativity simulations should be carried out.
Since they are so expensive, it is desirable that if further accuracy is
needed, then it should be possible to achieve it by adding more points (as
opposed to starting the process from scratch).

  5. 5.

Low complexity. That is, the cost of all steps 6-12 in Alg. 1 are independent
of $i$, leading to a total cost linear in $n$, the total number of basis
elements. This observation follows from the fact that

${\cal P}_{i+1}h_{\lambda}={\cal P}_{i}h_{\lambda}+\langle
e_{i+1},h_{\lambda}\rangle e_{i+1}\,.$ (68)

So, if the projection coefficients are stored while building the basis, at
each sweep only the dot product between the elements of the training set and
the last basis element has to be computed.

#### Further reading

A comprehensive review of reduced basis is [112]. For textbooks see [70] and
[111].

## Part II Predictive Models

###  9 Polynomial interpolation

There are several reasons why we find it useful to discuss classical
polynomial interpolation:

  * •

It is used in ROM for gravitational wave physics and other fields; we discuss
some advantages and disadvantages.

  * •

There is a close relationship between interpolation, quadrature rules,
spectral methods and Reduced Basis. It is easier to introduce the standard
theory first, followed by the application specific version at the heart of ROM
for parametrized systems.

The emphasis of this section is on approximation quality quantified through
convergence rates and error bounds. Except for some side comments pointing out
implementations known to be poorly conditioned, we shall not discuss numerical
algorithms or their implementation which can be found in standard references.

####  9.1 Representation versus prediction

In standard approximation theory, discussed in Sections 5 and 6, one expands a
function in some basis. Assuming that the approximation quality is measured by
some weighted $L_{2}$ norm, the optimal approximation is the orthogonal
projection onto the representation space. This projection is independent of
the basis, in the sense that it only depends on its span. Assuming these are
orthonormal for simplicity and conditioning purposes, the orthogonal
projection is given by

$h(x)\approx{\cal P}_{n}h(x)=\sum_{i=1}^{n}\langle e_{i},h\rangle e_{i}(x)\,.$
(69)

Approximating $h$ by its projection as in (69), while perhaps useful by
itself, is not a predictive tool: we need to know the values of $h(x)$ for all
$x$ in order to calculate the projection coefficients in Eq. (69). Here, by
prediction it is meant a new evaluation within the range used to build the
approximating model. This means, for example, within the same parameter and
physical (time/frequency for example) ranges in the case of GWs. We are not
dealing with a much more difficult problem: namely, that one of prediction of
time series, for example, based on previous history. Then, to get a predictive
power we need to first turn to an interpolation scheme.

Building an interpolant for $h(x)$ requires only the knowledge of $n$ values
of $h(x)$ whereas the unknown values of the function are predicted from this
interpolant. As previously mentioned, the applications are many, including
numerical schemes to solve partial differential equations, the numerical
computation of integrals and building predictive surrogate models. As we will
see in Section 9.4 these apparently radically different approaches for
approximation, such as interpolation and projection, are closely related, and
in some cases both are equivalent within certain frameworks.

####  9.2 Prediction through polynomial interpolation

Suppose we have a collection of $(n+1)$ different points
$\\{x_{i}\\}_{i=0}^{n}$555Note that we count from $i=0$, only in this Section,
in order to match our notation with the common definition of Lagrange
polynomials., and a set of associated function values
$\\{f_{i}:=f(x_{i})\\}_{i=0}^{n}$. We assume that $f$ is a real valued, one-
dimensional function. We will discuss the more challenging multivariate case
in Section 9.5.

The interpolation problem consists of using the partial function sampling
$\\{f_{i}\\}_{i=0}^{n}$ to approximate $f(x)$ by another function which agrees
with $f$ at the interpolating nodes $\\{x_{i}\\}_{i=0}^{n}$. To be concrete,
let’s focus on linear interpolation: we require the interpolant of $f$ to be
expressible as a linear combination of the $(n+1)$ independent basis functions
$\\{e_{i}(x)\\}_{i=0}^{n}$,

${\cal I}_{n}[f](x):=\sum_{i=0}^{n}C_{i}e_{i}(x)\,.$ (70)

We already discussed in section 6 that the optimal least-squares-sense
approximation of the above form is provided by orthogonal projection onto the
basis $\\{e_{i}(x)\\}_{i=0}^{n}$. In the interpolation scheme one trades the
least-squares criteria for one which requires less information about $f$.
Let’s define this criteria by making the interpolant being equal to the
function at the nodes,

${\cal I}_{n}[f](x_{j})=f(x_{j})\,,\quad j=0,\ldots,n\,.$ (71)

When the basis are linearly independent polynomials of degree $\leq n$ we call
this special case Polynomial interpolation. One can show by explicit
construction that there is a unique polynomial which satisfies the
interpolation problem (70,71). Before doing so we give an example.

Example 5.

Suppose we have two points, $\\{x_{0},x_{1}\\}$, sampling a function $f(x)$.
Let the sample values be $f_{0}=f(x_{0})$ and $f_{1}=f(x_{1})$. The functions
$\\{1,x\\}$ comprise a suitable basis for our $n=1$ polynomial interpolant.
Then the interpolating polynomial is of degree $\leq 1$, i.e. a straight line
joining the two points:

${\cal I}_{1}[f](x)=ax+b$

We explicitly solve for $(a,b)$ by writing the interpolation problem (71)

$\displaystyle{\cal I}_{1}[f](x_{0})$ $\displaystyle=$ $\displaystyle
ax_{0}+b=f_{0}$ (72) $\displaystyle{\cal I}_{1}[f](x_{1})$ $\displaystyle=$
$\displaystyle ax_{1}+b=f_{1}$ (73)

which has as solutions

$a=\frac{f_{0}-f_{1}}{x_{0}-x_{1}}\,\,,\,\,b=\frac{x_{0}f_{1}-x_{1}f_{0}}{x_{0}-x_{1}}\,.$
(74)

Thus,

${\cal
I}_{1}[f](x)=\left(\frac{f_{0}-f_{1}}{x_{0}-x_{1}}\right)x+\left(\frac{x_{0}f_{1}-x_{1}f_{0}}{x_{0}-x_{1}}\right)\,.$
(75)

General solution: Lagrange polynomials

Existence. The existence of a polynomial which solves the interpolation
problem (71) can be shown by explicit construction. Suppose we have access to
a basis $e_{i}(x)=\ell_{i}^{(n)}(x)$ which satisfies

${\cal I}_{n}[f](x)=\sum_{j=0}^{n}f_{j}l_{j}^{(n)}(x)$ (76)

where, for each $j=0\ldots n$, $l_{j}^{(n)}(x)$ is a polynomial of degree
$\leq n$ such that

$l_{j}^{(n)}(x_{i})=\delta_{ij}\,\,\,\mbox{for }i=0\ldots n$ (77)

and $\delta_{ij}$ is the Kronecker delta ($\delta_{ij}=1$ if $i=j$ and zero
otherwise). Then ${\cal I}_{n}[f](x)$ as given by Eq. (76) would be a solution
to the interpolation problem.

The Lagrange polynomials are precisely those with the property (77). They are
given by

$l_{j}^{(n)}(x)=\left(\prod_{k=0,k\neq
j}^{n}(x-x_{k})\right)/\left(\prod_{k=0,k\neq j}^{n}(x_{j}-x_{k})\right)\,.$
(78)

Example 6.

We look at the same case ($n=1$) of Example 9.2, but now we construct an
interpolant through Lagrange polynomials. By Eq. (78) these are

$\l_{0}^{(1)}(x)=\frac{x-x_{1}}{x_{0}-x_{1}}\;\;\;,\;\;\;\l_{1}^{(1)}(x)=\frac{x-x_{0}}{x_{1}-x_{0}}\,,$
(79)

and using (76) the Lagrange form of the interpolant becomes

$\displaystyle{\cal I}_{1}[f](x)$ $\displaystyle=$ $\displaystyle
f_{0}l_{0}^{(1)}(x)+f_{1}l_{1}^{(1)}(x)$ (80) $\displaystyle=$
$\displaystyle\left(\frac{f_{0}-f_{1}}{x_{0}-x_{1}}\right)x+\left(\frac{x_{0}f_{1}-x_{1}f_{0}}{x_{0}-x_{1}}\right)$

The result (80) is the same as (75) due to the uniqueness of the polynomial
interpolant, which we now prove.

Uniqueness. Suppose there are two polynomials $P_{n}$ and $R_{n}$ of degree
$\leq n$ satisfying the interpolation condition
$P_{n}(x_{j})=R_{n}(x_{j})=f_{j},\,j=0,\ldots,n$. The polynomial
$Q(x):=P_{n}(x)-R_{n}(x)$ then satisfies $Q(x_{j})=0,\,j=0,\ldots,n$. This
means that $Q(x)$ has $(n+1)$ zeros, in contradiction with the hypothesis that
it is at most of degree $n$. Ergo, $Q(x)=0$ and $P_{n}(x)=R_{n}(x)$.

The Vandermonde form

Yet another way of showing existence and uniqueness for the polynomial
interpolant in a constructive way is the following. We write down the
polynomial interpolant using monomial basis functions,

${\cal I}_{n}[f](x)=\sum_{j=0}^{n}a_{j}x^{j}\,.$

The interpolation conditions (71) become a set of linear equations for the
coefficients $a_{j}$

${\bf V}a=f\,,$ (81)

where ${\bf V}_{ij}=x_{i}^{j-1}$, $a=[a_{0},\ldots,a_{n}]^{\top}$ and
$f=[f(x_{0}),\ldots,f(x_{n})]^{\top}$. The matrix ${\bf V}$ is called the
Vandermonde matrix and its determinant, also called the Vandermonde
polynomial, takes the form

$\det({\bf V})=\prod_{0\leq i<j\leq n}(x_{j}-x_{i})\,.$

In order to have one and only one solution to Eq. (81) this determinant has to
be non-zero. In one spatial dimension this condition is automatically
satisfied whenever the interpolation points are unique. There are ways of
inverting ${\bf V}$ in order to solve Eq. (81) using Lagrange polynomials,
which we will not pursue here. The main reason for briefly mentioning the
Vandermonde matrix is because a generalization of it will play a prominent
role in the Empirical Interpolation Method, described in Section 10.

####  9.3 Convergence rates, collocation points and Runge’s phenomenon

Here we discuss the behavior of the error of the polynomial interpolant ${\cal
I}_{n}[f](x)$ of some function $f(x)$ at $(n+1)$ nodes, to illustrate the
difficulties in obtaining an accurate interpolation scheme on a unstructured
set of nodes which, from the point of view of ROM, one wishes to be sparse.

Leaving aside both sparsity and unstructured meshes, one might imagine that
for “nice enough” functions (for example, infinite differentiable), this error
decreases as the number of nodes and, therefore, degree of the polynomial
basis $n$ increase. This is not the case, the classical counterexample is
Runge’s one, here described, and the overall problem is referred to as Runge’s
phenomenon.

The distinction is between local approximations of fixed degree and global
ones. Because these two concepts play a key role in reduced order modeling we
discuss them in this context with some detail. Suppose that one has
$\\{x_{i}\\}_{i=0}^{n}$ nodes at which to interpolate a function.

Local interpolation

We consider local interpolants of degree one because we have already
explicitly discussed them in Example 75. One would partition the interval of
interest, where on each subinterval

$I_{j}:=[x_{j},x_{j+1}]\,,\quad j=0,\ldots,n-1$ (82)

one builds a local interpolant of degree one using the end points as
interpolating nodes,

${\cal
I}[f](x)=\left(\frac{f_{j}-f_{j+1}}{x_{j}-x_{j+1}}\right)x+\left(\frac{x_{j}f_{j+1}-x_{j+1}f_{j}}{x_{j}-x_{j+1}}\right)\,,\quad\text{for
}x\in I_{j}\,,\quad j=0,\ldots,n-1\,.$ (83)

Such fixed order interpolants, where their degrees are independent of the
number of nodes $n$, is – as explained later below – guaranteed to converge to
the interpolated function as the number of nodes in a given fixed physical
interval increases (the gridspacing decreases). Furthermore, such convergence
is guaranteed regardless of the structure/location of the interpolating nodes,
they can even be randomly located. Even more, to some degree, convergence is
also guaranteed regardless of the smoothness of the data; this is why fixed
order methods, especially very low order ones, are in general used for noisy
data.

This class of methods is usually very robust, but at the expense of accuracy.
The latter is exacerbated in the absence of a dense amount of data. On the
contrary, in Reduced Order Modeling one seeks very sparse representations of
very high accuracy. This combination of, and even trade-off between,
robustness (convergence) and accuracy is a delicate and challenging issue in
ROM.

Global interpolation

Instead, one could use a global interpolant, where all $(n+1)$ nodes and
function values are used to build a single interpolant,

${\cal I}_{n}[f](x)=\sum_{i=0}^{n}\ell_{i}^{(n)}(x)f_{i}\,.\quad\forall
x\in[x_{0},x_{n}]\,.$

As the density of points and $n$ increases, one would keep using all of them,
thereby increasing the degree of the interpolant. This might seem desirable
and more accurate than a local approach of fixed degree, since more
information is used. The classical counterexample is due to Runge [118], and
it consists of interpolating the function

$f(x)=\frac{1}{1+(5x)^{2}}\;\,,\quad x\in[-1,1]\,.$ (84)

This function is $C^{\infty}$ (has infinite derivatives) yet the error of a
global polynomial interpolant at equally spaced points diverges near the
boundaries [52] at higher orders (see Fig. 5). For unstructured points, such
as those usually appearing in reduced models, the lack of convergence of
global approximations is in general worse.

Figure 5: Runge’s phenomenon in polynomial interpolation. In this example the
approximation uses 19 equispaced nodes over the interval $[-1,1]$ and the
interpolant polynomial is of degree 18. The interpolation tends to diverge
near the boundaries. This issue cannot be solved by increasing the resolution
of the sampling but the opposite: divergencies get worse for equispaced nodes.

As discussed below, this can be resolved and fast convergence achieved if the
interpolation nodes can be appropriately chosen. However, standard choices for
such special nodes are in general not suitable for reduced order modeling. In
particular, they are neither necessarily sparse nor adapted to the problem of
interest. They are generic and in general (and most important) not
hierarchical, requiring recomputing from scratch the sought approximation if
higher accuracy is required.

Error analysis and convergence I

In what follows we assume that $f$ is smooth enough as needed from the
context. Denoting the interpolation error by

$E_{n}(x):=\left|f(x)-{\cal I}_{n}[f](x)\right|$

the following result holds:

$E_{n}(x)=\left|\frac{1}{(n+1)!}f^{(n+1)}(\xi)\omega_{n+1}(x)\right|$ (85)

where

$\omega_{n+1}(x):=\prod_{i=0}^{n}(x-x_{i})$

is called the nodal polynomial of order $(n+1)$, and $\xi$ is in the smallest
interval $I_{x}$ containing $x_{0}\ldots x_{n}$ and $x$. Equation (85) can be
proven as follows. Let’s define the function

$F(t)=f(t)-{\cal I}_{n}[f](t)-\alpha(x)\omega_{n+1}(t)$

where

$\alpha(x)=\frac{f(x)-{\cal I}_{n}[f](x)}{\omega_{n+1}(x)}$

This function is well defined since $x\neq x_{i}$, and it is as smooth as
$f(t)$. Observe that $F(x)=0$ and, in particular, $F(x_{i})=0$ ($i=1,...,n$).
This means that $F(t)$ has $(n+2)$ distinct zeros in the interval in which
$f(t)$ is smooth. Now, the Mean Value Theorem tells us that, between any
couple of zeros, there is one zero for the derivative $F^{\prime}(t)$. Then
the function $F^{\prime}(t)$ has $(n+1)$ zeros. Repeating the previous
argument, we conclude that $F^{(n+1)}(t)$ has exactly one zero. Let’s call it
$\xi$. Thus $F^{(n+1)}(\xi)=f^{n+1}(\xi)-\alpha(\xi)(n+1)!=0$. The statement
in (85) follows directly from this relation.

Let’s make a comment here. The point $x$ at which the interpolant is evaluated
when computing the error (85) does not need to be in the smallest interval
containing $x_{0}\ldots x_{n}$. Assuming an ordering $x_{0}<\ldots<x_{n}$,
sometimes the process of approximating $f(x)$ by ${\cal I}_{n}[f](x)$ is
called interpolation only if $x\in[x_{0},x_{n}]$ and extrapolation otherwise.

We get back to the error formula (85) and analyze its behavior as $n$
increases. In particular, we ask ourselves under what conditions does

$E_{n}(x)\rightarrow 0\;\;\;\mbox{as}\;\;\;n\rightarrow\infty\,.$

The error (85) can be decomposed in two terms, one related to the behavior of
the derivatives of $f$,

$\left|\frac{f^{(n+1)}(\xi)}{(n+1)!}\right|$

and another one related to the distribution of nodes through
$\omega_{n+1}(x)$.

We first discuss why fixed order interpolants are guaranteed to converge,
regardless of the distribution of nodes (assuming the interpolated function is
smooth enough), since the argument is simple.

Discussion 1. Local interpolants

For definiteness consider the case of equally spaced points (the
generalization to an arbitrary set is straightforward),

$x_{j}=x_{0}+j\Delta x\,,\quad j=0,\ldots n\,.$

Then the error of the interpolant (83) on each subinterval (82) is bounded by

$E\leq\frac{1}{(m+1)!}\max_{x\in I_{j}}\left|f^{(m+1)}(x)\right|(\Delta
x)^{m}\,,$ (86)

where $m$ is the (fixed) number of interpolation nodes and degree of the
interpolants (clearly, $m\leq n$). Since $m$ is fixed, the convergence rate is
determined by $\Delta x\rightarrow 0$, and is of order $m$.

Discussion 2. Global interpolants

In this case $m=n$. One can see from the interpolation error (85) that
convergence is guaranteed if the derivatives of $f$ decay with $n$
sufficiently fast. In fact, the interpolant converges in the infinity norm if
all derivatives $f^{(s)}$ are bounded by the same constant in the interval of
interest. However, it is not obvious from (85) to find out in a
straightforward way the allowed growth rate of the derivatives of $f$ such
that convergence/divergence of the error takes place.

In fact, it is not easy to see how such error formula explains Runge’s
phenomenon without resorting to complex calculus [52]. However, we can ask
ourselves the next related questions:

(Q1) What can be said about the maximum of $|\omega_{n+1}(x)|$ in the interval
of interest?

(Q2) If the physical problem allows such freedom, how can one choose the nodes
so as to minimize such maximum?

For simplicity and without loss of generality we focus on the interval
$x\in[-1,1]$. We can answer these questions as:

Answer to (Q1): For all choices of nodes

$\max_{[-1,1]}|\omega_{n+1}(x)|\geq 2^{-n}$ (87)

Answer to (Q2): If

$\omega_{n+1}(x)=2^{-n}T_{n+1}(x)\,,$

where $T_{j}$ denotes the Chebyshev polynomial of degree $n$, then

$\max_{[-1,1]}|\omega_{n+1}(x)|=2^{-n}\,.$

In other words, using the roots of the Chebyshev polynomial as interpolating
nodes minimizes the error associated with the location of such nodes. It is in
this sense that Chebyshev nodes are many times referred to as optimal points
for interpolation. They have the disadvantage, though, that they are note
hierarchical/nested, or adapted to the problem of interest. The former failure
means, from a ROM perspective, that any model relying on them, including
perhaps expensive numerical simulations, have to be completely built from
scratch if higher accuracy is required. This is solved by a nearly-optimal set
of nested nodes for application-specific interpolants designed with ROM
problems in mind, discussed in Section 10.

Error analysis and convergence II

In principle, direct approximation through projection requires the knowledge
of the continuum waveform, while interpolation only needs information at a
(specified by the user) set of nodes. We next address what is the relationship
between these two strategies. In terms of representation accuracy, we focus on
the $L_{2}$ norm, since it is directly related to the overlap error commonly
used in data analysis. Since, as discussed, given a basis and linear
approximation, projection is optimal in the least squares sense, one knows
that

$\|h-{{\cal P}}_{n}h\|^{2}\leq\|h-{\cal I}_{n}[h]\|^{2}\,.$ (88)

The question is how “sub-optimal” in the above sense interpolation is. The
answer obviously depends on the basis and interpolation nodes. A general
framework to study this is through the definition of the Lebesgue constant.

For any arbitrary basis and interpolation scheme one can derive the more
precise bound related to (88):

$\|h-{\cal I}_{n}[h]\|^{2}\leq\Lambda_{n}\|h-{{\cal P}}_{n}h\|^{2}\,,$ (89)

where

$\Lambda_{n}:=\|{\cal I}_{n}\|^{2}\,.$ (90)

Comment 4.

  1. 1.

The $\Lambda_{n}$’s are referred as the Lebesgue constants and depend, given a
vector norm, purely on the interpolation operator ${\cal I}_{n}$.

  2. 2.

In terms of accuracy, the strategy is to build an interpolant (by choosing
both the basis and nodes) such that the Lebesgue constant grows as slow as
possible with $n$. It is a computable quantity, so one can actually quantify
how much is “lost” with respect to projection through Eq. (89).

  3. 3.

Strictly speaking, the Lebesgue constant is usually referred to in the context
of the infinity norm,

$\|u\|_{\infty}:=\max_{i}|u_{i}|\,,$

but throughout this review, unless otherwise stated, we refer to it in the
$L_{2}$ norm.

As an example, and going back to Runge’s phenomenon at equally spaced versus
Chebyshev nodes for global polynomial interpolation, the following illustrates
the advantages of the latter and explains why it is such a popular choice –
when it is a feasible approach at all within the problem of interest, which is
not always the case in practice.

Example 7.

Here the $\|.\|_{\infty}$ norm is used. For large $n$ one has the following
behavior (see, for example, [71]).

  1. 1.

For equally spaced nodes, the Lebesgue constant grows exponentially, as

$\Lambda_{n}\sim\frac{2^{n+1}}{n\log n}\,.$ (91)

  2. 2.

For Chebyshev nodes, it only grows logarithmically,

$\Lambda_{n}\sim\log{(n)}\,.$ (92)

In the context of polynomial interpolation, it can be proved that the Lebesgue
constant grows at least logarithmically in the infinity norm for any selection
of nodes. In particular, Chebyshev nodes belong to a family that induce a
nearly-optimal behavior on the Lebesgue constant, since the growth (92) is
only logarithmic. For thorough discussions around this topic see [71].

####  9.4 Discrete expansions and interpolation

In approximation through projection one in principle needs full knowledge of
the function $f$ to be approximated through the projection coefficients
$\langle e_{i},h\rangle$ in Eq. (47). In interpolation, instead, one only
needs knowledge of the function at the interpolating nodes.

In practice, the integrals to compute $\langle e_{i},h\rangle$ are replaced by
a quadrature rule, i.e. a numerical approximation (we discuss quadrature rules
and specifically those built using ROM in Section 13). These are called then
discrete projection approximations. Here it suffices to say that, if the
quadrature nodes can be chosen, one can find an optimal family of quadrature
rules that maximizes the degree for which the quadrature is exact for
polynomials. These are called Gaussian quadratures and the nodes Gaussian
nodes; examples of the latter are Legendre and Chebyshev nodes. The relevant
result here relating approximation by discrete projection and interpolation is
the following:

Result 1. Discrete projection and interpolation

If the projection coefficients are approximated by Gaussian quadratures and
are used Gaussian nodes for interpolation, then the resulting approximation
${\cal P}_{n}h$ exactly equals the polynomial interpolant ${\cal I}[h]$ [71].

Through this result one can guarantee fast convergence of interpolation using
Gauss nodes. The problem is that in many parametrized problems of interest,
one cannot choose Gaussian nodes in parameter space as representative
solutions. Either because data is not available (for example, if taken from
experiments), or because it is not efficient to do so from a modeling
perspective. In Section 10 we discuss a generalization for ROM and
parametrized systems that attempts to generalize Result 9.4.

####  9.5 Multiple dimensions

Evaluation cost of multivariate polynomials

Even when counting with enough data for multi-dimensional fits, enough for
accurate local approximations which avoid Runge’s phenomenom, a perhaps not so
well known fact is that the evaluation cost of multivariate polynomials in
general grows exponentially with the dimensionality of the problem.

Naive evaluation of a $1$-dimensional polynomial of degree $n$, that is,
evaluating all the monomials in the standard format

$p_{n}(x)=a_{0}+a_{1}x+\ldots+a_{n}x^{n}$ (93)

carries an operation count of ${\cal O}(n^{2})$. Through a simple
factorization, Horner’s rule (see [132]), the operation count can be reduced
to $2n$. Suppose we want to evaluate $p_{n}(x)$. First note that $p_{n}(x)$
can be decomposed in a nested way

$p_{n}(x)=a_{0}+x(a_{1}+x(a_{2}+\cdots x(a_{n-1}+a_{n}x)\cdots)\,.$

Horner’s rule consists in computing iteratively the recursion

$\begin{split}b_{0}=&a_{n}\\\ b_{k}=&a_{n-k}+b_{k-1}x\quad
k=1,\cdots,n\end{split}$

reaching $p_{n}(x)=b_{n}$ in $2n$ operations (multiplication and addition, $n$
times). It can be shown that this is the minimum number of operations to
evaluate a one-dimensional problem unless some offline factorization is
carried out for multiple online evaluations, in which case the computational
cost can be slightly reduced.

Next, consider multivariate polynomial evaluation. Suppose one attempts to
evaluate the generalization of (93) to the $d$-dimensional case. For
definiteness, the form in the ${\tt dim}=2$ case would take the form

$p_{n}(x,y)=\sum_{i=0}^{n_{x}}\sum_{j=0}^{n_{y}}a_{ij}x^{i}y^{j}\,.$

The evaluation cost of such an approach in the $d$-dimensional case, assuming
for simplicity the same degree $n$ in each dimension, is of order ${\cal
O}(n^{2d})$. Ideally, one would wish a generalization of Horner’s rule with,
say, evaluation cost of order ${\cal O}(d\times n)$. Unfortunately not such
algorithm is known, neither a proof of which could be the optimal number of
operations needed to perform a multivariate polynomial evaluation. In recent
years there have been several proposals to face this problem, ranging from
rigorous mathematical attempts (see e.g. [106, 98, 141, 81, 18, 18, 77, 142])
to the proposal of heuristic greedy algorithms [37] aiming to reduce the
number of calculations.

Multivariate polynomial evaluations scale exponentially with the
dimensionality of the problem. This is a severe issue in most cases of
interest which cannot be underestimated. In the context of ROM, the
limitations of multivariate polynomial approximations appear in almost every
possible context. This point is discussed in more detail and emphasized in
Section 15.

#### Further reading

Reference [58] reviews multivariate polynomial interpolation, and [99]
presents a proposal for polynomial interpolation on unstructured grids in
multiple dimensions.

There are schemes for fast multiple evaluations of multi-variate polynomials
though. They are based on an offline factorization of the given polynomial
(which can be expensive, but done only once), with a total cost of ${{\cal
O}}(n^{{\tt dim}+1})$ for ${\cal O}(n^{\tt dim})$ evaluations, leading to an
average of ${\cal O}(n)$ for each evaluation [84]. To our knowledge this
problem of multi-dimensionality has not been yet tackled in ROM for
gravitational waves and is still largely an open problem in approximation
theory. Instead, in GW science local interpolants in the form of splines are
being used, which results in degradation of the accuracy of the surrogate
models being built. In terms of higher (than one) dimensions, sparse grids are
being used. One of their problems is that they are not necessarily
hierarchical, leading to a problem when building training sets from numerical
relativity to later apply the RB framework.

###  10 The Empirical Interpolation Method

The Empirical Interpolation Method (EIM) was proposed in 2004 [20] as a way of
identifying a good set of interpolation nodes on multi-dimensional
unstructured meshes and has since found numerous applications (for a very
short list, see [1, 39, 40, 51, 90]). When the basis is application-specific,
so is the EIM.

Here, perhaps contrary to what might be suspected, interpolation is in the
physical dimension(s), not on the parametric ones. For definiteness, in the
case of gravitational waves this is either the time or frequency domain.

Before getting into technical details, we highlight some of the properties of
the EIM:

  * •

It is hierarchical (nested). The “most relevant” points in the physical
dimension(s) and associated interpolants are selected and built, respectively.
This follows a process that is dual, in a precise sense, to the construction
of a reduced basis using a greedy algorithm.

  * •

It is designed to overcome the difficulties of stability, accuracy and
evaluation cost on sparse and scattered grids, especially in multiple
dimensions.

  * •

It is highly accurate. By design the algorithm attempts to control the
behavior of the Lebesgue constant (see Section 10.3) at each greedy EIM sweep.

  * •

It provides an affine parametrization which has multiple consequences. In
particular, it allows for the design of Reduced Order Quadratures, described
in Section 13, which lead in particular to fast likelihood evaluations, as
discussed in Section 14.

The material below relies on standard interpolation ideas which we have
summarized in Section 9. We will refer to the physical variable(s) as $x$, the
nodes selected by the EIM as $X_{i}\,(i=1\ldots n)$ 666Except for the GW case,
in which we will explicitly use $t$ and $T$ – or $f$ and $F$, in the case of
frequency., the (in general multidimensional) parameter as $\lambda$ and its
corresponding greedy points as $\Lambda_{i}\,(i=1\ldots n)$. It is not a
coincidence that there is the same number of EIM nodes and greedy points, in
fact they go hand by hand.

#### EIM: algorithm and properties

Within the EIM method one seeks to find an empirical (that is, problem-
dependent) global interpolant, which also provides an affine parametrization
of functions. That is, the EIM approximation of a function $h(x;\lambda)$ is
of the form

$h(x;\lambda)\approx{\cal
I}_{n}[h](x;\lambda)=\sum_{i=1}^{n}B_{i}(x)h(X_{i};\lambda)\,.$ (94)

The affine parametrization of ${\cal I}_{n}[h](x;\lambda)$ refers to the fact
that the r.h.s. of (94) is a sum of terms which depend only on $x$ (the
$B_{i}$ coefficients) multiplied by other ones which depend only on $\lambda$
(the $h(X_{i};\lambda)$); one might also want to refer to this as
“separability”.

Below in Section 10.1 we discuss how one arrives at (94), how to compute in
the offline stage the $B_{i}$ coefficients, and the quality of this
approximation.

In addition, even though it is perhaps not usually viewed this way, or not
emphasized enough compared to online evaluations to solutions of parametrized
PDEs, the EIM also provides an application-specific down- and upsampling that,
for practical purposes, beats Nyquist downsampling, with implications for
signal or more generally data processing. See for example Section 11.2 and
discussion around Figures 12 and 13.

####  10.1 From projection to interpolation

In Section 3 we discussed that the optimal linear representation in any
weighted $L_{2}$ norm, given a basis (assuming for simplicity and conditioning
that it is orthonormal) of cardinality $n$ is through an expression of the
form

$h(x;\lambda)\approx{\cal
P}_{n}h(x;\lambda):=\sum_{i=1}^{n}c_{i}(\lambda)e_{i}(x)\,,$ (95)

where

$c_{i}(\lambda)=\langle e_{i}(\cdot),h(\cdot;\lambda)\rangle\,.$ (96)

This is the standard approximation by projection, the most common one being
using polynomials or Fourier modes as basis. Now, approximation through
projection requires knowledge of the function to be represented at sufficient
values of $x$ so as to accurately compute the coefficients (96). In Section 9
we discussed how this is usually replaced by interpolation, and the strong
relationship between projection and interpolation in standard spectral theory.
Next we describe the EIM strategy, which mimics the spectral approach.

The approximation through projection onto a reduced basis, Eq. (95), is
replaced by interpolation (in the physical dimension(s), $x$) as follows.
First a rb $\\{e_{i}(x)\\}$ is chosen, for example through a POD or greedy
approach. The interpolant then is sought to have the form

${\cal I}_{n}[h](x;\lambda):=\sum_{i=1}^{n}C_{i}(\lambda)e_{i}(x)\,,$ (97)

where on purpose the expansion, which are not projection, coefficients $C_{i}$
are denoted by capital letters to distinguish them from the projection ones
(96). Instead, they are defined to be solutions of the interpolation problem;
namely that the interpolant exactly agrees with the function at the EIM nodes
(the construction of which we discuss below)

$\displaystyle{\cal
I}_{n}[h](X_{i};\lambda)=h(X_{i};\lambda),\qquad\forall\,i=1,\dots,n.$ (98)

For the moment, we shall assume that the EIM nodes $X_{i}$ are known and
proceed to describe how to use them to find the EIM interpolant. This can be
somewhat misleading, since it is not the way it works in practice, which is:
the first EIM node is found, its associated interpolant built, the second EIM
node found, the interpolant enriched to take it into account, and so on. They
are not disjoint processes, unlike (again) standard spectral methods, where
all the Gaussian nodes are found and afterwards the interpolant built. This is
not only a procedural difference, but highlights a big difference of the EIM:
namely that the approach (nodes and interpolant) is hierarchical. Hopefully
this gradual presentation is intuitive and pedagogical, later in this section
we will present the full, coupled algorithm together with its numerical
intricacies.

Equation (98) is equivalent to solving the $n$-by-$n$ system

$\sum_{i=1}^{n}{\bf V}_{ji}C_{i}(\lambda)=h(X_{j};\lambda)\,,\quad
j=1,\ldots,n$ (99)

for the coefficients $\\{C_{i}\\}_{i=1}^{n}$, where the interpolation matrix

$\bf{V}:=\left(\begin{array}[]{cccc}e_{1}(X_{1})&e_{2}(X_{1})&\cdots&e_{n}(X_{1})\\\
e_{1}(X_{2})&e_{2}(X_{2})&\cdots&e_{n}(X_{2})\\\
e_{1}(X_{3})&e_{2}(X_{3})&\cdots&e_{n}(X_{3})\\\
\vdots&\vdots&\ddots&\vdots\\\
e_{1}(X_{n})&e_{2}(X_{n})&\cdots&e_{n}(X_{n})\\\ \end{array}\right)$ (100)

is a generalization of the Vandermonde matrix – see Eq. (81) – when polynomial
basis are used. We have already discussed that in that case the Vandermonde
matrix can easily be very ill-conditioned if the nodes are chosen, for
example, equally spaced, not to mention if they are scattered. So one can
anticipate that the solution to the problem (97) is already non-trivial and
one of the goals of the EIM is to make sure that it does not lead to an ill
conditioned problem, as well as providing a high accuracy interpolant.

The choice of empirical nodes given by the EIM together with the linear
independence of the reduced basis ensure that $\bf{V}$ (Eq. (100)) is
invertible so that

$C_{i}=\sum_{j=1}^{n}\left({\bf V}^{-1}\right)_{ij}h(X_{j};\lambda)$ (101)

is the unique solution to (98). It then follows upon substituting (101) into
(97) that the empirical interpolant is

${\cal I}_{n}[h](x;\lambda)=\sum_{j=1}^{n}B_{j}(x)h(X_{j};\lambda)$ (102)

where

$B_{j}(x):=\sum_{i=1}^{n}e_{i}(x)\left({\bf V}^{-1}\right)_{ij}$ (103)

is independent of $\lambda$.

Note that (102) is a linear combination of the waveform itself evaluated at
the empirical nodes. The coefficients $\\{B_{i}\\}_{i=1}^{n}$ satisfy
$B_{i}(X_{j})=\delta_{ij}$ and are built directly from the reduced basis. They
provide a clean offline/online separation. Because of this the
$\\{B_{i}\\}_{i=1}^{n}$ functions can be pre-computed offline once the reduced
basis is generated while the (fast) interpolation is evaluated during the
online stage from (102) when the parameter $\lambda$ is specified by the user.
Evaluations of the waveform at the EIM nodes are still needed at the
arbitrarily chosen parameter $\lambda$ in order to construct the interpolant
in (102). One can wonder how this can be of any use. In Section 11 we explain
how to build surrogate predictive models for the evaluations of
$h(X_{j};\lambda)$ for any $\lambda$ (that is, not present in the original
training space).

Next we discuss how to compute the EIM nodes in a rather qualitative way,
before presenting the general algorithm later on.

Example 8.  Finding the EIM nodes and building the interpolant.

Consider for definiteness a time series. The algorithm takes as input the
basis set $\\{e_{i}\\}_{i=1}^{n}$ and an arbitrary number and choice of time
samples $\\{t_{i}\\}_{i=1}^{L}$ from which the empirical interpolation nodes
$\\{T_{i}\\}_{i=1}^{n}$ are to be selected. The EIM algorithm proceeds as
follows

  1. 1.

The first time node is chosen to maximize the value of $|e_{1}(t_{i})|$; that
is, $\left|e_{1}(T_{1})\right|\geq\left|e_{1}(t_{i})\right|$ for all time
samples $t_{i}$.

  2. 2.

Next, an empirical interpolant for the second basis function is built using
only the first basis function: From Eqs. (97,98) we have ${\cal
I}_{1}[e_{2}](t)=C_{1}e_{1}(t)$ where $C_{1}=e_{2}(T_{1})/e_{1}(T_{1})$ has
been found from Eq. (98) with $i=1$.

  3. 3.

The second empirical interpolation node is chosen to maximize the value of the
pointwise interpolation error of ${\cal I}_{1}[e_{2}](t)-e_{2}(t)$; that is,
$\left|{\cal I}_{1}[e_{2}](T_{2})-e_{2}(T_{2})\right|\geq\left|{\cal
I}_{1}[e_{2}](t_{i})-e_{2}(t_{i})\right|$ for all data samples.

  4. 4.

Steps $2$ and $3$ are then repeated to select the remaining $(n-2)$ nodes. The
full algorithm for generating the EIM nodes is shown in Algorithm 2.

Comment 5.  As described, the EIM follows a greedy approach, albeit somewhat
different from that one usually used to build a reduced basis. While a greedy
algorithm to build a RB selects the most relevant points in parameter space,
the EIM selects the most relevant points in the physical dimension(s). In
addition, the former uses a (possibly weighted) $L_{2}$ norm while the EIM
uses the infinity one.

Comment 6.  Notice also that even though constructing the EIM is also done
offline, it has a much smaller computational cost than constructing the basis,
since it operates only on the latter, the training set is not involved
anymore.

The following example graphically shows the first three iterations of the EIM
algorithm applied to Legendre polynomials.

Example 9.

Consider the set $\\{P_{i}(x)\\}_{i=0}^{n}$ of $(n+1)$ normalized Legendre
polynomials defined on $[-1,1]$. These form an orthonormal basis for the space