# Model reduction in acoustic inversion by artificial neural network

Janne Koponena,b, Timo Lähivaaraa, Jari Kaipioc, Marko Vauhkonena
aDepartment of Applied Physics, University of Eastern Finland, Kuopio, Finland
bSchool of Engineering and Technology, Savonia University of Applied Sciences,
Kuopio, Finland
cDepartment of Mathematics, University of Auckland, Auckland, New Zealand

## Abstract

In ultrasound tomography, the speed of sound inside an object is estimated
based on acoustic measurements carried out by sensors surrounding the object.
An accurate forward model is a prominent factor for high-quality image
reconstruction, but it can make computations far too time-consuming in many
applications. Using approximate forward models, it is possible to speed up the
computations, but the quality of the reconstruction may have to be
compromised. In this paper, a neural network -based approach is proposed, that
can compensate for modeling errors caused by the approximate forward models.
The approach is tested with various different imaging scenarios in a simulated
two-dimensional domain. The results show that with fairly small training data
sets, the proposed approach can be utilized to approximate the modelling
errors, and to significantly improve the image reconstruction quality in
ultrasound tomography, compared to commonly used inversion algorithms.

##  1 Introduction

In ultrasound tomography, speed of sound (SOS) distribution in an object is
estimated based on acoustic measurements. Mathematically, the reconstruction
of SOS is a nonlinear inverse problem. Many different techniques have been
developed to solve this problem, for example, algebraic reconstruction methods
[13, 22, 1], adjoint methods [40, 10, 11, 35], and artificial neural networks
(ANN) [45, 41, 9, 25]. Algebraic reconstruction methods are typically
computationally quite light methods to solve the inverse problem. On the other
hand, a fairly large number of sensors are required to achieve high-quality
reconstructions. Adjoint methods, on the other hand, can produce high-quality
reconstructions with a smaller number of measurements. However, the adjoint
methods are computationally more demanding, and thus they can be impractical
in some cases. The ANN-based inversion methods are typically computationally
less intensive, but require a large number of training samples, and therefore
these methods can be infeasible in some applications.

Solving the acoustic inverse problem with an adjoint method, a large number of
corresponding forward problems need to be solved. Computational requirements
of computing the reconstruction can be reduced by using computationally less
expensive forward models. However, these models are typically less accurate
compared to complete models, and hence, the quality of the reconstruction
tends to be lower. Various methods have been developed to reduce the accuracy
requirements of the forward model while preserving the accuracy of the inverse
solution. These methods include, for example, the Bayesian Approximation Error
method (BAE) [20, 21], correction methods [8, 31], and Hybrid Algorithm for
Robust Breast Ultrasound Tomography (HARBUT) [18]. For example in BAE, the
modeling errors are approximated in the statistical framework by Gaussian
distributed random variables, and the reconstruction is achieved by computing
a maximum a posteriori estimate of the parameters.

Artificial neural networks are a computational model that is inspired by
neurobiological systems. The ANN has been used in numerous applications, such
as image recognition [46, 36], natural language processing [12], economy [33,
2, 27], and medicine [3, 34]. Also, as mentioned earlier, the ANN has been
utilized to solve inverse problems[5]. In the SOS inverse problem, the
standard approach is to train the ANN with a large set of simulated SOS
distributions and the corresponding synthetic acoustic measurements. In the
image reconstruction, the input to the network is the acoustic measurement
data and the outcome is the desired SOS distribution.

A drawback of common ANN approaches for solving inverse problems is the
computationally expensive training phase, requiring tens of thousands of
training data sets. This may require a notably large amount of computational
resources already in the training phase. In this paper, we propose a novel
image reconstruction approach that utilizes the ANN to create a numerical
function that converts the measured signal to a corresponding low dimensional,
surrogate model, signal. This means, that the chosen ANN can be trained with a
(fairly small) set of SOS distributions with the corresponding accurate
acoustic data and the desired light (surrogate) model data. In the inverse
solution phase, the measured data is converted to the surrogate model data
with the trained ANN and the final image reconstruction is carried out using a
standard adjoint approach, with the chosen surrogate model. In the rest of the
paper, the proposed approach is called the Neural Network based Approximation
Error Estimation method (NNAEE).

The rest of the paper is organized as follows. The acoustic forward model is
presented in Section 2, modelling errors and the ANN approach of using a
surrogate model for solving the inverse problem is discussed in Sections 3 and
4. In Section 5, the proposed method is demonstrated by numerical examples,
and finally, in Section 6 the conclusions are given.

##  2 The forward model

In ultrasound tomography, $M_{s}$ transmitters and $M_{r}$ receivers are
located around the region of interest $\Omega$. A schematic picture of the
measurement setup is presented in Figure 1. One transmitter at a time sends an
acoustic signal, which propagates through the object. The propagated signal is
measured by all the receivers simultaneously. Wave propagation is modeled by
the wave equation in an unbounded domain, $x\in\mathbb{R}^{2}$ and $t>0$, as

$-\nabla^{2}u+\frac{1}{c^{2}}\frac{\partial^{2}u}{\partial
t^{2}}=f\delta(x-x_{t}),$ (1)

with the initial conditions

