# Randomized Reduced Basis Methods for Parameterized Fractional Elliptic PDEs

Harbir Antil Department of Mathematical Sciences and Center for Mathematics
and Artificial Intelligence, George Mason University Email<EMAIL_ADDRESS>Arvind K. Saibaba Department of Mathematics, North Carolina State University
Email<EMAIL_ADDRESS>
###### Abstract

This paper is interested in developing reduced order models (ROMs) for
repeated simulation of fractional elliptic partial differential equations
(PDEs) for multiple values of the parameters (e.g., diffusion coefficients or
fractional exponent) governing these models. These problems arise in many
applications including simulating Gaussian processes, geophysical
electromagnetics. The approach uses the Kato integral formula to express the
solution as an integral involving the solution of a parametrized elliptic PDE,
which is discretized using finite elements in space and sinc quadrature for
the fractional part. The offline stage of the ROM is accelerated using a
solver for shifted linear systems, MPGMRES-Sh, and using a randomized approach
for compressing the snapshot matrix. Our approach is both computational and
memory efficient. Numerical experiments on a range of model problems,
including an application to Gaussian processes, show the benefits of our
approach.

Keywords— Fractional elliptic PDEs, reduced order models, iterative methods,
randomization, Gaussian processes

##  1 Introduction

We want to derive efficient methods for the parameterized fractional partial
differential equation (PDE)

$(\mathcal{A}({\bf x};{\boldsymbol{\mu}}))^{\alpha}y({\bf
x};{\boldsymbol{\mu}})=b({\bf x};{\boldsymbol{\mu}}),\qquad{\bf x}\in\Omega$
(1)

where $\mathcal{A}$ is an elliptic differential operator, parameterized by the
set of parameters ${\boldsymbol{\mu}}\in\mathcal{P}\subset\mathbb{R}^{p}$.
Here $0<\alpha<1$ is the fractional exponent. Throughout this paper, unless
otherwise specified, we consider the boundary conditions to be homogeneous
Neumann, but the approach can be modified slightly to accommodate other
boundary conditions as well. We assume that the fractional power is defined in
the spectral sense. Solving the fractional elliptic PDE is a computationally
challenging problem in itself. Indeed, (1) is nonlocal and traditional ‘local
basis’ approaches such as finite element method require dealing with dense
numerical linear algebra. Therefore, the efficiency of traditional iterative
methods becomes questionable.

In this paper, we want to address the more challenging problem of efficiently
solving (1) for multiple instances of the parameter ${\boldsymbol{\mu}}$ and
the exponent $\alpha$. This is accomplished by the use of the reduced basis
approach. These kinds of problems arise in imaging [1], geophysical
electromagnetics [31], fractional Laplacian regularizers in optimization [4],
harmonic maps [2]. In [29], the authors investigate the well-posedness of a
Bayesian inverse problem with a fractional elliptic PDE forward model. To
estimate the diffusion coefficient from data, multiple instances of the
fractional PDE have to be solved. A similar situation occurs in optimal
control problems where the fractional exponent is considered as a random
variable [5], or optimization problems to identify the fractional exponent [6,
27]. Another application of interest is the stochastic PDE (SPDE) approach to
Gaussian processes, which is explored further in Section 4.1. In view of all
these applications, there is a clear need for an efficient reduced order model
(ROM).

In many applications of interest, high-fidelity models are governed by a
single PDE or systems of PDEs that are expensive to simulate. Reduced order
models refer to a collection of techniques for efficiently reducing the
complexity of high-fidelity models (also known as full-order models or FOMs)
with perhaps only a small sacrifice in accuracy compared to the FOMs. ROMs are
especially beneficial for problems in which the full-order model has to be
simulated for several parameter conditions as the “inner loop” for a variety
of applications such as optimization, control, inversion, uncertainty
quantification, design, or digital twins. This is known as parametric ROM and
this is the setting in the present paper.

While there are several possible techniques for ROMs, we follow the so-called
projection-based model reduction approach, in which there is a clear
separation into two stages. The first stage is called the offline stage in
which the FOM is simulated under a wide range of “representative” conditions
of the parameters to collect snapshots of the system. This is then followed by
a compression step to extract a basis that is used to then project the FOM
onto a lower-dimensional subspace yielding a ROM. In the online stage, the
lower-dimensional ROM is then simulated for a given parameter. The
computational tradeoff comes from the fact that the offline costs may be large
but upfront, and acceptable in certain situations, so long as the online costs
are small and the accuracy is good.

While ROM for elliptic PDE-based problems is a well-established field, ROM for
fractional elliptic PDEs is relatively less well-studied and is the main focus
of this paper. We briefly review the literature in this area and summarize the
main contributions of this paper.

#### Literature review

A recent paper [3] considers a reduced basis approach for the spectral
definition when realized using the Stinga-Torrea (Caffarelli-Silvestre)
extension [12, 28]. The articles [13, 14] employ a reduced basis approach by
interpolating operator norms. For another reduced ordering modeling approach
for problems with nonlocal integral kernels, we refer to [19]. The paper [18]
considers ROMs for time-fractional PDEs. The recent review paper [16] reviews
numerical methods, including reduced order models, for a different definition
of the nonlocal problem.

The following papers [15, 10] are closest to our setting, in which
$\mathcal{A}({\bf x};{\boldsymbol{\mu}})y({\bf x};{\boldsymbol{\mu}})=(-\Delta
y)({\bf x})$ is the standard Laplacian and $g({\bf
x};{\boldsymbol{\mu}})=g({\bf x})$. The reduced order model is constructed
over the exponent $\alpha$, but the parametric dependence of the PDE on other
parameters ${\boldsymbol{\mu}}$ is not considered. Notice that, these papers
build on the traditional deterministic reduced order method. Furthermore, in
our setting, we use randomized approaches, which offer additional benefits as
highlighted below. In recent work [7], we have used the MPGMRES-Sh approach
[8] for efficiently applying the fractional operator for a single parameter
value ${\boldsymbol{\mu}}$. It was observed in this paper that the same basis
can be used for multiple values of the exponent, which is important in our
approach as well.

We also briefly mention that the idea of using randomization in ROMs has
previously been explored; see, e.g., [25] and references within. The
randomized compression technique we use in this paper is applicable to the
streaming setting, in which several columns of the snapshot matrix are added
sequentially.

#### Contributions

In this paper, we derive reduced order models for solving (1). Our approach is
based on the following features:

  1. 1.

We use the spectral definition of the fractional elliptic operator and use the
Kato integral formula to express the solution of (2) as an integral involving
the solution of a parametric elliptic PDE with an additional parameter to
compute the integral. The operator is then discretized using standard Lagrange
finite elements, and the integral is discretized using the sinc quadrature
approach.

  2. 2.

We use a reduced basis approach for efficiently computing the solutions of the
parametric elliptic PDE. It contains an offline stage in which a reduced basis
is generated, and in the online stage, a projected system involving the
reduced basis is solved.

  3. 3.

The construction of the offline stage exploits the structure of the problem.
We use an efficient iterative method MPGMRES-Sh to accelerate the convergence
of shifted systems; this approach drastically reduces the number of PDE solves
needed for generating the basis. Furthermore, we use randomized solvers to
mitigate the cost of storing and compressing the snapshots, and an efficient
method involving generalized eigendecomposition to accelerate the costs of the
online stage.

As a part of future work, we plan to explore rigorous a posteriori error
analysis [20, 24] in the context of our randomized ROM approach.

##  2 Background

In this section, we discuss the integral approach to solving the fractional
elliptic PDE. We also discuss the ingredients required for the ROM, namely the
iterative solver MPGMRES-Sh (Section 2.2) and the randomized approach we use
for compressing the snapshots (Section 2.3).

###  2.1 Integral approach

Here we follow the approach in [15, Section 2.2] that combines an integral
formulation for the inverse fractional operator with a parameterized elliptic
PDE. Using the Kato integral formula [21, Theorem 2], with a simplification
$\lambda=0$, the solution to (1) can be expressed as

$y({\bf
x};{\boldsymbol{\mu}})=\frac{\sin(\alpha\pi)}{\pi}\int_{-\infty}^{\infty}e^{(1-\alpha)z}(e^{z}+\mathcal{A}({\bf
x};{\boldsymbol{\mu}}))^{-1}b({\bf x};{\boldsymbol{\mu}})dz.$

After setting ${\boldsymbol{\rho}}=({\boldsymbol{\mu}},z)$, we get the
parametric PDE

$(e^{z}+\mathcal{A}({\bf x};{\boldsymbol{\mu}}))u({\bf
x};{\boldsymbol{\rho}})=b({\bf x};{\boldsymbol{\mu}})\qquad{\bf x}\in\Omega,$
(2)

with appropriate boundary conditions, and arrive at the following integral
formula for $y$

$y({\bf
x};{\boldsymbol{\mu}})=\frac{\sin(\alpha\pi)}{\pi}\int_{-\infty}^{\infty}e^{(1-\alpha)z}u({\bf
x};{\boldsymbol{\rho}})dz.$ (3)

This integral approach is the foundation for the ROM that we construct.
Following [11], we use finite element discretization in space and sinc
quadrature for the integral, which we now discuss.

#### Discretization

