# An Autoencoder Compression Approach for Accelerating Large-scale Inverse
Problems

Jonathan Wittmer
Oden Institute for Computational Science and Engineering
University of Texas at Austin
Austin, TX 78712
<EMAIL_ADDRESS>
& Jacob Badger
Oden Institute for Computational Science and Engineering
University of Texas at Austin
Austin, TX 78712
<EMAIL_ADDRESS>
& Hari Sundar
Department of Computer Science
University of Utah
Salt Lake City, Utah 84112
<EMAIL_ADDRESS>& Tan Bui-Thanh
Department of Aerospace Engineering and Engineering Mechanics
Oden Institute for Computational Science and Engineering
University of Texas at Austin
Austin, TX 78712
<EMAIL_ADDRESS>
###### Abstract

PDE-constrained inverse problems are some of the most challenging and
computationally demanding problems in computational science today. Fine meshes
that are required to accurately compute the PDE solution introduce an enormous
number of parameters and require large scale computing resources such as more
processors and more memory to solve such systems in a reasonable time. For
inverse problems constrained by time dependent PDEs, the adjoint method that
is often employed to efficiently compute gradients and higher order
derivatives requires solving a time-reversed, so-called adjoint PDE that
depends on the forward PDE solution at each timestep. This necessitates the
storage of a high dimensional forward solution vector at every timestep. Such
a procedure quickly exhausts the available memory resources. Several
approaches that trade additional computation for reduced memory footprint have
been proposed to mitigate the memory bottleneck, including checkpointing and
compression strategies. In this work, we propose a close-to-ideal scalable
compression approach using autoencoders to eliminate the need for
checkpointing and substantial memory storage, thereby reducing both the time-
to-solution and memory requirements. We compare our approach with
checkpointing and an off-the-shelf compression approach on an earth-scale ill-
posed seismic inverse problem. The results verify the expected close-to-ideal
speedup for both the gradient and Hessian-vector product using the proposed
autoencoder compression approach. To highlight the usefulness of the proposed
approach, we combine the autoencoder compression with the data-informed active
subspace (DIAS) prior to show how the DIAS method can be affordably extended
to large scale problems without the need of checkpointing and large memory.

_Keywords_ Machine Learning $\displaystyle\cdot$ Autoencoder
$\displaystyle\cdot$ Inverse Problems $\displaystyle\cdot$ High Performance
Computing

##  1 Introduction

Modern supercomputers are an essential tool for today’s scientists and
engineers, enabling them to simulate larger and more complex systems than ever
before. The availability of large scale computing resources has enabled
numerous breakthroughs in scientific knowledge in areas such as quantum
computing [Liu et al., 2021c, Doi et al., 2019, Mandrà et al., 2021],
computational chemistry [De Jong et al., 2010, Kowalski et al., 2021], drug
discovery [Ge et al., 2013, Schmidt and Hildebrandt, 2017, Bharadwaj et al.,
2021, Sukumar et al., 2021], biology [McFarlane and Biktasheva, 2008, Bukowski
et al., 2010], and plasma physics [Bhattacharjee and Wells, 2021, Fedeli et
al., 2022]. While the increased number of floating-point operations per second
(FLOPS) has enabled larger and more challenging problems to be solved in
reasonable amounts of time, applications are increasingly becoming
bottlenecked by limited memory of computing systems — especially those
applications dealing with big data [Denis et al., 2022, Imani et al., 2019].

There are at least three common approaches to mitigating this memory
bottleneck: develop new algorithms that require less storage (for example,
using density functional theory to approximate electron interactions rather
than directly solving a many body problem in quantum mechanics [Bickelhaupt
and Baerends, 2000]), utilize additional storage devices such as hard drives
to increase the amount of data that can be stored, and trade extra computation
for reduced storage requirements. Algorithmic improvement is always desirable,
but is not always possible. Additionally, solving a closely related problem
that has computational or storage advantages may cause a loss in accuracy.
Serializing data to disk was previously only a viable solution for problems
with high compute intensity. The high cost of reading from and writing to disk
can be hidden when sufficient computation is performed. However, recent
advances in storage technology such as non-volatile memory (NVM) have made
serialization approaches more viable [Peng et al., 2020]. Indeed, NVM powers
the large-memory nodes on Frontera, the Texas Advanced Computing Center’s
largest supercomputer at the time of writing. This allows for traditional DRAM
to be used as an extra level of cache [Wu et al., 2017, Stanzione et al.,
2020]. The last mentioned approach, trading computation for storage, includes
checkpointing methods [Wang et al., 2009, Zhang and Constantinescu, 2023,
Griewank and Walther, 1997] and sparse decomposition [Zhao et al., 2020]. More
generally, this means storing the inputs required to regenerate the desired
output rather than storing the output directly [Akturk and Karpuzcu, 2018]. In
addition to saving time, it is shown in [Akturk and Karpuzcu, 2018] that
recomputation can have energy saving benefits over storing and retrieving
data.

One type of problem that calls for advanced big data management techniques is
full-waveform seismic inversion (FWI). As motivation, we give a high-level
description of why advanced data management techniques are required for FWI
here and give a detailed description in Section 2. FWI is a partial
differential equation (PDE) constrained inverse problem with dependencies on
both space and time derivatives. The dependency on time leads to the
generation of large amounts of data. With spatio-temporal measurements of the
velocity field, the inverse problem is to infer the underlying spatially
varying material properties, namely the acoustic wave speed. Typical inverse
solution methods reframe the problem as a constrained optimization problem and
use gradient-based methods to minimize an objective function. Employing the
adjoint approach [Giles and Pierce, 2000, Fichtner et al., 2006] to compute
the gradient necessitates solving an adjoint PDE that is solved backward in
time and depends on the corresponding forward solution at each timestep. As a
result, the entire solution history of the forward problem is needed in
reverse order to solve the adjoint equation. For large scale problems with
high spatial resolution and many timesteps, storing the entire forward
solution history is infeasible, requiring a different approach. As a quick
back-of-the-envelope calculation, consider a domain with 1 billion degrees of
freedom and 2 000 timesteps. Using a 4th order Runge-Kutta scheme to integrate
in time, there are 4 solutions to be stored per timestep. Each solution has 6
fields — 3 velocity fields and 3 strain fields (assuming a velocity-strain
formulation). This would require 384 TB to store in double precision. Perhaps
one of the most commonly used techniques for such a problem is to employ
checkpointing. Checkpointing is a technique whereby solutions are stored at
only a select few timesteps and the rest discarded. The solution at any
timestep can then be recreated by solving the forward problem from the nearest
previous checkpoint. While reducing memory requirements, checkpointing
effectively increases the computation requirements by 1 PDE solve per gradient
evaluation.

An alternative approach employs compression to reduce both the memory
requirements compared to full storage and computation requirements compared to
the checkpointing case. The idea to use compression to mitigate the effects of
limited memory or storage is not new. Compression techniques have been used
extensively in the audio and video processing community for decades [Blesser,
1969, Lewis and Knowles, 1990]. There are two kinds of compression: lossless
compression and lossy compression. Lossless compression is a class of
techniques that are able to exactly recover the original input data. The PNG
image format is an example of a format that employs lossless compression [Kaur
and Choudhary, 2016]. The checkpointing strategy can be viewed as a form of
lossless compression where information is encoded in the checkpoints and the
physics model provides a decompression algorithm. Lossy compression, on the
other hand, allows for errors to be made in reconstructing the input data. The
JPEG protocol is an example of lossy compression in the image processing
domain [Usevitch, 2001].

Various compression techniques have also been proposed to aid in scientific
computation. Specific to the seismic inversion community, there are two main
avenues of compression: compressing the source/receiver data [Habashy et al.,
2011] and compressing the forward solution to avoid checkpointing.
Source/receiver compression is useful in the case where large amounts of data
are available, such as when many shots are recorded. The gradient must be
evaluated at each shot; each requires the solution of the forward and adjoint
PDEs, implying that the cost scales linearly in the number of shots [Duarte et
al., 2020]. The goal of source/receiver compression is to find linear
combinations of simultaneous sources and receivers that are (nearly)
equivalent to the original problem. When the number of linear combinations
required to closely replicate the original problem is small, the number of
forward and adjoint solves required to compute a single gradient is much
smaller, reducing both the memory footprint and the computation time [Habashy
et al., 2011].

Several other works have already explored the idea of using compression to
mitigate high adjoint-induced memory requirements [Cyr et al., 2015, Kukreja
et al., 2019, Kukreja et al., 2022, Boehm et al., 2016]. Those works
demonstrate that data can be compressed without negatively affecting the
inverse solution and that compression can be faster than forward simulation.
Several techniques of compressing in space as well as a spline interpolation
technique for compression in time were proposed in [Boehm et al., 2016].
Results for a variety of off-the-shelf compression algorithms, including
extensive results using the ZFP compression algorithm [Lindstrom, 2014], were
given in [Kukreja et al., 2019]. Each of these papers demonstrates how
different compression techniques can be used to solve seismic inverse problems
on simple domains with large amounts of data. That is, the inverse problem is
not ill-posed and regularization techniques are not required. While
perturbations to the gradient due to compression may minimally affect the
inverse solution in the case of large data, it is unclear that such approaches
will be viable in solving ill-posed problems with limited data. Additionally,
these approaches have limited compression capabilities, limiting the maximum
problem size that can be solved.

Input layer Latent space Output layer EncoderDecoder Figure 1: Autoencoder
architecture with deep neural networks

Meanwhile, the scientific and technology communities have seen an explosion in
the applications of machine learning techniques. Deep learning in particular
has found fruitful application in enabling progress in speech recognition
[Kamath et al., 2019], computer vision [Bjerge et al., 2022, Esteva et al.,
2021], and in accelerating innovation in scientific computation [Reichstein et
al., 2019, Kates-Harbeck et al., 2019]. One approach that has emerged is to
use an autoencoder architecture to perform data compression [Cheng et al.,
2018, Liu et al., 2021b, Hinton and Salakhutdinov, 2006]. An autoencoder is
simply a deep neural network that approximates the identity map. The output at
some intermediate layer can be extracted as a representation of the input and
is referred to as the latent representation which resides in the latent space
[Dillon et al., 2021]. Autoencoders are typically decomposed into two separate
networks: an encoder which consists of the layers from the input to the latent
space and a decoder which is a network that maps the latent space back to the
original dimension. To achieve a compressed representation, the latent space
must have dimension strictly smaller than the input. Compared to traditional
dimension reduction methods, autoencoders can be viewed as a non-linear
version of principal component analysis (PCA) [Ladjal et al., 2019, Kneer et
al., 2021]. In the most simple case of a single layer encoder and single layer
decoder with linear activation, the autoencoder learns the same compression as
PCA. That is, the single layer linear autoencoder projects the input data onto
the same lower dimensional subspace that PCA identifies, though the
representation is different [Baldi and Hornik, 1989, Ladjal et al., 2019]. The
actual PCA basis vectors can be obtained from the weights of the trained
encoder, showing the two are equivalent [Plaut, 2018]. The advantage of a deep
autoencoder with nonlinear activations is that it is often able to achieve
lower reconstruction error for the same dimension latent space compared to PCA
[Liu et al., 2021a, Liu et al., 2021b], enabling higher feasible compression
ratios in application.