$\left\\{\begin{array}[]{lcll}u(x,0)&=&0\\\ \frac{\partial}{\partial
t}u(x,0)&=&0,\end{array}\right.$

where $u=u(x,t)$ is pressure, $c=c(x)$ is the speed of sound (SOS), $f=f(t)$
and $f\delta(x-x_{t})$ is a source term, where $\delta(x-x_{t})$ is Dirac
delta function, $t$ is time, and $x_{t}$ is a position of the transmitter. The
measured data is collected at fixed time instants over a certain measurement
period. When collecting data at $M_{t}$ time instants, a total of $M_{r}\times
M_{s}\times M_{t}$ data points are obtained in the measurement vector
$\hat{y}$. With a fixed source term $f\delta(x-x_{t})$, the measured data
depends only on the parameter $v(x):=\frac{1}{c^{2}(x)}$.

In practice, when numerical models are utilized, the unbounded computational
domain must be truncated. The arising boundaries generate reflections into the
modeled signal. There are several ways to reduce the reflections by
implementing Absorbing Boundary Conditions (ABC). These include methods such
as perfectly matching layers [29, 4], high-order boundary conditions [14, 15]
and absorbing layers near boundaries [19, 39]. In addition, there are numerous
methods to solve the wave equation, see for example [26, 37, 44, 32, 28, 23,
6, 7]. In this paper, a finite-difference time-domain (FDTD) solver is used to
solve the forward problem. In addition, absorbing layers near the boundaries
are utilized to implement the ABC. The absorbing layers are implemented by
adding a damping term, $\eta\frac{\partial u}{\partial t}$, into left hand
side of equation (1). Inside of area surrounded by the sensors, $\eta=0$, and
$\eta>0$ elsewhere.

Figure 1: Schematic figure of a measurement setup in ultrasound tomography.
The transmitter $s$ sends waveform, which propagates through the object. The
propagated signal is measured by sensors located on the dashed circle.

##  3 Model selections, inverse problem and modelling error

Let us write the true physical observation model as

$\hat{y}=\hat{A}(\hat{v})+\hat{e},$ (2)

where $\hat{y}$ is the measured signal, $\hat{A}$ is the true physical model,
$\hat{v}$ is the parameter vector of the model, and $\hat{e}$ is the
measurement noise. If synthetic measurements are used, this model and the
corresponding data vector are created with very accurate FDTD discretization.

In practice, the true physical model is rarely known accurately and is
therefore approximated by some high accuracy observation model as

$\bar{y}=\bar{A}(\bar{v})+\bar{e},$ (3)

where the model $\bar{A}$ and the corresponding parameter vector $\bar{v}$
approximate the true physical model as well as possible. This model, however,
is computationally demanding and the use of it in inverse problem solution is
not desirable. To overcome this, a less accurate forward model, also called as
surrogate model, is chosen and is written in the form

${y}={A}({v})+e,$ (4)

where ${A}$ is a very approximative forward model and ${y}$ is the
corresponding data vector. Here the parameter vector ${v}=P\hat{v}$, where the
operator $P$ projects the true $\hat{v}$ to the lower dimensional basis of the
surrogate model.

###  3.1 The proposed NNAEE approach

In solving the inverse problems, the accuracy of the forward model and its
numerical implementation play a crucial role. In general, accurate numerical
models are expensive to compute, and on the other hand, computationally
lighter models are less accurate. In the proposed NNAEE approach we proceed as
follows. Equation (2) can be written as

$\hat{y}={A}({v})+(\hat{A}(\hat{v})+\hat{e}-{A}({v}))={A}({v})-\nu(\hat{y}).$
(5)

This can further be written as

${A}({v})=\hat{y}+\nu(\hat{y})=:\phi(\hat{y}),$ (6)

where $\phi(\hat{y})$ is the unknown mapping that converts the real world
measurement signal $\hat{y}$ to a corresponding surrogate model signal
${A}({v})={y}$. Mapping $\phi$ is extremely complex, and finding a closed form
approximation of the mapping would be difficult. On the other hand, the
universal approximation theorem states [30, 16] that a large space of real-
valued functions can be approximated by ANN. Therefore, in the NNAEE approach,
the mapping $\phi$ is trained and approximated by an ANN. To compute the
reconstruction $v$ from the measured signal $\hat{y}$, the surrogate numerical
model (4) can be utilized without compromising the quality, if the mapping
$\phi$ is accurate enough. In the current paper, we utilize an adjoint method,
similar to [38], to solve the inverse problem. However, other methods could
also be used.

##  4 The structure of the proposed approach

In this section, the structure of the neural network including the
autoencoder/-decoder phase is described in detail. Also, the algorithm of the
proposed NNAEE approach is introduced.

###  4.1 Notations and parameters of the models used

In the numerical tests later on, we utilize the following three different
models

  1. 1.

Physically ”true” model $\hat{A}(\hat{v})$

  2. 2.

High quality numerical model $\bar{A}(\bar{v})$

  3. 3.

Computationally light surrogate model ${A}({v})$

The first, simulated physical model mimics the real physical world. In
simulations, the model utilizes finite-difference discretization in the
spatial domain and leapfrog integration in the time domain. In addition, the
measurement noise is simulated by Gaussian noise.

The second, accurate model is a high-quality numerical model of the physical
system. The model contains, also, measurement error. However, in real world,
the true measurement noise is rarely known, and hence we use a slightly
different distribution for the measurement noise compared to the noise of the
simulated physical model.

The third, the surrogate model, is a numerical model which is computationally
light and is used in the final image reconstruction. The model utilizes coarse
spatial discretization, and low-quality ABC.

###  4.2 Autoencoder and autodecoder

The measurement signal (2) contains measurement noise, and in addition, the
true information content of the measurement vector is much less than the
actual size of the signal. An autoencoder (AE) and autodecoder (AD) can be
utilized to suppress measurement noise [43, 42], and to reduce the measurement
vector dimension [24, 17].

In the proposed NNAEE approach for each transmitter-receiver pair, we utilize
a three-layer ANN for AE and a three-layer ANN for AD, presented in Figure 2.
This means that each measurement signal of length $M_{t}$ of one transmitter-
receiver pair is reduced to size $M_{p}$ with the AE unit. The total length of
the reduced data set is thus $\mathbb{R}^{M_{r}\times M_{s}\times M_{p}}$. The
AD unit is used to decode the signal of one transmitter-receiver pair back to
the original length $M_{t}$. In the simulation test cases, we have used
$M_{t}=2363$ and the width of the information bottleneck is chosen as
$M_{p}=50$.

Figure 2: Three layers implementing the autoencoder, mapping $\alpha$, on
left, and three layers implementing the autodecoder, mapping $\beta$, on
right. Layers 1 and 2 in both cases are dense layers using ReLU (rectified
linear unit) activation, and the last layers are dense layers with linear
activation. At the training phase, the autoencoder output is connected to the
autodecoder input.

We train the AE+AD by using both the surrogate model generated signals and
accurate model generated signals. We utilize both of the signals because they
are needed at a later stage for training the mapping $\phi$ and a large number
of training data sets decreases the risk of over-fitting. Because there are
$M_{s}\times M_{r}$ transmitter-receiver pairs, there are in total $2\times
M_{s}\times M_{r}\times N_{\mathrm{train}}$ samples for AE+AD training. Input
of the AE+AD, denoted as $X=A_{s,r}^{*}(v^{*})+e^{*}$, where
$A_{s,r}^{*}(v^{*})+e^{*}$ is either ${A}({v})+{e}$ or
$\bar{A}(\bar{v})+\bar{e}$, and $s,r$ is a transmitter-receiver pair, and
$e^{*}$ is simulated measurement noise. The desired output of the AE+AD is
$Y=A_{s,r}^{*}(v^{*})$, that is, the same signal as the input signal but
without the measurement noise. From now on, we denote AE as mapping $\alpha$,
and AD as mapping $\beta$.

Figure 3: A partial figure of the ANN approximating $\varphi$. Each input
contains all the measurements of one transmitter, compressed with the AE. The
input layer contains $M_{s}=16$ inputs, each size of $M_{r}\times
M_{p}=64\times 50$. Each ReLU-net at the first layer inputs data from 7
inputs, each net at the second layer utilizes data from 9 outputs of the first
layer etc. At the last layer, each input is added to the corresponding output
of the second last layer.

###  4.3 Mapping $\phi$

Let us define a composite mapping
$\phi(\hat{y})=\beta(\varphi(\alpha(\hat{y})))$, where $\varphi$ is
approximated by ANN, shown in Figure 3, and $\alpha$ and $\beta$ are the AE
and AD, respectively. To keep notation simple, we utilize here $M_{s}\times
M_{r}$ parallel AE and AD, one for each transmitter-receiver pair. The input
and output vectors of $\varphi$, both are in $\mathbb{R}^{M_{s}\times
M_{r}\times M_{p}}$. Using multiple layers of dense networks of such size is
infeasible. Decreasing the number of neuron connections it is possible to
create an ANN which does not require a large amount of GPU memory and can be
implemented in a standard desktop computer with, for example, 8 GB of GPU
memory. Decreasing the connections in the proposed NNAEE approach is
implemented by forming layers that contain parallel dense networks, with the
feature that all the networks are not connected to all neurons on the
consecutive layers, see Figure 3. Between each layer, selected outputs of
earlier level are merged for the input of the next level networks. Each subnet
input of the next level has unique set of earlier outputs. The network
contains three layers of dense networks using rectified linear units (ReLU),
and one layer with linear activation function. The last layer sums the input
and output of the linear units forming the output of $\varphi$.

###  4.4 The NNAEE image reconstruction algorithm

Collecting steps defined in earlier subsections, Algorithm 1 can be obtained.
The algorithm has two phases; one is done beforehand, only once, and the other
phase does the reconstruction and utilizes the actual measurement signal. As
the first step is done only once, the step can utilize large amounts of
computer resources. In practice, the step can be computationally quite
demanding, because the training of an ANN can require a large amount of
training data. The second step, on the other hand, utilizes the trained ANN,
which does not require an excessive amount of computing power, and most of the
computing power at this phase is utilized for solving the actual inverse
problem with the adjoint method.

A. Once, before any measurements

  1. 1.

Create training measurement signals

    1. (a)

Generate SOS fields representing the assumed true distributions

    2. (b)

Generate simulated measurement signals with the accurate model
$\bar{A}(\bar{v})$

    3. (c)

Generate simulated measurement signals with the surrogate model ${A}({v})$

  2. 2.

Train AE/AD

  3. 3.

Encode training signals by AE – Output is the training data set

  4. 4.

Train ANN representing the mapping $\varphi$ using encoded training data set
to form $\phi(\cdot)=\beta(\varphi(\alpha(\cdot)))$.

B. With the actual (or simulated) measurement data

  1. 1.

Perform a measurement (or simulate with $\hat{A}(\hat{v})$)

  2. 2.

Encode measurement signal

  3. 3.

Convert signal to corresponding surrogate model signal with the mapping $\phi$

  4. 4.

Solve the inverse problem (reconstruct the unknown $v$ field) using the
converted signal and the chosen surrogate model.

Algorithm 1 NNAEE

##  5 Numerical tests

In this section, we demonstrate the NNAEE method by investigating the
dependency of the quality of reconstruction on the size of the training
dataset. We also compare the proposed method to an accurate model
reconstruction and a conventional reconstruction method (CRM). In the accurate
model reconstruction, we use the model $\bar{A}(\bar{v})$ and in CRM the
reconstruction is obtained using the approximate model ${A}({v})$. In
addition, the method is compared to a direct ANN-based inversion. The ANN-
based inversion is implemented with a similar type of ANN as presented in
Figure 3. Comparison utilizes Root Mean Square Error,

${\rm
RMSE}(v,w)=\sqrt{\frac{1}{N_{\mathrm{grid}}^{2}}\sum_{i=1}^{N_{\mathrm{grid}}}\sum_{j=1}^{N_{\mathrm{grid}}}(v_{i,j}-w_{i,j})^{2}}\
\ ,$

where $v$ and $w$ are $N_{\mathrm{grid}}\times N_{\mathrm{grid}}$ SOS fields.

Training the ANN requires a large number of training data, and therefore using
real measurements to train the ANN is not feasible. The training dataset and
validation dataset are generated using a high quality numerical model
$\bar{A}(\bar{v})$, and the surrogate numerical model ${A}({v})$. The training
input data $X=\alpha(\bar{A}(\bar{v}))$, and the preferred output
$Y=\alpha(A(v))$, where $v=P\bar{v}$. In the simulations, the number of
transmitters and receivers were $M_{s}=16$ and $M_{r}=64$, respectively.

We use the following sizes of the training datasets:
$N_{\mathrm{\mathrm{train}}}\in\\{200,\ 400,\ 800,\ 1200,\ 1800,\ 2500,$
$5000,\ 10000,$ $25000,\ 50000\\}$ to see how the size of the training set
affects the accuracy of the reconstruction. Each training/validation data set
item is generated as shown in Algorithm 2. The size of a validation dataset
was 200 samples, and the size of a test dataset was $20$ samples. The
validation dataset was used to tune the hyperparameters of the network and
training. The test dataset was not used at the training phase. It was only
used to test the performance of the proposed method.

###  5.1 Model parameters

The parameters of the models are presented in Table 1. All models utilize FDTD
leapfrog solver with $\Delta t=10^{-7}\,{\textrm{s}}$ time step to solve the
partial differential equation. The spatial grid resolution $\Delta x$ and the
grid size $N_{\mathrm{grid}}$ change over models, as well as the quality of
the ABC. In AE+AD training phase an input signal contains Gaussian random
noise with variance of $\sigma^{2}_{e}$, which is added to the measurement
signals generated by the accurate and surrogate models. On the other hand, a
desired output signal does not contain noise. In training of the ANN, noisy
accurate simulated signal packed with AE is used along with the desired
noiseless surrogate model signal. In reconstruction and test phases, the
”true” simulated signal with added noise is utilized.

Table 1: Parameters of the models used in the study. The notations are explained in the text. Signals generated by accurate model and surrogate model are utilized, depending on the case, with or without simulated measurement noise. ABC layer thickness is the minimum thickness of the absorbing layer near the boundary, and $\lambda$ is the wavelength at the frequency of $50000\,\rm{Hz}$ with the SOS of $1500\,\rm{m/s}$. | Simulated physical | Accurate | Surrogate
---|---|---|---
| model | model | model
| $\hat{y}=\hat{A}(\hat{v})+\hat{e}$ | $\bar{y}=\bar{A}(\bar{v})+\bar{e}$ | ${y}={A}({v})+{e}$
$\Delta x$ | $0.2557\,{\textrm{mm}}$ | $0.2637\,{\textrm{mm}}$ | $0.8398\,{\textrm{mm}}$
$N_{\mathrm{grid}}$ | $1056\times 1056$ | $1024\times 1024$ | $256\times 256$
ABC layer | $35\,{\textrm{mm}}$ | $35\,{\textrm{mm}}$ | $7.5\,{\textrm{mm}}$
thickness | $1.17\lambda$ | $1.17\lambda$ | $0.25\lambda$
$\sigma^{2}_{e}$ | $0.0006^{2}$ | $0.0005^{2}\,/\,-$ | $0.0005^{2}\,/\,-$

  1. 1.

Generate a high quality SOS field $\bar{v}$ with the resolution of $1024\times
1024$.

     * •

Background $c(x)$ is Gaussian, $E[c(x)]=1500\,\rm{m/s}$, cutting the values of
$c(x)$ in $1100\,\rm{m/s}\dots\,1800\,\rm{m/s}$, standard deviation of
background SOS, $\sigma_{\mathrm{background}}=91\,\rm{m/s}$.

     * •

On the background, three circular shaped inclusions are inserted with constant
SOS values $c_{i},\ i=1,\ldots,3$, where $c_{i}$ is Gaussian with
$E[c_{i}]=1500\,\rm{m/s}$ and $\sigma_{c_{i}}=50\,\rm{m/s}$.

  2. 2.

Compute a low resolution SOS field. ${v}=P\bar{v}$ with the resolution of
$256\times 256$.

  3. 3.

Compute the simulated signal with the accurate model,
$\bar{y}=\bar{A}(\bar{v})+\bar{e}$ with the source term.

$f(t)=\exp\left(\xi(t-t_{0})^{2}\right)\,\sin(2\pi f_{0}t),$

where $\xi=-10^{10}\,$s-2, $f_{0}=50000\,$Hz, and $t_{0}=0.000025\,$s.

  4. 4.

Compute the simulated signal with the surrogate model, ${y}={A}({v})$, without
any noise.

  5. 5.

$\bar{y}$ is used as an input and ${y}$ is used as the desired output when
training the ANN approximating $\varphi$.

Algorithm 2 Generating a training/validation dataset item

###  5.2 Test dataset

Twenty test cases were generated to validate the method. Simulated measurement
signals were computed using the simulated physical model
$\hat{y}=\hat{A}(\hat{v})+\hat{e}$. The SOS fields used in the test cases are
significantly different from the SOS fields used at the training set, see
Figure 4. For example, SOS fields in the validation set contain star-shaped
inclusions while in the training set there are only circular inclusions. The
background variance of the sets are different, in test dataset
$\sigma_{\mathrm{background}}=23\,\rm{m/s}$. In addition, the variance of the
measurement error is different in the training dataset.

Figure 4: Examples of the used SOS fields. The upper row contains SOS fields
drawn from the training set, and the lower row contains SOS fields drawn from
the validation set.

###  5.3 Results

In the NNAEE method, we trained the ANN to obtain a mapping $\phi$ that
converts the (simulated) measurement signals to corresponding signals of the
surrogate model. Figure 5 shows two examples of the performance of the
mapping.The mean value of the error between the true signal and the surrogate
model signal is 17.75, whereas when using the ANN converted signal, the
average error goes down to 1.2.

Examples of the reconstructions and the corresponding error images are
presented in Figure 6 and the detailed results of the errors can be found in
Figure 7 and Table 2. The method is seen to improve the quality of the
reconstruction significantly, even with a relatively small size training
dataset. For example, with the training set of size 200, the reconstruction
RMS-error decreases from over $49.5$, corresponding to a conventional
reconstruction method (CRM), down to below $9.8$, which is about two times the
error of the reconstruction using the accurate model. On the other hand,
growing the training set does not improve the results after above $5000$
training samples. Using a training dataset of size $5000$ decreases the RMS-
error to less than $6.9$, however, even with $50000$ training samples the
error stays around $6.8$.

Using the method does not increase the reconstruction time compared to CRM.
The reconstruction time of the NNAEE and CRM is approximately $25\,\rm{min}$,
which is significantly lower than the reconstruction time when using the
accurate forward model, approximately $4\,\rm{h}\ 30\,\rm{min}$. On the other
hand, the proposed method can generate reconstructions of similar quality as
the standard ANN-based inversion with a significantly smaller training
dataset. For example, the RMS-error of the NNAEE reconstruction with
$N_{\mathrm{train}}=200$ is $9.76$, while the standard ANN inversion requires
almost $N_{\mathrm{train}}=10000$ to reach the same quality.

Figure 5: Simulated measurement signals using different forward models. Dotted-dashed lines represent the simulated measurement signal without measurement noise, and dotted lines show the signals generated by the surrogate model with the same SOS field. Solid lines, following almost exactly the surrogate model signal, show the ANN converted signal, which is used for computing the reconstruction. On the upper figure, the transmitter and receiver are at $90^{\circ}$ angle with respect to each other, and on the lower figure, the angle is $180^{\circ}$. Figure 6: Target is located on left, and reconstructions on two upper rows from left to right are computed with accurate model, NNAEE, with $N_{\mathrm{train}}=5000$, NNAEE, with $N_{\mathrm{train}}=200$, ANN based inversion, with $N_{\mathrm{train}}=50000$, and CRM. At lower rows there are the corresponding error fields. Table 2: Average RMS-errors of the inversions using NNAEE and ANN inversion as a function of $N_{\mathrm{train}}$. The error decreases significantly first, and converges to the value of approximately 6.80 when $N_{\mathrm{train}}$ increases. Error of the reconstruction using the accurate model, $\bar{A}(\bar{v})$, is presented here as a lower limit reference value, which is possible to reach with a perfect $\phi$. In addition, error of the reconstruction using the surrogate model $A(v)$, is presented as the upper limit reference value. $N_{\mathrm{train}}$ | NNAEE error | ANN error
---|---|---
200 | 9.76 | 21.63
400 | 9.55 | 17.97
800 | 8.61 | 15.77
1200 | 8.13 | 15.25
1800 | 7.32 | 13.66
2500 | 6.99 | 13.63
5000 | 6.89 | 12.09
10000 | 6.89 | 9.70
25000 | 6.80 | 7.92
50000 | 6.80 | 8.02
Accurate | 5.03 |
Surrogate | 49.54 |
Figure 7: Average RMS error of the NNAEE inversion, ANN based inversion, and
the inversion using the accurate numerical model as a function of
$N_{\mathrm{train}}$. The error of the NNAEE method is plotted with solid
line, the error using the ANN inversion with dashed line, and the error using
the accurate model with dotted line.

##  6 Conclusions

In the current paper, we propose a method utilizing an ANN to compensate for
the influence of modelling error, when using computationally light surrogate
model in the reconstruction of SOS field in acoustic inverse problems. The
method was tested with 20 simulated measurements using different number of
training data sets and the reconstructions were compared to an accurate
reconstruction, CRM as well as to a conventional ANN inversion. The results
show that the proposed NNAEE method can produce almost as good reconstructions
as with the accurate model, while reducing the required reconstruction times
significantly. In addition, its was shown that the proposed approach can work
with fairly small number of training data sets. Increasing the size of the
training set up to a certain limit, did not improve the quality of the
reconstruction further.

Further study is required for finding optimal hyper parameters of the ANN
presenting $\varphi$. In addition, comparing the proposed method to other
approximation error methods, like BAE, will be of interest.

Even though, the method was demonstrated over acoustic inverse problem, we
believe that it can be successfully utilized in other types of inverse
problems.

## Acknowledgments

This research was supported by the Academy of Finland (the Finnish Center of
Excellence of Inverse Modeling and Imaging, project 312344) and Academy of
Finland project 321761.

## References

  * [1] A. H. Andersen and A. C. Kak.  Simultaneous algebraic reconstruction technique (sart): A superior implementation of the art algorithm.  Ultrasonic Imaging, 6(1):81–94, 1984.
  * [2] E. Angelini, G. di Tollo, and A. Roli.  A neural network approach for credit risk evaluation.  The Quarterly Review of Economics and Finance, 48(4):733–755, 2008\.
  * [3] A. Bate, M. Lindquist, I. R. Edwards, S. Olsson, R. Orre, A. Lansner, and R. M. De Freitas.  A bayesian neural network method for adverse drug reaction signal generation.  European Journal of Clinical Pharmacology, 54(4):315–321, 1998\.
  * [4] J.-P. Berenger et al.  A perfectly matched layer for the absorption of electromagnetic waves.  Journal of Computational Physics, 114(2):185–200, 1994.
  * [5] M. J. Bianco, P. Gerstoft, J. Traer, E. Ozanich, M. A. Roch, S. Gannot, and C.-A. Deledalle.  Machine learning in acoustics: Theory and applications.  The Journal of the Acoustical Society of America, 146(5):3590–3628, 2019.
  * [6] N. N. Bojarski.  The k-space formulation of the scattering problem in the time domain.  The Journal of the Acoustical Society of America, 72(2):570–584, 1982.
  * [7] N. N. Bojarski.  The k-space formulation of the scattering problem in the time domain: An improved single propagator formulation.  The Journal of the Acoustical Society of America, 77(3):826–831, 1985.
  * [8] C. H. Chapman, J. W. D. Hobro, and J. O. A. Robertsson.  Correcting an acoustic wavefield for elastic effects.  Geophysical Journal International, 197(2):1196–1214, 03 2014.
  * [9] I. N. Daliakopoulos, P. Coulibaly, and I. K. Tsanis.  Groundwater level forecasting using artificial neural networks.  Journal of Hydrology, 309(1-4):229–240, 2005.
  * [10] A. Fichtner, H.-P. Bunge, and H. Igel.  The adjoint method in seismology: I. theory.  Physics of the Earth and Planetary Interiors, 157(1):86–104, 2006\.
  * [11] A. Fichtner, H.-P. Bunge, and H. Igel.  The adjoint method in seismology: Ii. applications: traveltimes and sensitivity functionals.  Physics of the Earth and Planetary Interiors, 157(1):105–123, 2006\.
  * [12] Y. Goldberg.  A primer on neural network models for natural language processing.  Journal of Artificial Intelligence Research, 57:345–420, 2016.
  * [13] R. Gordon, R. Bender, and G. T. Herman.  Algebraic reconstruction techniques (art) for three-dimensional electron microscopy and x-ray photography.  Journal of Theoretical Biology, 29(3):471–481, 1970.
  * [14] T. Hagstrom and T. Warburton.  A new auxiliary variable formulation of high-order local radiation boundary conditions: corner compatibility conditions and extensions to first-order systems.  Wave motion, 39(4):327–338, 2004.
  * [15] T. Hagstrom and T. Warburton.  Complete radiation boundary conditions: minimizing the long time error growth of local methods.  SIAM Journal on Numerical Analysis, 47(5):3678–3704, 2009.
  * [16] B. Hanin.  Universal function approximation by deep neural nets with bounded width and relu activations.  Mathematics, 7(10), 2019.
  * [17] G. E. Hinton and R. S. Zemel.  Autoencoders, minimum description length and helmholtz free energy.  In Advances in Neural Information Processing Systems, pages 3–10, 1994.
  * [18] P. Huthwaite and F. Simonetti.  High-resolution imaging without iteration: A fast and robust method for breast ultrasound tomography.  The Journal of the Acoustical Society of America, 130(3):1721–1734, 2011.
  * [19] M. Israeli and S. A. Orszag.  Approximation of radiation boundary conditions.  Journal of computational physics, 41(1):115–135, 1981.
  * [20] J. Kaipio and E. Somersalo.  Statistical and Computational Inverse Problems, volume 160.  Springer Science & Business Media, 2006.
  * [21] J. P. Kaipio and V. Kolehmainen.  Bayesian Theory and Applications, chapter 32 Approximate marginalization over modelling errors and uncertainties in inverse problems.  OUP Oxford, 2013.
  * [22] A. C. Kak.  Computerized tomography with x-ray, emission, and ultrasound sources.  Proceedings of the IEEE, 67(9):1245–1272, 1979.
  * [23] D. Komatitsch, D. Michéa, G. Erlebacher, and D. Göddeke.  Accelerating a 3D finite-difference wave propagation code by a factor of 50 and a spectral-element code by a factor of 25 using a cluster of GPU graphics cards.  In EGU General Assembly Conference Abstracts, EGU General Assembly Conference Abstracts, page 7462, May 2010.
  * [24] M. A. Kramer.  Nonlinear principal component analysis using autoassociative neural networks.  AIChE journal, 37(2):233–243, 1991.
  * [25] T. Lähivaara, L. Kärkkäinen, J. M. Huttunen, and J. S. Hesthaven.  Deep convolutional neural networks for estimating porous material parameters with ultrasound tomography.  The Journal of the Acoustical Society of America, 143(2):1148–1158, 2018.
  * [26] T. Lähivaara, M. Malinen, J. P. Kaipio, and T. Huttunen.  Computational aspects of the discontinuous Galerkin method for the wave equation.  Journal of Computational Acoustics, 16(04):507–530, 2008.
  * [27] Y. Li and W. Ma.  Applications of artificial neural networks in financial economics: A survey.  In 2010 International Symposium on Computational Intelligence and Design, volume 1, pages 211–214, 2010.
  * [28] L. R. Lines, R. Slawinski, and R. P. Bording.  A recipe for stability of finite-difference wave-equation computations.  Geophysics, 64(3):967–969, 1999.
  * [29] Q.-H. Liu and J. Tao.  The perfectly matched layer for acoustic waves in absorptive media.  The Journal of the Acoustical Society of America, 102(4):2072–2082, 1997.
  * [30] Z. Lu, H. Pu, F. Wang, Z. Hu, and L. Wang.  The expressive power of neural networks: A view from the width.  In I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, editors, Advances in Neural Information Processing Systems 30, pages 6231–6239. Curran Associates, Inc., 2017.
  * [31] S. Lunz, A. Hauptmann, T. Tarvainen, C.-B. Schönlieb, and S. Arridge.  On learned operator correction in inverse problems.  SIAM Journal on Imaging Sciences, 14(1):92–127, 2021.
  * [32] C. W. Manry Jr and S. L. Broschat Jr.  FDTD simulations for ultrasound propagation in a 2-D breast model.  Ultrasonic Imaging, 18(1):25–34, 1996.
  * [33] J. Moody.  Economic forecasting: Challenges and neural network solutions.  In Proceedings of the International Symposium on Artificial Neural Networks, pages 789–800, 1995.
  * [34] T. Ngo Trong, J. Mehtonen, G. González, R. Kramer, V. Hautamäki, and M. Heinäniemi.  Semisupervised generative autoencoder for single-cell data.  Journal of Computational Biology, 2019.
  * [35] R.-E. Plessix.  A review of the adjoint-state method for computing the gradient of a functional with geophysical applications.  Geophysical Journal International, 167(2):495–503, 2006.
  * [36] W. Rawat and Z. Wang.  Deep convolutional neural networks for image classification: A comprehensive review.  Neural Computation, 29(9):2352–2449, 2017.
  * [37] J. Rose and J. M. Richardson.  Time domain born approximation.  Journal of Nondestructive Evaluation, 3(1):45–53, 1982.
  * [38] O. Roy, I. Jovanovic, A. Hormati, R. Parhizkar, and M. Vetterli.  Sound speed estimation using wave-based ultrasound tomography: theory and gpu implementation.  In Proceedings of the SPIE Medical Imaging, number LCAV-CONF-2010-002. Spie-Int Soc Optical Engineering, Po Box 10, Bellingham, Wa 98227-0010 Usa, 2010.
  * [39] J.-F. Semblat, L. Lenti, and A. Gandomzadeh.  A simple multi-directional absorbing layer method to simulate elastic wave propagation in unbounded domains.  International Journal for Numerical Methods in Engineering, 85(12):1543–1563, 2011.
  * [40] A. Tarantola.  Theoretical background for the inversion of seismic waveforms, including elasticity and attenuation.  In Scattering and Attenuations of Seismic Waves, Part I, pages 365–399. Springer, 1988.
  * [41] B. B. Thompson, R. J. Marks, M. A. El-Sharkawi, W. J. Fox, and R. T. Miyamoto.  Inversion of neural network underwater acoustic model for estimation of bottom parameters using modified particle swarm optimizers.  In Proceedings of the International Joint Conference on Neural Networks, 2003., volume 2, pages 1301–1306. IEEE, 2003.
  * [42] P. Vincent.  A connection between score matching and denoising autoencoders.  Neural Computation, 23(7):1661–1674, 2011.
  * [43] P. Vincent, H. Larochelle, Y. Bengio, and P.-A. Manzagol.  Extracting and composing robust features with denoising autoencoders.  In Proceedings of the 25th International Conference on Machine Learning, pages 1096–1103, 2008.
  * [44] D. Yang, M. Lu, R. Wu, and J. Peng.  An optimal nearly analytic discrete method for 2d acoustic and elastic wave equations.  Bulletin of the Seismological Society of America, 94(5):1982–1992, 2004.
  * [45] M. E. Zakharia, P. Chevret, and F. Magand.  Estimation of shell characteristics using time-frequency patterns and neural network.  In 1996 IEEE Ultrasonics Symposium. Proceedings, volume 1, pages 713–716. IEEE, 1996.
  * [46] B. Zoph, V. Vasudevan, J. Shlens, and Q. V. Le.  Learning transferable architectures for scalable image recognition.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 8697–8710, 2018.


# Machine learning and predicting the time dependent dynamics of
local yielding in dry foams

Leevi Viitanen  Jonatan R. Mac Intyre  Juha Koivisto  Antti Puisto  Mikko
Alava  Aalto University, School of Science, Department of Applied Physics,
P.O.B 11100, 00076 Aalto, Finland

###### Abstract

The yielding of dry foams is enabled by small elementary yield events on the
bubble scale, “T1”s. We study the large scale detection of these in an
expanding 2D flow geometry using artificial intelligence (AI) and nearest
neighbour analysis. A good level of accuracy is reached by the AI approach
using only a single frame, with the maximum score for vertex centered images
highlighting the important role the vertices play in the local yielding of
foams. We study the predictability of T1s ahead of time and show that this is
possible on a timescale related to the waiting time statistics of T1s in local
neighborhoods. The local T1 event predictability development is asymmetric in
time, and measures the variation of the local property to yielding and
similarly the existence of a relaxation timescale post local yielding.

††preprint: APS/123-QED

##  I Introduction

Dry foams are assemblies of gas pockets separated by thin films of liquid
forming a connected polygonal film structure Weaire and Hutzler (1999) at a
configuration globally minimizing the surface energy Brakke (1992). A finite
external stress is needed to make foams flow and, due to the absence of
thixotropy, foams are usually rheologically characterized as simple yield
stress fluids Moller  _et al._ (2009). The steady-state, global flow curves of
foams are considered to be typical examples of Herschel-Bulkley fluid-like
behavior Ovarlez  _et al._ (2013); Barry  _et al._ (2010).

At bubble scale, the viscoplastic flow of dry foams is enabled by small
elementary topological yield events referred to as T1’s and T2’s, analogous
for instance to shear transformation zones (STZs) in amorphous solids Falk and
Langer (1998); Manning and Liu (2011). The T2 events involve the disappearance
of three-sided bubbles, while the T1 events refer to a neighbour swap between
four bubbles. Both events enable the system to jump from one metastable
surface energy minimum to another, including a local relaxation of the stored
elastic energy of the foam Petit  _et al._ (2015); Cohen-Addad  _et al._
(2013). T2 events require either film breakup or gas exchange to occur, making
them less frequent compared to T1s in systems under continuous deformation. A
quasi-2D setting or simplifying geometry, is often used to study the flow of
foams and T1s therein Dollet and Graner (2007); Van Hecke (2010); Chevalier
_et al._ (2017); Viitanen  _et al._ (2019a, b). The T1 events can be either
reversible or non-reversible depending on the geometric configuration and the
stress direction Lundberg  _et al._ (2008a); Keim and Arratia (2013).

Figure 1: Radially symmetric 2D Hele-Shaw cell creates an expanding flow
field. a) An angled photograph of the Hele-Shaw cell shows dry foam with
typically hexagonal bubbles. The 15 mm inlet hose is located at the bottom
center of the device. b) Schematic side view illustrates how the fluid and air
are foamed at the inlet. c) The schematic drawing of the top view of the
circular Hele-Shaw cell is overlaid by a single raw image from the video
stream. The dashed line illustrates the analysis area that is skeletonized to
a binary image. The inlet at the center and outlet holes have diameter
$d_{h}=15$ mm shown as dashed blue circles.