We pick a finite-dimensional basis for the subspace $\mathcal{V}_{h}\subset
L^{2}(\Omega)$ denoted by $\\{\phi_{j}\\}_{j=1}^{N_{h}}$ so that
$\textrm{dim}(\mathcal{V}_{h})=N_{h}$ and any function $v_{h}\in V_{h}$ can be
represented by its expansion in the basis
$v_{h}=\sum_{j=1}^{N_{h}}v_{h,j}\phi_{j}$. We denote the vector of its
coefficients as ${\bf
v}=\begin{bmatrix}v_{h,1}&\dots&v_{h,N_{h}}\end{bmatrix}^{\top}$. For the sinc
quadrature to discretize the integral representation (3), let $h$ represent
the mesh parameter and define

$Z_{+}=\lceil\frac{\pi^{2}}{4\alpha\zeta^{2}}\rceil,\qquad
Z_{-}=\lceil\frac{\pi^{2}}{4(1-\alpha)\zeta^{2}}\rceil,$ (4)

where $\zeta=1/\log(1/h)$. We also define the quadrature weights and nodes as

$w_{j}=\frac{\zeta\sin(\alpha\pi)}{\pi}e^{(1-\alpha)z_{j}}\qquad z_{j}=\zeta
j\qquad-Z_{-}\leq j\leq Z_{+}.$ (5)

Let $a(u,v;{\boldsymbol{\mu}})$ denote the weak form of the operator
$\mathcal{A}({\bf x};{\boldsymbol{\mu}})$ with the corresponding boundary
conditions; similarly, let $g(v;{\boldsymbol{\mu}})$ be the weak form of
$b({\bf x};{\boldsymbol{\mu}})$. We seek a solution $u_{h,j}\in V_{h}$ such
that

$a(u_{h,j}({\bf x};{\boldsymbol{\mu}}),v_{h}({\bf
x});{\boldsymbol{\mu}})+e^{z_{j}}\int_{\Omega}u_{h,j}({\bf
x};{\boldsymbol{\mu}})v_{h}({\bf x})d{\bf
x}=g(v_{h};{\boldsymbol{\mu}}),\qquad\text{for all }v_{h}\in V_{h}$ (6)

and for $-Z_{-}\leq j\leq Z_{+}$. Having obtained the solutions $u_{h,j}$, the
approximate solution $y_{h}$ to (3) can be obtained using the sinc-quadrature
formula as $y_{h}({\bf
x};{\boldsymbol{\mu}})=\sum_{j=-Z_{-}}^{Z_{+}}w_{j}u_{h,j}({\bf
x};{\boldsymbol{\mu}})$. However, as can be seen, computing this quadrature
formula requires the solution of many PDEs of the form (6). This can be
computationally intensive if this process has to be repeated for multiple
values of the exponent $\alpha$ and the parameter ${\boldsymbol{\mu}}$. This
motivates the need for a reduced-order approach. We first discuss an efficient
method for solving the linear systems arising from the system (6).

The reference [11] establishes exponential convergence in quadrature and also
derives appropriate error estimates with respect to spatial discretization.
Provided that the so-called Kolmogorov $n$-width in the parameter space decays
exponentially, which is indeed the case when the solution depends smoothly on
parameters (cf. [23]), one can establish that reduced basis error also decays
exponentially. In the case that the reduced bases are generated using a greedy
method, we refer to [9] for such results.

###  2.2 MPGMRES-Sh

MPGMRES-Sh [8] is an iterative method for solving a sequence of shifted
systems of the form

$({\bf C}_{1}+\sigma_{j}{\bf C}_{2}){\bf x}_{j}={\bf b}\qquad 1\leq j\leq
N_{\sigma},$ (7)

for a set of shifts $\\{\sigma_{j}\\}_{j=1}^{N_{\sigma}}$ and matrices ${\bf
C}_{1},{\bf C}_{2}\in\mathbb{R}^{n\times n}$. The main idea of this approach
is to use shift-and-invert preconditioners of the form $\\{{\bf
C}_{1}+\tau_{j}{\bf C}_{2}\\}_{j=1}^{n_{p}}$, where $n_{p}$ is the number of
preconditioners used. A simplified version of the algorithm is given in
Algorithm 1.

0: Matrices ${\bf C}_{1},{\bf C}_{2}$, Shifts $\\{\sigma_{j}\\}$ for $1\leq
j\leq{N_{\sigma}}$, Vector ${\bf b}$, Preconditioners ${\bf P}_{j}={\bf
C}_{1}+\tau_{j}{\bf C}_{2}$, $1\leq j\leq n_{p}$

1: Compute $\beta=\|{\bf b}\|_{2}$ and $\widehat{\bf V}^{(1)}={\bf b}/\beta$

2: for $k=1,\dots,$ until convergence do

3: $\widehat{{\bf v}}_{k}=\widehat{\bf V}^{(k)}{\bf e}_{n_{p}}$ and ${\bf
W}^{(k)}=\begin{bmatrix}{\bf P}_{1}^{-1}\widehat{{\bf v}}_{k}&\dots&{\bf
P}_{n_{p}}^{-1}\widehat{{\bf v}}_{k}\end{bmatrix}$

4: $\widehat{{\bf W}}={\bf C}_{2}{\bf W}^{(k)}$

5: for $j=1,\dots,k$ do

6: ${\bf H}^{(j,k)}=(\widehat{\bf V}^{(j)})^{\top}\widehat{\bf W}$

7: $\widehat{{\bf W}}=\widehat{{\bf W}}-{\bf V}^{(j)}{\bf H}^{(j,k)}$

8: end for

9: $\widehat{\bf W}={\bf Q}{\bf R}$ {Thin QR factorization}

10: Construct $\widehat{\bf V}^{(k+1)}={\bf Q}$ and ${\bf H}^{(k+1,k)}={\bf
R}$

11: Check for convergence for all $\\{\sigma_{j}\\}_{j=1}^{N_{\sigma}}$ and
terminate, if converged, at iteration $k_{m}$.

12: end for

13: Consolidate matrices ${\bf Z}_{k_{m}}=\begin{bmatrix}{\bf
W}^{(1)}&\dots&{\bf W}^{(k_{m})}\end{bmatrix}$

14: return Basis ${\bf Z}_{k_{m}}$, where $k_{m}$ is the number of iterations
upon termination

Algorithm 1 Simplified version of MPGMRES-Sh

The idea behind this approach is to form a basis for the search space ${\bf
Z}_{k_{m}}\in\mathbb{R}^{n\times k_{m}n_{p}}$ where $k_{m}$ is the number of
MPGMRES-Sh iterations and at each iteration, we increment the dimension of the
search space by $n_{p}$. Sometimes the basis fails to be incremented by
$n_{p}$ but this can be handled using rank-revealing QR factorizations; we do
not consider this in the present work but the reader is referred to [8]. This
search space is common to all the shifted systems and is used to solve a
projected subproblem for each shift $\sigma$ of size
$(k_{m}n_{p})\times(k_{m}n_{p})$ to obtain a solution ${\bf
p}(\sigma)\in\mathbb{R}^{(k_{m}n_{p})\times 1}$. In a bit more detail, the
approximate solution is computed as ${\bf x}_{k_{m}}(\sigma)={\bf
Z}_{k_{m}}{\bf p}_{k_{m}}(\sigma)$, where ${\bf p}_{k_{m}}$ solves the
optimization problem

$\min_{{\bf p}\in\mathbb{R}^{(k_{m}n_{p})\times 1}}\|{\bf b}-({\bf
C}_{1}+\sigma_{j}{\bf C}_{2}){\bf Z}_{k_{m}}{\bf p}\|_{2}^{2},$

for each shift $\sigma_{j}$, $1\leq j\leq N_{\sigma}$. In the application of
interest in this present paper, we will focus on computing the basis ${\bf
Z}_{k_{m}}$ rather than the solutions ${\bf x}(\sigma)$ as we now explain in
this remark.

######  Remark 1.

Since ${\bf Z}_{k_{m}}$ is the search space corresponding to the shifted
systems (7), we can write ${\bf x}_{j}\approx{\bf Z}_{k_{m}}{\bf
p}_{k_{m}}(\sigma_{j})$. Therefore, the matrix

${\bf X}=\begin{bmatrix}{\bf x}_{1}&\dots&{\bf
x}_{N_{\sigma}}\end{bmatrix}\approx{\bf Z}_{k_{m}}\begin{bmatrix}{\bf
p}_{k_{m}}(\sigma_{1})&\dots&{\bf
p}_{k_{m}}(\sigma_{N_{\sigma}})\end{bmatrix}.$

That is, the range of ${\bf X}$ is approximately captured by the matrix ${\bf
Z}_{k_{m}}$. Since the leading left singular vectors approximate the range of
${\bf X}$, it is sufficient to store ${\bf Z}_{k_{m}}$. This key insight is
important to reduce computational costs in the offline stage.

The number of iterations is determined when all the relative residual for the
shifted systems is smaller than $10^{-8}$. For more implementation details and
theoretical justifications, the reader is referred to [8]. If the number of
MPGMRES-Sh iterations, $k_{m}$, and the number of preconditioners, $n_{p}$,
are both small, then the cost of solving the projected problem across multiple
shifts is small and the cost of MPGMRES-Sh is dominated by the cost of forming
the basis for the search space. In our applications, the number of
preconditioners is $3$ and the number of iterations is smaller than $30$
yielding a basis ${\bf Z}_{k_{m}}$ of size at most $90$. We now discuss the
second main ingredient in our approach, the randomized low-rank
approximations.

###  2.3 Randomized low-rank approximation

A reduced-order model requires compression of the matrix containing snapshots
of the state. In this paper, we pursue a randomized approach since it has
reduced computational and storage requirements compared to the full SVD.
Furthermore, as we will explain, this approach is beneficial in our setting
where the matrix is available in a streaming fashion.