While autoencoders have empirically demonstrated impressive compression ratios
and low error rates, there is no theory that we are aware of at this time
bounding the error of an autoencoder. However, a practical method of bounding
the error of an autoencoder is to compute and store some form of the residual
between the desired output ($\displaystyle\boldsymbol{y}_{\text{true}}$) and
the actual output ($\displaystyle\boldsymbol{y}_{\text{pred}}$), i.e.
$\displaystyle\boldsymbol{r}=\boldsymbol{y}_{\text{true}}-\boldsymbol{y}_{\text{pred}}$.
Various techniques of approximately representing this residual vector have
been proposed such as compressing $\displaystyle\boldsymbol{r}$ with an off-
the-shelf algorithm when
$\displaystyle\left\|\boldsymbol{r}\right\|>\text{tol}$ for some tolerance,
tol [Lee et al., 2022, Liu et al., 2021a] or storing only the elements
$\displaystyle r_{i}$ of $\displaystyle\boldsymbol{r}$ where
$\displaystyle|r_{i}|>\text{tol}$ [Liu et al., 2021b, Abu Alsheikh et al.,
2014].

To address the problems of additional expensive PDE solves required by the
checkpointing approach and the limited compression capability of off-the-shelf
compression techniques, we propose a novel autoencoder approach for
compressing the forward solution. We show empirically that this compression
method is faster than the checkpointing method while achieving comparable
solution accuracy. In addition to showing that a trained autoencoder has
memory and computational advantages over existing methods, we develop an
efficient data generation and training procedure based on the Bayesian
inversion formulation that enables our approach to maintain its advantages,
even when including the cost of training the autoencoder. While variations of
autoencoders are often used for generative modeling [Kingma and Welling, 2013]
and others assign a physical interpretation to the latent space [Goh et al.,
2019], we use the auto-encoding capabilities of the autoencoder architecture
to allow the training process to find the most effective compressed latent
representation for seismic input data. We provide a mathematical description
of the FWI problem and the inverse problem in Section 2. Two autoencoder
compression variants are proposed in Section 3 along with detailed explanation
of autoencoder architecture, training data selection, and normalization. We
show in Section 4 that machine learning can be leveraged to accelerate the
solution of seismic inverse problems in very high dimensions, even in the ill-
posed case with few measurements. Through numerical experiment, we show that
this approach is faster than the checkpointing approach, and as fast as the
state-of-the-art off the shelf compression approach, ZFP [Lindstrom, 2014],
while achieving much higher compression ratios. We then extend the data
informed active subspace (DIAS) regularization method [Nguyen et al., 2022] to
nonlinear problems and show that autoencoder compression can be used to enable
the affordable application of the DIAS prior for large-scale inverse problems.

##  2 Full Waveform Inversion

As a simplified model of the earth’s seismic properties, we consider the
acoustic wave equation. In velocity-strain form, this gives rise to a system
of first order PDEs. We choose the velocity-strain form because of its close
relation to the acoustic-elastic wave equation [Wilcox et al., 2015]. This
enables simultaneous derivation of gradients and Hessians for both acoustic,
elastic, and acoustic-elastic formulations. Additional computational and
mathematical advantages can be found in [Wilcox et al., 2015]. Following the
setup given in [Bui-Thanh et al., 2013], the acoustic wave equation is given
by,

$\displaystyle\displaystyle\rho\frac{\partial\boldsymbol{v}}{\partial
t}-\nabla\left(\rho c^{2}e\right)=\boldsymbol{g},$ (1a)
$\displaystyle\displaystyle\frac{\partial e}{\partial
t}-\nabla\cdot\boldsymbol{v}=0,$ (1b)

where $\displaystyle\rho=\rho(\boldsymbol{x})$ is the density, $\displaystyle
c=c(\boldsymbol{x})$ is the acoustic wave speed for which we are inverting,
$\displaystyle\boldsymbol{g}=\boldsymbol{g}(\boldsymbol{x},t)$ is a forcing
function, $\displaystyle\boldsymbol{v}=\boldsymbol{v}(\boldsymbol{x},t)$ is
the velocity vector, $\displaystyle e=e(\boldsymbol{x},t)$ is the trace of the
strain tensor (dilatation), and $\displaystyle\boldsymbol{x}$ denotes the
spatial coordinate in the domain $\displaystyle\Omega$. We specify initial
conditions for the velocity and dilatation fields

$\boldsymbol{v}(\boldsymbol{x},0)=\boldsymbol{v}_{0}(\boldsymbol{x})\quad\text{and}\quad
e(\boldsymbol{x},0)=e_{0}(\boldsymbol{x}),\quad\boldsymbol{x}\in\Omega.$ (2)

Lastly, we impose traction-free boundary conditions

$e(\boldsymbol{x},t)=0,\quad\boldsymbol{x}\in\Gamma=\partial\Omega,\;t\in(0,T)$

with final time $\displaystyle T$. The 3 components of the velocity vector
field and the 3 diagonal components of the strain tensor field, used to
compute the dilatation, will be collectively referred to as the state
variables in the rest of this paper. It is the state variables that need to be
compressed and decompressed.

In this work, we consider the inverse parameter estimation problem — inferring
the wave speed $\displaystyle c\left(\boldsymbol{x}\right)$ from sparse
measurements of the velocity, $\displaystyle\boldsymbol{d}$. Since we consider
the setting with limited measurements, the inverse problem is ill-posed and
regularization is required. The inverse problem can be formulated as the
optimization problem

$\boldsymbol{u}^{*}:=\underset{\boldsymbol{u}}{\text{arg min
}}\frac{1}{2}\left\|\boldsymbol{\mathcal{F}}\left(\boldsymbol{u}\right)-\boldsymbol{d}\right\|_{{\Gamma}^{-1}_{\text{noise}}}^{2}+\frac{1}{2}\left\|\boldsymbol{u}-\boldsymbol{u}_{0}\right\|_{{\Gamma}^{-1}_{\text{prior}}}^{2}$
(3)

where $\displaystyle\boldsymbol{u}$ is referred to as the parameter of
interest (PoI), $\displaystyle\boldsymbol{\mathcal{F}}$ is the parameter to
observable map (PtO map), $\displaystyle\Gamma_{\text{noise}}$ is the noise
covariance matrix of the observations $\displaystyle\boldsymbol{d}$,
$\displaystyle\boldsymbol{u}_{0}$ is the initial guess and
$\displaystyle\Gamma_{\text{prior}}$ characterizes the regularization. In the
Bayesian setting [Kaipio and Somersalo, 2006],
$\displaystyle\boldsymbol{u}_{0}$ and $\displaystyle\Gamma_{\text{prior}}$ are
interpreted as the prior mean and prior covariance matrix, respectively. While
we consider full waveform inversion in the deterministic setting, it is
convenient to consider the problem in the Bayesian framework and only work
toward estimating the maximum a posterior (MAP) point. The statistical framing
of the MAP problem allows us to choose regularization via selection of a prior
distribution rather than ad-hoc and allows for trivial extension of our
proposed compression approach to quantifying uncertainty. In our case, the PoI
is the acoustic wave speed which is mapped to our observations,
$\displaystyle\boldsymbol{d}$, by solving the acoustic wave equation and
observing the velocity field only at the receivers. Let
$\displaystyle\boldsymbol{u}\approx c$ be the discretized acoustic wave speed
that we will numerically estimate. Then

$\boldsymbol{\mathcal{F}}:=\boldsymbol{B}\boldsymbol{\mathcal{A}}$

where $\displaystyle\boldsymbol{\mathcal{A}}$ is the solution of (1) and
$\displaystyle\boldsymbol{B}$ is an observation operator, extracting the
velocity field at receiver locations.

As in [Bui-Thanh et al., 2013], we choose the prior to be the PDE-based
BiLaplacian which encodes smoothness and anisotropy of the wave speed along
with our certainty about our initial guess. The BiLaplacian prior term
$\displaystyle\frac{1}{2}\left\|\boldsymbol{u}-\boldsymbol{u}_{0}\right\|_{{\Gamma}^{-1}_{\text{prior}}}^{2}$
is computed by solving the following elliptic PDE:

$\displaystyle\displaystyle-\alpha\nabla\cdot\left(\Theta\nabla\boldsymbol{u}\right)+\alpha\boldsymbol{u}=s\quad$
$\displaystyle\displaystyle\text{in }\Omega$ (4a)
$\displaystyle\displaystyle\alpha\left(\Theta\nabla\boldsymbol{u}\right)\cdot\boldsymbol{n}=0\quad$
$\displaystyle\displaystyle\text{on }\partial\Omega$ (4b)

where $\displaystyle\boldsymbol{n}$ is the outward facing unit normal on the
boundary, $\displaystyle\partial\Omega$.
$\displaystyle{\Gamma}^{-1}_{\text{noise}}$ is the square of this differential
operator. Properties of this prior and further discussion can be found in
[Bui-Thanh et al., 2013]. Finally, derivation of the discrete forward,
adjoint, and gradient expressions for a discontinuous Galerkin discretization
can be found in [Wilcox et al., 2015].

##  3 Compression methods

There are a variety of ways that the state variables can be compressed. Each
state variable is a function of time and 3 space dimensions, resulting in 4
dimensions across which the state can be compressed. Additionally, there might
be some relationship between the state variables that could be exploited to
achieve higher compression ratios. In this work, we propose two compression
techniques that are amenable for implementation with autoencoders, treating
each state variable independently:

  1. 1.

Compression in space, allowing the autoencoder to discover the 3D structure of
the state across multiple elements.

  2. 2.

Compression in time on an element-by-element basis.

As a brief justification for pursuing only these two compression setups, let
us make a few notes. First, there is potential for scale mismatch between each
of the state variables. Consider again equation (1b). The time derivative of
the dilatation is related to the divergence of the velocity field. From this
we can see that there could be very large velocities, but no dilatation if the
velocity field is divergence free. This makes it difficult to design a method
that can garner any correlation between the value of the velocity field at a
given timestep and the value of the strains.

Secondly, we decided against pursuing methods requiring explicit specification
of the fields in three-dimensional coordinates. That is, we treat all the
nodal coefficients equally, not explicitly taking advantage of the embedding
of these nodal coefficients in 3-dimensional space. While it is in principle
possible to use something like a convolutional neural network with
3-dimensional kernels [Ji et al., 2012], we found two major difficulties: 1)
convolutional neural networks in 3D are rather slow and 2) the nodes of our DG
mesh are not equally spaced. Although equally spaced nodes would enable simple
application of convolutional autoencoders, non-uniformly spaced nodal
coordinates such as Gauss-Lobatto-Legendre (GLL) enable more accurate
integration [Quarteroni et al., 2010] and are a natural choice for node
location in a DG scheme. In Sections 3.1 and 3.2, we will present the
implementations of the two proposed autoencoder compression approaches.

###  3.1 Compression on a uniform mesh