In the macroscopic flow and yielding of foams the plastic deformation arises
from a complex dynamics of T1 events with localization and clustering Gopal
and Durian (1995); Durian (1997); Abd el Kader and Earnshaw (1999). Here, we
study T1 events from large-scale statistics. We first detect tens of thousands
of events from bubble raft flow dynamical data by comparing two sequential
frames using a nearest neighbor detection algorithm Dollet and Graner (2007).
Second, we feed single frames to a Convolutional Neural Network (CNN) and
successfully predict the T1 events by the film structure and extract the
essential features. We use this tool to explore T1 dynamics and the role of
the film structure by varying the region of interest (ROI). The relation of
the T1 dynamics to macroscopic flow have been studied and interestingly, the
time dependent deformation tensor Dollet (2010); Marmottant  _et al._ (2008)
and the strain rate Wang  _et al._ (2007), both dependent on the local
velocity, were found to correlate with the T1 frequency while it has been
suggested there is no correlation between the snapshot of the film structure
and the T1 frequency Dollet and Graner (2007). This would imply T1 detection
based on a single structure snapshot should fail most of the time.

We find that T1 events can be predicted by the localized changes in film
vertices, while the bubble shape itself does not contain the same information.
The local predictability evolution is shown by time-dependent AI predictions
to be asymmetric in time, and we interpret the results as a measure of the
decrease in local yield stresses and as a manifestation of the inherent local
relaxation after a T1. The foam structure exhibits weaker and stronger spots,
be it due to the network structure or film properties Denkov  _et al._
(2009a). In amorphous solids soft spots have been studied using measures for
non-affine deformation Schoenholz  _et al._ (2014), Voronoi cell anisotropy
Rieser  _et al._ (2016), and machine learning tools Cubuk  _et al._ (2015);
Sharp  _et al._ (2018) and the distribution of local yield stresses is an
inherent property of plasticity models Nicolas  _et al._ (2018) as is the
local relaxation dynamics after a yield event Rodney  _et al._ (2011).