Suppose we wanted to compute a low-rank approximation to the matrix ${\bf
S}\in\mathbb{R}^{m\times n}$. We follow the approach from [30]. We draw two
random matrices ${\boldsymbol{\Omega}}\in\mathbb{R}^{n\times\ell_{1}}$ and
${\boldsymbol{\Psi}}\in\mathbb{R}^{m\times\ell_{2}}$; we take them to be
standard Gaussian random matrices, which means that the entries are
independent and drawn from the Gaussian distribution with zero mean and unit
variance. We form the sketches ${\bf Y}_{1}={\bf S}{\boldsymbol{\Omega}}$ and
${\bf Y}_{2}={\bf S}^{\top}{\boldsymbol{\Psi}}$, which approximate the column
and row spaces of ${\bf S}$ respectively. Then compute the thin-QR
factorization ${\bf Y}_{1}={\bf Q}{\bf R}$ and discard ${\bf R}$. The low-rank
approximation of ${\bf S}$ can be computed as

${\bf S}\approx{\bf S}_{\rm rand}\equiv{\bf Q}({\boldsymbol{\Psi}}^{\top}{\bf
Q})^{\dagger}{\bf Y}_{2}^{\top},$

where † is the Moore-Penrose inverse. The low-rank approximation can be
converted into the SVD format to extract singular vectors with relatively less
additional work. If a rank-$r$ approximation of the matrix ${\bf S}$ is
desired, assuming that the matrix is stored in dense form, the computational
cost of this algorithm is $\mathcal{O}(rmn+r^{2}(m+n)+r^{3})$ flops. The
choice of the sketch sizes $\ell_{1}$ and $\ell_{2}$ is discussed in [30,
Section 4.4] which we follow here as well. We choose $\ell_{1}=2r+1$ and
$\ell_{2}=2\ell_{1}+1$. Based on this sketch size, the error in expectation of
the low-rank approximation satisfies [30, Theorem 4.3]

$\mathbb{E}\|{\bf S}-{\bf S}_{\rm rand}\|_{F}^{2}\leq
4\left(\sum_{j>r}\sigma_{j}^{2}\right),$

where $\\{\sigma_{j}\\}_{j=1}^{\min\\{m,n\\}}$ are the singular values of the
snapshot matrix ${\bf S}$. The main message is that the error in the low-rank
approximation due to the randomization is, in expectation, within a small
factor of the “best” low-rank approximation (note that this approximation need
not be unique). In particular, if the singular values of the snapshot matrix
decay rapidly then the error in the low-rank approximation is very small.

Besides lower computational costs compared to the full SVD algorithm (i.e.,
$\mathcal{O}(mn\min\\{m,n\\})$), the proposed approach is amenable to the
streaming setting where the matrix is updated sequentially. That is, if the
matrix ${\bf S}$ is of the form

${\bf S}={\bf S}_{1}+\dots+{\bf S}_{M},$

where the summands ${\bf S}_{j}$ are available sequentially. As an example
that is relevant to our application, consider the case, where the columns (or
a block of columns) of a matrix are computed sequentially. It is easy to see
that the sketches ${\bf Y}_{1}$ and ${\bf Y}_{2}$ can also be updated
sequentially. Furthermore, once the sketches have been computed, the
appropriate columns can be discarded since all the randomized algorithm only
requires the sketches to compute the low-rank approximation. These insights
will be relevant to the proposed algorithm resulting in an algorithm with
lower storage and memory costs.

##  3 Reduced basis approaches

In this section, we introduce the weak form and discretization of the
parametric PDE that is the foundation of our ROM (Section 2). Then we present
two approaches for ROM: a naïve approach (Section 3.2), and an efficient
approach (Section 3.3). We conclude this section with a discussion on
computational costs.

###  3.1 Weak form and discretization

We aim to solve the parameterized PDE (2) for multiple instances of the
parameters ${\boldsymbol{\rho}}$. To this end, with an eye towards numerical
approximation using finite elements, we consider the weak form of the problem,
which reads as: given ${\boldsymbol{\rho}}\in\mathcal{P}\times\mathbb{R}_{+}$
find $u({\boldsymbol{\rho}})\in V$ such that

$a(u({\boldsymbol{\rho}}),v;{\boldsymbol{\rho}})=g(v;{\boldsymbol{\rho}})\qquad\forall
v\in V,$

where the exact form of $a$ depends on the form of $\mathcal{A}$ and the
boundary conditions. We assume that the weak form has an affine parameter
dependence as

$a(u,v;{\boldsymbol{\rho}}):=\sum_{t=1}^{n_{a}}f_{t}^{a}({\boldsymbol{\mu}})a_{t}(u,v)+e^{z}\int_{\Omega}uvd{\bf
x},$ (8)

where $a_{t}:V\times V\rightarrow\mathbb{R}$ are independent of the parameters
and the right-hand side also satisfies an affine parametric dependence

$g(v;{\boldsymbol{\rho}}):=\sum_{t=1}^{n_{g}}f_{t}^{g}({\boldsymbol{\mu}})g_{t}(v).$
(9)

If the problem does not admit an affine parametric dependence, one can use
techniques such as the empirical interpolation method [24, Chapter 10].

We proceed to discuss the matrix representation of this parametric PDE that
will be starting point of the ROMs. As before, let $V_{h}\subset V$ be a
finite-dimensional subspace of $V$ with $\dim(V_{h})=N_{h}$ with a basis
$\\{\phi_{1},\dots,\phi_{N_{h}}\\}$. We denote the matrix representation of
the bilinear forms as

$[{\bf A}_{t}]_{ij}=a_{t}(\phi_{i},\phi_{j})\qquad 1\leq i,j\leq N_{h},\quad
1\leq t\leq n_{a}.$ (10)

We denote by ${\bf M}$, the mass matrix ${\bf M}\in\mathbb{R}^{N_{h}\times
N_{h}}$ with entries $[{\bf M}]_{ij}=\int_{\Omega}\phi_{i}\phi_{j}d{\bf x}$
for $1\leq i,j\leq N_{h}$. We also define the vectors ${\bf g}_{t}$ with
entries $[{\bf g}_{t}]_{i}=g_{t}(\phi_{i})$ for $1\leq i\leq N_{h}$, and
$1\leq t\leq n_{g}$. In summary, the discretized parametric system we aim to
solve is

$\left(\sum_{t=1}^{n_{a}}f_{t}^{a}({{\boldsymbol{\mu}}}){{\bf
A}}_{t}+e^{z}{\bf M}\right){\bf
u}=\sum_{t=1}^{n_{g}}f_{t}^{g}({\boldsymbol{\mu}}){\bf g}_{t},$ (11)

for the set of parameters ${\boldsymbol{\rho}}=({\boldsymbol{\mu}},z)$. We now
discuss ROMs for the parametric PDEs.

###  3.2 A naïve first approach

Before we present our efficient approach, we first discuss a naïve approach
that will be the basis for comparison with our proposed approach.

We first draw a set of random instances of the parameters
$\mathcal{M}=\\{{\boldsymbol{\mu}}^{(j)}\\}$ for $1\leq j\leq M$; strategies
for sampling are discussed in [24, Section 6.6-6.7]. Similarly, using the
definitions in Section 2.1 we set the parameters $\mathcal{Z}=\\{z_{j}\\}$ for
$-Z_{-}\leq j\leq Z_{+}$ with $Z=Z_{-}+Z_{+}+1$ corresponding to the setting
$\alpha=0.5$. This is motivated by the observation in [7, Section 5.1] that
the exponents $e^{z_{k}}$ are completely independent of the exponent, except
for the limits $Z_{-}$ and $Z_{+}$. As we will see in the numerical
experiments, it is sufficient to generate the projection basis corresponding
to a single exponent $\alpha$, which is then applicable to a wider range of
the exponent.

Thus, we now have the set of parameters $\mathcal{M}\times\mathcal{Z}$ for
${\boldsymbol{\rho}}$. For each parameter value ${\boldsymbol{\mu}}^{(j)}$ we
first construct the matrices and vectors

${\bf K}_{j}=\sum_{t=1}^{n_{a}}f_{t}^{a}({{\boldsymbol{\mu}}}_{j}){\bf
A}_{t},\qquad{\bf
f}_{j}=\sum_{t=1}^{n_{g}}f_{t}^{g}({{\boldsymbol{\mu}}}_{j}){\bf g}_{t},\quad
1\leq j\leq M.$ (12)

Then we solve the sequence of shifted linear equations

$({\bf K}_{j}+e^{z_{k}}{\bf M}){\bf u}_{j,k}={\bf f}_{j}\qquad 1\leq j\leq
M,\quad-Z_{-}\leq k\leq Z_{+},$ (13)

to generate the snapshot matrix

${\bf S}=\begin{bmatrix}{\bf U}^{(1)}&\dots&{\bf
U}^{(M)}\end{bmatrix}\in\mathbb{R}^{N_{h}\times(MZ)},\qquad{\bf
U}^{(j)}=\begin{bmatrix}{\bf u}_{1,-Z_{-}}&\dots&{\bf
u}_{j,Z_{+}}\end{bmatrix}\qquad 1\leq j\leq M.$ (14)

The next step is to compute the dominant left singular vectors of ${\bf S}$
denoted ${\bf V}\in\mathbb{R}^{N_{h}\times K}$. We form the matrices
$\widehat{\bf A}_{t}={\bf V}^{\top}{\bf A}_{t}{\bf V}$ for $1\leq t\leq
n_{a}$, $\widehat{\bf M}={\bf V}^{\top}{\bf M}{\bf V}$, and vectors
$\widehat{\bf g}_{t}={\bf V}^{\top}{\bf g}_{t}$ for $1\leq t\leq n_{g}$.