Consider a box domain where $\displaystyle\Omega=\left[0,1\right]^{3}$. Such a
domain may arise in reservoir exploration, the example used in [Kukreja et
al., 2019]. Figure 2 shows the mesh on a box domain. While the nodes are
unevenly spaced, they occur in a regular pattern. This allows us to flatten a
state variable across multiple elements into a single vector for compression.
Care must be taken to partition the domain evenly across MPI ranks and to
ensure that each rank stores the nodes in the same spatial orientation, i.e.
the elements are stacked the same way in memory on each rank. The autoencoder
then has the freedom to identify the most relevant spatial correlation in
order to most accurately compress the data. Note that the order in which the
elements are stacked does not matter so long as it is the same across all
ranks since a permutation can be implicitly learned. With the same state
layout on every rank, a single autoencoder can be trained that compresses the
local degrees of freedom on a single rank. Once trained, each rank stores its
own copy of the autoencoder, eliminating the need for any communication. The
fact that no communication is required is key to the scalability of our
proposed approach.

(a)

(b)

Figure 2:  (a) Course mesh on box domain showing non-uniform spacing of grid.
Different colors correspond to different MPI ranks. (b) mesh of single MPI
rank consisting of 8 third order elements.

There are a few competing factors that must be balanced in designing an
autoencoder for compression. The input vector should be:

  1. 1.

Large enough so that a high compression ratio can be attained with low
reconstruction error.

  2. 2.

Small enough so that the autoencoder runs quickly.

The second item poses a particular challenge in the present case as the DG
solver used to solve the forward acoustic wave equation is highly tuned.
Special care must be taken in order to develop a method that is faster than
checkpointing while maintaining sufficient accuracy. We found empirically that
an input dimension of 4096 adequately balanced these two constraints. This
corresponds to 8 elements per MPI rank with 3${}^{\text{rd}}$ order
polynomials, resulting in 64 nodes per element ($\displaystyle 4^{3}$).

###  3.2 Compression on non-uniform mesh

While using an autoencoder to compress in space is intuitive and simple to
implement, there are severe limitations to the practical applications of such
an approach. First, each MPI rank must have the exact same number of degrees
of freedom that are arranged in exactly the same manner in order to avoid
padding or communicating ghost elements. This is because autoencoders that
employ dense hidden layers have a fixed input size. Unless extreme care is
taken, this requires that the mesh be uniformly refined and have no hanging
nodes. Such a restriction poses major challenges for problems with a non-
rectangular domain or when non-hexahedral elements are desired. To mitigate
these limitations, we propose to compress along the time dimension, treating
each element as an independent state. This method works seamlessly with
h-adaptive mesh refinement (refining the mesh geometry), though p-adaptivity
(changing the polynomial order of the element) remains a challenge. For this
work, we only consider algorithms with h-adaptivity.

Figure 3:  A cross-section of the earth domain with h-adaptive mesh
refinement. Notice the non-uniform character of the mesh, with very small
elements near the crust required to resolve the seismic waves.

To demonstrate the success of this approach, we consider solving a seismic
inverse problem on a spherical domain. We employ h-adaptive mesh refinement to
ensure that the seismic waves can be resolved on each element with at least 3
nodes per wavelength. Figure 3 shows a cross section of the mesh. Hanging
nodes that result from local mesh refinement pose no problem for our proposed
compression strategy since each node is considered independently.

The same constraints as discussed in Section 3.1 apply in this scenario.
Similar to the box domain case, we found that an input dimension to the
autoencoder of 4096 worked well. Since we use a 4${}^{\text{th}}$ order Runge-
Kutta (RK) explicit timestepping scheme [Quarteroni et al., 2010], a 4096
dimensional input vector results from using $\displaystyle 3^{\text{rd}}$
order polynomials (64 nodes per element) and 16 timesteps, each with 4 RK
stages.

###  3.3 Autoencoder architecture and training details

Any discussion of results obtained using deep learning would be incomplete
without a discussion of the architecture used to obtain the results. While the
order of the data differs between space compression and time compression
methods, the same autoencoder architecture can be used in either case since
the input and output dimensions are the same. The difference between the
trained autoencoders used for each compression method lies solely in the
training data. Table 1 shows the architecture design of both the encoder and
the decoder. All layers are “dense” layers, which corresponds to an affine
transformation followed by a nonlinear activation. We found the exponential
linear unit (ELU) to perform better than the rectified linear unit (ReLU),
hyperbolic tangent (tanh), and sigmoid activations. The ELU activation
function is given by

$\displaystyle\displaystyle\text{ELU}(y):=\begin{cases}e^{y}-1&\text{if
}y<0\\\ y&\text{else}.\end{cases}$

Table 2 shows the relevant hyperparameters of the training process. Training
took approximately 3 hours using two Nvidia 2080Ti GPUs.

| Encoder | Decoder
---|---|---
Input dimension | 4096 | 64
Output dimension | 64 | 4096
# hidden layers | 7 | 7
Neurons per hidden layer | [512, 256, 256, 256, 128, 64, 64] | [128, 128, 256, 256, 256, 512, 4096]
Activation | ELU | ELU
Table 1:  Autoencoder architecture and hyperparameters used to train network. Precision | float32
---|---
Optimizer | LAMB [You et al., 2019]
Initial learning rate | $\displaystyle 10^{-3}$
Learning rate decay | $\displaystyle lr=0.5lr$
applied every 5 epochs
Epochs | 20
Batch Size | 512
# GPUs | 2$\displaystyle\times$ 2080Ti
Table 2:  Training hyperparameters used to fit autoencoder weights and biases
to data.

As mentioned in Section 3.1, it is difficult to design an autoencoder
compression system that both performs well in terms of accuracy and
compression ratio while also being fast. Recall that the evaluation of a
single dense layer requires the computation of

$\sigma\left(\boldsymbol{y}\boldsymbol{W}+\boldsymbol{b}\right).$ (5)

Considering the first layer of the encoder,
$\displaystyle\boldsymbol{W}_{1}\in\mathbb{R}^{4096\times 512}$,
$\displaystyle\boldsymbol{y}\in\mathbb{R}^{4096}$ and
$\displaystyle\boldsymbol{b}_{1}\in R^{512}$. There are 2,097,152 entries in
the matrix $\displaystyle\boldsymbol{W}_{1}$. Compare this to the second layer
which has a weight matrix
$\displaystyle\boldsymbol{W}_{2}\in\mathbb{R}^{512\times 256}$ and 131,072
entries. The first layer then has 16 times more entries in the weight matrix
than the second layer. Further layers of the encoder have even smaller weight
matrices. It is clear that the dominant cost in compressing a vector is the
computation of the first matrix-vector multiply,
$\displaystyle\boldsymbol{y}\boldsymbol{W}_{1}$. To mitigate this cost, we
develop a hybrid sparse-dense architecture whereby the first layer of the
encoder and last layer of the decoder are trained to be 95% sparse.

As an example to show the improvement of the sparse-dense architecture over
the fully dense architecture, consider compressing a batch of 864 vectors each
with size 4096. Since evaluating the output of the neural network requires
nothing more than matrix-vector multiplications, vector additions, and
element-wise application of an activation function, all 864 vectors can be
processed simultaneously by stacking into rows. Then the input
$\displaystyle\boldsymbol{y}$ has shape $\displaystyle(864,4096)$. Timing
results for both the encoder and decoder on a single core of an Intel Xeon
Platinum 8280 CPU on Frontera are shown in Table 3. We would then expect
approximately a 50% speedup with the sparse-dense architecture over the fully
dense architecture.

| Encoder | Decoder
---|---|---
Fully dense | 102 ms | 140 ms
Hybrid sparse-dense | 40 ms | 79 ms
Table 3:  Comparison of timings of fully dense architectures vs hybrid sparse-
dense architecture for encoder and decoder for compressing/decompressing 864
vectors on a single core of an Intel Xeon Platinum 8280 CPU.

Unfortunately, the application setting is slightly different since all cores
are used rather than just a single core. Consider the same test case
replicated on each CPU core simultaneously. There are 56 cores across 2
sockets on a single node of Frontera. Since each task executes independently
with no communication, we might expect perfect scaling and for the test to
complete in the same amount of time. However, there are now more cores
contending for data with limited cache and memory bandwidth. Indeed, in the
first case, 864 single precision input vectors require only 14.2 MB to store.
This easily fits in the 38.5 MB of L3 cache of a single Intel Xeon Platinum
8280 [int,]. Furthermore, the weights of the fully dense encoder require only
17.8 MB to store. The input data and the entire encoder can fit in L3 cache!
Since L3 cache is shared among all cores of this CPU, bytes must be streamed
from main memory when all cores are running simultaneously. Table 4 shows the
timings for the encoder and decoder with the fully dense and hybrid sparse-
dense architectures when running on all cores.

| Encoder | Decoder
---|---|---
Fully dense | 110 ms | 160 ms
Hybrid sparse-dense | 120 ms | 117 ms
Table 4:  Comparison of timings of fully dense architectures vs hybrid sparse-
dense architecture for encoder and decoder for compressing/decompressing 864
vectors on 56 cores of 2 Intel Xeon Platinum 8280 CPUs on a single node of the
Frontera supercomputer.

Although there is still significant advantage in using the hybrid architecture
for the decoder, the evaluation of the encoder proves to be faster using the
fully dense architecture than the hybrid architecture when all 56 cores are
being utilized. We are unsure why this is the case at the time of writing and
further optimization is a topic for future investigation. However, measuring
the performance showed that we could achieve faster compression by using a
fully dense matrix for the weights of the first layer while optimizing the
evaluation time of the decoder using a sparse matrix for the final layer.
Because of the strict speed requirements, we do not utilize any error bounding
methods which require the computation and potential storage of a residual
vector [Lee et al., 2022, Liu et al., 2021a, Liu et al., 2021b, Abu Alsheikh
et al., 2014].

####  3.3.1 Training data

One common criticism of deep learning methods is that large training datasets
are often required for the DNN to perform well in practice. In the case of
compression, other algorithms such as those implemented in the ZFP package
[Lindstrom, 2014] do not require training data. Even when deep learning
methods may perform well on training data, they may perform poorly in practice
when the training dataset is not sufficiently large, a concept referred to as
generalization error. While this is true, we argue that the cost of training
an autoencoder is a one-time upfront cost. Further, data generation and
training require several orders of magnitude fewer compute node-hours than
solving the inverse problem. In this section, we detail the data generation
process.

Consider a statistical interpretation where the training data
$\displaystyle\boldsymbol{y}\in\mathbb{R}^{n}$, is drawn from some
distribution, $\displaystyle\pi\left(\boldsymbol{y}\right)$. Suppose that we
have a dataset of $\displaystyle N$ samples
$\displaystyle\mathcal{D}^{N}=\left\\{{\boldsymbol{{y}}}^{1},...,{\boldsymbol{{y}}}^{N}\right\\}$
where
$\displaystyle{\boldsymbol{{y}}}^{i}\sim\pi\left({\boldsymbol{{y}}}\right)$.
Let $\displaystyle\Psi_{\mathcal{D}^{N}}\left({\boldsymbol{{y}}}\right)$
denote the output of the machine learning algorithm after being trained with
dataset $\displaystyle\mathcal{D}^{N}$. A common metric for evaluating the
performance of any machine learning algorithm is the generalization error.
Given some loss function denoted
$\displaystyle\mathcal{L}\left(\Psi\left({\boldsymbol{{y}}}\right),{\boldsymbol{{y}}}\right)$,
such as mean squared error,