We compare two methods of identifying T1 events from a video of a two-
dimensional flow. First, we detect tens of thousands of events from dynamical
data by comparing two sequential frames using well-established detection
methods. Second, we feed single frames to neural network and successfully
predict the T1 events by the snapshot of the film structure. Our main findings
show that T1 events are correlated with the localized change in film vertices,
while the bubble perimeter itself does not contain similar information.

##  II Methods

Fig. 1 illustrates the circular Hele-Shaw cell with a bubble inlet at the
center used in the experiments. The flow created is symmetric with a radial
expansion that ensures T1 events will be present due to shear Dollet and
Bocher (2015). In addition, due to the decreasing flow velocity towards the
edges of the cell, a single experiment inherently probes the system at a range
of shear rates. The cell diameter is $d=250$ mm and the gap height is $h=0.75$
mm. The inlet has the diameter of $15$ mm. Foaming is facilitated by reducing
the water surface tension by mixing common household dishwashing liquid
(Fairy) with the weight fraction of 1 to 39. The solution viscosity is
indistinguishable from that of water, $\eta=1~{}\mathrm{mPa}\cdot\mathrm{s}$,
confirmed using Anton Paar 302 rheometer in a bob-and-cup geometry in the
shear rate range of $\dot{\gamma}=0.1$ 1/s to $\dot{\gamma}=100$ 1/s.

