# High resolution 3D ultrasonic breast imaging by time-domain full waveform inversion

Felix Lucka, Mailyn Pérez-Liva, Bradley E Treeby, Ben T Cox

► To cite this version:

Felix Lucka, Mailyn Pérez-Liva, Bradley E Treeby, Ben T Cox. High resolution 3D ultrasonic breast imaging by time-domain full waveform inversion. Inverse Problems, 2021, 38 (2), pp.025008. 10.1088/1361-6420/AC3B64 . hal-05481881

HAL Id: hal-05481881

https://u-paris.hal.science/hal-05481881v1

Submitted on 28 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

PAPER • OPEN ACCESS

High resolution 3D ultrasonic breast imaging by
time-domain full waveform inversion

To cite this article: Felix Lucka et al 2022 *Inverse Problems* **38** 025008

View the [article online](#) for updates and enhancements.

You may also like

- Supercomputer technologies in inverse problems of ultrasound tomography
  Alexander V Goncharsky and Sergey Y Romanov

- Pulsed laser diode based optoacoustic imaging of biological tissues
  Paul Kumar Upputuri and Manojit Pramanik

- Influence of source electrode metal work function on polar gate prompted source hole plasma in arsenide/antimonide tunneling interfaced junctionless TFET
  Samriti Sharma and Rishu Chaujar
# High resolution 3D ultrasonic breast imaging by time-domain full waveform inversion

Felix Lucka<fnref n="1" />,<fnref n="2" />,* , Mailyn Pérez-Liva<fnref n="3" />, , Bradley E Treeby<fnref n="4" /> and Ben T Cox<fnref n="4" />

E-mail: Felix.Lucka@cwi.nl

Received 4 February 2021, revised 31 August 2021

Accepted for publication 18 November 2021

Published 30 December 2021

<figure><img src="image_3.png" /></figure>

CrossMark

## Abstract

Ultrasound tomography (UST) scanners allow quantitative images of the human breast's acoustic properties to be derived with potential applications in screening, diagnosis and therapy planning. Time domain full waveform inversion (TD-FWI) is a promising UST image formation technique that fits the parameter fields of a wave physics model by gradient-based optimization. For high resolution 3D UST, it holds three key challenges: firstly, its central building block, the computation of the gradient for a single US measurement, has a restrictively large memory footprint. Secondly, this building block needs to be computed for each of the 10³–10⁴ measurements, resulting in a massive parallel computation usually performed on large computational clusters for days. Lastly, the structure of the underlying optimization problem may result in slow progression of the solver and convergence to a local minimum. In this work, we design and evaluate a comprehensive computational strategy to overcome these challenges: firstly, we exploit a gradient computation based on time reversal that dramatically reduces the memory footprint at the expense of one additional wave simulation per source. Secondly, we break the dependence on the number of measurements by using source encoding (SE) to compute stochastic gradient

<figure><img src="image_4.png" /></figure>

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹ Centrum Wiskunde & Informatica, Science Park 123, Amsterdam, The Netherlands</footnote>

<footnote n="2">² Centre for Medical Image Computing, University College of London, WC1E 6BT London, United Kingdom</footnote>

<footnote n="3">³ Université de Paris, PARCC, INSERM, 56 Rue Leblanc, 75015 Paris, France</footnote>

<footnote n="4">⁴ Department of Medical Physics and Biomedical Engineering, University College of London, WC1E 6BT London, United Kingdom</footnote>

<footnote>* Author to whom any correspondence should be addressed.</footnote>
estimates. Also we describe a more accurate, TD-specific SE technique with a finer variance control and use a state-of-the-art stochastic LBFGS method. Lastly, we design an efficient TD multi-grid scheme together with preconditioning to speed up the convergence while avoiding local minima. All components are evaluated in extensive numerical proof-of-concept studies simulating a bowl-shaped 3D UST breast scanner prototype. Finally, we demonstrate that their combination allows us to obtain an accurate $442 \times 442 \times 222$ voxel image with a resolution of 0.5 mm using Matlab on a single GPU within 24 h.

Keywords: full waveform inversion, ultrasound tomography, breast imaging, computational inverse problems, high resolution 3D

(Some figures may appear in colour only in the online journal)

# 1. Introduction

Screening using x-ray-based mammography has saved many lives by catching early-stage breast cancer. Nevertheless, there remain concerns related to tissue superposition, overdiagnosis, the effect of ionising radiation, the lower sensitivity in dense breasts, and the pain caused by breast compression [65]. Alternative breast imaging approaches have therefore been explored, including MRI, which is expensive and time-intensive [47], optical approaches such as photoacoustic tomography [63] and diffuse optical tomography [101] which are limited by the significant scattering of light by tissue, and ultrasound (US). In conventional US imaging, high frequency acoustic waves are transmitted into the breast from a hand-held linear array of typically 128 detection elements. The back-scattered waves are detected by the same array and used to form qualitative 2D images in real time. It is often used as an adjunct modality to supplement mammography [44, 76], but the operator-dependence means that hand-held US requires an expert user. It is also time-consuming when imaging a whole breast, making it unsuitable as a screening modality. In an attempt to overcome these drawbacks, automated breast US (ABUS) has been introduced, which uses a machine to scan a linear array perpendicular to its length over the breast to form a qualitative 3D reflection image [91] but it does not give a quantitative image of tissue properties. In contrast, a principal aim in *ultrasound tomography* (UST) of the breast [17, 23] is to obtain accurate high resolution *quantitative* images of the tissue's acoustic properties. To achieve this, a pulse of US is transmitted into a pendant breast and measurements are made of the resulting time-varying acoustic field at multiple detector positions around the breast. This is repeated for multiple source positions (or, more generally, source distributions). The set of acoustic time series thus measured will incorporate both the unscattered (directly transmitted) and scattered parts of the acoustic field. To facilitate practical measurement times, large arrays are typically used, consisting of many hundreds of elements. In recent years, a number of such systems have been developed for UST breast imaging, including ones based on ring arrays [1, 24, 70], rotating planar or linear arrays [3, 48, 62], and bowl-shaped detector arrays [2, 36]. Once the data has been measured, the challenge is to use it to form accurate quantitative images of the tissue's acoustic properties.

## 1.1. Image reconstruction for ultrasound tomography

Typically three types of images have been produced from UST data: quantitative images of the acoustic attenuation and the speed of sound, and qualitative images indicating the 'strength' of the scattering at each point in the image [20]. To obtain these images requires a nonlinear
inverse scattering problem to be tackled. The various approaches to solving this, which can
be characterised by the approximations and assumptions they make, represent different trade-
offs between computational efficiency and imaging accuracy. UST systems and measurement
protocols are therefore often designed and optimized for a particular image formation method,
unlike other imaging modalities where it is typically the other way around. There is an extensive
literature on UST reconstruction, and it shares much common ground with other areas of wave-
based imaging, such as in seismology and non-destructive testing. This is not the place for a
comprehensive account of that history, so we give just a few examples of typical reconstruction
approaches, focusing on quantitative reconstructions of sound speed as that is the topic of this
paper.

Travel-time or time-of-flight tomography, widely used in geophysics [19], assumes the wave propagation can be accurately described in terms of rays: either straight rays [12, 53, 64, 70, 74], thereby neglecting scattering, diffraction and refraction, or bent rays [5, 22, 52, 54, 58, 59], which can account for refraction to some extent. The measured travel times are linearly related to the integral of the slowness (reciprocal of the sound speed) along the rays, and the system matrix, once constructed, can be inverted using standard linear algebra approaches. Ray-based approaches can lead to robust and computationally efficient algorithms to recover smoothed sound speed images, but usually fine details are lost. Because the significance of diffraction depends on the relative sizes of the scatterer and wavelength, the ray approximation becomes more accurate the higher the frequency used. However, at higher frequencies the sources and detectors become more directional, reducing the possible ray paths across the imaged region, and the acoustic attenuation of biological tissues is greater.

Another approach, sometimes called *diffraction tomography*, considers the material heterogeneities to be weakly scattering perturbations over a background and uses Born or Rytov-type linearisations of the nonlinear scattering problem to arrive at linear methods of image formation. The free-space Green's function is used to estimate the effect of the perturbations on the field. In general for breast tissue, the weak scattering approximation does not hold, so adaptations of these linearisations are required to develop techniques suitable for quantitative breast imaging [49, 66].

To solve the full nonlinear inverse scattering tomography problem requires an iterative approach. One way is to iteratively update both the forward operator and the unknown parameter in the linearised approximations above. In *bent-ray reconstructions*, for example, both the rays and the unknown sound speed are updated at each iteration [52]. In the *distorted Born iterative method* [16, 46], the Green’s function and the unknown parameter are updated. This works well for low contrast media, but suffers from the failure of the Born series to converge for large contrast heterogeneities.

Before we introduce the particular approach we follow in this work in more detail, we remark that as in many fields, there has been a recent flurry of interest into the question of whether deep neural networks can be utilized for ultrasonic imaging [89], see, e.g., [4, 21, 28–32, 55, 56, 61, 89, 100] for approaches to improve the speed and accuracy of 2D UST reconstructions. While it seems likely that neural networks will find wide application in UST, there are reasons to think that, as in other imaging modalities, they will complement rather than supersede existing approaches [6]. In particular in high resolution 3D settings, neural network components that link the object to be imaged to the measured data can hardly be fully learned. Efficient and accurate implementations will always contain components derived from physical models accurately describing wave propagation in biological tissue according to the best of our understanding.
## 1.2. Full waveform inversion

An approach that has gained increasing attention in recent years has been to formulate the image formation as a optimization problem, essentially fitting the measured data to a forward model of the acoustic propagation. This way of tackling the inverse scattering tomography problem is known as *waveform tomography*, *model-based inversion* or *full waveform inversion* (FWI). Because the framework of FWI allows all the measured data to be used in the inversion—the whole time series, or all the frequencies—it can in principle produce quantitatively correct, high resolution images of the acoustic properties with many fewer measurements than is required for travel-time tomography (which essentially reduces each time series to a single value). There is a sizable literature on FWI in the seismic community, see [86] for a recent review.

FWI constitutes a large-scale, non-convex, PDE-constrained optimization problem. In addition to the intrinsic difficulties of such optimization problems, iterative optimization schemes need to run numerous numerical wave simulations to compare artificial UST measurements caused by the current guess of the acoustic parameters to the real US data obtained. Depending on the scan design and the spatial resolution desired, this can present a very significant computational challenge.

Dedicated hardware design can circumvent this bottleneck to some extent. For instance, when using a 2D ring-array focused in the plane, 2D FWI can be used to reconstruct the 3D volume slice-by-slice [24, 38, 58], but the risk of out-of-plane artefacts and the poor image resolution perpendicular to the plane makes this solution less than ideal [73]. The approach taken by Wiskin et al [93–96] has been to design the hardware with a planar geometry so that a computationally-efficient approximation to the wave simulation can be exploited.

There will be a highest usable frequency dictated by inherent attenuation through the breast and the sensitivity of the hardware, and this will limit the achievable image resolution. If this frequency is too high for the computations to be practicable then reducing the highest frequency, at the loss of the achievable image resolution, may be a compromise it is necessary to make. Goncharsky et al [37–42] follow this idea and consider using low-frequency transducers. They further combine with a ‘2.5’D slice-by-slice reconstruction approach.

However, there is not always the freedom to choose the hardware. For instance, in the case we are interested in, the scanner has been designed to be optimal for *photoacoustic tomography* (PAT) of the breast [2] and uses a hemispherical array with unfocused, broadband US transducers. UST, which will be an adjunct modality, will be performed using the same array and we aim to obtain the same, isotropic sub-millimeter resolution across the whole breast volume. Another example of a fully 3D UST system, also a bowl array, is the KIT 3D USCT system [35, 36], which was designed for obtaining both high-quality reflectivity images and quantitative tissue images. To realize FWI for such UST scenarios, large high-performance computing clusters are typically used [13, 38, 40] and/or the computation takes days, which severely limits the range of clinical applications for which these methods are viable. A notable exception is given by the work presented in [8], which we were made aware of during the writing of this article.

## 1.3. Paper scope and structure

The key focus and contribution of our work is to develop and demonstrate a comprehensive computational strategy that achieves accurate high-resolution, 3D FWI for breast UST with a hemispherical array using only moderate computational resources (at least one GPU) and a self-imposed time limit of 24 h. This would allow the method to fit into clinical trajectories for diagnosis and therapy planing at competitive costs compared to other modalities such as MRI.
<table><caption>**Table 1.** List of commonly occurring abbreviations.</caption><thead><tr><td>Abbreviation</td><td>Meaning</td><td>Reference</td></tr></thead><tbody><tr><td>FD</td><td>Frequency domain</td><td>Section 2.3</td></tr><tr><td>FWI</td><td>Full waveform inversion</td><td>Section 2.3</td></tr><tr><td>US</td><td>Ultrasound</td><td>Section 1</td></tr><tr><td>UST</td><td>Ultrasound tomography</td><td>Section 2.2</td></tr><tr><td>SGD</td><td>Stochastic gradient descent</td><td>Section 3.2</td></tr><tr><td>SLBFGS</td><td>Stochastic L-BFGS method</td><td>Section 3.2</td></tr><tr><td>TD</td><td>Time domain</td><td>Section 2.3</td></tr><tr><td>TR</td><td>Time reversal</td><td>Section 3.1</td></tr></tbody></table>

In section 2 we summarize the mathematical modeling of UST and FWI, and describe the state of the art. Section 3 illustrates several separate computational challenges of FWI for high-resolution 3D UST, and describes how novel and established solutions can be combined into a comprehensive inversion strategy. Extensive numerical studies in section 4 first validate solutions to the separate sub-problems before the results of the whole scheme are demonstrated in section 5. Finally, we discuss the results of our work and point to future directions of research in section 6 before closing with conclusions in section 7. Table 1 lists all commonly occurring abbreviations for reference.

## 2. Full waveform inversion for ultrasound tomography

This section presents the model used to describe ultrasound propagation in soft tissue (section 2.1) and the mathematical formulation of the inverse problem of UST (section 2.2). The FWI approach to tackling this problem, the progress that has been made, and the challenges that remain, are then described (sections 2.3 and 2.4).

### 2.1. Modelling ultrasound propagation in soft biological tissues

It is common to model soft biological tissues as compressible fluids and describe acoustic propagation therein by linearising the equations of fluid dynamics derived from conservation laws and combining them with an empirical equation of state. Following this approach, a sufficiently low amplitude ultrasonic wave can be modelled by the following system of first order equations [82]:

$$ \frac{\partial}{\partial t} v(x, t) = - \frac{1}{\rho_0(x)} \nabla p(x, t) \qquad (\text{momentum conservation}) \tag{1} $$

$$ \frac{\partial}{\partial t} \rho(x, t) = -\rho_0(x) \nabla \cdot v(x, t) - v(x, t) \nabla \rho_0(x) + m(x, t) \quad (\text{mass conservation}) \tag{2} $$

$$ p(x, t) = c_0^2(x) (\rho(x, t) + d(x, t) \cdot \nabla \rho_0(x) - L\rho(x, t)), \quad (\text{pressure-density relation}) \tag{3} $$

where $v$ is the acoustic particle velocity (the time derivative of the acoustic displacement $d$), $p$ and $\rho$ are the acoustic pressure and acoustic density fluctuations, $\rho_0$ and $c_0$ are the ambient density and sound speed, and $m$ is a mass source term. $L$ models acoustic absorption and
dispersion using fractional Laplacians:

$$L = \tau(x) \frac{\partial}{\partial t} (-\nabla^2)^{\frac{y}{2}-1} + \eta(x) (-\nabla^2)^{\frac{y+1}{2}-1}, \quad (4)$$

$$\tau(x) = -2\alpha_0(x)c_0(x)^{y-1}, \quad \eta(x) = 2\alpha_0(x)c_0(x)^y \tan(\pi y/2), \quad (5)$$

where $\tau$ and $\eta$ denote absorption and dispersion proportionality coefficients, $\alpha_0$ is the power law prefactor, and $y$ is the power law exponent [81]. Equations (1)–(3) can be combined into the following lossy second-order wave equation for a heterogeneous medium:

$$A p(x, t) := \left( \frac{1}{c_0^2(x)} \frac{\partial^2}{\partial t^2} - \rho_0(x) \nabla \cdot \left( \frac{1}{\rho_0(x)} \nabla \right) + L \nabla^2 \right) p(x, t) \\ = \frac{\partial}{\partial t} m(x, t) := s(x, t), \quad (6)$$

where we introduced a lossy d'Alembert operator $A$. As initial conditions we have $p = 0$ and $\partial_t p = 0$. Suitable data pre-processing typically allows us to model the propagation as if it occurs in an unbounded domain, i.e., no explicit boundary conditions are required (reflections from experimental equipment can be time-gated out and waves propagating deeper into the body are absorbed). See [7, 80, 82] for more detailed discussions on the mathematical modeling and [10, 43, 79] for recent extensions of FWI to UST involving hard tissues such as bone.

## 2.2. The inverse problem of ultrasound tomography

In most UST systems, movable arrays of US transducers are placed around the sample $\Gamma \subset \mathbb{R}^3$ and acoustically coupled with it. For a single recording, a subset of transducers emit an acoustic pressure signal of length $T_s$ while the others are receiving for a time $T > T_s$, chosen long enough such that the acoustic pressure remaining inside the scanner is not measurable anymore. A whole scan consists of $n_s$ such recordings, between which the transducer arrays may be moved. Here, we model a UST scan in the following way: we have $n_s$ temporal sources $s_i(x, t)$, with $\text{supp}(s_i) \subset \Gamma^c \times [0, T_s]$ (in practice, each source is only defined over a small region corresponding to the front face of the transducer). The measurement process is modelled by applying a linear operator $M_i$ to the pressure fields $p_i(x, t), t \in [0, T]$ resulting from (6), i.e., we have

$$A(u) p_i(x, t) = s_i(x, t), \quad f_i = M_i p_i, \quad i = 1, \dots, n_s. \quad (7)$$

Here, $u$ denotes the unknown acoustic material properties inside the sample $\Gamma$, in the extreme case $u = (c_0, \rho_0, \alpha_0, y)$, where $c_0, \rho_0$ and $\alpha_0$ depend on $x$. We assume that all other properties are sufficiently well-known to ignore their modeling error. The general inverse problem of UST is to recover $u$ (or features of it) given a noisy data set $\{f_i^\delta\}_{i=1}^{n_s}$. A general discussion of the uniqueness and stability of this problem can be found in [78]. As mentioned in the introduction, many different approaches have been used to tackle this problem.

## 2.3. Full waveform inversion

In FWI, we assume that for a given $u$, we can solve (7) to simulate data, i.e., $f_i(u) := M_i A(u)^{-1} s_i$. Then, we try to optimize $u$ such that the discrepancies between simulated and measured data
become small:

$$
\begin{align}
\min_{u \in \mathcal{U}} \mathcal{J}(u), \quad \text{with} \quad \mathcal{J}(u) &:= \sum_{i}^{n_s} \mathcal{D}_i(u) := \sum_{i}^{n_s} \mathcal{D}(f_i(u), f_i^{\delta}) \nonumber \\
&= \sum_{i}^{n_s} \mathcal{D}(M_i A^{-1}(u) s_i, f_i^{\delta}), \tag{8}
\end{align}
$$

where $\mathcal{D}(f,g)$ is a loss function (see [25, 98] for a discussion of suitable loss functions), and $\mathcal{U}$ is a set of constraints on $u$, e.g., bound constraints. First-order optimization schemes solve (8) using only the gradient $\nabla \mathcal{J}(u)$, which is given by the sum over terms of the form $\nabla_u \mathcal{D}(MA^{-1}(u)s, f^{\delta})$. Such expressions can be computed efficiently using the adjoint state method [71], which we will summarize here in a short, instructive way. One starts by differentiating both sides of (7) w.r.t. $u$:

$$
\begin{align}
\frac{\partial A}{\partial u} p + A \frac{\partial p}{\partial u} = 0 \quad &\Rightarrow \quad \frac{\partial p}{\partial u} = -A^{-1} \frac{\partial A}{\partial u} p \quad &\Rightarrow \quad \frac{\partial f}{\partial u} = -MA^{-1} \frac{\partial A}{\partial u} p \nonumber \\
&\Rightarrow \frac{\partial \mathcal{D}}{\partial u} = \left(\frac{\partial f}{\partial u}\right)^T \frac{\partial \mathcal{D}}{\partial f} = -\left(\frac{\partial A}{\partial u} p\right)^T A^{-T} M^T \frac{\partial \mathcal{D}}{\partial f}. \tag{9}
\end{align}
$$

Here, $A^{-T}$ is the operator solving the adjoint wave equation [71]. For instance, for $u = c_0$ and $L = 0$, we get with $\frac{\partial A}{\partial c_0} = -\frac{2}{c_0^2} \frac{\partial^2}{\partial t^2}$:

$$
\nabla_{c_0} \mathcal{D}(f(c_0), f^{\delta}) = \int_0^T \frac{2}{c_0(x)^3} \left( \frac{\partial^2 p(x, t)}{\partial t^2} \right) q^*(x, T-t) \mathrm{d}t, \quad (10)
$$

where $q^*$ is the adjoint pressure field obtained by solving $A(u)q^* = s^*$, and $s^*(x,t)$ is the time-reversed data discrepancy gradient $\frac{\partial D}{\partial f}$ [71]. We will use $\mathcal{D}(f,g) = \frac{1}{2}\|f-g\|^2_2$, for which we have $\frac{\partial D}{\partial f} = f-g$. Corresponding formulas for $\rho_0$, $\alpha_0$ and $y$ are listed in appendix A.

The focus of this work is the practical feasibility of computing a sufficiently accurate and
spatially resolved approximation to the solution of the core problem (8) for 3D breast UST
within a reasonable amount of time and with reasonable computing resources. For this rea-
son, we will only include bound constraints on u via U in the main studies. To embed more
sophisticated a-priori knowledge on u and account for noise models and model discrepancies
in f^\delta, one needs to extend (8) by adding regularization functionals (e.g., to penalize unwanted
spatial features of u) or additional constraints (see [26, 67] and references therein). The addi-
tional noise stability study in appendix F will demonstrate the use of a simple regularization
strategy.

We introduced FWI in the *time domain* (TD) here, which is typically also the domain the measurements are obtained in. One can also formulate FWI in the *frequency domain* (FD). Then, computing $\mathcal{J}(u)$ and its gradient requires solving a direct and an adjoint Helmholtz equation for each frequency. If the emitted waves are narrow band and/or the transducers are only sensitive within a narrow band (modelled by $M_i$), this can be very advantageous but typically leads to high memory requirements if high resolution in 3D is desired (see [73] and references therein).

Depending on the scanning geometry, transducer characteristics, and measurement protocol and whether TD-FWI or FD-FWI is pursued, different numerical schemes are advantageous for solving the wave or Helmholtz equation. Direct methods such as finite-difference, pseudospectral, finite/spectral element methods or discontinous Galerkin methods aim to solve the PDEs
directly while integral equation methods such as boundary element methods first transform them. Asymptotic or approximate methods such as geometrical optics or Gaussian beams only capture a part of the wave-matter interaction while typically being computationally much more efficient. See [33, 50] for overviews on computational wave propagation methods for seismic FWI imaging.

## 2.4. Challenges of high resolution 3D time domain full waveform inversion

Many UST systems are designed to capture 3D information and emit unfocused, broadband waves to obtain sub-millimeter isotropic spatial resolutions over the whole breast volume [2, 36]. For instance, within the PAMMOTH project [2], we built a scanner with a bowl-shaped transducer array and aim to reach a resolution of $\le 0.5$ mm. Even when using the adjoint state method, solving TD-FWI (8) for such scenarios is computationally challenging: for the PAMMOTH setup, even an efficient $k$-space pseudospectral method [80] takes at least 10 min to solve a single wave simulation on a recent GPU. As $\mathcal{J}(u)$ involves the sum over $n_s$ sources (8), computing $\nabla \mathcal{J}(u)$ requires $n_s$ parallel gradient computations (10) each involving two wave simulations. On a single GPU, this would take 14 days for $n_s = 1024$. While using a cluster with many computational nodes could reduce this [13, 40], a closer look at (10) reveals a problem of the TD adjoint state method: ideally, the field $p(x, t)$ needs to be kept in GPU memory or at least in the main memory of the computational node while the adjoint field $q^*(x, t)$ is computed. For the PAMMOTH setup, this would require up to 150 GB working memory. Current GPU cards are limited to 48 GB and typical nodes in computational clusters are not equipped with that much main memory, either. In addition to these difficulties, first-order optimization methods for FWI (8) suffer from slow convergence and may get stuck in suboptimal local minima of the non-convex function $\mathcal{J}(u)$.

In the following section, we describe techniques to circumvent each of these problems and combine them into a comprehensive computational strategy to achieve high resolution 3D TD-FWI for ultrasonic breast imaging.

# 3. Improving memory footprint, efficiency and convergence

## 3.1. Memory-efficient gradient computation using time reversal

First, we describe an approach to circumvent storing $p(x, t)$, $x \in \Gamma$, $t \in [0, T]$. Note that in UST, all sources and sensors are outside $\Gamma$ and $p(x, 0) = \partial_t p(x, 0) = 0$. We first consider the case of no absorption, i.e., $L = 0$. By Huygens' principle, the field inside $\Gamma$ is uniquely determined by its trace on the boundary $\partial\Gamma \times [0, \infty]$. The theory of time reversal (TR) describes under which conditions the time reversal symmetry of the wave equation (6) allows $p(x, t)$ to be reconstructed by simulating the following wave equation for the time-reversed pressure $p^<(x, t)$:

$$ 
\begin{align} 
A(u)p^<(x,t) &= 0 && \text{in } \Gamma \times [0, T^<], \nonumber \\ 
p^<(x,t) &= p(x, T^< - t) && \text{on } \partial\Gamma \times [0, T], \tag{11} \\ 
p^<(x, 0) &= \partial_t p^<(x, 0) = 0 && \text{in } \Gamma. \nonumber 
\end{align} 
$$

Namely, $p^<(x, t)$ is a good approximation of $p(x, T^< - t)$ if $\Gamma$ is convex, $u$ is non-trapping and $T^<)$ is chosen large enough such that $p(x, T^<)$ and $\partial_t p(x, T^<)$ are small and monotonically decaying inside $\Gamma$ [18, 34, 57, 97]. In 3D UST of the breast, these conditions are fulfilled: $\Gamma$ can be chosen as a convex set including the breast (e.g., in the PAMMOTH system, the breast is
placed in a plastic cup with convex shape). The acoustic parameter variations are such that no waves get trapped in the sample and after all acoustic energy has entered the breast, it decays exponentially fast, see appendix E. Choosing $T^\triangle = T$ is a safe choice as we chose $T$ large enough to ensure the pressure inside the scanner is not measurable any more (cf section 2.2). Note that (11) is not the same as the adjoint wave equation [7].