$\mathcal{L}_{\text{MSE}}\left(\Psi\left(\boldsymbol{y}_{\text{true}}\right),\boldsymbol{y}_{\text{true}}\right):=\frac{1}{n}\left\|\boldsymbol{y}_{\text{true}}-\Psi\left(\boldsymbol{y}_{\text{true}}\right)\right\|^{2}_{2},$

the generalization error for an autoencoder is
$\displaystyle\mathbb{E}_{\pi\left({\boldsymbol{{y}}}\right)}\left[\mathcal{L}\left(\Psi\left({\boldsymbol{{y}}}\right),{\boldsymbol{{y}}}\right)\right]$
[Nadeau and Bengio, 1999]. That is, the generalization error is the average
loss over the entire distribution
$\displaystyle\pi\left({\boldsymbol{{y}}}\right)$.

Generalization error plays an important role in designing a data generation
process. There are two ways one can reduce the generalization error:

  1. 1.

Increase the size of the training dataset by drawing more samples from
$\displaystyle\pi\left({\boldsymbol{{y}}}\right)$.

  2. 2.

Choose $\displaystyle\pi\left({\boldsymbol{{y}}}\right)$ so that it can be
well-approximated with few samples.

The second method can be interpreted as choosing the most narrow distribution
that is appropriate for the application at hand. We will explore how this
option can be used in the context of generating data for seismic inversion.
For typical machine learning applications, this may not be possible as we
don’t know $\displaystyle\pi\left({\boldsymbol{{y}}}\right)$. The Bayesian
inversion setting makes this possible for our application, as we discuss next.

In order to compress seismic data, we need to generate examples similar to the
states that will be seen in application. One way to do this is to solve the
acoustic wave equation and save some of the data for training. However, it is
still unclear which solutions of the wave equation should be included in the
training dataset. Should we try to train the autoencoder to compress and
decompress all possible solutions to the wave equation on a given domain?
There are clearly wave equation solutions that are not relevant for our
application. Mathematically, the challenge becomes defining the relevant
distribution $\displaystyle\pi\left({\boldsymbol{{y}}}\right)$. The Bayesian
framework [Kaipio and Somersalo, 2006] helps us to choose a good distribution
for generating data. Bayes’ formula tells us,

$\pi_{\text{{post}}}\left(\boldsymbol{u}|\boldsymbol{d}\right)\propto\pi_{\text{{like}}}\left(\boldsymbol{d}|\boldsymbol{u}\right)\pi_{\text{{prior}}}\left(\boldsymbol{u}\right).$
(6)

This can be interpreted as the likelihood,
$\displaystyle\pi_{\text{{like}}}\left(\boldsymbol{d}|\boldsymbol{u}\right)$,
updating the prior density,
$\displaystyle\pi_{\text{{prior}}}\left(\boldsymbol{u}\right)$, based on
observed data, $\displaystyle\boldsymbol{d}$, to produce the posterior
density,
$\displaystyle\pi_{\text{{post}}}\left(\boldsymbol{u}|\boldsymbol{d}\right)$,
[Scales and Tenorio, 2001]. The observations $\displaystyle\boldsymbol{d}$ are
related to the states that will be compressed, denoted
$\displaystyle{\boldsymbol{{y}}}$ here, through an observation operator, i.e.
$\displaystyle\boldsymbol{d}=\boldsymbol{\mathcal{B}}{\boldsymbol{{y}}}$.
Since the posterior density is the product of the likelihood and the prior,
the posterior density is non-zero only where the prior density is non-zero.
That is, the posterior lies within the support of the prior. Therefore, we are
only interested in whether the autoencoder performs well in compressing
states, $\displaystyle{\boldsymbol{{y}}}$, that can be generated by samples
from the prior. We can generate relevant examples of possible states by
solving the acoustic wave equation with the PoI (the wave speed) set to
samples from the prior. The autoencoder does not need to perform well on every
possible state that could be generated from the wave equation, it only needs
to work well on states generated by solving the wave equation with parameters
sampled from the prior.

With motivation from statistical learning theory and Bayesian inverse
problems, we now detail an algorithm for generating training data.

Algorithm 1 Data generation algorithm for training compression autoencoder

Input: number of samples, $\displaystyle n$; final time, $\displaystyle T$

1: for $\displaystyle\,i=1,...,n$ do

2: $\displaystyle c=$ random draw from prior

3: for $\displaystyle\,t=1,...,T$ do

4: Solve acoustic wave equation for one timestep

5: Consolidate states, $\displaystyle\boldsymbol{v}$, $\displaystyle e$, into
compression data structure

6: if Data structure is full then

7: Write data structure to file

8: end if

9: end for

10: end for

The compression data structure is a vector that stores the state variables in
the order corresponding to the compression scheme in use. This setup allows
any ordering to be used and for training data to automatically have the
correct ordering — reducing effort to set up the training pipeline. Algorithm
1 can generate large amounts of data very quickly. For the earth-scale problem
with an 11.2 million DoF mesh, several terabytes can be generated in only a
few minutes. To reduce the amount of data written to file, we randomly chose
whether to write or discard the data structure. On 32 nodes of Frontera,
solving the wave equation and writing to file takes on the order of 1 minute.
With 10 draws from the prior, keeping only $\displaystyle 0.1\%$ of the data,
we generated 50 GB of training data in less than 1 hour. The total cost of
data generation is then $\displaystyle\mathcal{O}\left(10\right)$ node-hours.

Finally, we discuss normalization of the training data. The states have scales
that vary between $\displaystyle 10^{-16}$ and $\displaystyle 10^{5}$ over the
solution of (1). It would be difficult for an autoencoder to be able to
accurately compress and decompress states with such widely varying scales,
especially in single precision. To reduce the burden on the autoencoder of
learning scale, we normalize the data between 0 and 1. This enables us to
store the scale and offset of the data in double precision while capturing the
relative variation of the data in single precision. Consider an arbitrary
state vector that will be compressed, denoted
$\displaystyle{\boldsymbol{{y}}}$ following the notation in (5). Let
$\displaystyle{\boldsymbol{{y}}}_{\text{true}}$ be a consolidated state vector
that will be compressed (state DoFs consolidated either in space or in time).
The input vector to the encoder, $\displaystyle{\boldsymbol{{y}}}$, is then
given by

$\displaystyle\displaystyle r$
$\displaystyle\displaystyle=\max(\max({\boldsymbol{{y}}}_{\text{true}})-\min({\boldsymbol{{y}}}_{\text{true}}),\beta)$
$\displaystyle\displaystyle{\boldsymbol{{y}}}$
$\displaystyle\displaystyle=\left({\boldsymbol{{y}}}_{\text{true}}-\min({\boldsymbol{{y}}}_{\text{true}})\right)/r.$

In order to avoid division by zero, $\displaystyle r$ is set to be at least
some tolerance, $\displaystyle\beta$. We choose $\displaystyle\beta=10^{-7}$.

##  4 Numerical results

In this section we present numerical results showing the efficacy of our
proposed autoencoder compression approach for both spatial compression on the
box domain problem (Section 3.1) and temporal compression on a spherical
domain (Section 3.2). For both problems, we solve the inverse problem using a
Newton conjugate gradient (NewtonCG) method [Yang, 2009, Epanomeritakis et
al., 2008]. In each conjugate gradient iteration, the action of the Hessian on
a vector is required which entails solving two additional linearized PDEs,
referred to as the incremental forward problem and incremental adjoint
problem. Expressions for these PDEs in the case of the acoustic wave equation
can be found in [Bui-Thanh et al., 2012]. The incremental forward problem
relies on the solution of the forward problem at each timestep and so requires
decompression as well. The incremental adjoint problem requires the solution
of both the forward and incremental forward problems, and thus requires
another decompression of the forward states. Empirically, we find that the
gains of decompression are mostly offset by the costs of compression. However,
the states only need to be compressed once per Newton iteration and are
potentially decompressed many times. That is, the cost of compression is
amortized over every evaluation of the adjoint, incremental forward, and
incremental adjoint problems within a single Newton iteration. While it is
also possible to compress and decompress the incremental forward problem, it
is not explored in this paper since each incremental forward solution is only
re-used once.

###  4.1 Results: compression in space

We demonstrate the compression in space approach on a cube domain with all
sides of length 1km. Five sources are located at 1m depth and 1089 receivers
are evenly spaced at the surface, observing the velocity field for 3 seconds.
The sources are Ricker wavelets in time, smoothed by convolving with a narrow
Gaussian in space. The source term for a source at
$\displaystyle\boldsymbol{x}_{0}$ is then given by

$\boldsymbol{g}\left(\boldsymbol{x},t;\boldsymbol{x}_{0}\right)=\rho\left(\boldsymbol{x}\right)\frac{1}{\sqrt{2\pi}\sigma_{\boldsymbol{x}}}e^{-\frac{\left\|\boldsymbol{x}-\boldsymbol{x}_{0}\right\|_{2}^{2}}{2\sigma_{\boldsymbol{x}}^{2}}}\left(1-\frac{(t-t_{c})^{2}}{\sigma_{t}^{2}}\right)e^{-\frac{(t-t_{c})^{2}}{2\sigma_{t}^{2}}}\left[0,0,-1\right]^{T}.$

In this case, we consider a wavelet centered at $\displaystyle t_{c}=0.6$
seconds with $\displaystyle\sigma_{t}=\frac{1}{\pi}$ and
$\displaystyle\sigma_{\boldsymbol{x}}=0.05$.

A snapshot of the velocity magnitude is shown in Figure 4(a). Following the
setup in Section 3.1, we show both speedup and absolute error results compared
to the checkpointing solution. While the synthetic true solution is our
target, our compression strategy aims to replace checkpointing as a means for
managing limited memory, not to change the inverse solution. Therefore, the
appropriate comparison is with the inverse solution computed using checkpoints
rather than the true solution. Similar to the earth problem, we invert for a
variation in the background wavespeed field. A box-shaped inclusion is
introduced into the domain, shown in Figure 4(b), and is referred to as the
anomaly.

(a)

(b)

Figure 4:  (a) Snapshot of velocity magnitude at 1.2 seconds (b) Inclusion in
background wave speed.