Both, the fluid and air, are driven into the cell using a constant pressure
set by two SMC ITV0010 pressure controllers and a manual flow limiter. The
inlet flow rates are set to give an average bubble velocity of $v=1.0$ mm/s.
The air is injected through a 26-gauge needle. The outer rim outlets are open
to the ambient pressure. The imaging is implemented using a Canon M3 camera,
which takes compressed video at $1920\times 1080$ pixels using the h.264
encoding resulting in image resolution corresponding to 0.09 mm/pixel. The
resulting exposure time of one frame is 40 ms, which is significantly shorter
than the temporal separation between two sequential T1 events _at the same
node_. Fig. 1c) displays an example of a grayscale snapshot overlaid to the
schematic top view of the device. The video is then interpolated to 13 000
grayscale images at 10 Hz frequency , resulting on average 0.1 s time between
subsequent frames, and eventually cropped to a square corresponding to the
region of interest. The cropped square can be then used as a grayscale image
as is or it can be processed to a skeletonized binary image with exactly one
pixel wide lines separating bubbles from each other. The bubbles are thus
identified from images as connected 2D regions with an area $A$. Their size
distribution using effective radius $R_{eff}=\sqrt{A/\pi}$ is illustrated in
Fig. 2 highlighting the monodispersity of our foam.