TR was developed for focusing ultrasound waves through inhomogenous media [18, 34, 97] and is used for image reconstruction in photoacoustic tomography (PAT) [9, 57]. Here, we use it as a numerical trick to approximately replay the forward field $p(x, t)$ backwards in time, in parallel to solving the adjoint wave equation. This way, we do not need to keep the pressure field $p(x, t)$ in the sample's volume $\Gamma$ in memory, only its trace on the boundary $\partial\Gamma$, which reduces the memory footprint considerably. Note that similar techniques have been developed in the seismic literature. See [68] and references therein for the idea to shift the storage onto the boundary for our type of boundary conditions.

The integral in (10) can directly be accumulated during the parallel time stepping scheme solving the both TR and adjoint wave equations: the state variables in each computation are exactly what is needed to compute the contribution to the integral. For instance, the first time step of the TR wave simulation results in $p^\triangle(x, 0)$, which is approximately $p(x, T)$, while the first time step of the adjoint wave simulation results in $q(x, 0)$. With these variables, the contribution to the integral (10) for time $t = T$ can be computed. The price we pay for the heavily reduced memory footprint is that we have to run three instead of two wave simulations to compute (10). However, the simulation domain of the TR wave simulation only needs to enclose $\Gamma$ and can therefore often be chosen smaller. Appendix B describes the implementation of the TR-based gradient computation using the $k$-space pseudospectral method implemented in the k-Wave Matlab toolbox [80] and section 4.3 validates it numerically.

In the case of absorption, i.e., $L \neq 0$ [85], describes how to modify TR to recover the pressure fields (essentially, the sign of the absorption needs to be switched). Note however, that compared to the application in PAT image reconstruction considered in [85], no additional regularization is needed in our case as the boundary time series are not contaminated with measurement noise.