This basis ${\bf V}$ can be used to approximate the solution space in the
online stage. Given a new parameter value $\overline{{\boldsymbol{\mu}}}$, and
an exponent $\alpha\in(0,1)$, we solve the projected system (see (4) for
definitions of $Z_{+},Z_{-}$)

$\left(\sum_{t=1}^{n_{a}}f_{t}^{a}(\overline{{\boldsymbol{\mu}}})\widehat{{\bf
A}}_{t}+e^{z_{k}}\widehat{{\bf M}}\right)\widehat{\bf
u}_{k}=\left(\sum_{t=1}^{n_{g}}f_{t}^{g}(\overline{{\boldsymbol{\mu}}})\widehat{{\bf
g}}_{t}\right)\qquad-Z_{-}\leq k\leq Z_{+}.$ (15)

The approximation solution is $\widehat{{\bf y}}={\bf
V}\left(\sum_{k=-Z_{-}}^{Z_{+}}w_{k}\widehat{\bf u}_{k}\right)$. The details
of this approach are given in Algorithm 2.

0: Sample parameter $\overline{{\boldsymbol{\mu}}}$, Precomputed matrices
$\\{\widehat{\bf A}_{t}\\}_{t=1}^{n_{a}}$ and vectors $\\{{\bf
g}_{t}\\}_{t=1}^{n_{g}}$, exponent $\alpha$

1: Compute quadrature nodes $\\{z_{k}\\}_{-Z_{-}}^{Z_{+}}$ and weights
$\\{w_{k}\\}_{-Z_{-}}^{Z_{+}}$ corresponding to the exponent $\alpha$; see (4)
and (5).

2: Form and solve the projected system (15)

3: Evaluate $\widehat{\bf w}={\bf
V}\left(\sum_{k=-Z_{-}}^{Z_{+}}w_{k}\widehat{\bf u}_{k}\right)$.

4: return Approximate solution $\widehat{\bf w}$

Algorithm 2 Reduced basis approach for fractional elliptic PDE: Online stage

However, the cost of generating and compressing the snapshot matrix can be
prohibitively expensive in practice as we now explain. This is because $MZ$
linear systems need to be solved in total to compute the snapshot matrix. Note
that $M$ can range from $10-100$ whereas $Z$ can range from $100-1000$,
depending on the value of the exponent $\alpha$ and the discretization of the
system. Furthermore, the size of ${\bf S}$ can be large so the cost of
compressing this system is large; the cost of SVD is
$\mathcal{O}(N_{h}MZ\min\\{N_{h},MZ\\})$. Even the cost of storing the matrix
may be expensive in practice. To address these issues we propose a new
approach for efficiently generating and compressing the basis.

###  3.3 Proposed approach

There are two components to the proposed approach in the offline stage: the
use of MPGMRES-Sh to solve the shifted systems and randomization to compress
the snapshot matrix. In the online stage, we use a generalized
eigendecomposition to speed up computational costs.

#### Offline stage

First, we solve the systems (13) for each value of ${\boldsymbol{\mu}}^{(j)}$
using the MPGMRES-Sh approach (Section 2.2). The algorithm is implemented with
${\bf C}_{1}={\bf M}$ and ${\bf C}_{2}={\bf A}_{j}$ and shifts
$\sigma_{k}=1/\exp(z_{k})$ for $-Z_{-}\leq k\leq Z_{+}$. We get a basis ${\bf
Z}^{(j)}\in\mathbb{R}^{N_{h}\times(n_{m})}$ where $n_{m}$ is the size of the
MPGMRES-Sh basis (obtained as $n_{m}=k_{m}n_{p}$, i.e., the number of MPGMRES-
Sh iterations $k_{m}$ times the number of preconditioners $n_{p}$). For
simplicity, we have assumed that the MPGMRES-Sh basis size is the same for
every $j$ but in practice, it may be different. We have also dropped the
dependence of the basis on the number of iterations $k_{m}$. In all the
experiments, for MPGMRES-Sh we take three preconditioners $n_{p}=3$ with the
corresponding values of $\tau\in\\{10^{-8},10^{-4},10^{-2}\\}$. This is the
same setup as in [7] and we refer to this paper for additional implementation
details. The preconditioners are factorized and their factorizations are
stored.

After solving the $j$th set of systems, we have the basis ${\bf Z}^{(j)}$ with
which we can approximate the solutions of ${\bf U}^{(j)}$; see (14).
Therefore, rather than compute the matrix ${\bf U}^{(j)}$ and then compress it
subsequently, following Remark 1, we directly compute the new “snapshot”
matrix as $\widehat{\bf S}=\begin{bmatrix}{\bf Z}^{(1)}&\dots&{\bf
Z}^{(M)}\end{bmatrix}\in\mathbb{R}^{N_{h}\times Mn_{m}}$. This is because the
range of ${\bf U}^{(j)}$ is approximately captured by the range of ${\bf
Z}^{(j)}$, since ${\bf Z}^{(j)}$ is the search space for the set of shifted
systems. Although it is not technically a matrix of snapshots, we still refer
to $\widehat{\bf S}$ as a snapshot matrix. Notice that the number of columns
of $\widehat{\bf S}$ is much smaller than that of ${\bf S}$ if $n_{m}\ll Z$,
as is the case in numerical experiments.

Even with this reduction, the cost of storing and computing the dominant
singular vectors can be expensive. To alleviate these costs, we use the single
view randomized approach as described in Section 2.3. To this end, we consider
two sketching matrices
${\boldsymbol{\Omega}}\in\mathbb{R}^{(Mn_{m})\times\ell_{1}}$ and
${\boldsymbol{\Psi}}\in\mathbb{R}^{N_{h}\times\ell_{2}}$ which are standard
Gaussian random matrices. Given the target basis size $K$, we choose
$\ell_{1}=2K+1$ and $\ell_{2}=2\ell_{1}+1.$ We observe that the matrix
$\widehat{{\bf S}}$ can be written as

$\widehat{{\bf S}}=\widehat{{\bf S}}_{1}+\dots+\widehat{{\bf
S}}_{M}\qquad\widehat{{\bf S}}_{j}={\bf Z}^{(j)}{\bf E}_{j}^{\top},$

where, using MATLAB notation, ${\bf E}_{j}={\bf
I}_{Mn_{m}}(:,(j-1)*n_{m}+1:j*n_{m})$ for $1\leq j\leq M$ contains columns
from the $(Mn_{m})\times(Mn_{m})$ identity matrix. Therefore, we can easily
update the sketches ${\bf Y}_{1}$ and ${\bf Y}_{2}$ for each $j$ and then
discard the bases ${\bf Z}^{(j)}$. In fact, the matrix ${\boldsymbol{\Omega}}$
does not need to be formed all at once; this is beneficial since the basis
size $n_{m}$ may not be known in advance. Consider
${\boldsymbol{\Omega}}=\begin{bmatrix}{\boldsymbol{\Omega}}_{1}^{\top}&\dots&{\boldsymbol{\Omega}}_{M}^{\top}\end{bmatrix}^{\top}$,
where ${\boldsymbol{\Omega}}_{j}\in\mathbb{R}^{n_{m}\times\ell_{1}}$. Then the
sketch ${\bf Y}_{1}$ can be formed as

${\bf Y}_{1}={\bf Z}^{(1)}{\boldsymbol{\Omega}}_{1}+\dots+{\bf
Z}^{(M)}{\boldsymbol{\Omega}}_{M}.$

The details of this approach are given in Algorithm 3.

#### Online stage

In the online stage, the steps involve solving the $K\times K$ systems for
$Z=Z_{+}+Z_{-}+1$, see (15) which can be expensive. Rather than factorizing
each linear system, we can use a generalized eigendecomposition. To this end,
define $\widehat{\bf
K}=\sum_{j=1}^{n_{a}}f_{t}^{a}(\overline{{\boldsymbol{\mu}}})\widehat{{\bf
A}}$, $\widehat{\bf
g}=\sum_{j=1}^{n_{g}}f_{t}^{g}(\overline{{\boldsymbol{\mu}}})\widehat{\bf g}$,
and consider a generalized eigendecomposition of the pair $(\widehat{\bf
K},\widehat{\bf M})$; that is, compute a Cholesky factorization $\widehat{\bf
M}={\bf L}{\bf L}^{\top}$ and compute an eigendecomposition ${\bf
L}^{-1}\widehat{\bf K}{\bf L}^{\top}={\bf U}{\boldsymbol{\Lambda}}{\bf
U}^{\top}$. We are guaranteed such an eigendecomposition since both ${\bf M}$
and $\widehat{\bf M}$ are symmetric positive definite. Then, we can obtain the
solution $\widehat{\bf y}$ as

$\widehat{\bf y}=({\bf V}{\bf L}^{-\top}{\bf
U})\left(\sum_{k=-Z_{-}}^{Z_{+}}w_{k}({\boldsymbol{\Lambda}}+e^{z_{k}}{\bf
I})^{-1}\right){\bf U}^{\top}{\bf L}^{-1}\widehat{\bf g}.$ (16)

The modification to Algorithm 2 is straightforward and therefore omitted.

0: Parameter samples
$\mathcal{M}=\\{{\boldsymbol{\mu}}^{(1)},\dots,{\boldsymbol{\mu}}^{(M)}\\}$,
basis size $K$