In the skeletonized image the nearest neighbors of a bubble are exactly 1
pixel apart. This information is further utilized in identifying the time and
location of the T1 events: If there is a change in nearest neighbors, there is
either one or more T1 or T2 events. Here, we did not observe any T2 events as
no significant coarsening occurs at the time-scales during which the bubbles
occupy the imaging area. The raw images and the identified locations of the T1
events are then divided into training and test sets as detailed in
Supplementary Tables 1 and 2 sup .

Figure 2: The foam studied is rather mono disperse having the average
effective radius $R_{eff}=2.0\pm 0.1$ mm, where the error is the standard
deviation. The histogram contains bubble sizes from three snapshots separated
by 200 s with data from 378 bubbles. Figure 3:  a) Architecture of neural
network (illustrated using a tool from LeNail ) consists of three subsequent
convolutional layers and subsampling layers and finishes with two fully
connected layers. b) Example of a raw grayscale image with one vertex that
produces a T1 event highlighted with a red square. The inset shows one of the
3 by 3 convolution filters of a fully trained CNN. The dark diagonal of the
filter captures the similarly aligned films of the raw images. c) The
convolution of the filter and the raw image illustrates the captured elements
of the raw image. The convolution is thresholded for clarity.

We pose the following question: given a single frame, are we able to predict
the occurrence of a T1 event in the subsequent frame or frames? We start with
a convolutional neural network (CNN) Bishop (2006); Salmenjoki  _et al._
(2018) and modify it to accept grayscale and skeletonized image as an input.
Fig. 3a) depicts a schematic illustration of the CNN application leading to a
single binary classification: either T1 event occurs after $\Delta t$ seconds
or not. The neural network convolves and coarse grains the input image in
three layers. The values are then run through fully connected neural network
with two hidden layers until one value is left that is thresholded to produce
the prediction. In the training, we use Adam optimizer Kingma and Ba (2014).
In total 10000 training steps are used with learning rate of $5\cdot 10^{-5}$,
and the prediction results were found to saturate around 2000 training steps.
The network was trained using sample batches at each step, rather than using
the whole data set at every step. The batch size was 1000 samples and it was
tested that increasing the batch size did not improve the prediction results
anymore. An example of a grayscale input and a trained filter is depicted in
Fig. 3b) and the first convolution layer is depicted in Fig. 3c). Here, the
filter is from a fully trained network. Initially the filters contains
randomly chosen values that converge to ones capturing essential features by
reinforcing the features producing correct prediction using back propagation
during training. The other filters capture different orientations of the
films, while the function of some filters are not clear.

The input data is split into two, roughly equally large, sets for training and
testing. In the training, a sequence of small regions of images are used. The
essential structural features for T1 detection and prediction may be studied
by limiting the amount of input data given to the network, here by limiting
the size and location of regions of interest (ROI) around the possible T1
event. First, we generate a training set consisting of small ROIs that precede
a T1 event in the next frame. Then, we balance the data set by adding randomly
picked non-overlapping regions which do not precede a T1 event until the data
set contains 45-55 % portions of each ROI type. This results in data sets that
contain roughly equal amount of both sample types so the AI does not benefit
by favoring either of the outcomes. This prescreening is necessary since T1
events are rare. Without it the algorithm scores high by always predicting “no
T1”. The quality of the prediction is measured with a score parameter $\xi$
that is the ratio of correct predictions (T1 or no T1) per all predictions.
Some examples of input images given for CNN are shown in Fig. 4 with the
outcomes of the prediction. In Fig. 4, the left column has samples that do not
precede a T1 event within 0.1 second time frame while the samples on the right
column precede a T1 event. The figure demonstrates that better accuracy of the
predictions can be reached by using the vertex centered sample set than by
using the bubble centered sample set. The CNN trained with bubble centered
samples highlights features such as bubble shapes and stable vertices that are
rather irrelevant for T1 prediction.

Figure 4: Examples of two different kind of samples inputted for the CNN and
the prediction results. On the left are negative samples that do not precede a
T1 event within 0.1 seconds and on the right are positive samples which
precede a T1 event. The sample types are starting from the first row: vertex
centered grayscale image and bubble centered grayscale image. The texts
“Correct” and “Incorrect” in the figure indicate the success of the CNN
prediction. Figure 5: Comparing the score parameter $\xi$ to the size of
regions of interest reveals that for vertex centered images the score
parameter is virtually independent of the ROI size. However, for bubble
centered images the accuracy of the algorithm increases with the size of the
region of interest ($L\times L$ mm) until the area is at least equivalent to
that of four bubbles. There the performance is similar for both the grayscale
(filled symbols) and skeletonized (unfilled symbols) images with the same
centering. For the largest ROI, the algorithm works best for skeletonized
network (unfilled circles) reaching almost 90 % accuracy. The result is worse
for grayscale images (blue) or if the data is reduced down to a single bubble
shape (inset in b), open diamonds).

##  III Prediction results

###  III.1 Machine learning

Next, we explore the T1 predictions produced by the CNN. For the purpose we
use the score parameter defined earlier. We study the CNN for five different
combinations of locations and grayscale or skeletonized frames for the full
range of ROIs. The combinations selected allows us to compare i) the effect of
film or vertex thicknesses by comparing the predictions produced using
grayscale (filled symbols) and skeletonized images (unfilled) or ii) the
effect of focusing on the vertices (triangles) or bubble centers (circles) or
iii) concenrtrating on the shape of the bubble perimeter without its neighbors
(diamonds). These choices of ROIs allow us to explore the image features, such
as bubble shape, liquid fraction, or neighbor bubbles, by which the CNN is
able to produce the best predictions. Fig. 5b) shows the score parameter $\xi$
for different sizes ($4\leq L\leq 13$ [mm]) and locations (bubble or vertex
centered) of ROIs.

The best score $\xi=98$ % is obtained using vertex centered ROIs of $4\times
4$ mm illustrated by the red triangle and a corresponding example input in
panel 5a). As a sanity check, we get the same score $\xi_{+}=98$ % if we
restrict to true positives by defining $\xi_{+}$ as correctly predicted T1
events per T1 predictions only. This indicates that there is no bias in the
classification. Increasing the size of the ROI does not improve the score
$\xi$ and thus we conclude that the information in the local surroundings of
the vertex provides a good indicator of a T1 event (filled triangles).

We considered that the change in the local liquid fraction is the reason for
lower yield point and T1 events, what would appear as darker nodes. Thus,
removing the information about the local liquid fraction by skeletonizing, we
expected that the score parameter would decrease its value significantly. We
performed the same CNN analysis on the skeletonized frames, which result
roughly the same scores with the gray scale images as witnessed in Fig. 5b)
(unfilled triangles). This is indicative that the CNN is not in fact capturing
the gray scale levels in the images, but rather predicts solely based on the
orientations of the films surrounding the vertices. In other words, local
liquid fraction or liquid motion in the films plays no significant role in the
local yielding. Similarly to the gray scale data, the skeletonized data shows
no significant dependence on the ROI size.

Moving the ROI center from the vertex to the bubble center of volume has a
significant impact on the score parameter. First, let us concentrate on the
skeletonized bubble shape (black diamonds and inset in Fig. 5b)) all the
information of the neighboring bubbles is removed from ROIs. The score is
close to $\xi=65$ % with only weak increase with ROI size $L$. This is
significantly lower compared to the ROIs that focus on the vertex indicating
that the bubble shape does not contain enough information for accurate
prediction of T1 events. The score calculated only for the positive
predictions is even lower $\xi_{+}=60$ % showing that the CNN does a slightly
better job at predicting samples without a T1 event than predicting the
events.

Now we include the films of neighboring bubbles but keep the ROI center at the
center of bubble. This data is plotted as a function of the bubble centered
ROI size for both grayscale (filled circles) and skeletonized (unfilled
circles) images in Fig. 5b). The predictions on the smallest ROIs for
skeletonized or grayscale images gives similar scores as the predictions based
on bubble shapes considered above. For both skeletonized and grayscale ROIs,
the score $\xi$ increases with the size of the ROI. As the ROI increases to
include the entire bubble and parts of its neighbors (green circle), the score
reaches over 80 %, lower than the one that focused on the vertices, yet
significantly higher than the one excluding all the bubbles neighbors. This
supports our previous interpretation that the essential information on the T1
event is encapsulated by the structures and locations of the vertices and not
by the bubble shape.