### 3.2. Tuneable stochastic gradient optimization using delayed source encoding

FWI (8) is a finite sum minimization problem. As such, one can approximate $\nabla \mathcal{J}(u)$ by

$$g_S(u) := \frac{n_s}{|S|} \sum_{j \in S} \nabla \mathcal{D}_j(u), \quad S \subset \{1, \dots, n_s\}. \qquad (12)$$

In incremental gradient or ordered sub-set methods [11], one chooses a different subset $S_k$ in each iteration $k$ in a predetermined way. If the subsets are chosen at random, $g_S(u)$ becomes a stochastic estimator of $\nabla \mathcal{J}(u)$ and stochastic gradient descent (SGD) schemes need to be used. Fortunately, supervised training of deep neural networks are also instances of finite sum minimization problems (empirical risk minimization) and the recent success of deep learning techniques has stimulated research into efficient SGD schemes [14]. In section 4.4, we will demonstrate the benefits of using a recently developed stochastic L-BFGS (SLBFGS) method (see [27] and appendix D) for UST over earlier schemes.

Besides being unbiased, i.e., $\mathbb{E}[g(u)] = \nabla \mathcal{J}(u)$, a desirable property of a gradient estimator $g(u)$ is that it is both computationally efficient (fast to compute) and stochastically efficient (small variance/error). Computing $g_S(u)$ with a subset of size one takes only two wave simulations (or three if TR is used). It turns out that due to the structure of (8), there are gradient estimators with the same computational but higher stochastic efficiency.
In this work, we assume $M_i := M$. Let $w_i$, $i = 1, \dots, n_s$ be random weights with $\mathbb{E}[w] = 0$, $\mathrm{Cov}[w] = I$ and define the synthetic source activations $\hat{s} := \sum_{i}^{n_s} w_i s_i$ and corresponding synthetic data $\hat{f}^{\delta} := \sum_{i}^{n_s} w_i f_i^{\delta}$. Then, $g_{\varepsilon}(u) := \nabla_u \mathcal{D}(MA(u)^{-1}\hat{s}, \hat{f}^{\delta})$ is an unbiased estimator of $\nabla \mathcal{J}(u)$. Practically, we run two 'super-wave' simulations where all sources are fired simultaneously. This technique exploits the linearity of the wave equation and is called source encoding, see [45] for the general theory, [92] for the form we use it here (but we apply it in 3D), and appendix E for a more detailed description. Choosing $w_i = \pm 1$ with equal probability (Rademacher distribution) minimizes the variance of $g_{\varepsilon}$ and will be used from now on (note that one could also define $w$ such that $g_{\varepsilon}$ is recovered). The case $M_i \neq M$ will be examined in forthcoming work that applies the techniques discussed here to concrete experimental data scenarios.

While SGD methods need to reduce the variance of the update progressively in order to converge at all, tailored variance reduction schemes can improve the convergence rate to that of deterministic methods (cf section 5 in [14]). If the variance of the gradient estimator can not be influenced in a direct way, variance reduction can only be achieved via reducing the step size or averaging iterates or gradients. Having fine control over the variance of $g(u)$ is more advantageous. For $g_{\varepsilon}(u)$, one can only increase the size of $\mathcal{S}$; for $g_{\varepsilon}$, one would average multiple realizations. In both cases, we would have a sudden increase of the computational effort which means a rather coarse control over the variance [88]. If multiple GPUs are available, one can perform these computations in parallel, which we will investigate in section 4.5. However, if different GPU architectures are used in this way, the slowest one limits the overall efficiency, cf section 4.6. To overcome these problems, we propose a novel source encoding scheme $g_{\varepsilon}^{\tau}$ with fine variance control that also exploits the time-invariance of the wave equation and the fact that all waves eventually leave the region of interest (which is essential for the sequential scanning of UST, cf section 2.2): let $d_i \sim U[0, 1]$ be random, $\tau \ge 0$ a fixed maximal delay and define

$$ \hat{s}(x, t) := \sum_{i}^{n_s} w_i s_i(x, t - d_i \tau), \quad \hat{f}^{\delta}(x, t) := \sum_{i}^{n_s} w_i f_i^{\delta}(x, t - d_i \tau). \quad (13) $$

For $\tau = 0$, we recover the conventional source encoding where all sources are fired simultaneously while for $\tau > 0$, we activate the sources with random delays. While this necessitates running the super-wave simulation for a longer time $T + \tau$, it dampens the cross-talk between the pressure fields coming from different sources, which is the origin of the error of conventional source encoding. For $n_s = 2$ and $\tau = 2T$, the average delay between the two sources is $T$ and we essentially fire sources with the delay we deemed necessary for regarding the measurements as completely separate (cf section 2.2). More generally, $g_{\varepsilon}^{\tau}$ converges to $\nabla \mathcal{J}(u)$ for growing $\tau$ for any realization of Rademacher weights $w$ and $d$. More details can be found in appendix E. Note that our approach extends the 2D FWI random time-delay approach described in [99] by combining it with the rigorous stochastic analysis from [45]. Recently, another interesting TD approach to reduce the cross-talk between different sources by running longer 'super-wave' simulations has been presented in [8, 87]: the encoding assigns a monochromatic time course to each source $i$, chosen from a set of regularly spaced frequencies. The forward and adjoint wave simulations are run long enough to produce the steady-state pressure fields from which the contributions of each source can be decoded by Fourier decomposition.

In section 4.4, we will validate the proposed gradient estimator numerically. On a single GPU, this technique allows one to increase the gradient estimator's precision in arbitrarily small amounts while on heterogeneous multi-GPU environments, it allows faster GPUs
to compute more precise gradient estimates using tailored $\tau > 0$ while the slowest GPUs use $\tau = 0$.

### 3.3. Improving convergence using coarse-to-fine schemes and preconditioning

FWI (8) is a non-convex optimization problem. Iterative optimization techniques based on local descent directions like gradient descent schemes may experience slow convergence or even get stuck in local minima that only explain parts of the data. In seismic imaging, the latter is known as *cycle skipping* [25, 86, 98]. A common way to avoid it in FD-FWI schemes is to restrict the frequency range to the lower frequencies first and increase it progressively [94]. In TD schemes, this corresponds to solving the FWI on a coarse spatio-temporal grid first, which is then progressively refined. In the most basic case, the interpolation of the solution computed on the coarse grid is used to initialize the optimization scheme on the finer grid but more sophisticated multi-grid schemes can be derived [51]. In the $k$-space pseudospectral method used here, all grids are regular and implementing correct up and down-sampling is easy, see appendix C. In 3D, a coarsening of the spatio-temporal grid by factor of $\eta = 2$ (i.e., $dx^c = \eta dx^f$, $dt^c = \eta dt^f$, ...) leads to a reduction of the computational operations by a factor of $\eta^4 = 16$. We will illustrate in section 4.7 that it is essential to exploit this advantage to obtain fast computational schemes.

Preconditioning techniques try to improve the convergence of numerical optimization schemes by reformulating the original problem. For instance, instead of using $u = c_0^2(x)$ in (8), we could solve for $u = 1/c_0(x)$ (*slowness*), $u = 1/c_0^2(x)$ (*squared slowness*) or a weighted version thereof, $u = w(x)/c_0^2(x)$. The weights $w(x)$ could be a function of the average distance of $x$ to the transducer locations. Due to the non-convex nature of $\mathcal{J}(u)$ and the way SLBFGS constructs low-rank quadratic approximations to $\mathcal{J}(u)$ based on stochastic function and gradient evaluations, this choice is not trivial. In section 4.8, we will examine different possibilities.

## 4. Numerical accuracy and efficiency studies

### 4.1. Setup: 3D digital phantom and transducer arrangement

The setup used in the numerical validation studies resembles the PAMMOTH scanner and is illustrated in figure 1. The computational domain is 221 × 221 × 111 mm wide and the ultrasound transducers are located on a half-spherical surface of radius 104.5 mm. We want to compute a FWI reconstruction at spatial resolution $dx = 0.5$ mm, which leads to a computational domain of 442 × 442 × 222 voxels. The speed of sound bounds enforced during the optimization are $c_0^{\min} = 1350$ m s⁻¹ and $c_0^{\max} = 1800$ m s⁻¹. Based on this and using a cfl number of 0.3, we set $\Delta t = \frac{cfl \cdot \Delta x}{c_0^{\max}} = 83.3$ ns. We set $T$ to twice the length of the computational domain divided by $c_0^{\min}$, which leads to $T/\Delta t = 3912$ time steps.

The anatomically realistic numerical breast phantom used is part of the OA-breast database<fnref n="5" />, which was designed for (photo-)acoustic simulations and is described in [60]. It is based on the tissue segmentation of a clinical contrast enhanced MRI of a breast in prone, free-hanging position. The original segmentation is fitted and interpolated into our setup. The tissue-dependent speed of sound values we use are the same as in [60] and can be found in the caption of figure 1. The mass density $\rho_0$ was assumed constant at 1000 kg m⁻³ for simplicity here.

<footnote n="5">⁵ https://anastasio.bioengineering.illinois.edu/downloadable-content/oa-breast-database/.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 1.** (a) Illustration of the setup used in the numerical studies: the surfaces of different breast tissues have been extracted from a high-res version of the voxel-grid-based breast phantom and rendered with different colors and transparencies. Shown here are blood vessels (red, in-transparent), fibro-glandular tissue (blue, transparent) and skin (yellow, transparent). Source and sensor positions are indicated by the blue and orange dots distributed on a half-sphere, respectively. (b) Illustration of the central x-slice of the 3D volume (corresponding to the green line in (a)). Vertical and horizontal axis correspond to y and z axis, respectively. Different materials are colour coded depending on their speed of sound values ($c_0$), ranging from fat (1470 m s⁻¹) over water/background (1500 m s⁻¹), fibro-glandular tissue (1515 m s⁻¹), blood vessels (1584 m s⁻¹) to skin (1650 m s⁻¹). The pink pixels indicate the half-sphere on which the transducers are located, the reconstruction region $\Gamma$ is composed by all pixels containing breast tissues.</figcaption></figure>

We distribute 2048 transducer locations over the half-sphere by the golden section method. Half of them will be emitting (sources, $n_s = 1024$) while the other half will be receiving (sensors). Simulated data of this artificial source-sensor setup is an array of size $1024 \times 1024 \times 3912$ which corresponds to 15.28 GB in single precision. We assume idealized spatio-temporal transducer characteristics here, i.e., we model them as point-like in space and their temporal impulse response is a delta function. The sources emit a single broadband pressure pulse ('click') smoothed in time to remove temporal frequencies not supported by the computational grid (the grid supports up to 1.5 MHz at $c_0 = 1500 \text{ m s}^{-1}$ [80]). More details can be found in appendix C.

Note that the setup we use here was primarily designed to validate the accuracy and efficiency of the computational solutions we described in the previous section in the best possible way. For this reason, it reproduces the key properties of 3D FWI for breast imaging but omits modeling features of real-world systems that would complicate the interpretation of our results. It is furthermore not designed to examine the condition of the inverse problem (7).

## 4.2. *Baseline reconstructions*

Before we can examine the different aspects discussed in section 3 in depth, we need to illustrate some basic features of the reconstructions. Figure 2(a) shows the speed of sound reconstruction of the data simulated at 0.5 mm resolution when reconstructed using a grid size
<figure><img src="image_6.png" /><figcaption>**Figure 2.** (a) Zoom into central x-slice of the FWI solution computed using a grid size of 2 mm, $c_0^{\text{ini}}$. Colormap in m s⁻¹. All reconstructions will be displayed in this way. (b) Error $c_0^{\text{ini}} - c_0^{\dagger}$ displayed with bidirectional colour map in m s⁻¹.</figcaption></figure>

of 2 mm. This reconstruction was computed using SLBFGS with source encoding gradient estimator $g_E$ for 128 iterations (computing time ~30 min) and could be the multi-grid based initialization of a FWI run on the finer resolution. For this reason, we will call it $c_0^{\text{ini}}$ and take it as a reference point in the following studies. Note that here and in the following, we only reconstruct the speed of sound inside the breast region $\Gamma$, while the speed of sound outside is known and fixed. Figure 2(b) plots the reconstruction error $c_0^{\text{ini}} - c_0^{\dagger}$ locally, where $c_0^{\dagger}$ is the ground-truth speed of sound inside the breast. We will measure the global error of a variable $a$ by their relative $\ell_2$ distances to a reference $b$ computed as $\text{rel } \ell_2(a, b) = \|a - b\|_2/\|b\|_2$, and expressed in percent, e.g., we have that $\text{rel } \ell_2(c_0^{\text{ini}}, c_0^{\dagger}) = 2.18\%$.

First, we run SLBFGS initialized at $c_0^{\text{ini}}$ with the source encoding gradient estimator $g_E$. Figure 3 shows the iterates after 32/64/128 gradient evaluations (eval), the corresponding $\text{rel } \ell_2$ errors are 1.53%/1.29%/1.07%. Note that we deliberately chose to not add measurement noise to the data (see above), the noisy appearance of the reconstructions are due errors of the gradient estimator (*gradient noise*), a purely numerical phenomena that we need to highlight. On the difference plots in figures 3(d)–(f), one can see that the error is larger around tissue interfaces and close to the chest wall. The former comes from the finite spatial grid spacing: the *k*-space pseudospectral method cannot propagate the highest spatial frequencies present for a given spatial grid which means that we can at best hope to reconstruct a slightly smoothed version of $c_0^{\dagger}$. The higher errors close to the chest wall are caused by the insufficient sensor coverage of the half-spherical sensor array (cf figure 1). To illustrate this effect in more detail, we computed error statistics for each depth slice. For slices deeper than 2 cm into the bowl, the errors are almost normally distributed, with zero mean and a standard deviation that decays as the SLBFGS advances as can be seen in figure 4(a). We gain an accuracy of ~5 m s⁻¹ with each doubling of number of iterations, i.e., the computational effort.
<figure><img src="image_7.png" /><figcaption>**Figure 3.** (a)–(c) $c_0^{\text{rec}}$ (m s⁻¹) computed by SLBFGS initialized at $c_0^{\text{init}}$ with the source encoding gradient estimator $g_{\varepsilon}$. (d)–(f) Corresponding error plots (m s⁻¹).</figcaption></figure>

### 4.3. Accuracy of time reversal based gradient computation