Assemble the matrices ${\bf A}_{t}$ for $1\leq t\leq n_{a}$ and vectors ${\bf
g}_{t}$ for $1\leq t\leq n_{g}$.

Compute the nodes $z_{k}={k\zeta}$ where $\zeta=1/\log(1/h)$ and $-Z_{-}\leq
k\leq Z_{+}$ where $Z_{+}=\lceil\frac{\pi^{2}}{2\zeta^{2}}\rceil$,
$Z_{-}=\lceil\frac{\pi^{2}}{2\zeta^{2}}\rceil$

Draw random matrix ${\boldsymbol{\Psi}}\in\mathbb{R}^{N_{h}\times\ell_{2}}$

Initialize sketches ${\bf Y}_{1}={\bf 0}\in\mathbb{R}^{N_{h}\times\ell_{1}}$,
${\bf Y}_{2}=[]$

{Stage 1: Constructing sketches of the snapshot matrix}

for $j\in\\{-Z_{-},\dots,Z_{+}\\}$ do

Form matrix ${\bf
K}_{j}=\sum_{t=1}^{n_{a}}f_{t}^{a}({{\boldsymbol{\mu}}}_{j}){\bf A}_{t}$ and
vector ${\bf f}_{j}=\sum_{t=1}^{n_{g}}f_{t}^{g}({{\boldsymbol{\mu}}}_{j}){\bf
g}_{t}$

Solve the systems of linear equations (13) using MPGMRES-Sh to generate a
basis ${\bf Z}^{(j)}\in\mathbb{R}^{N_{h}\times n_{m}}$.

Draw random matrix
${\boldsymbol{\Omega}}_{j}\in\mathbb{R}^{n_{m}\times\ell_{1}}$

Update sketches ${\bf Y}_{1}={\bf Y}_{1}+{\bf
Z}^{(j)}{\boldsymbol{\Omega}}_{j}$ and ${\bf Y}_{2}=[{\bf Y}_{2};({\bf
Z}^{(j)})^{\top}{\boldsymbol{\Psi}}]$

end for

{Stage 2: Compression of snapshot matrix}

Compute thin QR ${\bf Y}_{1}={\bf Q}{\bf R}$

Compute ${\bf W}=({\boldsymbol{\Psi}}^{\top}{\bf Q})^{\dagger}{\bf
Y}_{2}^{\top}$ and its thin SVD ${\bf U}_{W}{\boldsymbol{\Sigma}}_{W}{\bf
V}_{W}^{\top}$

Compute ${\bf V}={\bf Q}{\bf U}_{W}(:,1:K)$ where $K$ is the basis size

Precompute $\widehat{{\bf A}}_{t}={\bf V}^{\top}{\bf A}_{t}{\bf V}$ for $1\leq
t\leq n_{a}$, $\widehat{{\bf M}}={\bf V}^{\top}{\bf M}{\bf V}$, and vector
$\widehat{{\bf g}}_{t}={\bf V}^{\top}{\bf M}{\bf w}$

return Basis vectors ${\bf V}$, matrices $\\{\widehat{{\bf
A}}_{t}\\}_{t=1}^{n_{a}}$ and vectors $\\{\widehat{{\bf
g}}_{t}\\}_{t=1}^{n_{g}}$

Algorithm 3 Reduced basis approach for fractional elliptic PDE: Offline stage

###  3.4 Computational Cost

We now summarize the computational cost of the offline and online stages.

In the offline stage, there are two sources of computational cost: the
snapshot matrix construction and compression. The construction essentially
involves the cost of applying MPGMRES-Sh to $M$ set of $Z$ shifted systems.
Let $m_{p}$ be the number of preconditioners and $k_{m}$ be the number of
MPGMRES-Sh iterations and their product is $n_{m}=k_{m}n_{p}$; then the cost
of one set of MPGMRES-Sh iterations is

$T_{\rm MP}=n_{p}T_{\rm PC}+T_{\rm
PA}n_{m}+\mathcal{O}(N_{h}Zn_{m}^{2}+Zn_{m}^{3}))\quad\text{flops},$

where $T_{\rm PC}$ is the cost of constructing each preconditioner and $T_{\rm
PA}$ is the cost of applying the combination of the matrix ${\bf C}_{2}$ and
the preconditioner ${\bf P}_{j}$ (lines 3-4 of Algorithm 1). The total cost of
forming the snapshot matrix is $T_{\rm S}=MT_{\rm MP}$. The cost of the
randomized low-rank approximation is analyzed in Section 2.3; assuming the
basis size $K$, the cost is $T_{\rm
R}=\mathcal{O}(N_{h}n_{m}MK+K^{2}(N_{h}+n_{m})+K^{3})$ flops. The total cost
of the offline stage is then simply the sum of $T_{\rm S}$ and $T_{\rm R}$.

For the offline stage, the naïve approach requires us to solve $Z$ systems
with $\mathcal{O}(K^{3})$ cost per system, with a total cost of
$\mathcal{O}(ZK^{3}+N_{h}K^{2})$ flops. The proposed approach only requires
$\mathcal{O}(K^{3})$ cost for the eigendecomposition and $\mathcal{O}(K^{2})$
flops per system solve for a total cost $\mathcal{O}(K^{3}+K^{2}(Z+N_{h}))$
flops. If $K$ is large, computing the eigendecomposition is computationally
expensive or even prohibitive. Then one can use techniques from [17] to
further reduce the computational cost. Another option is to view (15) through
the lens of a generalized Sylvester equation; one can then use techniques from
[26] to reduce the computational costs.

##  4 Numerical Experiments

In this section, we explore the performance of the proposed ROMs on three
illustrative applications. The relative error is measured in the $L^{2}$ norm
with a reference finite element solution computed using MPGMRES-Sh. The finite
element discretization was carried out using continuous piecewise linear
elements. The MPGMRES-Sh solver has been validated in [7, Section 5.1]. The
timing results were performed on a Mac Mini 2020, with $16$ GB RAM and running
MATLAB 2021a.

###  4.1 Gaussian process

In this application, we consider simulation from the Gaussian process, by
solving the stochastic PDE on the domain $\Omega\in\mathbb{R}^{d}$

$(\kappa^{2}-\Delta)^{(\nu+d/2)/2}u({\bf x})=\xi({\bf x}),\qquad{\bf
x}\in\Omega$ (17)

where $\xi$ is the white noise process, and in the exponent $\nu$ is the
smoothness parameter and $d$ the spatial dimension. The parameter $\kappa^{2}$
is related to the correlation length of the process. We want to compute the
solution of the SPDE for multiple values of the parameters $\kappa^{2}$ and
$\alpha=(\nu+d/2)/2$. For this application, there is a single parameter
${\boldsymbol{\mu}}=(\kappa^{2})$. Therefore, $n_{a}=2$ with
$f_{1}^{a}({\boldsymbol{\mu}})=1$ and
$f_{2}^{a}({\boldsymbol{\mu}})=\kappa^{2}$; similarly, $n_{g}=1$ with
$f_{1}^{g}({\boldsymbol{\mu}})=1$. The weak forms are
$a_{1}=\int_{\Omega}\nabla u\cdot\nabla vd{\bf x}$,
$a_{2}=\int_{\Omega}uvd{\bf x}$, and $g_{1}=\int_{\Omega}\xi vd{\bf x}$.

Figure 1: (left) Singular values of the snapshot matrix $\widehat{\bf S}$.
(right) The relative error ($\log_{10}$) in the test set.

#### Experiment 1: Accuracy

We take the domain $\Omega=(0,1)^{2}$ and consider homogeneous Neumann
boundary conditions for the SPDE. The domain is discretized using a grid of
size $257\times 257=66,049$ points. We use a basis size of $K=100$ based on
the decay of the singular values. We generate a training set of parameters
$\kappa^{2}\in[10,200]$ in increments of $5$, giving a training set of size
$|\mathcal{M}|=39$. We investigate the accuracy of the approach for several
instances of $\kappa^{2}$ and $\alpha$. Specifically, we generate a test set
involving $20$ values of $\kappa^{2}\in[10,200]$ and $\alpha\in[0.1,0.9])$ in
increments of $0.1$ giving a test set of $180$ candidates. As can be seen from
Figure 1, the relative error is small throughout the test set with errors
ranging from $10^{-9}-10^{-7}$. The errors show a mild increase with the
exponent $\alpha$.

#### Experiment 2: Timings

In this next experiment, we investigate the reduction in computational costs.
We focus on the timing of the compression. The time of solving (13) for a
single parameter sample ${\boldsymbol{\mu}}^{(j)}$ using the naïve approach is
$23.33$ seconds whereas for MPGMRES-Sh the cost is $1.28$ seconds (this time
includes the time to factorize the preconditioners). Furthermore, in the naïve
approach, the size of the snapshot matrix ${\bf S}$ is $66,049\times 11895$,
whereas in the proposed approach the size of the snapshot matrix is
$66,049\times 1872$. The size of the MPGMRES-Sh basis for each
${\boldsymbol{\mu}}^{(j)}$ was $48$. Since the cost of forming and compressing
the snapshot matrix is enormous, we do not consider the naïve approach for
timing purposes anymore. For the baseline approach, we use the approach in
Section 3.3 but use SVD instead of the randomized approach.

| Baseline | Proposed
---|---|---
Snapshot | $47.03$ | $46.66$
Compression | $18.69$ | $0.59$
Table 1: Wall-clock time (in seconds) of the time for computing the snapshot
matrix $\widehat{\bf S}$ and the compression stages.