Comparison of the pure skeletonized data with and without the neighboring
bubbles offers a visual confirmation to these observations. These structures
are shown in the inset of Fig. 5b) and highlighted by green color in Fig. 5c).
In the frame in question, a T1 event will take place in the lowest vertex, a
location that is not obvious by only looking at the bubble shape. Thus, if the
information on the neighboring bubbles is removed, the score parameter
understandably drops dramatically.

###  III.2 Comparison to established methods

In the simplest view, the CNN could be predicting the T1 events simply by
applying the Plateau’s laws. This seems to be the case in the sense, that the
CNN obtains the best predictions using the vertex centred samples and the
Plateau’s rules apply to the vertices here. As shown in Fig. 6 the number of
neighbors, on the other hand, only correlates weakly with the probability of
the occurrence of a T1 event. As the Figure shows, the most probable number of
neighbors for bubble going through a T1 event is actually six. This
observation may explain why bubble centred samples yield worse predictability
compared to vertex centered. Here, however, one must note that we have a
finite accuracy in determining if the Plateau’s rules are obeyed.
Supplementary Video sup  highlights the lifetimes of these metastable states
apparently violating the Plateau’s rule lasting up to several seconds. The
actual violation of Plateau’s rule lasts only a few milliseconds, observed in
a 3D case with high magnification Petit  _et al._ (2015). Therefore, more
rigorous statement is that the CNN bases the predictions on apparent violation
of Plateau’s rules within the spatial and temporal resolution of the
measurement.

To further confirm the irrelevance of the bubble shape, we checked for any
correlation with the textural tensor

$\hat{\ell}=\left\langle\frac{\vec{\ell}\otimes\vec{\ell}}{\ell}\right\rangle$
(1)

capturing the bubble shape Marmottant  _et al._ (2008) and the T1 event rate
which has already been used for various bubble monolayers with similar results
Dollet and Graner (2007); Dollet (2010). Here $\vec{\ell}$ represents the
vector joining two neighboring Plateau borders within a single bubble. Using a
rectangular mesh of $8\times 8$ boxes dividing each image into 64 sub-regions,
Fig. 7 shows no correlation between T1 events and the bubble shape for single
frames.

Figure 6: Histogram of number of nearest neighbors for all bubbles and bubbles
right before T1 event. Clearly the distribution is wider for the bubbles about
to have T1 event indicating that local configurations where bubble has exactly
six neighbors are more stable than configurations where bubbles have more or
less than six neighbors.

Supplementary Tab. 3 shows the score parameter $\xi$ along with other
benchmark parameters of our CNN with various sanity checks sup . These
included a virgin data set without any reduction of vertices, comparison to a
hand made algorithm and prescreening the input data to contain only vertices
involving four films, those appearing to violate Plateau’s 120 degree rule.
The last method has the best performance, reaching $\xi=99$ % score and
capturing 34 % of all the T1 events.

Figure 7: Textural tensor and T1 event rate averaged over $8\times 8$ sub-
regions of imaged area and time. Although, both quantities show variation
between different locations the quantities are not correlated.

Based on our extensive testing we conclude that there exist a significant
subset of T1 events triggered by the changes of local microstructure and film
orientations in the vertex that can be captured with a CNN without any
information on the time dependence. Thus, already a small region of interest
enables the CNN to predict the T1 events. In the future, to further improve
the predictions we will train other neural networks using parameters, which
describe features of the local bubble geometry. We suspect adding some
information about the magnitude and direction of the velocity or recent T1
history could improve the predictability.

##  IV Time-dependent properties of elementary yield events

We next study to which extent in time the T1 events are predictable. To answer
this, we study the temporal development of T1 events using the data plotted in
Fig. 8.

We identify the score parameter $\xi$ as a measure of predictability.
Essentially, predictability here means the existence of a heterogeneous
feature or “a defect” - since it leads to yielding - different from a
featureless material, allowing the CNN algorithm to make a prediction. We
perform the analysis for vertex and bubble centered grayscale samples
separately for the ROI size $9\times 9$ mm2, corresponding to sample image
highlighted green in Fig. 5 a) and c). The size is chosen as the smallest area
where the score parameter is saturated to its limiting value. Also, it is not
necessary to include analysis with skeletonized images since the
skeletonization only has small effect on the prediction score ($<5$ %) if any.

Figure 8: a) Temporal evolution of the score parameter $\xi$ as a measure for
structural inhomogeneity for different lag times between the CNN input frame
and the T1 event $\Delta t=t_{frame}-t_{event}$, for the vertex (brown
triangles) and bubble (blue circles) centered cases. The solid curves are the
exponential fits to the data before the T1 event. The same fits are also
plotted at the positive side (dashed). This shows the time reversal asymmetry
with the vertex centered case highlighted with arrows. b) The characteristic
time scale of subsequent T1 events is of the order of one second from the
waiting time distributions $P(\tau)$.

Fig. 8a) depicts predictability $\xi$ against the time difference $\Delta
t=t_{frame}-t_{event}$ between the T1 event and the single input frame used by
the CNN. This allows to examine the structural signatures of a T1 event before
($\Delta t<0$), or even after ($\Delta t>0$) the event. The brown triangles
refer to data set with vertex centered ROI, where the CNN is retrained and
evaluated for each measurement point. Here, the vertex is always at the center
of the ROI even if the lag $\Delta t$ is increased.

We find that the prediction probability follows an exponential fit (solid
line) and drops to one half at $\Delta t=-1.9$ s (before) the event. Using the
ROIs from frames after the T1 event the probability drops to one half faster
at $\Delta t=1.5$ s indicating that there is an asymmetry in predictability.
Interestingly, a similar breakdown of time reversibility related to local
geometry of T1 events has been reported previously in bubble raft shearing
experiment measuring orientation of T1 events Wang  _et al._ (2007). The
asymmetry seems reasonable as shear rate drives the film shrinkage while force
balance at the liquid-air interfaces drives the film growth for small shear
rates (see supplemental video sup ) Biance  _et al._ (2009). Although, this
picture is only crude simplification as evident based on the video and the
previous studies T1 events also show clustering due to redistribution of
stress Desmond and Weeks (2015) causing non-local deformation in the foam
Goyon  _et al._ (2008); Elias  _et al._ (1999).

The bubble centered predictions are symmetric in $\Delta t$, although the
prediction probability is not as good as in the vertex centered case. This,
might be due to the CNN focusing more on other parts than the film created
during a T1 event. Therefore, while the film instability triggers the event,
the local neighborhood still plays a significant role. This can be understood
in the context of slow energy dissipation, where a relaxation of the energy
landscape enables reversible T1 events Lundberg  _et al._ (2008a). Thus, the
configuration does not essentially change in the scale of few bubbles, even if
the local film geometry can become completely different.

Fig. 8b) shows the T1 event waiting time distributions using the inter-event
arrival time in a restricted area of $50\times 50$ mm. The data based on all
the CNN prediction algorithms agree with that obtained using the nearest
neighbor detection algorithm capturing all the events (red triangles).
Moreover, the data shows that the typical T1 time scales are short, of the
order of one second in agreement with Fig. 8a). To compare with other dry
foams with synthetic surfactants that do not affect the viscosity of the
carrier liquid Denkov  _et al._ (2009b) the time scales for stress jumps was
measured $t_{\mathrm{jump}}=1.5\pm 0.5$ s Cantat and Pitois (2006) and the
average time to reach 90 % of the final film length was measured
$t_{\mathrm{90}}=0.5$ s Durand and Stone (2006), that is in both cases of the
same order of magnitude as here suggesting that film growth is driven by
interfacial tension. With the average bubble velocity of $v=1.0$ mm/s the
timescales are in the ballpark of the average bubble diameter of $d=4$ mm
making T1 events highly localized.

##  V Conclusions

We have successfully trained a general purpose CNN to recognize the neighbor
swap T1 events in radial 2D foam flow using only snapshots of the structure
with no time dependent information. We capture the essential features of the
images, namely vertices and film orientations. Using these features we show
that typically T1 events initiate from the unstable vertices that appear to
violate Plateau’s rules while the bubble shape is a less relevant quantity
(Fig. 5). This highlights the importance of local film geometry and
microstructure in rearrangements similarly to bubble coalescence Forel  _et
al._ (2019). The development of the shape or perimeter of a bubble and its
neighborhood is symmetric in time for the bubble participating to a T1 event.
In contrast, the changes in films are asymmetric in time. Here, the emphasis
is on the neighbor swap aspect, the bubble is still a relevant unit for other
processes such as coarsening Khakalo  _et al._ (2018) and recoverable elastic
response Lundberg  _et al._ (2008b). The elastic energy stored in the system
does not vanish instantly; major part of it is stored in a different film
configuration making it possible to have reversible local rearrangements.

On average, a vertex about to yield looks very different before than after the
event, as seen for example in Supplementary Video 1, that is, one can see
whether the configuration under scrutiny is close to the local yield threshold
or not sup . This manifests as temporal asymmetry in the score $\xi$ that
drops relatively fast to the baseline after the event. In other amorphous
solids e.g. granular pillars Cubuk  _et al._ (2015) the local configuration,
namely the particle contacts correlate with the yield stress and exhibit local
variations (soft spots).