First, we validate the TR based gradient computation described in section 3.1 when implemented using a k-space pseudospectral method (cf appendix B). In principle, this schemes relies on regular spatial grids and is not well suited for representing arbitrarily shaped boundaries $∂Γ$, cf (11). To record or impose pressure values on $∂Γ$, we need to instead record and impose them on a layer of grid points lying on the inner surface of $Γ$. Determining the thickness of this layer is a trade-off: a thin layer leads to a small memory footprint when storing $p(x, t)$ on them, which was our original motivation to introduce the TR-based gradient computation. A large layer will lead to a more accurate reproduction of the forward field, i.e., $p^{\triangle}(x, t)$ is a better approximation of $p(x, T^{\triangle} - t)$. To examine the effects of this, we first compute the gradient (10) at $c_0^{\text{init}}$ for a single source storing the complete field $p(x, t)$ in $Γ \times [0, T]$. Then, we compute it using the TR formulation using boundary layers of thickness 1/2/4/8 voxels. Table 2 shows the average relative errors of the TR gradients computed over 16 randomly picked sources and the working memory requirements in our setup. The errors decay very fast as a function of boundary layer thickness. However, we need to check that the errors in the TR gradients do not accumulate during an iterative reconstruction: table 2 shows the reconstruction errors when running SLBFGS with 32 gradient evaluations initialized at $c_0^{\text{init}}$. Interestingly,
<figure><img src="image_8.png" /><figcaption>**Figure 4.** (a) Mean and standard deviation (std) of the reconstruction error $c_0^{\text{rec}} - c_0^{\dagger}$ vs depth (=distance from chest wall) for the reconstructions shown in figure 3. (b) Difference between reconstruction using full gradient computation and TR based gradient with voxel layer of size 1.</figcaption></figure>

<table><caption>**Table 2.** Comparison of gradient computation storing the full forward field (full) with TR based gradient computation using different number of boundary layer voxels.</caption><thead><tr><th></th><th>full</th><th>TR 1</th><th>TR 2</th><th>TR 4</th><th>TR 8</th></tr></thead><tbody><tr><td>rel ℓ<sub>2</sub> gradients at c<sub>0</sub><sup>ini</sup></td><td>0</td><td>22.8085%</td><td>5.6131%</td><td>0.8953%</td><td>0.1376%</td></tr><tr><td>Memory requirements</td><td>98.00 GB</td><td>2.57 GB</td><td>5.95 GB</td><td>11.34 GB</td><td>21.75 GB</td></tr><tr><td>rel ℓ<sub>2</sub>(c<sub>0</sub><sup>rec</sup>, c<sub>0</sub><sup>†</sup>)</td><td>1.5328%</td><td>1.5492%</td><td>1.5401%</td><td>1.5319%</td><td>1.5328%</td></tr><tr><td>rel ℓ<sub>2</sub>(c<sub>0</sub><sup>rec,TR</sup>, c<sub>0</sub><sup>rec,full</sup>)</td><td>0%</td><td>0.2389%</td><td>0.0681%</td><td>0.0178%</td><td>0.0038%</td></tr></tbody></table>

even using a boundary layer of just one voxel, the error between reconstructions with normal gradient and TR-based gradient computations are small. Furthermore, figure 4(b) shows that the error appears noise-like with no clear pattern visible. This suggests that the errors rather cancel than accumulate over the course of the iteration. The difference in computational time will be examined in section 4.6. For the rest of the studies, a TR based gradient computation will be used with a boundary layer of eight voxels.

### 4.4. Accuracy of stochastic gradient optimization

First, we demonstrate the impact the choice of the stochastic gradient method has and compare SLBFGS to plain SGD and SGD with inertia/momentum [14] in figure 5. All methods use source encoding as a gradient estimator and a maximum of 100 gradient evaluations. Furthermore, once an increase in the estimated energy is detected, they do not return the iterates $u^k$ directly but a weighted average $\bar{u}^k := (\sum_l^k l^3 u^l)/\sum_l^k l^3$ to reduce variance, cf appendix D. SLBFGS reaches the final rel $\ell_2$ error of plain SGD already after 55 evaluations, so twice as fast. This number increases to 76 when we add inertia of $\beta = 0.5$ [14] to SGD. However, we
<figure><img src="image_9.png" /><figcaption>**Figure 5.** Comparison of different stochastic gradient methods. (a) Estimated energy $\mathcal{J}(c_0^{\text{rec}})$ is plotted vs the number of gradient evaluations. (b) $\text{rel } \ell_2(c_0^{\text{rec}}, c_0^\dagger)$ is plotted vs the number of gradient evaluations.</figcaption></figure>

see that while SGD with inertia converges as fast as SLBFGS in the beginning, the decay levels off after ~60 evaluations.

Second, we examine the impact of using different gradient estimators: figure 6(a) shows the full gradient computed at $c_0^{\text{ini}}$, which took a week of computing time on a server with 4 GPUs (cf section 4.6). Then, we computed the different stochastic estimates of this gradient, $g_S$, $g_\varepsilon$, $g_\varepsilon^\tau$ (cf section 3.2) and examined how their precision increases by investing more computational effort. Figure 6(b) shows the results averaged over 16 i.i.d. gradient realizations. One can see that using source encoding ($g_\varepsilon$, $g_\varepsilon^\tau$) indeed leads to a more accurate gradient estimation at the same computational effort. We will use source encoding for the rest of the studies. Also, we can see that using the novel delayed source encoding scheme $g_\varepsilon^\tau$ not only gives a finer variance control compared to increasing the precision of $g_\varepsilon$ via averaging, it also provides a more accurate gradient estimate with the same computational effort. Within a stochastic gradient method, the question is whether to do more steps with a less accurate but computationally cheaper gradient estimator vs less steps with a more accurate but computationally more costly estimator. The former is typically more beneficial in the beginning of the iteration, while the latter becomes more beneficial towards the end of it [14]. To illustrate this, we first initialized SLBFGS with $c_0^{\text{ini}}$ and then with the result shown in figure 3(c), which corresponds to 128 gradient evaluations (roughly 4 days computing time on a TITAN RTX GPU). In both cases, we ran it for another 72 h using $g_\varepsilon^\tau$ with $\tau = b(i - 1)$ for different values of $b$ and with and without starting iterate averaging after an increase in the energy estimate in the case $b = 0$ (so no delay). Figure 7(a) shows that in the low accuracy regime, all methods perform rather similar for the first 24 h. In the case $b = 0$, iterate averaging is only activated after 66 h (the point where the first two plots diverge). The delayed source encoding methods $b > 0$ eventually perform better than $b = 0$, but are slower at the start. Figure 7(b) shows the results for the high accuracy regime: while using iterate averaging stabilizes the iteration for $b = 0$, using increasing values of $b$ leads to faster convergence in terms of computing time.
<figure><img src="image_10.png" /><figcaption>**Figure 6.** (a) Gradient $\nabla \mathcal{J}(c_0^{\text{init}})$, normalized to [-1, 1]. (b) Error of different gradient estimators and different ways to increase their precision by investing more computational effort, averaged over 16 independent realizations. A computational effort of 1 corresponds to computing the gradient for a single source.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 7.** Comparison of different delayed source encoding strategies. rel $\ell_2(c_0^{\text{rec}}, c_0^\dagger)$ is plotted vs computing time on a TITAN RTX GPU. All schemes are initialized with (a) $c_0^{\text{init}}$, cf figure 2(a) and (b) the solution shown in figure 3(c).</figcaption></figure>

### 4.5. Multi-GPU acceleration

Computing platforms with multiple GPUs can be used in different ways to accelerate the image reconstruction. In the TR based gradient computation described in section 3.1, the adjoint and TR wave simulation can be run on two separate GPUs in parallel, which would lower the
<figure><img src="image_12.png" /><figcaption>**Figure 8.** Effect of utilizing multiple GPUs to decrease the variance of the SE gradient estimator by averaging independent realizations. (a) rel $ℓ_2(c_0^{\text{rec}}, c_0^{\dagger})$ is plotted vs the number of gradient evaluations. (b) Mean and standard deviation (std) of the reconstruction error $c_0^{\text{rec}} - c_0^{\dagger}$ vs depth, cf figure 4(a).</figcaption></figure>

computational cost to that of the standard gradient computation again. To reconstruct 3D images on higher resolutions, one can implement 3D domain decomposition methods to distribute the computations over several GPUs [69, 83, 84]. A straightforward way which does not need sophisticated implementation is to average statistically independent gradient estimates each computed on a different GPU. This reduces the variance of the gradient estimator at least by a factor corresponding to the number of GPUs, cf section 3.2: figure 8 shows the reconstruction errors of SLBFGS initialized at $c_0^{\text{init}}$ when using a growing number of GPUs in this way (using source encoding without time shifting). One can see that the convergence is in general faster. However, figure 8(a) highlights that due to the non-convexity of $\mathcal{J}(u)$, there is not a simple relationship between the accuracy of the gradient estimator and the convergence curve. The reconstruction error obtained after running SLBFGS using 2/4/8 GPUs for 32 gradient evaluations is the same as running it on a single GPU for 51/89/103 gradient evaluations. Comparing figures 4(a) and 8(b) shows also the decrease in standard deviation saturates with a growing number of GPUs.

## 4.6. Comparison of different GPU architectures

Table 3 compares the performance of different GPUs. One can see that the run times of normal and TR-based gradient computations are more alike than expected based on the number of computational operations (flops). One reason is that the huge memory footprint of the normal gradient computation requires a lot of working memory operations, which are no longer negligible compared to the computational operations.

## 4.7. Acceleration by coarse-to-fine multigrid initialization

Next, we illustrate that a coarse-to-fine initialization strategy can speed up the convergence and avoid local minima as described in section 3.3. We use SLBFGS with source encoding without time shifting and initialize it with the sound speed value of water (1500 m s⁻¹) everywhere. The
<table><caption>**Table 3.** Comparison of different Nvidia GPUs. Shown are the computation times for simulating data generated by a single source (fwd), computing a gradient for it via the standard approach (10) (grad) and the TR based computation (cf section 3.1) (TR grad). N/A indicates that the computer the GPU was installed in did not have insufficient working memory (RAM).</caption><thead><tr><th></th><th>Tesla P40</th><th>TITAN RTX</th><th>GeForce RTX 2080 Ti</th><th>GeForce GTX 1070</th></tr></thead><tbody><tr><th>Memory</th><td>24 GB</td><td>24 GB</td><td>11 GB</td><td>8 GB</td></tr><tr><th>fwd</th><td>19 m 54 s</td><td>10 m 59 s</td><td>11 m 42 s</td><td>30 m 27 s</td></tr><tr><th>grad</th><td>46 m 9 s</td><td>26 m 44 s</td><td>27 m 34 s</td><td>n/a</td></tr><tr><th>TR grad</th><td>61 m 48 s</td><td>34 m 22 s</td><td>36 m 15 s</td><td>1 h 34 m</td></tr></tbody></table>

<figure><img src="image_13.png" /><figcaption>**Figure 9.** (a) Comparison of using multi-grid initializations: all schemes are initialized with a uniform background ($c_0 = 1500 \text{ m s}^{-1}$). A different number of gradient evaluations (evals) are run on each level, the TR-based approach is only used on the finest level ($dx = 0.5 \text{ mm}$). In principle, the vertical axis displays flops which scale as $\eta^4$ (cf section 3.3) where $\eta$ is the coarsening factor between levels. For a better interpretability, we re-scaled flops to GPU hours by using the TR-gradient computation time of the TITAN RTX on the finest level (cf table 3). The black crosses mark where a scheme transitions between levels. The rel $\ell_2$ error of solutions on coarser grids than $dx = 0.5 \text{ mm}$ is measured by interpolating them first. Note that each level transition involves a smoothing which can lead to an increase in rel $\ell_2$ error. (b) rel $\ell_2(c_0^{\text{rec}}, c_0^{\dagger})$ vs the number of gradient evaluations for SLBFGS initialized at $c_0^{\text{ini}}$ and different preconditioning techniques.</figcaption></figure>

first scheme runs on the finest grid ($dx = 0.5 \text{ mm}$) only. The second scheme has two levels: it starts on $dx = 2 \text{ mm}$ before switching to $dx = 0.5 \text{ mm}$, so $\eta = 4$. The third scheme has three levels $dx = 2, 1, 0.5 \text{ mm}$ ($\eta = 2$), and the fourth $dx = 2, 1.41, 1, 0.70, 0.5 \text{ mm}$ ($\eta = \sqrt{2}$). Figure 9(a) shows the error plots vs computational effort and shows how essential it is to use multi-grid schemes in UST: when the two-level scheme switches to the finest level after 23 min, it has already reached a lower error than the single-level scheme will reach after 16 h.
4.8. Acceleration by preconditioning

Figure 9(b) shows that using the non-linear slowness and squared-slowness preconditioning
transformations discussed in section 3.3 do indeed speed up the convergence and provide
clear advantages over the standard parameterization. All spatial error plots clearly indicate
that the errors increase with decreasing depth, i.e., closer to the chest wall where the sensor
coverage is getting worse, cf figures 1–4. This observation suggests that introducing a local
weighting w(x) that depends on depth or sensor coverage could be beneficial. We tried differ-
ent variants of such weightings but could not find significant differences. One explanation for
this is that quasi-Newton schemes like SLBFGS have low sensitivity to linear preconditioning
transformations.

5. Full inversion scheme demonstration

In this section, we combine all the techniques described in the previous section to demonstrate
the results one can obtain within a computing time limit of 24 h. While this limit may seem
arbitrary at this point, it is motivated by the perspective of adding UST to clinical work-flows
for breast-cancer diagnosis and treatment planning: after a positive screening result, several
examinations have to be scheduled and carried out and before a treatment decision can be
made. While health care providers try to shorten the overall time for this as much as possible,
it can rarely be done in a single day [75]. Within this context, a computing time of 24 h would
allow UST to fit into clinical trajectories. Similar considerations hold for the computational
resources needed: a main advantage of UST over MRI could be lower operating costs for a
single examination. Having to use or rent a large computational cluster over several days to
compute the results of a single examination may diminish this advantage.

To appreciate the computational advantages described in this work, we first show the recon-
struction result of applying a computational technique similar to the one used in [72] in
figure 10(a). One can see that even after 24 h the result can not be used to reliably identify
anatomic structures. Then, we use

* TR based gradient computation with a eight voxel boundary layer as examined in section 4.3.

* SLBFGS with source encoding and weighted iterate averaging as examined in section 4.4.