We report both the time to construct the snapshot matrices $\widehat{\bf S}$
and the time for compression (SVD vs randomization). The results are reported
in Table 1. Both methods give a comparable time for constructing the matrix.
However, the baseline approach takes slightly more time (since the total
number of MPGMRES-Sh iterations is not known in advance, and memory has to be
allocated on-the-fly which makes it inefficient). Note that the proposed
approach does not form the matrix $\widehat{\bf S}$ and the reported time
includes the time to form the sketches ${\bf Y}_{1}$ and ${\bf Y}_{2}$.
However, when it comes to the compression costs, our approach shows a
$30\times$ speedup over the SVD. This speedup will be more significant for
larger-scale problems.

Finally, we time for evaluating the ROM for $\alpha=0.5$ is $0.01$ seconds,
which is much more efficient than both the MPGMRES-Sh (1.28 seconds) and the
naïve approaches (23.33 seconds).

###  4.2 Fractional cookies

In this experiment, inspired by the cookies test problem from [22], we
consider the fractional parameterized PDE (2) with the domain
$\Omega=(-1,1)^{2}$ with Dirichlet boundary conditions, and the elliptic
operator is given by $\mathcal{A}({\bf
x};{\boldsymbol{\mu}})=-\nabla\cdot(D({\bf
x};{\boldsymbol{\mu}})\nabla(\cdot))$ where

$D({\bf x};{\boldsymbol{\mu}})=1+\sum_{t=1}^{p}\mu_{j}\chi_{\Omega_{t}}({\bf
x})\qquad$

where $\chi_{A}$ is the characteristic function (that takes the value $1$ if
${\bf x}\in A\subset\Omega$ and $0$ otherwise), and $p$ is the number of
regions (“cookies”). We consider two different cases for the discs:

Figure 2: (left) Visualization of the diffusion coefficient $D({\bf
x};{\boldsymbol{\mu}})$; (right) first $100$ singular values of the snapshot
matrix $\widehat{\bf S}$.

  1. 1.

Case A: the number of discs $p=1$ and $\Omega_{1}$ is a disc of radius $0.5$
and center $(0,0)$.

  2. 2.

Case B: the number of discs $p=4$ and $\Omega_{j}$ are discs of radii $0.3$
with centers $(-0.5,-0.5)$, $(-0.5,0.5)$, $(0.5,0.5)$, and $(0.5,-0.5)$.

See Figures 2 and 3 for visualizations of the diffusion coefficient $D$.

We take the values of ${\boldsymbol{\mu}}\in[0,1]^{p}$. The fractional PDE is
augmented with zero Dirichlet boundary conditions. For this problem,
$n_{a}=p+1$ with $a_{t}(u,v)=\int_{\Omega_{t}}\nabla u\cdot\nabla vd{\bf x}$
for $1\leq t\leq p$ and $a_{5}=\int_{\Omega}\nabla u\cdot\nabla vd{\bf x}$;
similarly, $n_{g}=1$ and $g_{1}(v)=\int_{\Omega}vd{\bf x}$.

Figure 3: (left) Visualization of the diffusion coefficient $D({\bf
x};{\boldsymbol{\mu}})$ for Case $B$; (right) first $700$ singular values of
the snapshot matrix $\widehat{\bf S}$.

We discretize the domain using a grid of size $257^{2}$. To generate the
training set, we pick $M=100$ parameters ${\boldsymbol{\mu}}^{(j)}$ from
$[0,1]^{p}$ sampled using Latin hypercube sampling. Based on the decay of the
singular values, we take the reduced basis size to be $K=100$ (Case A) and
$K=700$ (Case B). Compared to the test problem in Section 4.1, the singular
values of the snapshot matrix do not decay as sharply in Case B (see Figure
3). For the test set, similarly, we consider $10$ parameters uniformly chosen
at random from the range ${\boldsymbol{\mu}}\in[0,1]^{p}$ and for values of
$\alpha\in[0.1,0.9]$. Figure 4 shows the distribution of relative errors for
Case A and Case B.

Figure 4: Relative error in the test set (left) Case A, (right) Case B. The
mean of the error over $10$ random samples is plotted for different values of
$\alpha$, with the error bar denoting one standard deviation.

In Table 2, we display the computational time for the offline and online
stages. As can be seen, for both cases, the online computational time
outperforms the Direct and the MPGMRES-Sh algorithms. However, the offline and
online times of Case B are higher since the reduced basis has a higher rank
$K$. Note that when the reduced basis rank $K$ is high, the ROM may offer
computational benefits over the conventional approaches only if $N_{h}$ is
sufficiently high.

| Snapshot | Compression | Direct | MPGMRES-Sh | Online
---|---|---|---|---|---
Case A | $143.72$ | $0.67$ | $71.83$ | $1.48$ | $0.006$
Case B | $236.49$ | $20.51$ | $71.71$ | $1.93$ | $0.37$
Table 2: Wall-clock time (in seconds) for various components of the cookies
test problem.

###  4.3 Anisotropic diffusion

In this example, we consider the fractional PDE (1) with the anisotropic
diffusion model $\mathcal{A}({\bf
x};{\boldsymbol{\mu}})=-\nabla\cdot({\boldsymbol{\Theta}}({\boldsymbol{\mu}})\nabla(\cdot))$,
where ${\boldsymbol{\Theta}}$ takes the form

${\boldsymbol{\Theta}}({\boldsymbol{\mu}})=\begin{bmatrix}\cos\theta&-\sin\theta\\\
\sin\theta&\cos\theta\end{bmatrix}\begin{bmatrix}D_{1}&\\\
&D_{2}\end{bmatrix}\begin{bmatrix}\cos\theta&-\sin\theta\\\
\sin\theta&\cos\theta\end{bmatrix}^{\top}.$

The fractional PDE is posed on the domain $\Omega=(0,1)^{2}$ and is augmented
with zero Dirichlet boundary conditions. The parameters
${\boldsymbol{\mu}}=(D_{1},D_{2},\theta)$ include the two diffusion
coefficients and the angle $\theta$.

We can write the weak form as
$a(u,v;{\boldsymbol{\rho}})=\sum_{t=1}^{3}f_{t}({\boldsymbol{\mu}})a_{t}(u,v)$,
where $a_{1}(u,v)=\int_{\Omega}{\bf e}_{1}{\bf e}_{1}\nabla u\cdot\nabla
vd{\bf x}$, $a_{2}(u,v)=\int_{\Omega}({\bf e}_{1}{\bf e}_{2}^{\top}+{\bf
e}_{2}{\bf e}_{1}^{\top})\nabla u\cdot\nabla vd{\bf x}$, and
$a_{3}(u,v)=\int_{\Omega}{\bf e}_{2}{\bf e}_{2}^{\top}\nabla u\cdot\nabla
vd{\bf x}$. Here, ${\bf e}_{1},{\bf e}_{2}$ are the two columns from the
$2\times 2$ identity matrix. Furthermore,
$f_{1}^{a}({\boldsymbol{\mu}})=(\cos^{2}\theta)D_{1}+(\sin^{2}\theta)D_{2}$,
$f_{2}^{a}({\boldsymbol{\mu}})=\sin\theta\cos\theta(D_{2}-D_{1})$, and
$f_{3}^{a}({\boldsymbol{\mu}})=(\cos^{2}\theta)D_{2}+(\sin^{2}\theta)D_{1}$.
We take $b({\bf x};{\boldsymbol{\mu}})=1$ so that
$f_{1}^{g}({\boldsymbol{\mu}})=1$ and $g_{1}(v)=\int_{\Omega}vd{\bf x}$.

Figure 5: (left) Singular values of the snapshot matrix $\widehat{\bf S}$,
(right) the relative error in the test set. The mean of the error over $10$
random samples is plotted for different values of $\alpha$, with the error bar
denoting one standard deviation.

For the training set, we generate $100$ samples generated uniformly at random
from the ranges $(D_{1},D_{2},\theta)\in[0.5,4.5]^{2}\times[0,\pi/2]$. The
basis size $K$ is taken as $300$. The test set is drawn randomly similarly to
the training set and the maximum error over $90$ samples is $\sim 2.5\times
10^{-4}$; the distribution of the error for different values of the exponent
$\alpha$ is given in Figure 5.

##  5 Conclusions

This paper proposes ROMs for parametric fractional elliptic PDEs. The
recommended approach combines efficient solvers for shifted linear systems and
randomization to dramatically reduce the computational and storage costs over
a naïve version of the algorithm. Numerical results demonstrate the
computational benefits and viability of the ROMs. Future work involves
adaptive construction of the ROM to guarantee a certain user-specified error
tolerance and the development of error indicators and metrics necessary to
develop stopping criteria for the adaptive algorithms. Another line of future
inquiry involves adaptively estimating the size of the reduced basis.

##  6 Acknowledgements

AKS was supported in part by the National Science Foundation through the grant
DMS-1845406. HA was partially supported by NSF grant DMS-2110263 and the Air
Force Office of Scientific Research under Award NO: FA9550-22-1-0248. AKS
would like to acknowledge Hussam Al Daas for helpful discussions.