To demonstrate the weak scaling of this method, we show results for several
refinement levels, each refinement increasing the total number of degrees of
freedom by 8 times. Table 5 shows the speedup compared to the checkpointing
solution. While the overall speedup may not be very impressive, it is not the
only objective of our proposed compression approach. As we have discussed, the
other objectives are to reduce memory footprint and to eradicate
checkpointing. Furthermore, there is a maximum speedup of
1.3$\displaystyle\times$ for the gradient since the time to solve the forward
problem, the adjoint problem, and evaluate the gradient are approximately
equal. That is, the cost of computing the gradient in the ideal case is 3 PDE
solves while the checkpointing case requires 4 PDE solves. Thus, the speedup
of 1.22 for the largest mesh in Table 5 is close to ideal. Similarly, the
Hessian-vector product computation has the maximum speedup of 1.5 (requiring 4
PDE solves in the ideal case vs 6 PDE solves with checkpointing), and the
speedup of 1.27, without compressing the incremental forward solve, for our
compression approach is already an achievement. Furthermore, the acoustic wave
equation, (1), is linear and not difficult to solve with explicit time
stepping. Our proposed approach would show closer to ideal speedup when used
to replace implicit solvers or when used on nonlinear problems with much more
time consuming forward solves since the cost of autoencoder compression only
depends on the number of parameters of the neural network. Table 5 shows that
the speedup approaches the ideal one and the relative error decreases as the
number of DoFs increases. A cutaway of the inverse solution with the largest
mesh for both the checkpointing case and the autoencoder compression case is
shown in Figure 5. The fact that the anomaly reconstructions are visibly
indistinguishable is due to 0.4% difference.

Mesh DoFs | Cores |  $\displaystyle\nabla$ Speedup |  $\displaystyle\mathcal{H}\boldsymbol{v}$ Speedup | Relative $\displaystyle l_{2}$ Error %
---|---|---|---|---
4,096 | 1 | 1.06 | 1.20 | 1.7
32,768 | 8 | 1.19 | 1.26 | 0.7
262,144 | 64 | 1.17 | 1.22 | 0.6
2,097,152 | 512 | 1.17 | 1.23 | 0.7
16,777,216 | 4096 | 1.22 | 1.27 | 0.4
Table 5:  Speedup of the autoencoder compression algorithm compared to the
checkpointing solution. $\displaystyle\nabla$ represents the gradient and
$\displaystyle\mathcal{H}\boldsymbol{v}$ represents a Hessian-vector product.
The speedup approaches the ideal speedup as the number of DoFs increases while
the relative error decreases.

(a) Checkpoints

(b) Autoencoder compression

(c) $\displaystyle c$, km/s

Figure 5:  Inverse solutions on box domain using (a) checkpoints and (b)
autoencoder compression in space for problem with 16,777,216 DoF mesh. The
solutions have the same color scale with the inclusion rendered as an
isosurface at 1.67 km/s.

###  4.2 Results: compression in time

While it is intuitive that seismic waves exhibit strong spatial redundancy and
thus have high compressibility, it is not immediately obvious that temporal
compression schemes will yield acceptable results. We test the temporal
compression scheme presented in Section 3.2 on the earth-scale seismic inverse
problem. We build upon the framework established in [Bui-Thanh et al., 2012,
Bui-Thanh et al., 2013] to show that our proposed compression scheme is a
viable method for large-scale seismic inverse problems. The results shown here
follow the setup of Case II from [Bui-Thanh et al., 2013]. There are 3 source
locations in the Northern Hemisphere, one at the North Pole and two placed at
90 degree increments along the equator. Here, we consider sources modeled as
Dirac delta functions in space centered at the source location,
$\displaystyle\boldsymbol{x}_{0}$, and a Gaussian in time given by

$\boldsymbol{g}\left(\boldsymbol{x},t;\boldsymbol{x}_{0}\right)=\frac{1}{\sqrt{2\pi}\sigma}\delta\left(\boldsymbol{x}-\boldsymbol{x}_{0}\right)e^{-\frac{\left(t-t_{c}\right)^{2}}{2\sigma^{2}}}$

where $\displaystyle t_{c}=60$ seconds and $\displaystyle\sigma=20$ seconds.
We place 130 receivers spaced at 7.5 degree increments throughout North
America as show in Figure 6, measuring the velocity field. Additive noise with
standard deviation 0.002 is introduced into the data. The earth is modeled as
a sphere with radius 6,371km. Both the sources and receivers are buried at
10km depth.

Figure 6:  Three sources shown in green at the North Pole and along the
equator. 130 receivers shown in red distributed in North America spaced 7.5
degrees apart. Figure from [Bui-Thanh et al., 2013], Copyright ©2013 Society
for Industrial and Applied Mathematics. Reprinted with permission. All rights
reserved.

Course meshes cannot resolve the sufficiently high frequency waves required to
effectively image the earth and thus may not be used. We found that a maximum
frequency resolution of 0.05Hz was sufficient to obtain good inversion
results. Thus, the smallest domain size that we consider has a mesh with 11.2
million degrees of freedom. This corresponds to a discontinuous Galerkin
discretization with 174,481 elements using 3${}^{\text{rd}}$ order
polynomials. We solve this problem on 64 nodes (3,584 CPU cores) of Frontera
running for 10 hours. Speedup results compared to the checkpointing case are
shown in Table 6. Here we report the average relative $\displaystyle l_{2}$
error of the solution obtained using autoencoder compression to the solution
obtained using checkpointing. Similar to the box domain, there is greater
speedup in the $\displaystyle\mathcal{H}\boldsymbol{v}$ product computation
than the gradient computation due to paying the cost of compression during the
gradient solve. The difference in speedup between the box problem and the
earth problem is likely due to variable per-process problem size induced by
the non-uniform spherical mesh and related cache contention issues. As
mentioned in Section 3, cache effects can greatly impact the performance of
the compression algorithm.

Mesh DoFs |  $\displaystyle\nabla$ Speedup |  $\displaystyle\mathcal{H}\boldsymbol{v}$ Speedup | Relative $\displaystyle l_{2}$ error %
---|---|---|---
11,166,784 | 1.05 | 1.17 | 0.03
Table 6:  Speedup of the autoencoder compression algorithm compared to the
checkpointing solution. $\displaystyle\nabla$ represents the gradient and
$\displaystyle\mathcal{H}\boldsymbol{v}$ represents a Hessian-vector product.
Timing results were obtained by running on 64 nodes (3,584 CPU cores) of the
Frontera supercomputer.

To demonstrate scalability of the autoencoder compression approach, we solve
the same problem on a refined mesh with 139,227,136 DoFs. The full state
vector (3 velocities and 3 strains) then has 835,362,816 DoFs. This problem
was solved on 256 nodes of Frontera (14,336 CPU cores). Table 7 compares the
autoencoder compression results to the checkpoint solution. We find that there
is a slight reduction in terms of gradient speedup and a slight increase in
speedup of the Hessian-vector product compared to the course mesh solution
detailed in Table 8. One possible reason for the decrease in speedup of the
gradient computation is the increased per-process work of the larger problem.
The smaller problem with 11 million mesh DoFs has an average of 4,681 mesh
DoFs per process, while the larger problem has an average of 9,711 mesh DoFs
per process. We explored the effects of limited cache size on the timing of
autoencoder compression in Section 3 and found that the compression is more
sensitive to cache effects than decompression, at least for our
implementation. It is therefore not surprising that the speedup of the
gradient computation, which includes the cost of compression, is reduced when
each process is responsible for compressing a larger amount of data. While
speedup is desired in all computations, there is an order of magnitude more
Hessian-vector products than gradient computations required to compute the MAP
solution, with
$\displaystyle\mathcal{O}\left(100\right)\mathcal{H}\boldsymbol{v}$ products
and $\displaystyle\mathcal{O}\left(10\right)$ gradient evaluations. The
overall cost is dominated by $\displaystyle\mathcal{H}\boldsymbol{v}$
products, and so the observation that the gradient computation sees no speedup
from compression has little impact on the total computation time.

Mesh DoFs |  $\displaystyle\nabla$ Speedup |  $\displaystyle\mathcal{H}\boldsymbol{v}$ Speedup |  Relative error %
---|---|---|---
139,227,136 |  1.00 |  1.18 |  0.06
Table 7:  Inverse results for earth-scale seismic inverse problem with higher
mesh refinement, solved on 14,336 CPU cores. The gains from decompression over
resolving from checkpoints for computation of the gradient is completely
offset by the cost of compression. However, this cost is paid once and there
is similar speedup of the Hessian-vector product as seen on a smaller problem
in Table 8

Having shown that autoencoders are a viable compression technique for
accelerating seismic inverse problems, we now compare our results to the
state-of-the-art off-the-shelf compression technique proposed for seismic
inversion in [Kukreja et al., 2019]. While other off-the-shelf packages have
been proposed, the most comprehensive results in the literature are given for
the ZFP compression package [Lindstrom, 2014], especially in the seismic
inversion literature [Kukreja et al., 2019, Kukreja et al., 2022]. Therefore,
we compare the speed, accuracy, and compression ratio of our autoencoder
approach to the ZFP compression package. While the ZFP package does not
require normalization, we tried both with and without normalization and found
that normalizing the states before compressing was necessary to obtain
accurate results. Table 8 shows the inversion results using ZFP compression
for several different required accuracies. The relative errors for all methods
are still quite small, but ZFP results have qualitatively different solutions
than the checkpointing solution when the compression tolerance
$\displaystyle\eta=10^{-2}$ and $\displaystyle\eta=10^{-3}$ as shown in Figure
7. Thus, $\displaystyle\eta\leq 10^{-4}$ is required in order to obtain
comparable quality compression results to the autoencoder solution. While
autoencoder compression does not show a significant advantage over ZFP in
terms of speed, there is a substantial difference in the compression ratio,
and hence memory footprint. With this tolerance, the compression ratio of ZFP
is 8.4. On the other hand, the autoencoder compression algorithm achieved a
compression ratio of $\displaystyle 128$, 15 times higher than ZFP while
maintaining comparable accuracy.

(a) Checkpoint

(b) Autoencoder

(c) zfp, $\displaystyle\eta=10^{-5}$

(d) zfp, $\displaystyle\eta=10^{-4}$

(e) zfp, $\displaystyle\eta=10^{-3}$

(f) zfp, $\displaystyle\eta=10^{-2}$

Figure 7:  Comparison of inverse solutions obtained using autoencoder compression and the ZFP compression package with various compression tolerances, $\displaystyle\eta$. Results are shown at 130km depth on the same color scale.  Method |  $\displaystyle\nabla$ Speedup |  $\displaystyle\mathcal{H}\boldsymbol{v}$ Speedup |  Relative $\displaystyle l_{2}$ error % |  Compression Ratio
---|---|---|---|---
Autoencoder |  1.05 |  1.17 |  0.03 |  128
ZFP, $\displaystyle\eta=10^{-5}$ |  1.0 |  1.09 |  0.004 |  7.2
ZFP, $\displaystyle\eta=10^{-4}$ |  1.05 |  1.16 |  0.02 |  8.4
ZFP, $\displaystyle\eta=10^{-3}$ |  1.05 |  1.16 |  0.10 |  11.2
ZFP, $\displaystyle\eta=10^{-2}$ |  1.06 |  1.17 |  0.12 |  17.6
Table 8:  Comparison of ZFP compression to autoencoder compression for various
compression tolerances. There is not a significant difference in the speedup
of the autoencoder compression vs the ZFP compression, but the autoencoder is
able to achieve a much higher compression ratio (and hence much smaller memory
footprint) while still maintaining solution accuracy. Relative $\displaystyle
l_{2}$ error is computed with respect to the checkpoint solution.

###  4.3 DIAS regularization with autoencoder compression