* Multigrid initialization with dx = 2, 1, 0.5 mm (η = 2) and 128/128/32 gradient evaluations as examined in section 4.7.

* Squared slowness preconditioning as examined in section 4.8.

The results of computing on a single GeForce RTX 2080 Ti installed in a conventional desktop PC (total price ~ 5K€) are shown in figures 11(a) and (d). While one can still perceive stochastic gradient noise and structured errors in the reconstruction, it already provides a useful high-resolution estimate of the underlying anatomy. Figures 11(b)/(e) and (c)/(f) show how these results improve using a server with 4 GeForce RTX 2080 Ti (total price ~ 20K€) or a cluster with 16 TITAN RTX (total price ~ 100K€). Figure 10(b) shows the corresponding depth distribution of the errors. At this point, we remind the reader that the images shown only illustrate a sub-set of the whole reconstruction domain, cf figure 1.

As discussed earlier, we did not add measurement noise to the data in our studies to allow for
an easier visual assessment of the convergence of the different computational schemes as the
stochastic gradient error also manifests as noise-like image features. To round off the numerical
studies, appendix F examines the sensitivity to noise.
<figure><img src="image_14.png" /><figcaption>**Figure 10.** (a) Reconstruction of a reconstruction method similar to the one used [72] within 24 h computing time on a single GPU. (b) Mean and standard deviation (std) of the reconstruction error $c_0^{\text{rec}} - c_0^{\text{std}}$ vs depth, cf figure 4(a).</figcaption></figure>

# 6. Discussion and outlook

Our work was motivated by the specifications of the PAMMOTH scanner that is currently being developed and tested [2]. The main contribution was to combine a number of techniques to design a comprehensive computational strategy to realize high resolution 3D TD-FWI for UST of the human breast with moderate computational resources within a day of computing time. The results were validated with extensive numerical proof-of-concept studies on different computational platforms.

**TR-based gradient computation.** For the biggest breast cup designed for PAMMOTH, storing $p(x, t)$ throughout the breast would require 146.81 GB at 0.5 mm spatial resolution. This memory footprint rules out keeping it in the memory of any currently available GPU. If kept in CPU memory, it severely limits the number of parallel source simulations one can perform, e.g., running parallel simulations on a server with 16 GPUs like showcased in section 5 would require more than 2 TB CPU memory. In addition, memory operations start to dominate the total computation time (cf table 3).

The time-reversal-based gradient computation presented in section 3.1 can overcome this bottleneck by shifting the storage from the breast volume to its boundary: with a moderate increase of computation time (cf table 3), the memory requirement can be lowered by a factor 28.79 to as little as 5.10 GB without compromising accuracy of the final solution by more than 0.25%, cf section 4.3. This allows us to utilize multiple GPUs without memory restrictions, which is examined in detail in section 4.5 and demonstrated in section 5, cf figures 8 and 11.

**Stochastic optimization.** The second crucial ingredient of our strategy is to combine novel stochastic optimization techniques with efficient stochastic gradient estimators and preconditioning. In section 4.4, we demonstrate the benefits of SLBFGS over plain SGD and show that source encoding schemes lead to much more accurate gradient estimators at the same computational cost as sub-sampling methods, cf figure 6. Furthermore, we presented a novel
<figure><img src="image_15.png" /><figcaption>**Figure 11.** Demonstration full inversion scheme. Shown are results computed within 24 h using (a) a desktop PC with one GeForce RTX 2080 Ti GPU, (b) a server with 4 GeForce RTX 2080 Ti GPUs, (c) a cluster with 16 TITAN RTX GPUs. (d)–(f) corresponding error plots. While reconstruction (a) still shows some stochastic gradient noise, all main anatomical features are visible and there is no systematic quantitative error.</figcaption></figure>

source encoding scheme that exploits the time-invariance of the wave equation by introduc-
ing random time delays to increase the accuracy of the gradient estimate at the expense of
longer computation times. This offers a much finer variance control than estimate averag-
ing as well as being more accurate, cf figure 6(b). Embedded in SLBFGS, we showed that
using delayed source encoding will be beneficial over iterate averaging in the high accuracy
regime, beyond the fixed computation time budget of 24 h we set ourselves in this work, cf
figure 7. Future work will examine this topic more carefully, e.g., by developing adaptive delay
choice rules and by examining how to balance delays to optimally utilize heterogeneous multi-
GPU computing environments. The encoding scheme presented here relies on the transient
nature of the wave fields to reduce cross-talk. It would be interesting to compare it to the
approach presented in [8], which obtains impressive results by encoding a harmonic steady
state instead. Finally, one can observe that the error of the source encoding gradient estimates
appears to be noise-like (cf figures 3 and 11). As such, it is worth examining whether coupling
source encoding with preconditioning that performs image denoising will further accelerate
the convergence.

**Multigrid schemes.** The results presented in section 4.7 reveal that using coarse-to-fine multi-grid schemes are another essential component of our strategy: even the simple multi-grid initialization technique presented here accelerates the overall convergence massively. In the future, we will examine more sophisticated multi-grid schemes.

**Further extensions.**

* The methods presented here were exclusively implemented in Matlab, using the GPU support offered by the parallel computing toolbox. Implementing it using optimized CUDA code will lead to further performance gains [69, 83, 84].

* While we presented the material in terms of a general acoustic parameter *u*, we only showcased the methods for *u* = *c*₀. The extension to acoustic density *ρ*₀ and attenuation (*α*₀, *y*) is described in appendix A and will be examined numerically in future work.

* In this numerical proof-of-concept study, we used idealized assumptions on the transducer locations and properties. In forthcoming work, we will address realistic transducer and scan modeling.

* Finally, we will demonstrate the application of this work to real data from experimental phantoms, healthy volunteers and breast-cancer patients acquired by the PAMMOTH scanner in the future.

**7. Conclusions**

In this work, we described and evaluated a comprehensive computational strategy for high res-
olution 3D TD-FWI for UST of the human breast. With moderate computational resources,
accurate results can be obtained within less than a day of computation time. The extent of
the improvements presented here can best be appreciated by comparing figure 10(a), which
corresponds to our starting point [72] with figure 11, which showcases our current results.
Similarly promising results were presented in [8], which was developed in parallel to the
work described here. While the computing time and costs of our scheme are still too high
for screening applications, it may well already fit into clinical trajectories for diagnosis and
therapy planing.

**Acknowledgments**

This work was supported by the Netherlands Organisation for Scientific Research (NWO)
Project Number 613.009.106 and by the European Union's Horizon 2020 Research and Inno-
vation Program H2020 ICT 2016–2017 under Grant Agreement No. 732411, which is an
initiative of the Photonics Public Private Partnership.

**Data availability statement**

No new data were created or analysed in this study.
## Appendix A. Extension to density and absorption

This section extends the formulas in section 2.3 and lists the remaining derivatives. First, including absorption with $L \neq 0$, we get

$$ \frac{\partial L}{\partial c_0} = -2\alpha_0(y-1)c_0^y \frac{\partial}{\partial t}(-\nabla^2)^{\frac{y}{2}-1} + 2\alpha_0 y c_0^{y-1} \frac{\partial}{\partial t}(-\nabla^2)^{\frac{y+1}{2}-1} \quad (14) $$

as additional terms for the derivative with respect to $c_0$. See [72] for a discussion of their relative importance. For $u = \rho_0$ we get

$$ \frac{\partial A}{\partial \rho_0} = \left( \nabla \cdot \frac{1}{\rho_0} - (\rho_0 \nabla \cdot) \frac{1}{\rho_0^2} \right) \nabla. \quad (15) $$

For $u = \alpha_0$ we get

$$ \frac{\partial A}{\partial \alpha_0} = \frac{\partial L}{\partial \alpha_0} = -2c_0^{y-1} \frac{\partial}{\partial t}(-\nabla^2)^{\frac{y}{2}-1} + 2c_0^y \tan(\pi y/2) (-\nabla^2)^{\frac{y+1}{2}-1}. \quad (16) $$

For $u = y$, we remark that functions of the self-adjoint operator $-\nabla^2$, are understood in terms of its spectral decomposition, i.e., in Fourier space (cf [81]) and get

$$ \begin{aligned} \frac{\partial A}{\partial y} = \frac{\partial L}{\partial y} ={}& \tau \log(c_0) \frac{\partial}{\partial t}(-\nabla^2)^{\frac{y}{2}-1} + \frac{\tau}{2} \frac{\partial}{\partial t}(-\nabla^2)^{\frac{y}{2}-1} \log(-\nabla^2) \\ & + \alpha_0 c_0^y \left( 2 \log(c_0) \tan\left(\frac{\pi y}{2}\right) + \pi \sec^2\left(\frac{\pi y}{2}\right) \right) (-\nabla^2)^{\frac{y+1}{2}-1} \\ & + \frac{\eta}{2} (-\nabla^2)^{\frac{y+1}{2}-1} \log(-\nabla^2). \end{aligned} \quad (17) $$

## Appendix B. k-space pseudospectral time domain implementation of time-reversal-based gradient computation

Here we sketch how to implement the proposed memory efficient gradient computation (11) using a k-space TD scheme to solve the underlying first order system of equations (1)–(3). It follows the same approach implemented in k-Wave [80]. For a more detailed description, we refer to the k-Wave manual<fnref n="6" />. The k-space TD method is a collocation scheme that interpolates between the collocation points using a truncated Fourier series. This allows gradients to be calculated using the fast Fourier transform (FFT) which leads to fast implementations on GPU architectures. The ‘k-space’ in the name refers to a correction $\kappa$ applied in the Fourier domain to account for the finite difference approximation of the time derivative. This correction is exact for a homogeneous medium [77] and reduces the errors for acoustically heterogeneous media. The spatial domain $\Omega$ has to be embedded into a rectangular region, which is then discretised by a regular grid of N collocation points with grid size $\Delta x$. To mimic free-space propagation, a perfectly matched layer (PML) absorbing boundary is wrapped around the box to damp outgoing waves without reflecting them. In this work, the size of the PML layer is automatically chosen to be between 10–20 voxels by minimizing the highest prime factor of the total grid size (cf function getOptimalPMLSize.m in k-Wave). This will lead to efficient fft executions.

<footnote n="6">⁶ www.k-wave.org.</footnote>
We discretize the measurement time interval $[0, T]$ by $t_n = n\Delta t, n = 0, \dots, N_t, \Delta t = T/N_t$. The discrete pressure is denoted by $p$ and the $\xi$-component of the discrete particle velocity by $v_\xi$ ($\xi$ stands for one of the spatial components $x, y, z$ of these vector fields). As the terms containing $\nabla\rho_0(x)$ in (1)–(3) cancel out, they are omitted [82]. The scalar density $\rho$ is split into three parts, $\rho_\xi$. This is nonphysical but allows anisotropic absorption to be included within the PML. The $k$-space derivative and the $k$-space operator $\kappa$ are defined as

$$ \frac{\partial}{\partial \xi} g := \mathcal{F}^{-1} \{ik_{\xi} \kappa \mathcal{F}\{g\}\}; \quad \kappa = \text{sinc} \left( \frac{c_{\text{ref}} \Delta t}{2} \sqrt{k_x^2 + k_y^2 + k_z^2} \right), \quad (18) $$

where $k_\xi \in \mathbb{R}^N$ is the discrete wavevector in $\xi$ direction and all multiplications between $N$-dimensional vectors are understood as componentwise. $c_{\text{ref}}$ is a reference sound speed, chosen to ensure stability. For consistency during the FWI inversion, we fix $c_{\text{ref}} = c_0^{\text{max}}$. Grid-staggering is incorporated into the calculation of the gradients. This means that a spatial translation of $\pm \Delta \xi/2$ is introduced as

$$ \frac{\partial^{\pm}}{\partial \xi} g := \mathcal{F}^{-1} \{ik_{\xi} \kappa e^{\pm ik_{\xi} \Delta \xi/2} \mathcal{F}\{g\} \}. \quad (19) $$

Staggering in time will be included by interleaving the gradient and updates steps. Multiplication operators $\Lambda_\xi$ and $\Lambda_\xi^s$ implement the PML on the normal and staggered grid, respectively. The measurement operator $M$, i.e., the receiving US transducers, is modeled by a linear mapping between spatial grid and measurement channels and is also denoted by $M$. This assumes that it acts instantaneously, i.e., in iteration $n$, we approximate $f(t_n) = Mp$ (otherwise, the pressure fields at the sensor locations need to be buffered). In the same way, a source operator $S$ will embed the discrete pressure source $s(t_n)$ into the spatial grid. After this, $Ss(t_n)$ is divided by $c_0$ pointwise and scaled by $2/(3\Delta x)$ to convert it into a mass source, cf section 2.4 in the k-Wave manual. To implement the TR-based gradient (11), we need to extract the pressure on the boundary of sample, $\partial\Gamma$, during the forward wave simulation. In our collocation scheme, we will simply extract the pressure in a layer of grid points with a defined thickness as examined in section 4.3, which we will denote as $p(\partial\Gamma)$ for simplicity. The scheme to compute the TR-based gradient consists of three wave simulations, forward, adjoint, and time-reversal of which the latter two are interleaved. For simplicity, we only consider $u = c_0$ here (as in the numerical studies) and any scaling factors like $\Delta t$ have been omitted if they cancel out in the end. Algorithm 1 describes the whole iteration. Step 13 corrects for the fact that the first order scheme effectively uses the time derivative of the mass source term provided, not the term itself, cf (6). Code will be released as part of a more general Matlab toolbox for ultrasonic and photoacoustic image reconstruction.

## Appendix C. Multigrid for $k$-space pseudospectral schemes

In the Fourier collocation scheme described above, the discrete, distributed field parameters $(p, v, \rho)$ implicitly represent band-limited interpolants of the corresponding infinite dimensional variables (cf section 2.8 in the k-Wave manual). As such, Fourier/trigonometric interpolation is the natural way to transfer variables from one spatial grid into another. For regular spatial grids, this can be implemented efficiently using FFTs (cf function interpftn.m in k-Wave), and can be combined with smoothing in Fourier space. To implement the coarse-to-fine initialization strategy described in section 3.3 and examined in section 4.7, we first need to down-sample the temporal dimension of data $f^\delta$ and our model of the corresponding (mass) source $m$ to a coarser temporal grid. This is done via a combination of Fourier interpolation and
<caption>**Algorithm 1.** Time reversal based speed-of-sound gradient computation.</caption>