## References

  * [1] H. Antil and S. Bartels.  Spectral Approximation of Fractional PDEs in Image Processing and Phase Field Modeling.  Comput. Methods Appl. Math., 17(4):661–678, 2017.
  * [2] H. Antil, S. Bartels, and A. Schikorra.  Approximation of fractional harmonic maps.  IMA J. Numer. Anal., 43(3):1291–1323, 07 2022.
  * [3] H. Antil, Y. Chen, and A. Narayan.  Reduced basis methods for fractional Laplace equations via extension.  SIAM J. Sci. Comput., 41(6):A3552–A3575, 2019.
  * [4] H. Antil, Z. W. Di, and R. Khatri.  Bilevel optimization, deep learning and fractional Laplacian regularization with applications in tomography.  Inverse Problems, 36(6):064001, May 2020.
  * [5] H. Antil, D. P. Kouri, and J. Pfefferer.  Risk-Averse Control of Fractional Diffusion with Uncertain Exponent.  SIAM J. Control Optim., 59(2):1161–1187, 2021.
  * [6] H. Antil, E. Otárola, and A. Salgado.  Optimization with Respect to Order in a Fractional Diffusion Model: Analysis, Approximation and Algorithmic Aspects.  J. Sci. Comput., 77(1):204–224, 2018.
  * [7] H. Antil and A. K. Saibaba.  Efficient algorithms for Bayesian inverse problems with Whittle–Mat$\backslash$’ern priors.  arXiv preprint arXiv:2205.04417, 2022.
  * [8] T. Bakhos, P. K. Kitanidis, S. Ladenheim, A. K. Saibaba, and D. B. Szyld.  Multipreconditioned GMRES for shifted systems.  SIAM J. Sci. Comput., 39(5):S222–S247, 2017.
  * [9] P. Binev, A. Cohen, W. Dahmen, R. DeVore, G. Petrova, and P. Wojtaszczyk.  Convergence rates for greedy algorithms in reduced basis methods.  SIAM J. Math. Anal., 43(3):1457–1472, 2011.
  * [10] A. Bonito, D. Guignard, and A. R. Zhang.  Reduced basis approximations of the solutions to spectral fractional diffusion problems.  Journal of Numerical Mathematics, 28(3):147–160, 2020.
  * [11] A. Bonito and J. Pasciak.  Numerical approximation of fractional powers of elliptic operators.  Math. Comput., 84(295):2083–2110, 2015.
  * [12] L. Caffarelli and L. Silvestre.  An extension problem related to the fractional Laplacian.  Comm. Part. Diff. Eqs., 32(7-9):1245–1260, 2007.
  * [13] T. Danczul and J. Schöberl.  A reduced basis method for fractional diffusion operators II.  J. Numer. Math., 29(4):269–287, 2021.
  * [14] T. Danczul and J. Schöberl.  A reduced basis method for fractional diffusion operators I.  Numerische Mathematik, 151(2):369–404, 2022.
  * [15] H. Dinh, H. Antil, Y. Chen, E. Cherkaev, and A. Narayan.  Model reduction for fractional elliptic problems using Kato’s formula.  Math. Control Relat. Fields, 12(1):115–146, 2022.
  * [16] M. D’Elia, Q. Du, C. Glusa, M. Gunzburger, X. Tian, and Z. Zhou.  Numerical methods for nonlocal and fractional models.  Acta Numer., 29:1–124, 2020.
  * [17] H. C. Elman and V. Forstall.  Preconditioning techniques for reduced basis methods for parameterized elliptic partial differential equations.  SIAM J. Sci. Comput., 37(5):S177–S194, 2015.
  * [18] H. Fu, H. Wang, and Z. Wang.  POD/DEIM reduced-order modeling of time-fractional partial differential equations with applications in parameter identification.  J. Sci. Comput., 74(1):220–243, 2018.
  * [19] Q. Guan, M. Gunzburger, C. G. Webster, and G. Zhang.  Reduced basis methods for nonlocal diffusion problems with random input data.  Comput. Methods Appl. Mech. Engrg., 317:746–770, 2017.
  * [20] J. S. Hesthaven, G. Rozza, and B. Stamm.  Certified reduced basis methods for parametrized partial differential equations.  SpringerBriefs in Mathematics. Springer, Cham; BCAM Basque Center for Applied Mathematics, Bilbao, 2016.  BCAM SpringerBriefs.
  * [21] T. Kato.  Note on fractional powers of linear operators.  Proc. Japan Acad., 36(3):94 – 96, 1960.
  * [22] D. Kressner and C. Tobler.  Krylov subspace methods for linear systems with tensor product structure.  SIAM J. Matrix Anal. Appl., 31(4):1688–1714, 2010.
  * [23] A. Pinkus.  $n$-widths in approximation theory, volume 7 of Ergebnisse der Mathematik und ihrer Grenzgebiete (3) [Results in Mathematics and Related Areas (3)].  Springer-Verlag, Berlin, 1985.
  * [24] A. Quarteroni, A. Manzoni, and F. Negri.  Reduced basis methods for partial differential equations, volume 92 of Unitext.  Springer, Cham, 2016.  An introduction, La Matematica per il 3+2.
  * [25] A. K. Saibaba.  Randomized discrete empirical interpolation method for nonlinear model reduction.  SIAM J. Sci. Comput., 42(3):A1582–A1608, 2020.
  * [26] V. Simoncini.  Computational methods for linear matrix equations.  SIAM Rev., 58(3):377–441, 2016.
  * [27] J. Sprekels and E. Valdinoci.  A new type of identification problems: optimizing the fractional order in a nonlocal evolution equation.  SIAM J. Control Optim., 55(1):70–93, 2017.
  * [28] P. Stinga and J. Torrea.  Extension problem and Harnack’s inequality for some fractional operators.  Comm. Part. Diff. Eqs., 35(11):2092–2122, 2010.
  * [29] N. G. Trillos and D. Sanz-Alonso.  The Bayesian formulation and well-posedness of fractional elliptic inverse problems.  Inverse Probl., 33(6):065006, 2017.
  * [30] J. A. Tropp, A. Yurtsever, M. Udell, and V. Cevher.  Practical sketching algorithms for low-rank matrix approximation.  SIAM J. Matrix Anal. Appl., 38(4):1454–1485, 2017.
  * [31] C. J. Weiss, B. G. van Bloemen Waanders, and H. Antil.  Fractional operators applied to geophysical electromagnetics.  Geophys. J. Int., 220(2):1242–1259, 2020.


SNSN-323-63

$t\overline{t}$ \+ $Z$ / $W$ / $t\overline{t}$ at ATLAS

Clara Nellist, on behalf of the ATLAS Collaboration

Georg-August-Universität Göttingen
Göttingen, Germany

> The newest results from the ATLAS Collaboration for the production of a top-
> quark pair in association with a $Z$ or $W$ boson, and for the production of
> four top quarks, are summarised in these proceedings. The measurements were
> performed with 36.1 fb-1 of proton-proton collision data from the Large
> Hadron Collider at a centre-of-mass energy of 13 TeV.

> PRESENTED AT
>
>
>
>
> $11^{\mathrm{th}}$ International Workshop on Top Quark Physics
> Bad Neuenahr, Germany, September 16–21, 2018

©2019 CERN for the benefit of the ATLAS Collaboration.
Reproduction of this article or parts of it is allowed as specified in the CC-
BY-4.0 license.

##  1 Introduction

The study of the production of a pair of top quarks in association with a $W$
or $Z$ boson is referred to as $t\overline{t}{V}$. The $t\overline{t}{Z}$
process provides a direct probe of the weak couplings of the top quark. In the
Standard Model of particle physics (SM), the $t\overline{t}{W}$ process is
only possible through initial state radiation from the incoming light quarks,
and therefore an enhancement of the $t\overline{t}{W}$ process could be an
indication of Beyond the Standard Model (BSM) physics. It is also an important
background for BSM physics searches with various final states, such as stop
searches and electroweak Super Symmetry productions. Deviations from the SM
can be probed with Effective Field Theory. The $t\overline{t}{V}$ process is
also an important background for $t\overline{t}{H}$ multi-lepton searches and
four-top production.

The measurement of the SM four-top cross section is now becoming possible with
current LHC statistics. Possible enhancements of SM cross sections from new
physics through the production of heavy objects in association with a top-
quark pair are also possible. The theoretical SM cross section for four-top
production at a centre-of-mass (CoM) of 13 TeV is 9.2 fb at NLO precision [1],
which means that $10^{5}$ top-quark pairs are produced at the LHC at
$\sqrt{s}=13$ TeV for every four-top event.

The ATLAS experiment [2] at the Large Hadron Collider at CERN, is a multi-
purpose particle detector with a forward-backward symmetric cylindrical
geometry and a near $4\pi$ coverage in solid angle.

##  2 $t\overline{t}{V}$

Previous measurements of the $t\overline{t}{V}$ cross section with the full
ATLAS 8 TeV dataset [3] measured a significance of 5.0$\sigma$ (4.2$\sigma$)
over the background-only hypothesis for $t\overline{t}{W}$
($t\overline{t}{Z}$) production. There were a number of limiting factors for
this measurement, including its statistical uncertainty. For the
$t\overline{t}{W}$ fit, the dominant systematic uncertainty was the modelling
of fake leptons and background processes with misidentified charge. For the
$t\overline{t}{Z}$ fit, the dominant systematic uncertainty source was the
modelling of backgrounds from simulation.

Performing the analysis at 13 TeV, instead of 8 TeV, was more advantageous for
$t\overline{t}{Z}$ than for $t\overline{t}{W}$ due to the rapidly increasing
$t\overline{t}$ background for $t\overline{t}{W}$ with respect to the more
moderate increase in signal cross section.

###  2.1 $t\overline{t}{V}$ at 13 TeV

A simultaneous measurement of the $t\overline{t}{Z}$ and $t\overline{t}{W}$
production cross sections in final states with two, three or four isolated
electrons or muons, using 36.1 fb-1 of data from 2015 and 2016 was performed
[4]. Each channel is further divided into multiple regions to maximise the
sensitivity of the measurement. The signal regions for the channel with three
charged leptons in the final state can be seen as an example in Figure 1
(left), and for regions targeting the $t\overline{t}{W}$ process in Figure 1
(right).