As a capstone result, we show how our proposed autoencoder compression
approach can be combined with the DIAS regularization [Nguyen et al., 2022]
approach on the earth-scale seismic inverse problem. While the DIAS algorithm
has the benefit of only applying regularization in the inactive subspace,
speed is not one of it’s strengths. Simply estimating the active subspace for
a nonlinear inverse problem requires numerous samples of the gradient — each
evaluation of which requires 2 PDE solves. In the case of large-scale seismic
inversion, each gradient evaluation can take several minutes on dozens of
compute nodes. We propose to combine the DIAS regularization algorithm with
autoencoder compression to alleviate the high cost and memory footprint of
checkpointing in estimating the active subspace. Since we have already
discussed the speedup and memory footprint gained by the autoencoder
compression at length above, we will focus on the solution quality using DIAS
plus autoencoder compression in the following.

Let us briefly recap the DIAS regularization. We will focus on the DIAS-F (see
[Nguyen et al., 2022, Section 4] variant, which uses the full data misfit and
only applies regularization in the inactive subspace. First, the active
subspace is estimated via,

$\begin{bmatrix}\boldsymbol{W}_{1}&\boldsymbol{W}_{2}\end{bmatrix}\begin{bmatrix}\Lambda_{1}&0\\\
0&\Lambda_{2}\end{bmatrix}\begin{bmatrix}\boldsymbol{W}_{1}^{T}\\\
\boldsymbol{W}_{2}^{T}\end{bmatrix}={\boldsymbol{C}}:=\int\nabla_{\boldsymbol{u}}f\left(\boldsymbol{u}\right)\nabla_{\boldsymbol{u}}f\left(\boldsymbol{u}\right)^{T}\pi_{\text{{\text{prior}}}}\left(\boldsymbol{u}\right)d\boldsymbol{u},$
(7)

where $\displaystyle
f:=\frac{1}{2}\left\|\boldsymbol{\mathcal{F}}\left(\boldsymbol{u}\right)-\boldsymbol{d}\right\|^{2}_{{\Gamma}^{-1}_{\text{noise}}}$,
$\displaystyle\boldsymbol{W}_{1}$ are the eigenvectors of
$\displaystyle{\boldsymbol{C}}$ corresponding to the active subspace, and
$\displaystyle\boldsymbol{W}_{2}$ are the eigenvectors corresponding to the
inactive subspace. As in Section 3.3.1,
$\displaystyle\pi_{\text{{prior}}}\left(\boldsymbol{u}\right)$ is the prior
probability density from the Bayesian formulation of the inverse problem. The
DIAS solution is then defined to be

$\boldsymbol{u}_{\text{DIAS}}=\underset{\boldsymbol{u}}{\text{argmin
}}\frac{1}{2}\left\|\boldsymbol{\mathcal{F}}\left(\boldsymbol{u}\right)-\boldsymbol{d}\right\|^{2}_{\Gamma_{\text{noise}}^{-1}}+\frac{1}{2}\left\|\boldsymbol{W}_{2}^{T}\left(\boldsymbol{u}-\boldsymbol{u}_{0}\right)\right\|^{2}_{\left(\boldsymbol{W}_{2}^{T}\Gamma_{\text{prior}}\boldsymbol{W}_{2}\right)^{-1}},$
(8)

that is, it is a Tikhonov solution with regularization acting only on the
inactive space.

For large scale problems, it is infeasible to compute and store the large
matrix $\displaystyle\boldsymbol{W}_{2}$, since the active subspace tends to
be much smaller than the inactive subspace. To avoid explicitly computing and
storing the potentially large matrix $\displaystyle\boldsymbol{W}_{2}$, it can
be shown using [Schott, 2016, Theorem 5.8] that

$\left\|\boldsymbol{W}_{2}^{T}\left(\boldsymbol{u}-\boldsymbol{u}_{0}\right)\right\|^{2}_{\left(\boldsymbol{W}_{2}^{T}\Gamma_{\text{prior}}\boldsymbol{W}_{2}\right)^{-1}}=\left\|\boldsymbol{P}_{2}\left(\boldsymbol{u}-\boldsymbol{u}_{0}\right)\right\|^{2}_{\left(\boldsymbol{P}_{2}\Gamma_{\text{prior}}\boldsymbol{P}_{2}\right)^{\dagger}}$
(9)

where
$\displaystyle\boldsymbol{P}_{2}:=\boldsymbol{W}_{2}\boldsymbol{W}_{2}^{T}=\mathcal{I}-\boldsymbol{W}_{1}\boldsymbol{W}_{1}^{T}$
is a projection matrix onto the inactive subspace and
$\displaystyle\left(\boldsymbol{P}_{2}\Gamma_{\text{prior}}\boldsymbol{P}_{2}\right)^{\dagger}$
denotes the pseudoinverse.

Lest we conclude that we have a free lunch in computing the regularization
term arising from the active subspace prior via
$\displaystyle{\boldsymbol{W}}_{1}$, there is one more challenge that remains.
Note that we need to apply not the prior covariance, but its (pseudo-)inverse.
That is, we need to compute
$\displaystyle\left(\boldsymbol{P}_{2}\Gamma_{\text{prior}}\boldsymbol{P}_{2}\right)^{\dagger}$.
There are several factors at play here.

  1. 1.

$\displaystyle\Gamma_{\text{prior}}$ is often a sparse matrix, while
$\displaystyle\boldsymbol{P}_{2}\Gamma_{\text{prior}}\boldsymbol{P}_{2}$ is
very likely a dense matrix.

  2. 2.

We need to invert the dense matrix
$\displaystyle\boldsymbol{P}_{2}\Gamma_{\text{prior}}\boldsymbol{P}_{2}$ using
the pseudo-inverse.

  3. 3.

This matrix is in the dimension of the parameter, $\displaystyle n$.

For cases where it is feasible to form and invert
$\displaystyle\boldsymbol{P}_{2}\Gamma_{\text{prior}}\boldsymbol{P}_{2}$, it
is also feasible to simply compute $\displaystyle{\boldsymbol{W}}_{2}$
explicitly. Then we can work with the lower dimensional form given in (8)
which still involves inversion, but of a lower dimensional, full rank matrix.
For cases where this is not feasible, we need to develop more sophisticated
techniques for computing or estimating
$\displaystyle\boldsymbol{P}_{2}\Gamma_{\text{prior}}\boldsymbol{P}_{2}$.

The first thing one might try is to find an identity that allows us to break
the projection out of the pseudo-inverse, leaving
$\displaystyle\Gamma_{\text{prior}}$ by itself. We might hope that

$\left(\boldsymbol{P}_{2}\Gamma_{\text{prior}}\boldsymbol{P}_{2}\right)^{\dagger}=\boldsymbol{P}_{2}\Gamma_{\text{prior}}^{-1}\boldsymbol{P}_{2},$

but this is unfortunately not the case, as we show using Schur complements.
Decompose the prior covariance as

$\displaystyle\displaystyle\Gamma_{\text{prior}}=\left[\begin{array}[]{c|c}{\boldsymbol{W}}_{1}^{T}\Gamma_{\text{prior}}{\boldsymbol{W}}_{1}&{\boldsymbol{W}}_{1}^{T}\Gamma_{\text{prior}}{\boldsymbol{W}}_{2}\\\
\hline\cr{\boldsymbol{W}}_{2}^{T}\Gamma_{\text{prior}}{\boldsymbol{W}}_{1}&{\boldsymbol{W}}_{2}^{T}\Gamma_{\text{prior}}{\boldsymbol{W}}_{2}\end{array}\right]$

and the inverse prior covariance as

$\displaystyle\displaystyle\Gamma_{\text{prior}}^{-1}=\left[\begin{array}[]{c|c}{\boldsymbol{W}}_{1}^{T}\Gamma_{\text{prior}}^{-1}{\boldsymbol{W}}_{1}&{\boldsymbol{W}}_{1}^{T}\Gamma_{\text{prior}}^{-1}{\boldsymbol{W}}_{2}\\\
\hline\cr{\boldsymbol{W}}_{2}^{T}\Gamma_{\text{prior}}^{-1}{\boldsymbol{W}}_{1}&{\boldsymbol{W}}_{2}^{T}\Gamma_{\text{prior}}^{-1}{\boldsymbol{W}}_{2}\end{array}\right].$

Using Schur complements gives

$\displaystyle\displaystyle{\boldsymbol{W}}_{2}^{T}\Gamma_{\text{prior}}^{-1}{\boldsymbol{W}}_{2}$
$\displaystyle\displaystyle=\left[{\boldsymbol{W}}_{2}^{T}\Gamma_{\text{prior}}{\boldsymbol{W}}_{2}-{\boldsymbol{W}}_{2}^{T}\Gamma_{\text{prior}}{\boldsymbol{W}}_{1}\left({\boldsymbol{W}}_{1}^{T}\Gamma_{\text{prior}}{\boldsymbol{W}}_{1}\right)^{-1}{\boldsymbol{W}}_{1}^{T}\Gamma_{\text{prior}}{\boldsymbol{W}}_{2}\right]^{-1}.$
(10)

While the inverse can be expanded using the Sherman-Morrison-Woodbury formula,
the important observation is that the second term,

${\boldsymbol{W}}_{2}^{T}\Gamma_{\text{prior}}{\boldsymbol{W}}_{1}\left({\boldsymbol{W}}_{1}^{T}\Gamma_{\text{prior}}{\boldsymbol{W}}_{1}\right)^{-1}{\boldsymbol{W}}_{1}^{T}\Gamma_{\text{prior}}{\boldsymbol{W}}_{2}$

is not necessarily 0, except when $\displaystyle\Gamma_{\text{prior}}$ is
scaled identity. Thus, making the approximation
$\displaystyle\left(\boldsymbol{P}_{2}\Gamma_{\text{prior}}\boldsymbol{P}_{2}\right)^{\dagger}=\boldsymbol{P}_{2}\Gamma_{\text{prior}}^{-1}\boldsymbol{P}_{2}$
incurs some error in the general case. Even though there is some error, it
becomes necessary from a practical perspective to use the approximate form
$\displaystyle\boldsymbol{P}_{2}\Gamma_{\text{prior}}^{-1}\boldsymbol{P}_{2}$
for large scale problems.

Before the DIAS algorithm can be applied to the seismic inverse problem, we
must first define an algorithm for applying DIAS to nonlinear inverse
problems. One advantage of using the active subspace to define the data-
informed subspace is that the active subspace is formed from the global
average of the outer product of the gradient. This would indicate that the
DIAS algorithm can be applied naively to nonlinear inverse problems:

$\boldsymbol{u}_{\text{DIAS}}=\underset{\boldsymbol{u}}{\text{argmin
}}\frac{1}{2}\left\|\boldsymbol{\mathcal{F}}\left(\boldsymbol{u}\right)-\boldsymbol{d}\right\|^{2}_{\Gamma_{\text{noise}}^{-1}}+\frac{1}{2}\left\|\boldsymbol{P}_{2}\left(\boldsymbol{u}-\boldsymbol{u}_{0}\right)\right\|^{2}_{\left(\boldsymbol{P}_{2}\Gamma_{\text{prior}}\boldsymbol{P}_{2}\right)^{-1}}.$
(11)

However, the fact that the gradient samples are computed from draws of the
prior biases the active subspace toward the directions in which the misfit
function $\displaystyle f$ is most sensitive near the prior mean,
$\displaystyle\boldsymbol{u}_{0}$.

Recalling that the goal of the DIAS algorithm is to remove regularization in
directions that the data are more informative and that the DIAS solutions are
likely closer to the usual Bayesian solution than to the prior mean for a
well-chosen prior, we propose a two-step algorithm:

  1. 1.

Approximately solve the usual MAP problem for
$\displaystyle\boldsymbol{u}_{\text{MAP}}$ using the full prior:

$\boldsymbol{u}_{\text{MAP}}=\underset{\boldsymbol{u}}{\text{argmin
}}\frac{1}{2}\left\|\boldsymbol{\mathcal{F}}\left(\boldsymbol{u}\right)-\boldsymbol{d}\right\|^{2}_{{\Gamma}^{-1}_{\text{noise}}}+\frac{1}{2}\left\|\boldsymbol{u}-\boldsymbol{u}_{0}\right\|^{2}_{{\Gamma}^{-1}_{\text{prior}}}.$
(12)

  2. 2.

Compute the active subspace centered at
$\displaystyle\boldsymbol{u}_{\text{MAP}}$ and, with initial guess
$\displaystyle\boldsymbol{u}_{\text{MAP}}$, solve DIAS problem:

$\boldsymbol{u}_{\text{DIAS}}=\underset{\boldsymbol{u}}{\text{argmin
}}\frac{1}{2}\left\|\boldsymbol{\mathcal{F}}\left(\boldsymbol{u}\right)-\boldsymbol{d}\right\|^{2}_{\Gamma_{\text{noise}}^{-1}}+\frac{1}{2}\left\|\boldsymbol{P}_{2}\left(\boldsymbol{u}-\boldsymbol{u}_{0}\right)\right\|^{2}_{\boldsymbol{P}_{2}{\Gamma}^{-1}_{\text{prior}}\boldsymbol{P}_{2}}.$

Even with acceleration via autoencoder compression, the cost of gradient
evaluations is still high so we estimate the active subspace using 30 gradient
samples and take the active subspace to have dimension 5. The difference
between the DIAS and the MAP estimate are small in this case. Compared to the
true solution, the MAP estimate has average relative error of 0.3293% while
the DIAS refined solution has average relative error of 0.3291%. This is
reflected visually in Figure 8.

(a) MAP

(b) DIAS

(c) True

---
$\displaystyle c,$ km/s

Figure 8:  The DIAS solution does not differ noticeably from the MAP solution
for the seismic inverse problem. Results shown for mesh with 11,166,784 DoFs.

##  5 Conclusions

In this paper, we proposed an autoencoder compression algorithm to mitigate
the high storage requirements for solving large-scale time-dependent PDE
constrained inverse problems. We show the feasibility and scalability of this
approach against the popular checkpointing strategy for a seismic inverse
problem on both a simple box domain with uniform refinement and a complex
spherical domain with adaptive mesh refinement and a nonconforming mesh. In
order to effectively train the autoencoder, we present a data generation
procedure based on the Bayesian formulation that provides a problem-focused
data generation scheme. We proposed two separate autoencoder compression
algorithms for spatial and temporal domains. The scalability of each algorithm
was shown, providing a clear close-to-ideal speed advantage on large scale
problems compared to the traditional checkpointing approach and a substantial
improvement in the compression ratio of the state-of-the-art floating point
compression package, ZFP. While there are clear advantages of using
autoencoders to perform compression, care must be taken in their
implementation in order for there to be an advantage over the checkpointing
approach. As an important application of the proposed autoencoder algorithm,
we combined the proposed autoencoder compression approach with the DIAS prior
to show how the DIAS method can be affordably extended to large scale problems
without the need of checkpointing and large memory.

##  6 Acknowledgements

We would like to thank Georg Stadler for his many correspondences and helpful
advice on stably solving the forward and inverse problem in large scale runs;
Jau-Uei Chen for sharing his insight on the discontinuous Galerkin method;
Sheroze Sheriffdeen and Hwan Goh for their many discussions on autoencoders,
compression, and HPC; and the Texas Advanced Computing Center (TACC) at the
University of Texas at Austin for providing HPC resources that contributed to
the results presented in this work. URL: http://www.tacc.utexas.edu. This
research is partially funded by the National Science Foundation awards NSF-
OAC-2212442, NSF-2108320, NSF-1808576 and NSF-CAREER-1845799; and by the
Department of Energy awards DE-SC0018147 and DE-SC0022211.

## References

  * [int, ] Intel® xeon® platinum 8280 processor (38.5m cache, 2.70 ghz) - product specifications.
  * [Abu Alsheikh et al., 2014] Abu Alsheikh, M., Poh, P. K., Lin, S., Tan, H.-P., and Niyato, D. (2014).  Efficient data compression with error bound guarantee in wireless sensor networks.  In Proceedings of the 17th ACM international conference on Modeling, analysis and simulation of wireless and mobile systems, pages 307–311.
  * [Akturk and Karpuzcu, 2018] Akturk, I. and Karpuzcu, U. R. (2018).  Trading computation for communication: A taxonomy of data recomputation techniques.  IEEE Transactions on Emerging Topics in Computing, 9(1):496–506.
  * [Baldi and Hornik, 1989] Baldi, P. and Hornik, K. (1989).  Neural networks and principal component analysis: Learning from examples without local minima.  Neural networks, 2(1):53–58.
  * [Bharadwaj et al., 2021] Bharadwaj, K. K., Srivastava, A., Panda, M. K., Singh, Y. D., Maharana, R., Mandal, K., Singh, M., Singh, D., Das, M., Murmu, D., et al. (2021).  Computational intelligence in vaccine design against covid-19.  In Computational intelligence methods in COVID-19: surveillance, prevention, prediction and diagnosis, pages 311–329. Springer.
  * [Bhattacharjee and Wells, 2021] Bhattacharjee, A. and Wells, J. (2021).  Preface to special topic: Building the bridge to the exascale—applications and opportunities for plasma physics.
  * [Bickelhaupt and Baerends, 2000] Bickelhaupt, F. M. and Baerends, E. J. (2000).  Kohn-sham density functional theory: predicting and understanding chemistry.  Reviews in computational chemistry, pages 1–86.
  * [Bjerge et al., 2022] Bjerge, K., Mann, H. M., and Høye, T. T. (2022).  Real-time insect tracking and monitoring with computer vision and deep learning.  Remote Sensing in Ecology and Conservation, 8(3):315–327.
  * [Blesser, 1969] Blesser, B. (1969).  Audio dynamic range compression for minimum perceived distortion.  IEEE Transactions on Audio and Electroacoustics, 17(1):22–32.
  * [Boehm et al., 2016] Boehm, C., Hanzich, M., de la Puente, J., and Fichtner, A. (2016).  Wavefield compression for adjoint methods in full-waveform inversion.  Geophysics, 81(6):R385–R397.
  * [Bui-Thanh et al., 2012] Bui-Thanh, T., Burstedde, C., Ghattas, O., Martin, J., Stadler, G., and Wilcox, L. C. (2012).  Extreme-scale uq for bayesian inverse problems governed by pdes.  In SC’12: Proceedings of the international conference on high performance computing, networking, storage and analysis, pages 1–11. IEEE.
  * [Bui-Thanh et al., 2013] Bui-Thanh, T., Ghattas, O., Martin, J., and Stadler, G. (2013).  A computational framework for infinite-dimensional bayesian inverse problems part i: The linearized case, with application to global seismic inversion.  SIAM Journal on Scientific Computing, 35(6):A2494–A2523.
  * [Bukowski et al., 2010] Bukowski, R., Sun, Q., Howard, M., and Pillardy, J. (2010).  Biohpc: computational biology application suite for high performance computing.  Journal of Biomolecular Techniques: JBT, 21(3 Suppl):S23.
  * [Cheng et al., 2018] Cheng, Z., Sun, H., Takeuchi, M., and Katto, J. (2018).  Deep convolutional autoencoder-based lossy image compression.  In 2018 Picture Coding Symposium (PCS), pages 253–257. IEEE.
  * [Cyr et al., 2015] Cyr, E. C., Shadid, J., and Wildey, T. (2015).  Towards efficient backward-in-time adjoint computations using data compression techniques.  Computer Methods in Applied Mechanics and Engineering, 288:24–44.
  * [De Jong et al., 2010] De Jong, W. A., Bylaska, E., Govind, N., Janssen, C. L., Kowalski, K., Müller, T., Nielsen, I. M., van Dam, H. J., Veryazov, V., and Lindh, R. (2010).  Utilizing high performance computing for chemistry: parallel computational chemistry.  Physical Chemistry Chemical Physics, 12(26):6896–6920.
  * [Denis et al., 2022] Denis, A., Jeannot, E., and Swartvagher, P. (2022).  Modeling memory contention between communications and computations in distributed hpc systems.  In IPDPS-2022-IEEE International Parallel and Distributed Processing Symposium Workshops, page 10.
  * [Dillon et al., 2021] Dillon, B., Plehn, T., Sauer, C., and Sorrenson, P. (2021).  Better latent spaces for better autoencoders.  SciPost Physics, 11(3):061.
  * [Doi et al., 2019] Doi, J., Takahashi, H., Raymond, R., Imamichi, T., and Horii, H. (2019).  Quantum computing simulator on a heterogenous hpc system.  In Proceedings of the 16th ACM International Conference on Computing Frontiers, pages 85–93.
  * [Duarte et al., 2020] Duarte, E. F., da Costa, C. A., de Araújo, J. M., Wang, Y., and Rao, Y. (2020).  Seismic shot-encoding schemes for waveform inversion.  Journal of Geophysics and Engineering, 17(5):906–913.
  * [Epanomeritakis et al., 2008] Epanomeritakis, I., Akçelik, V., Ghattas, O., and Bielak, J. (2008).  A newton-cg method for large-scale three-dimensional elastic full-waveform seismic inversion.  Inverse Problems, 24(3):034015.
  * [Esteva et al., 2021] Esteva, A., Chou, K., Yeung, S., Naik, N., Madani, A., Mottaghi, A., Liu, Y., Topol, E., Dean, J., and Socher, R. (2021).  Deep learning-enabled medical computer vision.  NPJ digital medicine, 4(1):1–9.
  * [Fedeli et al., 2022] Fedeli, L., Huebl, A., Boillod-Cerneux, F., Clark, T., Gott, K., Hillairet, C., Jaure, S., Leblanc, A., Lehe, R., Myers, A., et al. (2022).  Pushing the frontier in the design of laser-based electron accelerators with groundbreaking mesh-refined particle-in-cell simulations on exascale-class supercomputers.  In 2022 SC22: International Conference for High Performance Computing, Networking, Storage and Analysis (SC), pages 25–36. IEEE Computer Society.
  * [Fichtner et al., 2006] Fichtner, A., Bunge, H.-P., and Igel, H. (2006).  The adjoint method in seismology: I. theory.  Physics of the Earth and Planetary Interiors, 157(1-2):86–104.
  * [Ge et al., 2013] Ge, H., Wang, Y., Li, C., Chen, N., Xie, Y., Xu, M., He, Y., Gu, X., Wu, R., Gu, Q., et al. (2013).  Molecular dynamics-based virtual screening: accelerating the drug discovery process by high-performance computing.  Journal of chemical information and modeling, 53(10):2757–2764.
  * [Giles and Pierce, 2000] Giles, M. B. and Pierce, N. A. (2000).  An introduction to the adjoint approach to design.  Flow, turbulence and combustion, 65(3):393–415.
  * [Goh et al., 2019] Goh, H., Sheriffdeen, S., Wittmer, J., and Bui-Thanh, T. (2019).  Solving bayesian inverse problems via variational autoencoders.  arXiv preprint arXiv:1912.04212.
  * [Griewank and Walther, 1997] Griewank, A. and Walther, A. (1997).  Treeverse: An implementation of checkpointing for the reverse or adjoint mode of computational differentiation.  Preprint IOKOMO-04-1997, Technische Universitat, Dresden.
  * [Habashy et al., 2011] Habashy, T., Abubakar, A., Pan, G., and Belani, A. (2011).  Source-receiver compression scheme for full-waveform seismic inversion.  Geophysics, 76(4):R95–R108.
  * [Hinton and Salakhutdinov, 2006] Hinton, G. E. and Salakhutdinov, R. R. (2006).  Reducing the dimensionality of data with neural networks.  science, 313(5786):504–507.
  * [Imani et al., 2019] Imani, M., Gupta, S., Kim, Y., Zhou, M., and Rosing, T. (2019).  Digitalpim: Digital-based processing in-memory for big data acceleration.  In Proceedings of the 2019 on Great Lakes Symposium on VLSI, pages 429–434.
  * [Ji et al., 2012] Ji, S., Xu, W., Yang, M., and Yu, K. (2012).  3d convolutional neural networks for human action recognition.  IEEE transactions on pattern analysis and machine intelligence, 35(1):221–231.
  * [Kaipio and Somersalo, 2006] Kaipio, J. and Somersalo, E. (2006).  Statistical and computational inverse problems, volume 160.  Springer Science & Business Media.
  * [Kamath et al., 2019] Kamath, U., Liu, J., and Whitaker, J. (2019).  Deep learning for NLP and speech recognition, volume 84.  Springer.
  * [Kates-Harbeck et al., 2019] Kates-Harbeck, J., Svyatkovskiy, A., and Tang, W. (2019).  Predicting disruptive instabilities in controlled fusion plasmas through deep learning.  Nature, 568(7753):526–531.
  * [Kaur and Choudhary, 2016] Kaur, R. and Choudhary, P. (2016).  A review of image compression techniques.  Int. J. Comput. Appl, 142(1):8–11.
  * [Kingma and Welling, 2013] Kingma, D. P. and Welling, M. (2013).  Auto-encoding variational bayes.  arXiv preprint arXiv:1312.6114.
  * [Kneer et al., 2021] Kneer, S., Sayadi, T., Sipp, D., Schmid, P., and Rigas, G. (2021).  Symmetry-aware autoencoders: s-pca and s-nlpca.  arXiv preprint arXiv:2111.02893.
  * [Kowalski et al., 2021] Kowalski, K., Bair, R., Bauman, N. P., Boschen, J. S., Bylaska, E. J., Daily, J., de Jong, W. A., Dunning Jr, T., Govind, N., Harrison, R. J., et al. (2021).  From nwchem to nwchemex: Evolving with the computational chemistry landscape.  Chemical reviews, 121(8):4962–4998.
  * [Kukreja et al., 2019] Kukreja, N., Hückelheim, J., Louboutin, M., Hovland, P., and Gorman, G. (2019).  Combining checkpointing and data compression to accelerate adjoint-based optimization problems.  In European Conference on Parallel Processing, pages 87–100. Springer.
  * [Kukreja et al., 2022] Kukreja, N., Hückelheim, J., Louboutin, M., Washbourne, J., Kelly, P. H., and Gorman, G. J. (2022).  Lossy checkpoint compression in full waveform inversion: a case study with zfpv0. 5.5 and the overthrust model.  Geoscientific Model Development, 15(9):3815–3829.
  * [Ladjal et al., 2019] Ladjal, S., Newson, A., and Pham, C.-H. (2019).  A pca-like autoencoder.  arXiv preprint arXiv:1904.01277.
  * [Lee et al., 2022] Lee, J., Gong, Q., Choi, J., Banerjee, T., Klasky, S., Ranka, S., and Rangarajan, A. (2022).  Error-bounded learned scientific data compression with preservation of derived quantities.  Applied Sciences, 12(13):6718.
  * [Lewis and Knowles, 1990] Lewis, A. and Knowles, G. (1990).  Video compression using 3d wavelet transforms.  Electronics Letters, 26(6):396–398.
  * [Lindstrom, 2014] Lindstrom, P. (2014).  Fixed-rate compressed floating-point arrays.  IEEE transactions on visualization and computer graphics, 20(12):2674–2683.
  * [Liu et al., 2021a] Liu, J., Di, S., Zhao, K., Jin, S., Tao, D., Liang, X., Chen, Z., and Cappello, F. (2021a).  Exploring autoencoder-based error-bounded compression for scientific data.  In 2021 IEEE International Conference on Cluster Computing (CLUSTER), pages 294–306. IEEE.
  * [Liu et al., 2021b] Liu, T., Wang, J., Liu, Q., Alibhai, S., Lu, T., and He, X. (2021b).  High-ratio lossy compression: Exploring the autoencoder to compress scientific data.  IEEE Transactions on Big Data.
  * [Liu et al., 2021c] Liu, Y., Liu, X., Li, F., Fu, H., Yang, Y., Song, J., Zhao, P., Wang, Z., Peng, D., Chen, H., et al. (2021c).  Closing the" quantum supremacy" gap: achieving real-time simulation of a random quantum circuit using a new sunway supercomputer.  In Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis, pages 1–12.
  * [Mandrà et al., 2021] Mandrà, S., Marshall, J., Rieffel, E. G., and Biswas, R. (2021).  Hybridq: A hybrid simulator for quantum circuits.  In 2021 IEEE/ACM Second International Workshop on Quantum Computing Software (QCS), pages 99–109. IEEE.
  * [McFarlane and Biktasheva, 2008] McFarlane, R. and Biktasheva, I. V. (2008).  Beatbox-a computer simulation environment for computational biology of the heart.  In Visisions of Computer Science—BCS International Academic Conference, pages 99–109. British Computer Society.
  * [Nadeau and Bengio, 1999] Nadeau, C. and Bengio, Y. (1999).  Inference for the generalization error.  Advances in neural information processing systems, 12.
  * [Nguyen et al., 2022] Nguyen, H., Wittmer, J., and Bui-Thanh, T. (2022).  Dias: A data-informed active subspace regularization framework for inverse problems.  Computation, 10(3):38.
  * [Peng et al., 2020] Peng, I., Wu, K., Ren, J., Li, D., and Gokhale, M. (2020).  Demystifying the performance of hpc scientific applications on nvm-based memory systems.  In 2020 IEEE International Parallel and Distributed Processing Symposium (IPDPS), pages 916–925. IEEE.
  * [Plaut, 2018] Plaut, E. (2018).  From principal subspaces to principal components with linear autoencoders.  arXiv preprint arXiv:1804.10253.
  * [Quarteroni et al., 2010] Quarteroni, A., Sacco, R., and Saleri, F. (2010).  Numerical mathematics, volume 37.  Springer Science & Business Media.
  * [Reichstein et al., 2019] Reichstein, M., Camps-Valls, G., Stevens, B., Jung, M., Denzler, J., Carvalhais, N., et al. (2019).  Deep learning and process understanding for data-driven earth system science.  Nature, 566(7743):195–204.
  * [Scales and Tenorio, 2001] Scales, J. A. and Tenorio, L. (2001).  Prior information and uncertainty in inverse problems.  Geophysics, 66(2):389–397.
  * [Schmidt and Hildebrandt, 2017] Schmidt, B. and Hildebrandt, A. (2017).  Next-generation sequencing: big data meets high performance computing.  Drug discovery today, 22(4):712–717.
  * [Schott, 2016] Schott, J. R. (2016).  Matrix analysis for statistics.  John Wiley & Sons.
  * [Stanzione et al., 2020] Stanzione, D., West, J., Evans, R. T., Minyard, T., Ghattas, O., and Panda, D. K. (2020).  Frontera: The evolution of leadership computing at the national science foundation.  In Practice and Experience in Advanced Research Computing, pages 106–111.
  * [Sukumar et al., 2021] Sukumar, S. R., Balma, J. A., Rickett, C. D., Maschhoff, K. J., Landman, J., Yates, C. R., Chittiboyina, A. G., Peterson, Y. K., Vose, A., Byler, K., et al. (2021).  The convergence of hpc, ai and big data in rapid-response to the covid-19 pandemic.  In Smoky Mountains Computational Sciences and Engineering Conference, pages 157–172. Springer.
  * [Usevitch, 2001] Usevitch, B. E. (2001).  A tutorial on modern lossy wavelet image compression: foundations of jpeg 2000.  IEEE signal processing magazine, 18(5):22–35.
  * [Wang et al., 2009] Wang, Q., Moin, P., and Iaccarino, G. (2009).  Minimal repetition dynamic checkpointing algorithm for unsteady adjoint calculation.  SIAM Journal on Scientific Computing, 31(4):2549–2567.
  * [Wilcox et al., 2015] Wilcox, L. C., Stadler, G., Bui-Thanh, T., and Ghattas, O. (2015).  Discretely exact derivatives for hyperbolic pde-constrained optimization problems discretized by the discontinuous galerkin method.  Journal of Scientific Computing, 63(1):138–162.
  * [Wu et al., 2017] Wu, K., Ober, F., Hamlin, S., and Li, D. (2017).  Early evaluation of intel optane non-volatile memory with hpc i/o workloads.  arXiv preprint arXiv:1708.02199.
  * [Yang, 2009] Yang, J. (2009).  Newton-conjugate-gradient methods for solitary wave computations.  Journal of Computational Physics, 228(18):7007–7024.
  * [You et al., 2019] You, Y., Li, J., Reddi, S., Hseu, J., Kumar, S., Bhojanapalli, S., Song, X., Demmel, J., Keutzer, K., and Hsieh, C.-J. (2019).  Large batch optimization for deep learning: Training bert in 76 minutes.  arXiv preprint arXiv:1904.00962.
  * [Zhang and Constantinescu, 2023] Zhang, H. and Constantinescu, E. M. (2023).  Optimal checkpointing for adjoint multistage time-stepping schemes.  Journal of Computational Science, 66:101913.
  * [Zhao et al., 2020] Zhao, Y., Chen, X., Wang, Y., Li, C., You, H., Fu, Y., Xie, Y., Wang, Z., and Lin, Y. (2020).  Smartexchange: Trading higher-cost memory storage/access for lower-cost computation.  In 2020 ACM/IEEE 47th Annual International Symposium on Computer Architecture (ISCA), pages 954–967. IEEE.