1: **Input:** *s*
   **Forward iteration**
2: *p*, *v*<sub>ξ</sub>, *ρ*<sub>ξ</sub> ← 0 for all *ξ* ∈ [*x*, *y*, *z*]
3: **for** *n* = 0, . . . , *N*<sub>*t*</sub> **do**
4: &nbsp;&nbsp;&nbsp;&nbsp;*v*<sub>ξ</sub> ← Λ<sub>ξ</sub><sup>s</sup> (Λ<sub>ξ</sub><sup>s</sup>*v*<sub>ξ</sub> − Δ<i>t</i>/ρ<sub>0</sub> ∂<sup>+</sup>/∂ξ *p*)
5: &nbsp;&nbsp;&nbsp;&nbsp;*ρ*<sub>ξ</sub> ← Λ<sub>ξ</sub> (Λ<sub>ξ</sub>*ρ*<sub>ξ</sub> − Δ<i>t</i>˜ρ<sub>0</sub> ∂<sup>-</sup>/∂ξ *v*<sub>ξ</sub>) + <sup>2Δ<i>t</i></sup>⁄<sub>3Δ<sub>x</sub>c<sub>0</sub></sub> *Ss*(*t*<sub>*n*</sub>)
6: &nbsp;&nbsp;&nbsp;&nbsp;*p* ← *c*<sub>0</sub><sup>2</sup> (*ρ*<sub>*x*</sub> + *ρ*<sub>*y*</sub> + *ρ*<sub>*z*</sub>)
7: &nbsp;&nbsp;&nbsp;&nbsp;*f*<sup>*n*</sup> ← *Mp*
8: &nbsp;&nbsp;&nbsp;&nbsp;*g*<sup>*n*</sup> ← *p*(*∂*Γ)
9: **end for**
**setup interleaved adjoint and time reversal iterations**
10: ∇*D*, *p*, *v*<sub>ξ</sub>, *ρ*<sub>ξ</sub>, *p*<sup>⊤</sup>, *v*<sub>ξ</sub><sup>⊤</sup>, *ρ*<sub>ξ</sub><sup>⊤</sup> ← 0 for all *ξ* ∈ [*x*, *y*, *z*]
11: *q*<sup>*n*</sup> ← *f*<sup>*N*<sub>*t*</sub>−*n*</sup> − *f*<sup>*s*</sup>(*T* − *t*<sub>*n*</sub>) for all *n*
12: *D* ← <sup>1</sup>⁄<sub>2</sub> ‖*q*‖<sub>2</sub><sup>2</sup>
13: *q* ← *cumsum*(*q*)
14: *p*<sup>⊤</sup>(*∂*Γ) ← *g*<sup>*N*<sub>*t*</sub></sup>
15: *g*<sup>*n*</sup> ← *g*<sup>*n*</sup>/(3c<sub>0</sub>(*∂*Γ)<sup>2</sup>) for all *n*
16: *p*<sub>rec</sub><sup>0</sup> ← *p*<sup>⊤</sup>(*Γ*)
**Advance time reversal iteration by one step to get ahead of adjoint iteration:**
17: *v*<sub>ξ</sub><sup>⊤</sup> ← Λ<sub>ξ</sub><sup>s</sup> (Λ<sub>ξ</sub><sup>s</sup>*v*<sub>ξ</sub><sup>⊤</sup> − Δ<i>t</i>/ρ<sub>0</sub> ∂<sup>+</sup>/∂ξ *p*<sup>⊤</sup>)
18: *ρ*<sub>ξ</sub><sup>⊤</sup> ← Λ<sub>ξ</sub> (Λ<sub>ξ</sub>*ρ*<sub>ξ</sub><sup>⊤</sup> − Δ<i>t</i>˜ρ<sub>0</sub> ∂<sup>-</sup>/∂ξ *v*<sub>ξ</sub><sup>⊤</sup>)
19: *ρ*<sub>ξ</sub><sup>⊤</sup>(*∂*Γ) ← *g*<sup>*N*<sub>*t*</sub>−1</sup>
20: *p*<sup>⊤</sup> ← *c*<sub>0</sub><sup>2</sup> (*ρ*<sub>*x*</sub><sup>⊤</sup> + *ρ*<sub>*y*</sub><sup>⊤</sup> + *ρ*<sub>*z*</sub><sup>⊤</sup>)
21: *p*<sub>rec</sub><sup>+</sup> ← *p*<sup>⊤</sup>(*Γ*)
**Advance adjoint and time reversal iterations simultaneously:**
22: **for** *n* = 1, . . . , *N*<sub>*t*</sub> − 1 **do**
23: &nbsp;&nbsp;&nbsp;&nbsp;*v*<sub>ξ</sub> ← Λ<sub>ξ</sub><sup>s</sup> (Λ<sub>ξ</sub><sup>s</sup>*v*<sub>ξ</sub> − Δ<i>t</i>/ρ<sub>0</sub> ∂<sup>+</sup>/∂ξ *p*)
24: &nbsp;&nbsp;&nbsp;&nbsp;*ρ*<sub>ξ</sub> ← Λ<sub>ξ</sub> (Λ<sub>ξ</sub>*ρ*<sub>ξ</sub> − Δ<i>t</i>˜ρ<sub>0</sub> ∂<sup>-</sup>/∂ξ *v*<sub>ξ</sub>) + <sup>2Δ<i>t</i></sup>⁄<sub>3Δ<sub>x</sub>c<sub>0</sub></sub> *M*<sup>*T*</sup> *q*<sup>*n*</sup>
**Algorithm 1.** Continue.

25: $p \leftarrow c_0^2 (\rho_x + \rho_y + \rho_z)$

26: $v_{\xi}^{\triangleleft} \leftarrow \Lambda_{\xi}^{s} \left( \Lambda_{\xi}^{s} v_{\xi}^{\triangleleft} - \frac{\Delta t}{\rho_0} \frac{\partial^{+}}{\partial \xi} p^{\triangleleft} \right)$

27: $\rho_{\xi}^{\triangleleft} \leftarrow \Lambda_{\xi} \left( \Lambda_{\xi} \rho_{\xi}^{\triangleleft} - \Delta t \bar{\rho}_0 \frac{\partial^{\leftarrow}}{\partial \xi} v_{\xi}^{\triangleleft} \right)$

28: $\rho_{\xi}^{\triangleleft}(\partial\Gamma) \leftarrow g^{N_1 - n - 1}$

29: $p^{\triangleleft} \leftarrow c_0^2 (\rho_x^{\triangleleft} + \rho_y^{\triangleleft} + \rho_z^{\triangleleft})$

30: $\bar{p}_{\text{rec}} \leftarrow p_{\text{rec}}^0$

31: $p_{\text{rec}}^0 \leftarrow p_{\text{rec}}^+$

32: $p_{\text{rec}}^{+} \leftarrow p^{\triangleleft}(\Gamma)$

33: $\nabla D \leftarrow \nabla D + \frac{2}{\Delta t c_0^2} (p_{\text{rec}}^{+} - 2p_{\text{rec}}^{0} + \bar{p}_{\text{rec}}) p(\Gamma) \quad \triangleright \text{ Assemble gradient, cf (10), } p \text{ corresponds to } q^*$

34: end for

35: return: $D, \nabla D$
<figure><img src="image_16.png" /><figcaption>**Figure 12.** Mass source signals $m(t)$ and corresponding simulated measurements $f(t)$ on different spatio-temporal grids ($\eta = 2$) for a random source-receiver pair. The signals on coarser grids result from Fourier interpolation of the fine grid signals followed by low-pass filtering. The signals have been rescaled for comparison. On each level, the temporal filtering used to remove the highest frequencies not supported by the grid introduces a different phase shift. However, as the shift is the same for source and measured time series, it does not affect the FWI results.</figcaption></figure>

low-pass filtering with a Kaiser window (cf function `filterTimeSeries`.min in k-Wave) to remove frequencies not supported by the numerical scheme. The data additionally needs to be scaled by $\eta^{d-1}$. Figure 12 illustrates the results. After the FWI solution $u^c$ is computed on the coarse grid, we need to up-sample it to initialize the FWI on the next, finer grid. For this, we use a combination of Fourier interpolation and smoothing with a Blackman window (cf section 2.8 in the k-Wave manual): $u^f \leftarrow interpftn(u^c, [N_x^f, N_y^f, N_z^f], 'Blackman')$, where $[N_x^f, N_y^f, N_z^f]$ are the dimensions of the fine spatial grid.

## Appendix D. L-BFGS with stochastic gradient estimates

We want to minimize an energy $\mathcal{J}(u)$ as in (8) while only having access to an unbiased estimator $g(u, \theta, n)$ of $(\mathcal{J}(u), \nabla \mathcal{J}(u))$. Here, $\theta \in [0, 1]$ determines the approximation strength of $g$ where $\theta = 0$ corresponds to no approximation, e.g., $g(u, 0, \cdot) = \mathcal{J}(u)$, and $\theta = 1$ to maximal approximation, e.g., by using a data sub-set $S$ in $g_S$ of size one. In our implementation, $\theta$ can be a function of the iteration count $i$. The last input $n \in \mathbb{N}$ in $g(u, \theta, n)$ controls the stochastic approximation performed by $g$. For instance, it could seed the random generator used to select $S$ in $g_S$. Important for the algorithm is only that it leads to the same type of approximation for different $u$, e.g., for a given $\theta, n, g_S(\cdot, \theta, n)$ always uses the same subset of data $S$.

Algorithm 2 describes the stochastic L-BFGS (SLBFGS) scheme used in this work, which is based on the one presented in [27]. See section 6.2 in [14] and references therein for a detailed discussion of stochastic quasi-Newton schemes. The key idea is that for computing the vectors used to construct the low-rank approximation of the inverse Hessian, $H_k$, only iterate pairs evaluated with the same stochastic approximation should be used. This increases the number of gradient computations to at least two per iteration (more if line-search is used) but the second gradient evaluation can be utilized to compute a second update. In [27], no bound constraints
$u_{\min} \le u \le u_{\max}$ were considered. Here, we include them but assume that they are not active at the vicinity of the solution. This means that they do not stabilize the image reconstruction problem. Rather, they merely safeguard against running the numerical wave simulation scheme with values of the parameter fields $u$ that would lead to instability. For this reason, we enforce them by a simple projection of the updated variable $u^k$ at the end of each iteration (step 23). We use progressive iterate averaging with $w(i) = i^3$, but only start averaging after the estimated energy increases for the first time (step 31). We use history size $m_h = 64$ throughout the paper. Code will be released as part of a more general Matlab toolbox for deterministic and stochastic reconstruction.

## Appendix E. Delayed source encoding for time-invariant systems

Here, we examine the conventional and delayed source encoding described in section 3.2 in more detail. It is crucial that we assumed that the acoustic properties are known everywhere outside the sample $\Gamma$ (e.g., $\text{supp}(u) \subset \Gamma$), and that $p(x, t)$ is not practically measurable for $t > T$ anymore. Let $\hat{p}$ and $\hat{q}^*$ be the forward and adjoint wave fields corresponding to the delayed encoded sources $\hat{s}$ and measurements $\hat{f}^\delta$, cf (13). By the linearity and time invariance of the wave equation, we have $\hat{p}(x, t) = \sum_i w_i p_i(x, t - d_i \tau)$ and $\hat{q}^*(x, t) = \sum_i w_i q_i^*(x, t - d_i \tau)$. With this, we get (cf section 2.3 and (10)):

$$ 
\begin{aligned}
\nabla_u \mathcal{D}(\hat{f}(u), \hat{f}^\delta) &= \int_0^{T+\tau} \left(\frac{\partial A}{\partial u} \hat{p}(x, t)\right) \hat{q}^*(x, t) \, dt \\
&= \sum_{i,j}^{n_s} w_i w_j \int_0^{T+\tau} \left(\frac{\partial A}{\partial u} p_i(x, t - d_i \tau)\right) q_j^*(x, t - d_j \tau) \, dt \\
&= \sum_i^{n_s} w_i^2 \int_0^{T+\tau} \left(\frac{\partial A}{\partial u} p_i(x, t - d_i \tau)\right) q_i^*(x, t - d_i \tau) \, dt \\
&\quad + \sum_{i \neq j}^{n_s} w_i w_j \int_0^{T+\tau} \left(\frac{\partial A}{\partial u} p_i(x, t - d_i \tau)\right) q_j^*(x, t - d_j \tau) \, dt \\
&= \sum_i^{n_s} w_i^2 \int_0^T \left(\frac{\partial A}{\partial u} p_i(x, t)\right) q_i^*(x, t) \, dt \\
&\quad + \sum_{i \neq j}^{n_s} w_i w_j \int_0^{T+\tau} \left(\frac{\partial A}{\partial u} p_i(x, t - d_i \tau)\right) q_j^*(x, t - d_j \tau) \, dt \\
&= \sum_i^{n_s} w_i^2 \nabla_u \mathcal{D}(f_i(u), f_i^\delta) \\
&\quad + \sum_{i \neq j}^{n_s} w_i w_j \int_0^{T+\tau} \left(\frac{\partial A}{\partial u} p_i(x, t - d_i \tau)\right) q_j^*(x, t - d_j \tau) \, dt.
\end{aligned}
\tag{20}
$$

With $\mathbb{E}[w] = 0$, $\text{Cov}[w] = I$, it is easy to see now that $\mathbb{E}[\nabla_u \mathcal{D}(\hat{f}(u), \hat{f}^\delta)] = \nabla \mathcal{J}(u)$. If we choose the Rademacher distribution ($w_i = \pm 1$ with equal probability), we have that $w_i^2 = 1$
<caption>**Algorithm 2.** Stochastic limited memory BFGS (SLBFGS).</caption>