Figure 1: Event yields in data compared with the fit results in, left, the
trilepton signal regions targeting the $t\overline{t}{Z}$ process and, right,
regions targeting the $t\overline{t}{W}$ process [4].

Background events containing prompt leptons were estimated from the Monte
Carlo (MC) simulation. Normalisation corrections were obtained from the
control regions included in the log likelihood fit. A data-driven approach was
taken for the $t\overline{t}$ background in the channel with two charged
leptons of opposite sign in the final state (2LOS). The validation region was
selected by inverting the requirement that the charged leptons have the same
flavour. Contributions from charge-flip events, significant for events with
electrons in the final state, were estimated from data. This background in
events with two muons in the final state was deemed to be negligible since the
probability of misidentifying the charge of a muon in the transverse momentum
range is very small. Backgrounds with greater than or equal to one fake
charged lepton were modelled using data in dedicated control regions.

Figure 2: The result of the simultaneous fit to $\sigma_{t\overline{t}{Z}}$
and $\sigma_{t\overline{t}{W}}$ with the 68% and 95% confidence level contours
[4].

In the $t\overline{t}{W}$ analysis and the $t\overline{t}{Z}$ channel with
three charged leptons in the final state (3L), the fake lepton background was
estimated using the matrix method. For the region with four charged leptons in
the final state, the semi data-driven fake factor method was used. Correction
factors accounting for potential data / MC simulation discrepancies were
extracted in a dedicated control region and are enriched with processes that
contain greater than one fake electron or muon.

####  2.1.1 $t\overline{t}{V}$ Simultaneous fit

The result of the combined fit, shown in Figure 2, measured
$\sigma_{t\overline{t}{Z}}$ = 0.95 $\pm$ 0.13 pb and
$\sigma_{t\overline{t}{W}}$ = 0.87 $\pm$ 0.19 pb. This yielded a significance
of $t\overline{t}{W}$: 4.3$\sigma$ (3.4$\sigma$) observed (expected). The
$t\overline{t}{Z}$ result has a significance of greater than 5$\sigma$.

A summary of the systematic and statistical uncertainties for the combined fit
can be found in Table 1.

Uncertainty | $\sigma_{{t\overline{t}}Z}$ | $\sigma_{{t\overline{t}}W}$
---|---|---
Total systematic | 10% | 16%
Statistical | 8.4% | 15%
Total | 13% | 22%
Table 1: Summary of relative systematic and statistical uncertainties in the
measured $\sigma_{t\overline{t}{Z}}$ and $\sigma_{t\overline{t}{Z}}$ processes
from the fit [4].

####  2.1.2 Constraints on BSM

Interpretations of the inclusive cross section measurement in terms of
Effective Field Theory (EFT) were performed. Constraints were set on the five
operators which modify the $t\overline{t}{Z}$ vertex. The first two enter the
$t\overline{t}{Z}$ vertex as a linear combination, as the measurement is
sensitive to the difference. However, only one operator was considered at a
time.

##  3 Four-tops

The four top signal is very busy, since there are many jets in the final
signature, including those that originated from a $b$-quark.

The previous ATLAS measurement was performed with 3.2 fb-1 of data taken at
$\sqrt{s}=13$ TeV [5]. Limits on the observed (expected) upper limit of 21
(16) times the four-top standard model cross section at 95% confidence level
were obtained. This measurement was performed in the single lepton channel
only.

The measurements presented here are with 36.1 fb-1 data taken in 2015 and 2016
with $\sqrt{s}=13$ TeV. The analysis was split into different channels
depending on the decay of the $W$ boson (either hadronically or leptonically)
from each top quark and combined afterwards. Broadly the separation is to one
charged lepton or two charged leptons of opposite sign (1L & 2LOS,
respectively) [6] and two charged leptons with the same sign or three charged
leptons (2LSS & 3L, respectively) [7]. Again, charged leptons are only
referring to electrons or muons.

###  3.1 1L and 2LOS

Events were categorised according to jet, $b$-tagged jet and mass-tagged
reclustered large-$R$ jet multiplicities. Both channels were dominated by the
$t\overline{t}b\overline{b}$ background, which is not expected to be modelled
well, therefore a data-driven approach was used for its estimation.

For the 1L channel, the background from events with a fake or non-prompt
lepton was estimated directly from data using the matrix method. The method
called $t\overline{t}$ Tag Rate Function (ttTRF) was used to estimate
$t\overline{t}$ \+ jets and gives the probability that a jet is $b$-tagged.

In the 2LOS channel, the background was estimated from the MC simulation.

A simultaneous fit was performed in the twenty signal regions and resulted in
an observed (expected) $\sigma$($t\overline{t}$ $t\overline{t}$) of 47 fb (33
fb) with 95% CL upper limit. Consequently, the upper limit on
$\sigma_{t\overline{t}t\overline{t}}$ was measured to be 5.1 (3.6) times the
SM prediction.

For four-top production via an EFT model with a four-top-quark contact
interaction, an observed (expected) upper limit on $\sigma$ of 21 fb (22 fb)
at 95% CL was measured. In this case the SM four-top production was considered
as a background. This result is lower than SM as the contact interaction tends
to result in final state objects with larger momenta.

###  3.2 2LSS and 3L

These channels were dominated by lepton fakes and $t\overline{t}{V}$
background. Data driven techniques were used for fake and non-prompt lepton
backgrounds, which was estimated with the matrix method, and charge
misidentification in the 2LSS channel. Irreducible backgrounds (e.g.
$t\overline{t}{V}$) were modelled by the MC simulation. This measurement used
a cut and count analysis in eight signal regions and six validation regions.
The SM $\sigma_{t\overline{t}t\overline{t}}$ upper limit was 69 fb (29 fb) for
the observed (expected).

###  3.3 Combination

For the combination, all experimental systematic uncertainties were treated as
fully correlated between channels. Background modelling uncertainties were
treated as uncorrelated. An excess of events over the SM background
prediction, excluding the SM $\sigma_{t\overline{t}t\overline{t}}$ production,
was observed (expected) at 2.8 $\sigma$ (1.0 $\sigma$). Observed (expected) at
a 95% CL upper limit on $\sigma_{t\overline{t}t\overline{t}}$ of 49 fb (19 b).
Corresponding to an upper limit on $\sigma_{t\overline{t}t\overline{t}}$ of
5.3 (2.1) times the SM prediction.

Figure 3: Left: The 95% CL upper limits on
$\sigma_{t\overline{t}t\overline{t}}$ for the individual measurements and the
combined result, relative to the SM prediction [6]. Right: The signal-strength
measurements, $\mu$ in the individual channels and for the combination [6].

##  4 Summary

Measurement with 36 fb-1 13 TeV data. Combined fit measured
$\sigma_{t\overline{t}{Z}}$ = 0.95 $\pm$ 0.13 pb $\sigma_{t\overline{t}{W}}$ =
0.87 $\pm$ 0.19 pb and yields significance of 4.3$\sigma$ (4.3$\sigma$)
observed (expected) for $t\overline{t}{W}$. The significance of
$t\overline{t}{Z}$ was measured to be greater than 5$\sigma$. No significant
excess above the SM expectation was observed. Limits were set on EFT
coefficients, which matched or exceeded previous measurements.

The four-tops combined result with all channels using 36 fb-1 of 13 TeV data
was presented. An observed (expected) upper limit of 5.3 (2.1) times the four-
top SM cross section at 95% confidence limit was measured. A slight excess
above the SM expectation is observed (expected) at 2.8 $\sigma$ (1.0
$\sigma$), which is driven by the 2LSS / 3L channel. The kinematic properties
of the 2LSS / 3L channel events were compared with expectations from BSM four-
tops and found to agree poorly. The analysis is still statistically dominated;
consequently more data is expected to improve the result.

## References

  * [1] G. Bevilacqua and M. Worek,“Constraining BSM physics at the LHC: four top final states with NLO accuracy in perturbative QCD”, JHEP 07 (2012) 111
  * [2] ATLAS Collaboration, “The ATLAS Experiment at the CERN Large Hadron Collider”, JINST 3 (2008) S08003
  * [3] ATLAS Collaboration, “Measurement of the $t\overline{t}{W}$ and $t\overline{t}{Z}$ production cross sections in pp collisions at $\sqrt{s}=8~{}$TeV with the ATLAS detector”, JHEP 11 (2015) 172
  * [4] ATLAS Collaboration, “Measurement of the $t\overline{t}{Z}$ and $t\overline{t}{W}$ cross sections in proton-proton collisions at $\sqrt{s}=13$ TeV with the ATLAS detector”, arXiv:1901.03584
  * [5] ATLAS Collaboration, “Search for four-top-quark production in final states with one charged lepton and multiple jets using 3.2 fb-1 of proton-proton collisions at $\sqrt{s}=13$ TeV with the ATLAS detector at the LHC”, ATLAS-CONF-2016-020, https://cds.cern.ch/record/2144537
  * [6] ATLAS Collaboration, “Search for four-top-quark production in the single-lepton and opposite-sign dilepton final states in pp collisions at $\sqrt{s}=13$ TeV with the ATLAS detector”, arXiv:1811.02305
  * [7] ATLAS Collaboration, “Search for new phenomena in events with same-charge leptons and b-jets in pp collisions at $\sqrt{s}=13$ TeV with the ATLAS detector”, JHEP 1812 (2018) 039