Our work has focused on the experimentally amenable 2D case  considering the
case of constant driving pressure and liquid fraction. It would now be quite
interesting to investigate how one can change the foam yielding in terms of
the local predictability of T1 events and their spatio-temporal correlations.
This could be achieved altering the physical properties of the sample foam,
such as liquid fraction, polydispersity, geometry or shear rate. Although,
changing any of these parameters may impose difficulties in maintaining a
stable flow and keeping coarsening negligible. One particular direction would
be particle-laden foams where the reinforcing particles would induce other
dynamical timescales. This leads to a wide variety of industrial applications
with tunable orientation dependent properties such as tensile strength or heat
conductivity Wicklein  _et al._ (2014); Yang  _et al._ (2017).

###### Acknowledgements.

JK and AP acknowledge the funding from Academy of Finland (308235 and 278367),
Business Finland (211715) and Aalto University (974109903) as well as Aalto
Science IT project for computational resources. LV acknowledges the funding
from the Vilho, Yrjö and Kalle Väisälä Foundation via personal grant and
Academy of Finland (278367). We acknowledge Henri Salmenjoki for sharing CNN
code and helpful discussions.

## References

  * Weaire and Hutzler (1999) D. L. Weaire and S. Hutzler,  _Foam collapse_ (Oxford University Press, New York, 1999).
  * Brakke (1992) K. A. Brakke, Experimental Mathematics , 141 (1992).
  * Moller  _et al._ (2009) P. Moller, A. Fall, V. Chikkadi, D. Derks, and D. Bonn, Philos. Trans. A Math. Phys. Eng. Sci. 367, 5139 (2009).
  * Ovarlez  _et al._ (2013) G. Ovarlez, S. Cohen-Addad, K. Krishan, J. Goyon, and P. Coussot, J. Non-Newtonian Fluid Mech. 193, 68 (2013).
  * Barry  _et al._ (2010) J. D. Barry, D. Weaire, and S. Hutzler, Rheologica acta 49, 687 (2010).
  * Falk and Langer (1998) M. L. Falk and J. S. Langer, Phys. Rev. E 57, 7192 (1998).
  * Manning and Liu (2011) M. L. Manning and A. J. Liu, Phys. Rev. Lett. 107, 108302 (2011).
  * Petit  _et al._ (2015) P. Petit, J. Seiwert, I. Cantat, and A.-L. Biance, Journal of Fluid Mechanics 763, 286 (2015).
  * Cohen-Addad  _et al._ (2013) S. Cohen-Addad, R. Höhler, and O. Pitois, Annual Review of Fluid Mechanics 45, 241 (2013).
  * Dollet and Graner (2007) B. Dollet and F. Graner, Journal of Fluid Mechanics 585, 181 (2007).
  * Van Hecke (2010) M. Van Hecke, Journal of Physics Condensed Matter 22, 033101 (2010).
  * Chevalier  _et al._ (2017) T. Chevalier, J. Koivisto, N. Shmakova, M. J. Alava, A. Puisto, C. Raufaste, and S. Santucci, Journal of Physics: Conference Series (2017).
  * Viitanen  _et al._ (2019a) L. Viitanen, J. Koivisto, A. Puisto, M. Alava, and S. Santucci, European Physical Journal B 14 (2019a).
  * Viitanen  _et al._ (2019b) L. Viitanen, J. Koivisto, M. Korhonen, A. Puisto, and M. Alava, Bioresources , 5716 (2019b).
  * Lundberg  _et al._ (2008a) M. Lundberg, K. Krishan, N. Xu, C. S. O’Hern, and M. Dennin, Phys. Rev. E 77, 041505 (2008a).
  * Keim and Arratia (2013) N. C. Keim and P. E. Arratia, Soft Matter 9, 6222 (2013).
  * Gopal and Durian (1995) A. D. Gopal and D. J. Durian, Phys. Rev. Lett. 75, 2610 (1995).
  * Durian (1997) D. J. Durian, Phys. Rev. E 55, 1739 (1997).
  * Abd el Kader and Earnshaw (1999) A. Abd el Kader and J. C. Earnshaw, Phys. Rev. Lett. 82, 2610 (1999).
  * Dollet (2010) B. Dollet, Journal of Rheology 54, 741 (2010).
  * Marmottant  _et al._ (2008) P. Marmottant, C. Raufaste, and F. Graner, The European Physical Journal E 25, 371 (2008).
  * Wang  _et al._ (2007) Y. Wang, K. Krishan, and M. Dennin, Philosophical Magazine Letters 87, 125 (2007).
  * Denkov  _et al._ (2009a) N. D. Denkov, S. Tcholakova, K. Golemanov, and A. Lips, Phys. Rev. Lett. 103, 118302 (2009a).
  * Schoenholz  _et al._ (2014) S. S. Schoenholz, A. J. Liu, R. A. Riggleman, and J. Rottler, Physical Review X 4, 031014 (2014).
  * Rieser  _et al._ (2016) J. M. Rieser, C. P. Goodrich, A. J. Liu, and D. J. Durian, Phys. Rev. Lett. 116, 088001 (2016).
  * Cubuk  _et al._ (2015) E. D. Cubuk, S. S. Schoenholz, J. M. Rieser, B. D. Malone, J. Rottler, D. J. Durian, E. Kaxiras, and A. J. Liu, Phys. Rev. Lett. 114, 108001 (2015).
  * Sharp  _et al._ (2018) T. A. Sharp, S. L. Thomas, E. D. Cubuk, S. S. Schoenholz, D. J. Srolovitz, and A. J. Liu, Proceedings of the National Academy of Sciences 115, 10943 (2018).
  * Nicolas  _et al._ (2018) A. Nicolas, E. E. Ferrero, K. Martens, and J.-L. Barrat, Rev. Mod. Phys. 90, 045006 (2018).
  * Rodney  _et al._ (2011) D. Rodney, A. Tanguy, and D. Vandembroucq, Modelling and Simulation in Materials Science and Engineering 19, 083001 (2011).
  * Dollet and Bocher (2015) B. Dollet and C. Bocher, The European Physical Journal E 38, 123 (2015).
  * (31) See Supplemental Material at [URL will be inserted by publisher] for a video of T1 events and further details about machine learning data sets and prediction results.
  * (32) A. LeNail, “Publication-ready nn-architecture schematics,” http://alexlenail.me/NN-SVG/LeNet.html, accessed: 7.8.2019.
  * Bishop (2006) C. M. Bishop,  _Pattern recognition and machine learning_ (springer, 2006).
  * Salmenjoki  _et al._ (2018) H. Salmenjoki, M. J. Alava, and L. Laurson, Nature Communications 9, 5307 (2018).
  * Kingma and Ba (2014) D. P. Kingma and J. Ba, arXiv preprint arXiv:1412.6980 (2014).
  * Biance  _et al._ (2009) A.-L. Biance, S. Cohen-Addad, and R. Höhler, Soft Matter 5, 4672 (2009).
  * Desmond and Weeks (2015) K. W. Desmond and E. R. Weeks, Physical review letters 115, 098302 (2015).
  * Goyon  _et al._ (2008) J. Goyon, A. Colin, G. Ovarlez, A. Ajdari, and L. Bocquet, Nature 454, 84 (2008).
  * Elias  _et al._ (1999) F. Elias, C. Flament, J. Glazier, F. Graneri, and Y. Jiang, Philosophical Magazine B 79, 729 (1999).
  * Denkov  _et al._ (2009b) N. D. Denkov, S. Tcholakova, K. Golemanov, K. P. Ananthpadmanabhan, and A. Lips, Soft Matter 5, 3389 (2009b).
  * Cantat and Pitois (2006) I. Cantat and O. Pitois, Physics of Fluids 18, 083302 (2006).
  * Durand and Stone (2006) M. Durand and H. A. Stone, Phys. Rev. Lett. 97, 226101 (2006).
  * Forel  _et al._ (2019) E. Forel, B. Dollet, D. Langevin, and E. Rio, Phys. Rev. Lett. 122, 088002 (2019).
  * Khakalo  _et al._ (2018) K. Khakalo, K. Baumgarten, B. P. Tighe, and A. Puisto, Physical Review E 98, 012607 (2018).
  * Lundberg  _et al._ (2008b) M. Lundberg, K. Krishan, N. Xu, C. S. O’Hern, and M. Dennin, Phys. Rev. E 77, 041505 (2008b).
  * Wicklein  _et al._ (2014) B. Wicklein, A. Kocjan, G. Salazar-Alvarez, F. Carosio, G. Camino, M. Antonietti, and L. Bergström, Nature Nanotechnology 10, 277 (2014).
  * Yang  _et al._ (2017) M. Yang, N. Zhao, Y. Cui, W. Gao, Q. Zhao, C. Gao, H. Bai, and T. Xie, ACS Nano 11, 6817 (2017).