1: **Inputs:** $u, H_0, w, m_h, N_{evl}, \nu, do_linesearch$
2: $\bar{u}, \bar{w} \leftarrow 0$
3: $J_0 \leftarrow \infty$
4: $S, Y \leftarrow []$
5: **average_iterates** $\leftarrow$ **false**
6: $i \leftarrow 0$
7: **while** $n_{evl} < N_{evl}$ **do**
8: &nbsp;&nbsp;&nbsp;&nbsp;$i \leftarrow i + 1$
9: &nbsp;&nbsp;&nbsp;&nbsp;$F_u, G_u \leftarrow g(u, \theta(i), i)$
10: &nbsp;&nbsp;&nbsp;&nbsp;$n_{evl} \leftarrow n_{evl} + 1$
11: &nbsp;&nbsp;&nbsp;&nbsp;$z \leftarrow$ -twoLoopRecursion($G_u, S, Y, H_0$)
12: &nbsp;&nbsp;&nbsp;&nbsp;**if** $do_linesearch$ **then**
13: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$F_z, \bar{G}_z, \nu, m_{evl} \leftarrow$ lineSearch($u, z, g(\cdot, \theta(i), i), \nu$)
14: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$F_z, \bar{G}_z \leftarrow g(u + \nu z, \theta(i), i)$
15: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$m_{evl} \leftarrow 1$
16: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$m_{evl} \leftarrow 1$
17: &nbsp;&nbsp;&nbsp;&nbsp;**end if**
18: &nbsp;&nbsp;&nbsp;&nbsp;$n_{evl} \leftarrow n_{evl} + m_{evl}$
19: &nbsp;&nbsp;&nbsp;&nbsp;$S \leftarrow$ updateBuffer($\nu z, S, m_h$)
20: &nbsp;&nbsp;&nbsp;&nbsp;$Y \leftarrow$ updateBuffer($G_z - G_u, Y, m_h$)
21: &nbsp;&nbsp;&nbsp;&nbsp;$z \leftarrow z -$ twoLoopRecursion($G_z, S, Y, H_0$)
22: &nbsp;&nbsp;&nbsp;&nbsp;$u \leftarrow u + \nu z$
23: &nbsp;&nbsp;&nbsp;&nbsp;$u \leftarrow \max(\min(u, u_{max}), u_{min})$
24: &nbsp;&nbsp;&nbsp;&nbsp;$J_i \leftarrow \min(F_u, F_z)$
25: &nbsp;&nbsp;&nbsp;&nbsp;**if** $average_iterates$ **then**
26: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\bar{u} \leftarrow \bar{u} + w(i)u$
27: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\bar{w} \leftarrow \bar{w} + w(i)$
28: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\hat{u} \leftarrow \bar{u}/\bar{w}$
29: &nbsp;&nbsp;&nbsp;&nbsp;**else**
30: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**if** $J_i > J_{i-1}$ **then**
31: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\bar{u} \leftarrow u$
32: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\bar{u} \leftarrow u$
33: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\hat{u} \leftarrow u$
34: &nbsp;&nbsp;&nbsp;&nbsp;**end if**
35: **end while**
36: **return:** $\hat{u}, J$
<figure><img src="image_17.png" /><figcaption>Figure 13. Temporal evolution of the pressure field power $P_i(t)$ for eight different sources.</figcaption></figure>

and

$$ \nabla_u \mathcal{D}(\hat{f}(u), \hat{f}^\delta) = \nabla \mathcal{J}(u) + \sum_{i \neq j}^{n_s} w_i w_j X_{ij}, \qquad (21) $$

$$ \text{where } X_{ij} := \int_{0}^{T+\tau} \left( \frac{\partial A}{\partial u} p_i(x, t - d_i \tau) \right) q_j^*(x, t - d_j \tau) dt $$

is the cross talk between sources i and j. To examine the impact of delays ($\tau > 0$) vs conventional source encoding ($\tau = 0$), we assume $\bar{d} = d_j - d_i > 0$ without loss of generality and shift the time in the integral by $d_i \tau$ to get

$$ X_{ij}(\tau) = \int_{0}^{T+\tau} \left( \frac{\partial A}{\partial u} p_i(x, t) \right) q_j^*(x, t - \bar{d}\tau) dt. \qquad (22) $$

For $\tau \ge T/\bar{d}$, the delay between both sources is bigger than $T$, and from our assumptions, it follows that $X_{i,j}(\tau) = 0$. To obtain generic estimates for $\tau < T/\bar{d}$, one needs to estimate how fast the field power

$$ P_i(t) := \int_{\Gamma} p_i^2(x, t) dx \qquad (23) $$

decays inside $\Gamma$. For such decay estimates, see section ‘local energy decay estimates’ in [57] and references therein. In 3D, there is a $T_d$ such that $P_i(t) < a e^{-bt}$ with $a, b > 0$ for all $t > T_d$. In figure 13, we plot $P_i(t)$ for eight different sources in the scenario used in the numerical studies to show that practically, the decay is quite rapid, which is also confirmed by the results in section 4.4, cf figure 6.

## Appendix F. Noise sensitivity

In this section, we demonstrate the impact of adding measurement noise to the simulated data. The noise characteristics of UST systems can differ quite a lot from conventional linear
<figure><img src="image_18.png" /><figcaption>Figure 14. (a) Comparison of clean data $f(t)$ and noisy data $f(t) + \varepsilon(t)$ for SNR 5: the clean data is the same time trace shown as 'f(t) fine level' in figure 12. The noise signal $\varepsilon(t)$ for SNR 10, 15, 20, 25, 30 is the same just multiplied with $10^{\frac{5-SNR}{20}} = 0.56, 0.32, 0.18, 0.10, 0.056$. (b) Quantitative results of the noise sensitivity study.</figcaption></figure>

array US systems, e.g., both the US transducers and the data acquisition (DAQ) system of the PAMMOTH scanner were custom designed and a detailed noise characterization will be included in forthcoming work. Here, we will assume that an additive Gaussian noise model is adequate to summarize the different noise contributions and further that a decorrelation transform will render its covariance matrix sufficiently close to a scaled identity:

$$f^{\delta} = A(u^{\dagger}) + \varepsilon, \quad \text{with } \varepsilon \sim \mathcal{N}(0, \sigma^2 I). \qquad (24)$$

In the following, we will use values of $\sigma$ determined by the signal-to-noise (SNR) level defined in terms of the root mean square (RMS) of the clean signal $f$ as

$$\sigma = \text{RMS}(f) 10^{-\frac{\text{SNR}}{20}}. \qquad (25)$$

Figure 14(a) shows an example of a clean and noisy time trace for SNR 5.

Without measurement noise ($\varepsilon = 0$), all noise-like features in the reconstructed images were a result of using stochastic gradient estimators in the optimization schemes to solve the FWI problem (8), i.e., a purely numerical phenomena. This 'gradient-noise' vanishes as the iteration progresses (cf figure 3) or if more accurate gradient estimators are used (cf figure 11). From the way that $\nabla \mathcal{J}(u)$ is computed as presented in section 2.3 is clear that any measurement noise $\varepsilon > 0$ will be back-propagated into the image alongside with the data residual. As the iteration progresses, the norm of this type of image noise will grow while its spectral content changes. This phenomena arises from the fact that in principle, we use a Landweber-type iteration to solve an ill-conditioned inverse problem (7). See [15, 90] for a general introduction into this topic. Here, we illustrate these two opposing effects by first using exactly the same settings as in section 5 on a single GPU for different SNR values. The first column of figures 15 and 16 show that while SNR 30 is visually indistinguishable from clean data result (figure 11(a)) the image noise increases towards SNR 5. Then we change the number of gradient evaluations on finest spatial level from 32 to 64. The second column of figures 15 and 16 and the first
<figure><img src="image_19.png" /><figcaption>**Figure 15.** Qualitative results of the noise sensitivity study. Shown are different combinations of noise level (SNR), maximal number of gradient evaluations (32 eval vs 64 eval) and of using additional total variation regularization.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>**Figure 16.** Error plots corresponding to figure 15.</figcaption></figure>
two plots in figure 14(b) show the results and confirm that for the highest noise level (SNR = 5), doing more iterations now indeed leads to a worse result. To prevent this, we will need to add a regularization functional to the data fidelity $\sum_{i}^{n_s} D_i(u)$ in (8). For simplicity, we use a smoothed version of the popular total variation (TV) functional as described in [90] here:

$$ \min_{u \in \mathcal{U}} \sum_{i}^{n_s} D_i(u) + \alpha \mathrm{TV}(u) \quad \text{with} \\ \mathrm{TV}(u) := \sum_{i} \sqrt{(\partial_x u)_i^2 + (\partial_y u)_i^2 + (\partial_z u)_i^2 + \beta^2} - \beta, \quad (26) $$

where we set $\beta$ to 1 m s$^{-1}$ and fix $\alpha > 0$ to an illustrative value that ensures observable impact of the regularization (for optimal results, one would need to choose it in dependence of $\sigma$). The third and fourth columns of figures 15 and 16 and the second two plots in figure 14(b) show that now the errors do not increase with increasing iteration count. Note however that it is not straight forward to predict the overall impact of adding regularization to FWI and in particular the question which type of functionals are advantageous for breast UST needs further research.

## ORCID iDs

Felix Lucka https://orcid.org/0000-0002-8763-5177

Mailyn Pérez-Liva https://orcid.org/0000-0002-4985-247X

Bradley E Treeby https://orcid.org/0000-0001-7782-011X

Ben T Cox https://orcid.org/0000-0001-7296-4093

## References

[1] Delphinus Medical Technologies 2020 http://delphinusmt.com/technology/ (accessed 12 December 2020)

[2] Pammoth 2020 Photoacoustic mammoscopy for evaluating screening-detected lesions in the breast https://pammoth-2020.eu/ (accessed 12 December 2020)

[3] QT Ultrasound 2020 https://qtultrasound.com/research/ (accessed 12 December 2020)

[4] Almansouri H, Venkatakrishnan S V, Buzzard G T, Bouman C A and Santos-Villalobos H 2018 Deep neural networks for non-linear model-based ultrasound reconstruction 2018 IEEE Global Conf. on Signal and Information Processing (GlobalSIP) (IEEE) pp 6–10

[5] Andersen A 1990 A ray tracing approach to restoration and resolution enhancement in experimental ultrasound tomography Ultrason. Imaging 12 268–91

[6] Arridge S, Maass P, Öktem O and Schönlieb C-B 2019 Solving inverse problems using data-driven models Acta Numer. 28 1–174

[7] Arridge S R, Betcke M M, Cox B T, Lucka F and Treeby B E 2016 On the adjoint operator in photoacoustic tomography Inverse Problems 32 115012

[8] Bachmann E and Tromp J 2020 Source encoding for viscoacoustic ultrasound computed tomography J. Acoust. Soc. Am. 147 3221–35

[9] Beard P 2011 Biomedical photoacoustic imaging Interface Focus 1 602–31

[10] Bernard S, Monteiller V, Komatitsch D and Lasaygues P 2017 Ultrasonic computed tomography based on full-waveform inversion for bone quantitative imaging Phys. Med. Biol. 62 7011–35

[11] Bertsekas D P 2015 Incremental gradient, subgradient, and proximal methods for convex optimization: a survey (arXiv:1507.01030)

[12] Birk M, Dapp R, Ruiter N V and Becker J 2014 GPU-based iterative transmission reconstruction in 3D ultrasound computer tomography J. Parallel Distrib. Comput. 74 1730–43
[13] Boehm C, Martiartu N K, Vinard N, Balic I J and Fichtner A 2018 Time-domain spectral-element ultrasound waveform tomography using a stochastic quasi-Newton method *Medical Imaging 2018: Ultrasonic Imaging and Tomography* ed N Duric and B C Byram (International Society for Optics and Photonics, SPIE) vol 10580 pp 92–100

[14] Bottou L, Curtis F E and Nocedal J 2018 Optimization methods for large-scale machine learning *SIAM Rev.* **60** 223–311

[15] Burger M, Kaltenbacher B and Neubauer A 2011 Iterative solution methods *Handbook of Mathematical Methods in Imaging* ed O Scherzer (New York: Springer) pp 345–84

[16] Carević A, Yun X and Almekkawy M 2019 Adaptive truncated total least square on distorted born iterative method in ultrasound inverse scattering problem *Medical Imaging 2019: Ultrasonic Imaging and Tomography* (International Society for Optics and Photonics) vol 10955 p 1095515

[17] Carson P L, Meyer C R, Scherzinger A L and Oughton T V 1981 Breast imaging in coronal planes with simultaneous pulse echo and transmission ultrasound *Science* **214** 1141–3

[18] Cassereau D and Fink M 1992 Time-reversal of ultrasonic fields: III. Theory of the closed time-reversal cavity *IEEE Trans. Ultrason. Ferroelectr. Freq. Control* **39** 579–92

[19] Cerveny V 2001 *Seismic Ray Theory* (Cambridge: Cambridge University Press)

[20] Cobbold R S C 2007 *Foundations of Biomedical Ultrasound* (Oxford: Oxford University Press)

[21] Conrath B C, Daft C M W and O'Brien W D 1989 Applications of neural networks to ultrasound tomography *Proc. IEEE Ultrasonics Symp.* (IEEE) pp 1007–10

[22] Dapp R, Gemmeke H and Ruiter N 2012 3D refraction-corrected transmission reconstruction for 3D ultrasound computer tomography *Medical Imaging 2012: Ultrasonic Imaging, Tomography, and Therapy* (International Society for Optics and Photonics) vol 8320 p 832014

[23] Duric N and Littrup P 2018 *Breast Imaging, Chapter Breast Ultrasound Tomography* (London: IntechOpen)

[24] Duric N et al 2013 Breast imaging with the SoftVue imaging system: first results *Medical Imaging 2013: Ultrasonic Imaging, Tomography, and Therapy* (International Society for Optics and Photonics) vol 8675 p 86750K

[25] Engquist B, Ren K and Yang Y 2020 The quadratic Wasserstein metric for inverse data matching *Inverse Problems* **36** 055001

[26] Esser E, Guasch L, van Leeuwen T, Aravkin A Y and Herrmann F J 2018 Total variation regularization strategies in full-waveform inversion *SIAM J. Imaging Sci.* **11** 376–406

[27] Fabien-Ouellet G, Gloaguen E and Giroux B 2017 A stochastic L-BFGS approach for full-waveform inversion *SEG Technical Program Expanded Abstracts* vol 22 pp 1622–7

[28] Fan Y, Wang H, Gemmeke H and Hesser J 2020 MI-net: a deep network for non-linear ultrasound computed tomography reconstruction *2020 IEEE Int. Ultrasonics Symp. (IUS)* (IEEE) pp 1–3

[29] Fan Y and Ying L 2019 Solving traveltime tomography with deep learning (arXiv:1911.11636)

[30] Feigin M, Freedman D and Anthony B W 2019 A deep learning framework for single-sided sound speed inversion in medical ultrasound *IEEE Trans. Biomed. Eng.* **67** 1142–51

[31] Feigin M, Makovsky Y, Freedman D and Anthony B W 2020 High-frequency full-waveform inversion with deep learning for seismic and medical ultrasound imaging *SEG Technical Program Expanded Abstracts 2020* (Society of Exploration Geophysicists) pp 3492–6

[32] Feigin M, Zwecker M, Freedman D and Anthony B W 2020 Detecting muscle activation using ultrasound speed of sound inversion with deep learning *2020 42nd Annual Int. Conf. of the IEEE Engineering in Medicine and Biology Society (EMBC)* (IEEE) pp 2092–5

[33] Fichtner A 2011 *Full Seismic Waveform Modelling and Inversion* (Berlin: Springer)

[34] Fink M 1992 Time reversal of ultrasonic fields: I. Basic principles *IEEE Trans. Ultrason. Ferroelectr. Freq. Control* **39** 555–66

[35] Gemmeke H, Berger L, Hopp T, Zapf M, Tan W, Blanco R, Leys R, Peric I and Ruiter N V 2018 The new generation of the kit 3D USCT *Proc. Int. Workshop on Medical Ultrasound Tomography* (Speyer, Germany 1–3 November 2017) ed T Hopp (KIT Scientific Publishing) pp 271–82

[36] Gemmeke H, Hopp T, Zapf M, Kaiser C and Ruiter N V 2017 3D ultrasound computer tomography: hardware setup, reconstruction methods and first clinical results *Nucl. Instrum. Methods Phys. Res.* **A** **873** 59–65

[37] Goncharsky A V, Romanov S Y and Seryozhnikov S Y 2018 Low-frequency 3D ultrasound tomography: dual-frequency method *Numer. Methods Program.* **19** 479–95

[38] Goncharsky A V and Romanov S Y 2013 Supercomputer technologies in inverse problems of ultrasound tomography *Inverse Problems* **29** 075004
[39] Goncharsky A V, Romanov S Y and Seryozhnikov S Y 2016 A computer simulation study of soft tissue characterization using low-frequency ultrasonic tomography *Ultrasonics* **67** 136–50

[40] Goncharsky A V, Romanov S Y and Seryozhnikov S Y 2019 Comparison of the capabilities of GPU clusters and general-purpose supercomputers for solving 3D inverse problems of ultrasound tomography *J. Parallel Distrib. Comput.* **133** 77–92

[41] Goncharsky A V, Romanov S Y and Seryozhnikov S Y 2019 Low-frequency ultrasonic tomography: mathematical methods and experimental results *Mosc. Univ. Phys.* **74** 43–51

[42] Goncharsky A V and Seryozhnikov S Y 2020 Capabilities of layered ultrasound tomography *Lobachevskii J. Math.* **41** 1454–62

[43] Guasch L, Calderón Agudo O, Tang M-X, Nachev P and Warner M 2020 Full-waveform inversion imaging of the human brain *npj Digit. Med.* **3** 28

[44] Guo R, Lu G, Qin B and Fei B 2018 Ultrasound imaging technologies for breast cancer detection and management: a review *Ultrasound Med. Biol.* **44** 37–70

[45] Haber E, Chung M and Herrmann F 2012 An effective method for parameter estimation with PDE constraints with multiple right-hand sides *SIAM J. Optim.* **22** 739–57

[46] Haddadin O S and Ebbini E S 1998 Imaging strongly scattering media using a multiple frequency distorted Born iterative method *IEEE Trans. Ultrason. Ferroelectr. Freq. Control* **45** 1485–96

[47] Heller S L and Moy L 2019 MRI breast screening revisited *J. Magn. Reson. Imaging* **49** 1212–21

[48] Huang L, Shin J, Chen T, Lin Y, Gao K, Intrator M and Hanson K 2016 Breast ultrasound tomography with two parallel transducer arrays *Medical Imaging 2016: Physics of Medical Imaging* (International Society for Optics and Photonics) vol 9783 p 97830C

[49] Huthwaite P and Simonetti F 2011 High-resolution imaging without iteration: a fast and robust method for breast ultrasound tomography *J. Acoust. Soc. Am.* **130** 1721–34

[50] Igel H 2016 *Computational Seismology* (Oxford: Oxford University Press)

[51] Javaherian A and Holman S 2017 A multi-grid iterative method for photoacoustic tomography *IEEE Trans. Med. Imaging* **36** 696–706

[52] Javaherian A, Lucka F and Cox B 2020 Refraction-corrected ray-based inversion for three-dimensional ultrasound tomography of the breast *Inverse Problems* **36** 125010

[53] Jirik R, Peterlik I, Ruiter N, Fousek J, Dapp R, Zapf M and Jan J 2012 Sound-speed image reconstruction in sparse-aperture 3D ultrasound transmission tomography *IEEE Trans. Ultrason. Ferroelectr. Freq. Control* **59** 254–64

[54] Johnson S A, Greenleaf J F, Samayoa W A, Duck F A and Sjostrand J 1975 Reconstruction of three-dimensional velocity fields and other parameters by acoustic ray tracing *Ultrason Symp. Proc.* (IEEE) pp 46–51 (75 CHO 994-4SU)

[55] Jush F K, Biele M, Dueppenbecker P M, Schmidt O and Maier A 2020 DNN-based speed-of-sound reconstruction for automated breast ultrasound *2020 IEEE Int. Ultrasonics Symp. (IUS)* (IEEE) pp 1–7

[56] Kothari K, de Hoop M and Dokmanić I 2020 Learning the geometry of wave-based imaging *Advances in Neural Information Processing Systems* vol 33

[57] Kuchment P and Kunyansky L 2011 Mathematics of photoacoustic and thermoacoustic tomography *Handbook of Mathematical Methods in Imaging* ed O Scherzer (New York: Springer) pp 817–65

[58] Li C, Duric N, Littrup P and Huang L 2009 *In vivo* breast sound-speed imaging with ultrasound tomography *Ultrasound Med. Biol.* **35** 1615–28

[59] Li S, Jackowski M, Dione D P, Varslot T, Staib L H and Mueller K 2010 Refraction corrected transmission ultrasound computed tomography for application in breast imaging *Med. Phys.* **37** 2233–46

[60] Lou Y, Zhou W, Matthews T P, Appleton C M and Anastasio M A 2017 Generation of anatomically realistic numerical phantoms for photoacoustic and ultrasonic breast imaging *J. Biomed. Opt.* **22** 041015

[61] Lyu J, Tan C and Dong F 2020 Image reconstruction for ultrasonic tomography using multi-channel convolutional neural network *2020 39th Chinese Control Conf. (CCC)* (IEEE) pp 7304–9

[62] Malik B, Terry R, Wiskin J and Lenox M 2018 Quantitative transmission ultrasound tomography: imaging and performance characteristics *Med. Phys.* **45** 3063–75

[63] Manohar S and Dantuma M 2019 Current and future trends in photoacoustic breast imaging *Photoacoustics* **16** 100134

[64] Manohar S, Willemink R G H, van der Heijden F, Slump C H and van Leeuwen T G 2007 Concomitant speed-of-sound tomography in photoacoustic imaging *Appl. Phys. Lett.* **91** 131911
[65] Marmot M G, Altman D G, Altman D G, Cameron D A, Dewar J A, Thompson S G and Wilcox M 2013 The benefits and harms of breast cancer screening: an independent review *Br. J. Cancer* **108** 2205–40

[66] Mast T D 2002 Aberration correction for time-domain ultrasound diffraction tomography *J. Acoust. Soc. Am.* **112** 55–64

[67] Matthews T P, Wang K, Li C, Duric N and Anastasio M A 2017 Regularized dual averaging image reconstruction for full-wave ultrasound computed tomography *IEEE Trans. Ultrason. Ferroelectr. Freq. Control* **64** 811–25

[68] Mulder W A 2018 Higher-order source-wavefield reconstruction for reverse time migration from stored values in a boundary strip just one point wide *Geophysics* **83** T31–8

[69] Nandapalan N, Jaros J, Rendell A P and Treeby B 2012 Implementation of 3D FFTs across multiple GPUs in shared memory environments *2012 13th Int. Conf. on Parallel and Distributed Computing, Applications and Technologies* pp 167–72

[70] Opieliński K J, Pruchnicki P, Szymanowski P, Szepieniec W K, Szweda H, Świś E, Jóźwik M, Tenderenda M and Bułkowski M 2018 Multimodal ultrasound computer-assisted tomography: an approach to the recognition of breast lesions *Comput. Med. Imaging Graph.* **65** 102–14

[71] Plessix R-E 2006 A review of the adjoint-state method for computing the gradient of a functional with geophysical applications *Geophys. J. Int.* **167** 495–503

[72] Pérez-Liva M, Herraiz J L, Udías J M, Miller E, Cox B T and Treeby B E 2017 Time domain reconstruction of sound speed and attenuation in ultrasound computed tomography using full wave inversion *J. Acoust. Soc. Am.* **141** 1595–604

[73] Sandhu G Y, West E, Li C, Roy O and Duric N 2017 3D frequency-domain ultrasound waveform tomography breast imaging *Medical Imaging 2017: Ultrasonic Imaging and Tomography* ed N Duric and B Heyde (International Society for Optics and Photonics, SPIE) vol 10139 pp 56–69

[74] Schmitt R M 2002 Ultrasound computed tomography: from the past to the future *Acoustical Imaging* (Berlin: Springer) pp 25–35

[75] Siesling S et al 2020 Diagnostic work-up in women suspect for breast cancer in The Netherlands *Eur. J. Cancer* **138** S87

[76] Sood R, Rositch A F, Shakoor D, Ambinder E, Pool K-L, Pollack E, Mollura D J, Mullen L A and Harvey S C 2019 Ultrasound for breast cancer detection globally: a systematic review and meta-analysis *J. Glob. Oncol.* **5** 1–17

[77] Tabei M, Mast T D and Waag R C 2002 A *k*-space method for coupled first-order acoustic propagation equations *J. Acoust. Soc. Am.* **111** 53–63

[78] Tarantola A 2005 *Inverse Problem Theory and Methods for Model Parameter Estimation* (Philadelphia, PA: SIAM)

[79] Taskin U, Solberg Eikrem K, Nævdal G, Jakobsen M, Verschuur D J and van Dongen K W A 2020 Ultrasound imaging of the brain using full-waveform inversion *2020 IEEE Int. Ultrasonics Symp. (IUS)* pp 1–4

[80] Treeby B E and Cox B T 2010 k-Wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wave fields *J. Biomed. Opt.* **15** 021314

[81] Treeby B E and Cox B T 2010 Modeling power law absorption and dispersion for acoustic propagation using the fractional Laplacian *J. Acoust. Soc. Am.* **127** 2741–8

[82] Treeby B E, Jaros J, Rendell A P and Cox B T 2012 Modeling nonlinear ultrasound propagation in heterogeneous media with power law absorption using a *k*-space pseudospectral method *J. Acoust. Soc. Am.* **131** 4324–36

[83] Treeby B E, Vaverka F and Jaros J 2016 Spectral domain decomposition using local Fourier basis: application to ultrasound simulation on a cluster of GPUs *Supercomput. Front. Innov.* **3** 40–55

[84] Treeby B E, Vaverka F and Jaros J 2018 Performance and accuracy analysis of nonlinear k-wave simulations using local domain decomposition with an 8-GPU server *Proc. Meetings on Acoustics* vol 34 p 022002

[85] Treeby B E, Zhang E Z and Cox B T 2010 Photoacoustic tomography in absorbing acoustic media using time reversal *Inverse Problems* **26** 115003

[86] Tromp J 2020 Seismic wavefield imaging of Earth’s interior across scales *Nat. Rev. Earth Environ.* **1** 40–53

[87] Tromp J and Bachmann E 2019 Source encoding for adjoint tomography *Geophys. J. Int.* **218** 2019–44
[88] van Herwaarden D P, Boehm C, Afanasiev M, Thrastarson S, Krischer L, Trampert J and Fichtner A 2020 Accelerated full-waveform inversion using dynamic mini-batches *Geophys. J. Int.* **221** 1427–38

[89] van Sloun R J G, Cohen R and Eldar Y C 2020 Deep learning in ultrasound imaging *Proc. IEEE 108* 11–29

[90] Vogel C R 2002 *Computational Methods for Inverse Problems* (Philadelphia, PA: SIAM)

[91] Vourtsis A 2019 Three-dimensional automated breast ultrasound: technical aspects and first results *Diagn. Intervent. Radiol.* **100** 579–92

[92] Wang K, Matthews T, Anis F, Li C, Duric N and Anastasio M 2015 Waveform inversion with source encoding for breast sound speed reconstruction in ultrasound computed tomography *IEEE Trans. Ultrason. Ferroelectr. Freq. Control* **62** 475–93

[93] Wiskin J, Borup D and Johnson S 2011 Inverse scattering theory *Acoustical Imaging* (Berlin: Springer) pp 53–9

[94] Wiskin J, Borup D T, Johnson S A, Berggren M, Abbott T and Hanover R 2007 Full-wave, non-linear, inverse scattering *Acoustical Imaging* (Berlin: Springer) pp 183–93

[95] Wiskin J, Borup D T, Johnson S A and Berggren M 2012 Non-linear inverse scattering: high resolution quantitative breast tissue tomography *J. Acoust. Soc. Am.* **131** 3802–13

[96] Wiskin J W, Borup D T, Iuanow E, Klock J and Lenox M W 2017 3D nonlinear acoustic inverse scattering: algorithm and quantitative results *IEEE Trans. Ultrason. Ferroelectr. Freq. Control* **64** 1161–74

[97] Wu F, Thomas J-L and Fink M 1992 Time reversal of ultrasonic fields: II. Experimental results *IEEE Trans. Ultrason. Ferroelectr. Freq. Control* **39** 567–78

[98] Yang Y 2018 Optimal transport for seismic inverse problems *PhD Thesis* University of Texas at Austin

[99] Zhang Z and Huang L 2014 Efficient implementation of ultrasound waveform tomography using data blending *Medical Imaging 2014: Ultrasonic Imaging and Tomography* ed J G Bosch and M M Doyley vol 9040 (International Society for Optics and Photonics, SPIE) pp 436–44

[100] Zhao W, Wang H, Gemmeke H, van Dongen K W A, Hopp T and Hesser J W 2020 Ultrasound transmission tomography image reconstruction with fully convolutional neural network *Phys. Med. Biol.* **65** 235021

[101] Zhu Q and Poplack S 2020 A review of optical breast imaging: multi-modality systems for breast cancer diagnosis *Eur. J. Radiol.* **129** 109067