# Dynamic image reconstruction in MPI with RESESOP-Kaczmarz

Marius Nitzsche  Bernadette N Hahn

(Department of Mathematics, University of Stuttgart, D-70569 Stuttgart,
Germany)

##  1 Introduction

Magnetic Particle Imaging (MPI), initially proposed by Gleich and Weizenecker
[1], is an emerging imaging modality that enables insights into a specimen by
measuring the response of superparamagnetic nanoparticles to an applied
magnetic field. To this end particles are injected with a certain tracer
concentration $c:\Omega\mapsto\mathbb{R}_{+}\cup{0}$ into the area of
interest, which is called the field of view $\Omega\subset\mathbb{R}^{3}$.
Through layering two magnetic fields, i.e. the drive field and the selction
field, a field-free point moves through the field of view along a Lissajous
trajectory. The time interval of one Lissajous trajectory is called _frame_.
The dynamic field excitation caused by the moving field-free point results in
a time-varying magnetization $M$ of the injected particles, which induces a
measurable voltage $v^{P}$ over a time interval $I:=[0,\text{T}]$ in the
receive coils. To reconstruct $c$ it is necessary to separate it from $M$ by
$M=c\bar{m}$ with the particles’ mean magnetic moment $\bar{m}:\Omega\times
I\mapsto\mathbb{R}^{3}$. Knowing this, at the core of MPI lies the fundamental
equation

$\displaystyle v^{P}(t)=$
$\displaystyle\int_{\Omega}c(x)\underbrace{(-a*\mu_{0}p^{R}(x)^{T}\frac{\partial}{\partial
t}\bar{m}(x,t))}_{=:S(x,t)}\text{d}x,\;t\in I$ (1)

where $\mu_{0}$ is the constant magnetic permeability of vacuum and
$p^{R}:\Omega\mapsto\mathbb{R}^{3}$ is the receive coil sensitivity. The
analog filter $a:[-\text{T},\text{T}]\to\mathbb{R}$ is applied to mitigate the
direct feedthrough, a voltage induced by the applied magnetic field.

In MPI, two key challenges arise. The determination of $\bar{m}(x,t)$ from
given tupels $(c,v^{P})$ is called the _calibration problem_. The standard
procedure relies on a time consuming measurement-based approach [2].
Alternatively, model-based techniques are explored which offer a less
resource-intensive path, although they may not deliver the same level of
precision, see e.g. [3].

In this research, our primary focus is on the _image reconstruction problem_.
In other words, the goal is to determine the tracer concentration
$c(x),\,x\in\Omega$ from observed measurements $v^{P}(t),\,t\in I$ with a
given system function. This process necessitates the discretization of (1),
yielding a linear system

$Ac=v,\quad A\in\mathbb{R}^{N\times M},c\in\mathbb{R}^{M},v\in\mathbb{R}^{N},$
(2)

with given system matrix $A$, measurement vector $v$ and searched-for discrete
concentration $c$. The dimensions are given by $N$, which is the number of
considered temporal points multiplied by the number of receive coils being
usually $3$, as well as the number of spatial points $M$ of the discretized
field of view.

The most common solver in MPI is the Kaczmarz method with Tikhonov
regularization, also called regularized Kaczmarz, which mainly consists of a
fixed point iteration. By iteratively projecting the approximate solution
orthogonally onto affine subspaces, the method solves the linear system in
(1). More insights can be found in the literature, see e.g. [2].

###  1.1 Dynamics in MPI

The regularized Kaczmarz method works well if the searched-for concentration
is stationary. Be that as it may, many potential clinical applications rely on
the reconstruction from dynamic MPI data, e.g. imaging blood flow [4],
tracking medical instruments [5] or monitoring of strokes [6, 7]. If the
concentration is time-dependent and changes during the data acquisition,
motion artefacts will arise in the reconstructions unless the dynamics is
compensated for within the reconstruction step. This issue is exacerbated in
the so-called multi-frame scenario, where the concentration is reconstructed
from multiple data sets that are blockwise averaged over time to improve the
signal-to-noise ratio and the spatial resolution of the image.

The regularized Kaczmarz method does not take dynamics of the concentration
into account. Figure 1 illustrates motion artefacts arising in single- as well
as multi-frame reconstructions. As example, we considered simulated noise free
data for a fast rotating cylinder, more precisely it performs one full
rotation during the acquisition of seven frames. A 2D slice of the 3D ground
truth at time instance $t=0$ is depicted in subfigure (a). Reconstruction
results from a single-frame data set with regularized Kaczmarz are shown in
subfigures (c) and (d) for different regularization parameters $\lambda$ with
the one in (d) being the optimal one. In both cases, the computed solution is
not sharp and includes motion artefacts. Even with the optimal parameter the
method is not able to provide a good approximation to the ground truth while
the averaging in the multi-frame scenario further increases the motion
artifacts, see subfigure (b).

(a) Phantom (b) $\lambda=8.5$: Averaging (c) $\lambda=2$: One frame (d)
$\lambda=8.5$: One frame

Figure 1: Reconstruction results with regularized Kaczmarz from noise free
data of a fast rotating cylinder.

Consequently, it becomes imperative to account for dynamics, i.e. to address
the image reconstruction task as a dynamic inverse problem. A universal,
standardized regularization technique for dynamic inverse problems does not
exist. While a multitude of approaches have been proposed for individual
traditional tomographic applications, only a few have been presented in the
field of MPI in response to this challenge.

For instance, a strategy to account for periodic motion by pre-processing the
measured data is developed in [8]. By estimating the motion frequency
associated with periodic motion, the authors effectively grouped data samples
and thus generated static single-frame data sets for the individual states of
the object. Consequently, an image of each state can be reconstructed with
standard static algorithms. This approach has also been extended to multi-
patch MPI [9].

A dynamic forward model within MPI is offered in [10], which is combined with
a reconstruction method using temporal splines to increase temporal resolution
in [11].

Another approach is based on estimating the motion using optical flow, a
computationally expensive method that has been examined in applications like
CT or MRI [12]. Its application in a dynamic MPI scenario was first proposed
in [13] and in more detail in [14] using a stochastic primal-dual algorithm to
jointly estimate the motion and reconstruct the image.

In this article we propose to apply the Regularized Sequential Subspace
Optimization method (RESESOP) for inexact model operators [15] to the dynamic
MPI problem. With this method we take motion implicitly as model inexactness
into account and are able to reconstruct images from motion corrupted data,
which are of higher quality than those computed by regularized Kaczmarz,
without the need of explicitly estimating the motion. We will see that the
only required a priori information can be gained from the measured MPI data
itself. Furthermore, no special pre-processing of the data beyond the standard
is required.

The article is organized as follows. In section 2 we first introduce the
mathematical model of time-dependent concentration reconstruction in MPI and
argue that dynamics induce model inexactness. We then present the general
concept of the RESESOP and transfer it to the dynamic reconstruction problem
in MPI. In particular, we discuss how the required a priori information on the
model inexactness can be extracted directly from the measured MPI data.
Finally, we provide a detailed numerical evaluation of the method. The
respective test cases, including simulated as well as real data sets, are
introduced in Section 3. In Section 4, we present respective reconstruction
results and analyze the performance of the proposed algorithm.

##  2 Methods

As illustrated in the introduction, reconstruction algorithms need to account
for dynamics of the studied concentration. In order to develop such
algorithms, the time-dependency of the concentration first needs to be
incorporated into the mathematical model of MPI, i.e. in equation (1) and (2),
respectively.

###  2.1 The mathematical model of MPI with time-dependent concentrations

Imaging a time-dependent particle concentration, i.e. $c:\Omega\times
I\mapsto\mathbb{R}_{+}\cup{0}$, implies that it is changing during the
measuring process. In the most extreme scenario, the time-scale of the
concentration dynamics coincides with the one of the data acquisition, i.e.

$v^{P}(t)=\int_{\Omega}c(x,t)S(x,t)\text{d}x,\quad t\in I$ (3)

or in the discretized setting

$A_{t}c_{t}=v_{t}\,\forall t\in[0,T],\quad A_{t}\in\mathbb{R}^{1\times
M},c_{t}\in\mathbb{R}^{M},v_{t}\in\mathbb{R}^{1}.$ (4)

In other words, each single measure $v^{P}(t)$ (respectively $v_{t}$), $t$
being one particular time instance, captures a different state of the
concentration. Even with strong a priori assumptions a dynamic concentration
reconstruction on the same temporal granularity as defined by the sampling
rate of the measurement process is hardly possible.

However, an important characteristic of MPI are the very fast measurement
times which imply that the temporal resolution of the concentration is
typically several orders coarser than the data sampling rate. Thus, it is
reasonable to assume that the searched-for concentration is piecewise constant
in time. This motivates to formally split the time interval $I$ into
$N_{\tau}$ pairwise disjoint time subintervals $I_{i}:=[\tau_{i},\tau_{i+1})$
with $i=0,1,\dots,N_{\tau}-1$ as well as $\tau_{0}:=0$ and
$\tau_{N_{\tau}}:=T$. As a consequence it holds

$\displaystyle I=\dot{\cup}_{i=0,1,\dots,N_{\tau}-1}I_{i},\quad
c(x,t)=c_{\tau_{i}}(x)\text{ for }t\in I_{i}.$

Both the time scale of the measurement and the coarser time scale of the
dynamic concentration can then be coupled by a function

$\displaystyle\gamma:I\rightarrow\\{\tau_{i},\ i=0,1,\dots,N_{\tau}-1\\},\
\gamma(t):=\tau_{i}\text{ if }t\in I_{i}.$

With this notation, the general dynamic reconstruction problem can be
formulated as

$v^{P}(t)=\int_{\Omega}c(x,\gamma(t))S(x,t)\text{d}x,\quad t\in I,$ (5)

or in discretized form

$A_{i}c_{\tau_{i}}=v_{i}\quad\forall\ i\in\\{0,1,\dots,N_{\tau}-1\\}$ (6)

with given $A_{i}\in\mathbb{R}^{(N_{i}\cdot 3)\times M}$, searched-for
$c_{\tau_{i}}\in\mathbb{R}^{M}$ and measured $v_{i}\in\mathbb{R}^{(N_{i}\cdot
3)}$, where $N_{i}$ denotes the number of time samples per sub-interval
$I_{i}$.

In Section 2.2, we will focus on solving the discretized dynamic inverse
problem. Considering each _subproblem_ $A_{i}c_{\tau_{i}}=v_{i}$ individually
would formally allow the application of a classic static reconstruction
method. However, in analogy to the multi-frame approach, the goal is to
improve the SNR and the spatial resolution of the images. In order to achieve
this in the presence of motion, correlations between the states have to be
exploited within the joint reconstruction step. Furthermore, if the length of
a subinterval $I_{i}$ is smaller than the time required to measure one
complete frame, the solution of the respective subproblem has to deal with a
limited-data problem. The actual decomposition of $I$ into subintervals will
depend on the specific application. Suitable choices are discussed in the
following subsection.

####  2.1.1 Choosing suitable subintervals

In the multi-frame scenario, an intuitive idea is to let each subproblem
coincide with one frame. Indeed, many concentration dynamics might be assumed
stationary within one frame. The actual time to acquire data for one frame is
very small, e.g. for a 3D field of view, it takes 21.54 ms. Thus, for
velocities lower than 10 centimeters per second, for instance, the travelled
distance within the 21.54 ms is smaller than the diameter of a voxel size
[16]. Additionally, the usage of frames as subproblems simplifies the
implementation within a MPI framework. Real data are usually pre-processed and
stored frame-by-frame. Thus, a dynamic reconstruction algorithm working on
frames can be embedded directly in existing MPI software pipelines.

However, if the examined structure changes faster than the acquisition of one
frame, e.g. due to a pulsating motion or a fast rotation, this rigid coupling
of time scales is not feasible and would still lead to motion artifacts. A way
to handle these scenarios is to choose smaller subproblems which are then
underdetermined. Looking at it from a practical perspective, one timestep $t$
is equivalent to $0.8\mu\text{s}$ in MPI. Thus, choosing subproblems
corresponding to either one frame, half a frame or a quarter of a frame seem
adequate to represent and capture most relevant concentration dynamics.

####  2.1.2 Interpreting the dynamics as model inexactness

A common strategy in dynamic imaging problems is the incorporation of motion
models to relate the different states of the object to each other, see e.g
[[17],, [18]]. In our case, this means that instead of considering a series of
independent concentrations $c_{\tau_{i}}$, it is assumed that they are all
linked together by an underlying motion model $\Gamma$, for instance

$\displaystyle c(x,t)=c_{0}(\Gamma_{t}(x))$

with a reference configuration $c_{0}$, e.g. the concentration at the inital
time of the data acquisition. Incorporating such a motion model into the
equation (5) and using an appropriate change of coordinates yields an inverse
problem for the reference concentration $c_{0}$ with a forward operator
depending on the motion information $\Gamma$, i.e. in the discretized setting

$A_{i,\Gamma}c_{0}=v_{i}\quad\forall\ i\in\\{0,1,\dots,N_{\tau}-1\\}$ (7)

with $A_{i,\Gamma}\in\mathbb{R}^{(N_{i}\cdot 3)\times M}$.

A detailed derivation in the general context of inverse problems can be found
in [17]. Thus, extracting the time-dependent concentrations $c(\cdot,t)$ is
equivalent to extracting the motion information $\Gamma$ and the reference
concentration $c_{0}$.

In general, the exact motion $\Gamma$ will be unknown, i.e. in practice, one
can only use an approximate forward operator to determine $c_{0}$, e.g. the
static model $A_{i}$ if no further information on the dynamics are available.
Therefore, we treat the dynamics in the following as inexactness in the
forward model and account for this inexactness while recovering $c_{0}$ from
the measured dynamic data $v_{i},\ i=0,1,\dots,N_{\tau}-1$.

###  2.2  RESESOP method for solving inverse problems with inexact forward
operator

In this work we propose to solve the dynamic reconstruction problem in MPI (7)
with the RESESOP-Kaczmarz method developed in [15] for solving general inverse
problems with inexact forward operator.

This method is based on the Sequential Subspace Optimization method (SESOP)
[19]. It is an iterative method that approximates the searched-for solution
$f$ of a given linear inverse problem $Bf=g$ by the following scheme: In the
$n$-th step, compute

$f_{(n+1)}=P_{H}(f_{(n)})$

with the metric projection $P_{H}$ onto an intersection $H=\cap_{j\in
J_{n}}H(u_{j},\alpha_{j})$ with hyperplanes

$H(u,\alpha)=\\{x\in X:\langle u,x\rangle=\alpha\\}.$

In particular, each iteration has the following choices of freedom:

  * •

a finite index set $J_{n}$, denoting the number of search directions $u_{j}$
and thus the number of hyperplanes used in the $n$-th iteration,

  * •

search-directions $u_{j}$ and parameters $\alpha_{j}$ which are chosen in
accordance to the right-hand side $g$ and the forward model $B$, for instance,
in case $J=\\{n\\}$ a suitable choice is $u_{n}=B^{*}(Bf_{(n)}-g)$ as well as
$\alpha_{n}=\langle Bf_{(n)}-g,g\rangle$.

This method relies on knowing the exact forward model $B$ and on exactly known
right-hand side $g$. In practice however, only noisy data $g^{\delta}$ with
noise level $\|g-g^{\delta}\|\leq\delta$ can be measured and, as in our case,
only an inexact forward operator $B^{\eta}$ with inexactness level
$\|B-B^{\eta}\|\leq\eta$ might be available for the reconstruction. Thus, the
idea of RESESOP (regularized SESOP) is to replace the hyperplanes
$H(u,\alpha)$ by stripes

$H(u,\alpha,\xi):=\\{x\in X:|\langle u,x\rangle-\alpha|\leq\xi\\},$

whose width $\xi$ is chosen in dependence on the noise and inexactness levels
$\delta$ and $\eta$, see Figure 2 for an illustration. This idea was
originally suggested in [20] for the noisy data case and was expanded to model
imperfections in [15]. The introduction of Morozov’s discrepancy principle
offers a stopping criterion; as a result, it can be proven under certain
conditions on the choice of $J_{n},u_{i},\alpha_{i}$, that the method provides
a regularized solution to the original inverse problem [20].

$H(u,\alpha)$$H(u,\alpha+\xi)$$H(u,\alpha-\xi)$$\xi$$H(u,\alpha,\xi)$ Figure
2: Illustration of a stripe

When solving an inverse problem that shows a subproblem structure like (7), it
is beneficial to include individual inexactness and noise levels
$\eta_{i},\,\delta_{i}$, in our case

$\|v_{i}-v_{i}^{\delta}\|\leq\delta_{i},\quad\|A_{i}-A_{i}^{\eta}\|\leq\eta_{i},\quad
i\in\\{0,1,\dots,N_{\tau}-1\\}.$

In particular, with such estimates, we obtain

$\|A_{i}^{\eta}c-v^{\delta}_{i}\|\leq\eta_{i}\,\rho+\delta_{i}$

as upper bound for the residual on the set of solutions whose norm is bounded
by $\rho>0$. For this case, Blanke et al proposed the so-called RESESOP-
Kaczmarz method: Each iteration corresponds to a RESESOP iteration with
subproblem specific inexactness and noise levels applied to every subproblem
that does not yet satisfy the discrepancy principle. In this way, the
algorithm combines all subproblems and incorporates information from all time
intervalls $I_{i}$ simultaneously.

A general form of the RESESOP-Kaczmarz method is stated in algorithm 1. With
$n$ being the current iteration index, the notation $[n]:=n\text{ mod
}N_{\tau}$ identifies the corresponding subproblem (6). One full iteration
consists of $N_{\tau}$ subiterations, each including a metric projection onto
the intersection of stripes $H_{(n)}^{\eta,\delta}$.

Algorithm 1 RESESOP-Kaczmarz

Choose $c_{(0)}$ and constant $\rho>0$.Let $n$ be the iteration index.

while $c_{(n)}\neq c_{(n-N_{\tau})}$ for $[n]=0$ do

if $||A_{[n]}c_{(n)}-v_{[n]}^{\eta,\delta}||\leq
1.001\cdot(\eta_{[n]}\rho+\delta_{[n]})$ then

$c_{(n+1)}=c_{(n)}$

else

Choose a finite index set $J_{n}\subset\\{0,1,\dots,n\\}$

Choose search directions $u_{n,j}$ for all $j\in J_{n}$

Define $H_{(n)}^{\eta,\delta}:=\underset{j\in
J_{n}}{\cap}H(u_{n,j},\alpha_{n,j},\xi_{n,j})$

with parameters $\alpha_{n,j}=\langle u_{n,j},c_{(n)}\rangle$ and
$\xi_{n,j}=(\eta_{[j]}\rho+\delta_{[j]})\|A_{[j]}c_{(n)}-v_{[j]}^{\eta,\delta}\|$

Compute $c_{(n+1)}=P_{H_{(n)}^{\eta,\delta}}(c_{(n)})$

end if

end while

RESESOP-Kaczmarz does not require extensive or unrealistic a priori
information to solve inverse problems with inexact forward model such as
dynamic reconstruction problems. The estimates of the uncertainty levels are
the only additional information used compared to standard static approaches.
Consequently, when applying the method to dynamic problems, the motion
information will be implicitly included in the inexactness levels.

However, a sufficiently good estimate is crucial for the performance of
RESESOP. If the estimates are too large the convergence of the method can be
fast but lead to a computed solution with a high error rate. If the levels are
estimated too low the noise and the model imperfections might not be
compensated for sufficiently resulting in a noisy reconstruction with
artefacts. The extreme case, where the levels are set to zero, corresponds to
a SESOP-Kaczmarz reconstruction without regularization and accounting for
model inexactness, i.e. in particular to a static reconstruction where the
Kaczmarz loop results in an averaging of the data.

Next, we discuss further adaptions of this general framework to the
application in dynamic MPI, in particular how the required estimates can be
obtained and how the search directions can be chosen.

###  2.3 Application in MPI

####  2.3.1 Inexactness levels in dynamic MPI

We first discuss how to determine estimates for the noise and inexactness
levels. In a static multi-frame setting, the noise levels $\delta_{i}$ can be
determined by comparing the data vectors measured per frame.

The inexactness levels depend on the selected model for the reconstruction
algorithm. Due to the time consuming process of calibration in real MPI and
since we assume no further prior information on the unknown motion $\Gamma$,
we propose to use the static system matrix $A_{i}$ as inexact forward model.
Without knowledge of $\Gamma$, it is difficult to compute a viable upper bound
$\eta_{i}$ with $\|A_{i}-A_{i,\Gamma}\|\leq\eta_{i}$.

However, if we compare again data sets for different frames in a dynamic
scenario, their difference captures noise as well as motion. For this reason,
we propose to jointly estimate both uncertainty levels together in form of
levels $\zeta_{i}\in\mathbb{R}$ characterizing the total inexactness of the
system, more precisely for each subproblem

$\zeta_{i}\approx\eta_{i}\rho+\delta_{i},\quad
i\in\\{0,1,\dots,N_{\tau}-1\\}.$

These values are computed by comparing the data $v_{0}$ of the initial
subproblem with data $v_{i}$ of the $i$-th subproblem. If the size of each
subproblem corresponds to one frame, it is simply $\zeta_{i}=\|v_{0}-v_{i}\|$.

However, in case of subframe-sized subproblems, this direct comparison is not
possible but requires an additional interpolation step beforehand. Each
subproblem corresponds to a number of time points $t$ of a frame.
Consequently, the data fragment of the initial subproblem is comparable to the
data acquired at the same time points of each frame, e.g. the uncertainty
levels between the first subproblem of the first frame and the first
subproblem of every other frame can be evaluated. However, to compute an
estimate for the uncertainty level of another subproblem of the same frame
with the same method is impossible. Therefore, we apply a cubic interpolation
between comparable data fragments to acquire uncertainty levels for each
subproblem.

An alternative method to estimate the uncertainty levels $\zeta_{i}$ is to
make use of preliminary regularized Kaczmarz reconstructions. To quantify the
error between Kaczmarz’s reconstructions, the mean-squared error MSE is used.
However, this approach depends on a priori reconstructions which will
typically contain artefacts and thus can result in a low MSE. Thus, an
extraction directly from measured data seems favorable.

####  2.3.2 Choice of search directions

The convergence properties of SESOP, and hence RESESOP-Kaczmarz, depend on the
number and choice of the search-directions. On the one hand, a higher number
of search directions results in more complex and computationally expensive
iterations. On the other hand, the algorithm requires in total less iterations
to achieve the desired accuracy. In this work we will focus on two search
directions which is a compromise between computational expense for one
iteration as well as total number of iterations.

More precisely, we choose in the $n$-th iteration step the finite index set
$J_{n}=\\{n_{-},n\\}$ with $n_{-}$ representing the last iteration in which
the discrepancy principle did not yet hold. The search directions $u_{n,j}$
with $j\in J_{n}$ are then chosen as

$\displaystyle u_{n,n_{-}}$
$\displaystyle=A_{{[n_{-}]}}^{*}(A_{{[n_{-}]}}c_{(n_{-})}-v_{{[n_{-}]}}),$
$\displaystyle u_{n,n}$
$\displaystyle=A_{{[n]}}^{*}(A_{{[n]}}c_{(n)}-v_{{[n]}}).$

In particular, this choice meets the criteria formulated in [15] that
guarantee the convergence of RESESOP-Kaczmarz to a regularized solution of the
underlying dynamic inverse problem.

####  2.3.3 Positivity constraints

In the MPI setting the reconstructed image approximates the concentration $c$,
which is a vector consisting of real, positive numbers. Therefore, a
positivity constraint is a common technique to improve reconstruction
algorithms [2]. To that end, after each iteration every negative part is set
to zero.

##  3 Data

Our test cases for the numerical evaluation comprise both real and simulated
data, each offering valuable insights into the dynamic reconstruction problem
in MPI and the performance of the proposed method. This dual test set strategy
allows for a robust evaluation under controlled conditions (simulated data) as
well as real-world scenarios (real data).

###  3.1 Simulated Data

We first conduct a fully simulated numerical experiment. This simulated
dataset was generated using the model B3 presented by Kluth, Szwargulski, and
Knopp [21].

The goal was to construct a similar setup as the available real data. To that
end, the phantom is a cylinder rotating through the field of view. There are
two different rotation speeds, depending on the number of frames per rotation
(fpr). They are equivalent to seven and 44 fpr. Figure 3 illustrates the
underlying rotating motion of the object in the simulated test data case.

(a) Frame 1

(b) Frame 2

(c) Frame 3

(d) Frame 4

(e) Frame 5

(f) Frame 6

(g) Frame 7

(h) Frame 8

Figure 3: The course of the motion illustrated by the ground truth phantom at
the beginning of each trajectory of the fast rotation (seven fpr).

The simulated system matrix is derived from dynamic simulations of Néel-type
particle magnetization dynamics. Here, the uniaxial particle anisotropy varies
with position within the field of view. The anisotropy is minimal at the
center but intensifies toward the boundaries. This choice aligns with the
assumptions in MPI, where the static component of the applied magnetic field
is believed to induce physical rotation of particles and influence their
combined anisotropy energy landscape.

The simulated data were generated on a $3$D grid of $61\times 61\times 5$
pixels, each with a resolution of $0.5$ mm. A $2$D Lissajous excitation was
used, with anisotropy constants reaching values as high as
$1250,\text{J/m}^{3}$.

To create the system matrix, we utilized a dedicated toolbox presented in
Albers, Kluth, and Knopp [22]. To prevent inverse crime, the data were
calculated using a slightly shifted, larger matrix, ensuring the integrity of
our simulation experiments. Furthermore, we added white Gaussian noise to the
data corresponding to a signal-to-noise ratio of 10.

###  3.2 Real Data

The real data were provided by Brandt from the University of Hamburg and
collected by a team at UKE Hamburg led by Knopp. These data encode a glass
capillary, which rotates at an average speed of seven Hz and was initially
introduced in the study by Gdaniec et al. [8]. This experiment represents a
very fast motion.

Working with real data presents some additional challenges compared to
simulated data. Typically, one lacks a ground truth which hinders a direct
assessment of the accuracy of the reconstruction. Moreover, real measurements
usually have high degree of noise, adding complexity to the reconstruction
task in the dynamic case where noise levels cannot be reduced by averaging
over multiple frames. This is further amplified in MPI due to the fact that
the system function is measured as well.

These experiments were conducted using a pre-clinical MPI scanner from Bruker
Biospin in Ettlingen. The imaging setup involved 3D MPI measurements with a
Lissajous excitation driven by three sinusoidal drive fields in the $x$-,
$y$-, and $z$-directions. These fields operated at frequencies of
$f_{x}=\frac{2.5\text{MHz}}{102}$, $f_{y}=\frac{2.5\text{MHz}}{96}$, and
$f_{z}=\frac{2.5\text{MHz}}{99}$, with an amplitude of
$14\frac{\text{mT}}{\mu_{0}}$. Each cycle of data acquisition took $21.54$ ms,
and the induced signal was sampled at intervals of $0.8\mu\text{s}$. The
measurements covered $15.625$ positions in a $25\times 25\times 25$ grid. More
detailed information can be found in [8].

###  3.3 Pre-processing

To the real data, the standard pre-processing steps were applied. After
applying the Fourier transform, a frequency selection was performed, more
precisely all frequencies below 80 k Hz and above 625 k Hz were eliminated.
This was followed by a SNR threshholding. Furthermore, a singular value
decompostion as proposed in [23] was performed to reduce the computational
cost of reconstructing a three dimensional image. With weigths computed from
background measurements, the randomized singular value decomposition trimmed
the system matrix to 5000 rows.

However, it is important to note that these conventional pre-processing steps
are tailored to full-frame data. Thus, they may encounter limitations when
applied to dynamic cases, where the dynamics requires a division of the data
into sub-frames. In order to realize this, the raw data first have to be
transformed into the time domain with an inverse Fourier transform. Only then
the frame-sized problems can be split into smaller subproblems. To this end,
access to the original raw-data is crucial for the solution of reconstruction
problems with fast dynamics. If only pre-processed data are provided, one is
restricted to decompositions, where the time intervalls $I_{i}$ correlate to
individual frames.

When working with simulated data, most of the pre-processing steps are not
necessary. Due to the 2D-Lissajous excitation, the system matrix is naturally
smaller compared to the system matrix of real data. Therefore, pre-processing
steps to reduce computational work are not critical and some steps e.g. SNR
threshholding not applicable to simulated data.

##  4 Results and Discussions

This section is dedicated to the numerical evaluation of the RESESOP-Kaczmarz
method on both simulated and real dynamic MPI data. In particular, we
investigate the influence of critical parameters of our algorithm on the
result.

All experiments are conducted with 3D data, but for clarity, each result is
shown with a 2D slice. In each experiment, the same frame and slice are
depicted, ensuring a comparable view across all scenarios.

###  4.1 Simulated Data

We first present reconstruction results from simulated noisy data, cf. Section
3.1, with two different rotation speeds - a very fast scenario of seven fpr as
well as the scenario of 44 fpr. Regarding the coupling of time scales, we
consider here each frame as one subproblem, i.e. we are choosing the least
intrusive version in the MPI pipeline.

To illustrate the motion compensation properties of the RESESOP-Kaczmarz
method, we compare the results with those generated by the regularized
Kaczmarz-algorithm which is the most commonly used algorithm within MPI. Since
averaging data is not recommendable for dynamic data, all reconstructions with
regularized Kaczmarz use one frame of data. As recommended in the literature,
e.g. [24], the algorithm is stopped early after five iterations since the
quality of the reconstruction deteriorates with more iterations. Furthermore,
at each iteration positivity of the approximated solution is enforced. The
initial value was set to zero. The regularization parameter $\lambda$ is
optimized by performing reconstructions for a wide range of parameters and
choosing the one that leads to the best result.

(a) Ground truth

(b) Ground truth

(c) Regularized Kaczmarz

(d) Regularized Kaczmarz

(e) RESESOP-Kaczmarz

(f) RESESOP-Kaczmarz

Figure 4: Dynamic Reconstructions of frame four from noisy simulated data with
two different rotation speeds. The left column shows results for a rotation
speed of 44 fpr, the right column for seven fpr.

The reconstruction results are illustrated in Figure 4. The first row depicts
the ground truth at frame four for both rotational velocities. The second row
shows the images reconstructed by regularized Kaczmarz. In both dynamic
scenarios significant artifacts are visible. Since the algorithm cannot rely
on averaging over multiple frames, the noise in the data is not satisfactorily
reduced, resulting in artifacts throughout the entire field-of-view.
Furthermore, the circular shape of the object is blurred, distorted and hence
no more recognizable. As to be expected, this effect is worse the faster the
motion. These examples motivate again that the dynamic behaviour of the
concentration needs to be taken into account within the reconstruction step.

The respective results of RESESOP-Kaczmarz are presented in the third row of
figure 4. In both dynamic scenarios, the actual shape of the object is
recognizable in the reconstructions, i.e. the algorithm does indeed compensate
for the motion. However, due to the strict restriction to frame-sized
subproblems, we can still observe some small distortions in case of the very
fast motion. Also the location of the object is not correctly reconstructed
but depicted in the middle of the traveled path during frame four and not in
the beginning. However, we further note that in both dynamic cases, the method
is able to eliminate the artifacts caused by the noise throughout the field of
view.

Altogether, Figure 4 demonstrates the advantage of applying RESESOP-Kaczmarz
to dynamic MPI problems.

###  4.2 Parameters for RESESOP-Kaczmarz

In this section, we want to study the robustness of the RESESOP-Kaczmarz
algorithm regarding estimates of the total inexactness levels $\zeta_{i}$,
number of iterations and size of subproblems. All experiments in this
subsection have been conducted with simulated data for the fast rotating
object with ten percent of added noise. The reconstruction in Figure 4 (f) was
obtained by determining $\zeta_{i}$ as described in Section II.III.1 with the
Euclidean norm and by performing ten full RESESOP-Kaczmarz iterations while
the size of all subproblems corresponds to one full frame.

To test stability regarding the estimates of the inexactness levels, we
performed the reconstruction with scaled versions of these inexactness levels
(more precisely with factors $10^{-1},\,10^{1}$ and $10^{5}$ respectively)
representing an under-, respectively over-estimations. In addition, we also
considered 15% noise added to the norm estimates $\zeta_{i}$ to test
robustness with respect to computation errors. Lastly, we further use the
alternative approach based on prior reconstruction discussed in section 2.3.1.

| Norm | 15% Noise | RegKacz
---|---|---|---
$\zeta\cdot 10^{-1}$ |  |  |
$\zeta\cdot 10^{0}$ |  |  |
$\zeta\cdot 10^{1}$ |  |  |
$\zeta\cdot 10^{5}$ |  |  |
Table 1: Dynamic Reconstructions with RESESOP-Kaczmarz for different estimates
of the total inexactness levels $\zeta_{i}$. Column 1: computed as distance of
two frames via euclidean norm, Column 2 : as in Column 1 with 15 % added noise
, Column 3: computed as MSE from regularized Kaczmarz reconstructions. Each
row corresponds to one additional scaling factor.

Table 1 illustrates the influence of different inexactness levels on the
reconstruction result of RESESOP-Kaczmarz. Both the proposed norm estimation
directly from the measured data as well as the estimation from prior static
reconstructions with regularized Kaczmarz result in images of similar quality
(unless in the underestimated case, see first row, where the computation from
the data is beneficial). Nevertheless, an underestimation of the inexactness
levels leads to an increase in motion artefacts, while the motion is well
compensated for even if largely overestimated levels are used. The algorithm
is also stable regarding (random) errors in the computed inexactness levels,
see Column 2. Despite strong deviations of overall 15%, the algorithm still
provides a reconstruction of comparable image quality as with the exact
values.

Next, we study how the reconstruction quality develops in dependence on the
number of iterations. As was shown in [15], the algorithm converges monotone
as long as adequate search directions are chosen. This is clearly visible in
figure 5 which shows the evolution of the mean-squared error of the RESESOP-
Kaczmarz solution compared to the ground truth with increasing number of
iterations. Furthermore, we observe that the error reduces rapidly during the
first iteration and then flattens out.

Figure 5: Mean squared error in dependence on the amount of RESESOP-Kaczmarz
loops.

Lastly, we want to examine the influence of the size of the subproblems.
Depending on the speed of the motion, assuming the concentration to be static
during a complete Lissajous trajectory can be too much simplified. Indeed,
Figure 4 shows that under this assumption in the very fast scenario, the
reconstructed object is not correctly located. This can be improved by
considering subproblems of smaller size, see Figure 6.

(a) Frame

(b) $\frac{1}{2}$-frame

(c) $\frac{1}{4}$-frame

(d) $\frac{1}{8}$-frame

(e) $\frac{1}{16}$-frame

(f) $\frac{1}{32}$-frame

Figure 6: Dynamic Reconstructions from noisy simulated data in the fast
scenario where different sizes of subproblems are considered within the
RESESOP-Kaczmarz algorithm.

Choosing the sizes of the subproblems as a quarter of a frame or smaller
improves the reconstructed location of the object, see e.g. 6(c). However,
choosing the size too small can reintroduce a low amount of noise in the
reconstruction, see Figures 6(e) and 6(f). Thus, in case of very fast motions,
one has to balance motion compensation and noise reduction when choosing
appropriate sizes of subproblems.

###  4.3 Real Data

In this section we evaluate the performance of RESESOP-Kaczmarz on the real
data introduced in Section 3.2.

| Frame |  $\frac{1}{4}$-frame
---|---|---
Frame 1 |  |
Frame 2 |  |
Frame 3 |  |
Frame 4 |  |
Table 2: Dynamic Reconstructions with RESESOP-Kaczmarz of four consecutive
frames of real data for a 7 Hz rotation with different sized subproblems. The
original numeration of the frames starts at frame 750.

Table 2 depicts reconstructions of four consecutive frames of real data. The
left column shows images computed with RESESOP-Kaczmarz using frame-sized
subproblems. From the reconstructed images, we can clearly deduce the temporal
evolution of the object, namely the rotation of the glass capillary. Also its
spatial location per time step is overall well captured and the noise well
suppressed. However, examining frames one, two and four the circular shape of
the reconstructed glass capillary is slightly distorted due to the fast
motion.

The right column of Table 2 shows the respective reconstruction results using
$\frac{1}{4}$-frames as subproblems. Compared to the frame-wise case, they
portray the object with a more constant round shape but are also more affected
by noise. This illustrates again the trade-off between noise reduction and
motion compensation. In a future step, it would be interesting to study
whether algorithms (or reconstructions) for different subproblem sizes could
be merged in order to combine the good motion compensation property of small
subproblems with the better noise reduction of large subproblems.

##  5 Conclusion

In this article, we propose the RESESOP-Kaczmarz algorithm to reconstruct
dynamic objects from MPI data. The method takes motion into account as a model
inexactness and, within the MPI framework, the only required a priori
information can be computed directly from the measured data. The potential of
the method was demonstrated on real and simulated data. In particular, our
detailed experiments on simulated data show the robustness of the method
regarding its various parameters and that it outperforms the regularized
Kaczmarz algorithm, a common solver in MPI based on a stationary assumption on
the object. Since the method further allows to consider subproblems smaller
than one complete frame, it is suitable for a variety of dynamic MPI problems,
including scenarios with rapid particle movements.

## Acknowledgments

The authors acknowledge the support by ’Deutsche Forschungsgemeinschaft’ under
grant HA 8176/2-1. The authors further thank Hannes Albers and Tobias Kluth
for providing us with simulated data as well as Christina Brandt and Tobias
Knopp for real data.

## References

  * [1] B. Gleich and J. Weizenecker  “Tomographic imaging using the nonlinear response of magnetic particles”  In _Nature_ 435.7046, 2005, pp. 1214–1217  DOI: 10.1038/nature03808
  * [2] T. Knopp and T.. Buzug  “Magnetic Particle Imaging: An Introduction to Imaging Principles and Scanner Instrumentation”  Berlin/Heidelberg: Springer, 2012  DOI: 10.1007/978-3-642-04199-0
  * [3] H. Albers et al.  “Modeling the magnetization dynamics for large ensembles of immobilized magnetic nanoparticles in multi-dimensional magnetic particle imaging”  In _Journal of Magnetism and Magnetic Materials_ , 2022
  * [4] J. Weizenecker et al.  “Three-dimensional real-time in vivo magnetic particle imaging”  In _Physics in medicine and biology_ 54, 2009, pp. L1–L10  DOI: 10.1088/0031-9155/54/5/L01
  * [5] J. Haegele et al.  “Visualization of Instruments for Cardiovascular Intervention Using MPI”  In _Magnetic Particle Imaging_ Berlin, Heidelberg: Springer Berlin Heidelberg, 2012, pp. 211–215  DOI: 10.1148/radiol.12120424
  * [6] M. Graeser et al.  “Human-sized magnetic particle imaging for brain applications”  In _Nature Communications_ 10.1, 2019, pp. 1936  DOI: 10.1038/s41467-019-09704-x
  * [7] P. Ludewig et al.  “Magnetic particle imaging for assessment of cerebral perfusion and ischemia”  In _WIREs Nanomedicine and Nanobiotechnology_ 14.1, 2022, pp. e1757  DOI: https://doi.org/10.1002/wnan.1757
  * [8] N. Gdaniec et al.  “Detection and Compensation of Periodic Motion in Magnetic Particle Imaging”  In _IEEE Transactions on Medical Imaging_ 36.7, 2017, pp. 1511–1521  DOI: 10.1109/TMI.2017.2666740
  * [9] N. Gdaniec et al.  “Suppression of Motion Artifacts Caused by Temporally Recurring Tracer Distributions in Multi-Patch Magnetic Particle Imaging”  In _IEEE Transactions on Medical Imaging_ 39.11, 2020, pp. 3548–3558  DOI: 10.1109/TMI.2020.2998910
  * [10] C. Brandt and C. Schmidt  “Modeling Magnetic Particle Imaging for Dynamic Tracer Distributions”  In _Sens Imaging_ 22.45, 2021  DOI: https://doi.org/10.1007/s11220-021-00368-w
  * [11] C. Brandt and C. Schmidt  “Motion compensation for non-periodic dynamic tracer distributions in multi-patch magnetic particle imaging”  In _Physics in Medicine and Biology_ 67.8, 2022, pp. 085005  DOI: 10.1088/1361-6560/ac5ce6
  * [12] M. Burger, H. Dirks and L. Frerking  “7. On optical flow models for variational motion estimation”  In _Variational Methods_ Berlin, Boston: De Gruyter, 2017, pp. 225–251  DOI: 10.1515/9783110430394-007
  * [13] T. Kluth, B.N. Hahn and C. Brandt  “Spatio-temporal concentration reconstruction using motion priors in magnetic particle imaging”  In _IWMPI_ , 2019
  * [14] C. Brandt, T. Kluth, T. Knopp and L. Westen  “Dynamic image reconstruction with motion priors in application to 3D magnetic particle imaging”  In _arXiv_ , 2023  DOI: https://doi.org/10.48550/arXiv.2306.11625
  * [15] S.E. Blanke, B.N. Hahn and A. Wald  “Inverse Problems with inexact forward operator: Iterative regularization and application in dynamic imaging”  In _Inverse Problems_ 36, 2020  DOI: 10.1088/1361-6420/abb5e1
  * [16] M. Kaul et al.  “Magnetic particle imaging for in vivo blood flow velocity measurements in mice”  In _Physics in Medicine and Biology_ 63.6, 2018
  * [17] B.. Hahn  “Efficient algorithms for linear dynamic inverse problems with known motion”  In _Inverse Problems_ 30.3 IOP Publishing, 2014, pp. 035008  DOI: 10.1088/0266-5611/30/3/035008
  * [18] L. Desbat, S. Roux and P. Grangeat  “Compensation of Some Time Dependent Deformations in Tomography”  In _IEEE Transactions on Medical Imaging_ 26.2, 2007, pp. 261–269  DOI: 10.1109/TMI.2006.889743
  * [19] F. Schöpfer, T. Schuster and A. Louis  “Metric and Bregman projections onto affine subspaces and their computation via sequential subspace optimization methods”  In _J. Inv. Ill-Posed Problems_ 15, 2007, pp. 1–29  DOI: https://doi.org/10.1515/JIIP.2008.026
  * [20] F. Schöpfer and T. Schuster  “Fast regularizing sequential subspace optimization in Banach spaces”  In _Inverse Problems_ 25, 2009, pp. 015013  DOI: 10.1088/0266-5611/25/1/015013
  * [21] T. Kluth, P. Szwargulski and T. Knopp  “Towards accurate modeling of the multidimensional magnetic particle imaging physics”  In _New Journal of Physics_ 21.10 IOP Publishing, 2019, pp. 103032  DOI: 10.1088/1367-2630/ab4938
  * [22] H. Albers, T. Kluth and T. Knopp  “Simulating magnetization dynamics of large ensembles of single domain nanoparticles: Numerical study of Brown/Néel dynamics and parameter identification problems in magnetic particle imaging”  In _Journal of Magnetism and Magnetic Materials_ 541, 2022, pp. 168508  DOI: https://doi.org/10.1016/j.jmmm.2021.168508
  * [23] T. Kluth and B. Jin  “Enhanced reconstruction in magnetic particle imaging by whitening and randomized SVD approximation”  In _Physics in Medicine and Biology_ 64.12 IOP Publishing, 2019, pp. 125026  DOI: 10.1088/1361-6560/ab1a4f
  * [24] T. Knopp et al.  “Weighted iterative reconstruction for magnetic particle imaging”  In _Physics in Medicine and Biology_ 55.6, 2010, pp. 1577  DOI: 10.1088/0031-9155/55/6/003


# The hierarchical generalized linear model and the bootstrap estimator of the
error of prediction of loss reserves in a non-life insurance company

Alicja Wolny-Dominiak

###### Abstract

This paper presents the hierarchical generalized linear model (HGLM) for loss
reserving in a non-life insurance company. Because in this case the error of
prediction is expressed by a complex analytical formula, the error bootstrap
estimator is proposed instead. Moreover, the bootstrap procedure is used to
obtain full information about the error by applying quantiles of the absolute
prediction error. The full R code is available on the Github
https://github.com/woali/BootErrorLossReserveHGLM.

Keywords. Loss Reserves; Hierarchical Generalized Linear Model; Error of
prediction; MSE; Parametric Bootstrap

JEL: G22, C21, C53

##  1 Introduction

This analysis concerns a non-life insurance company. The largest item on an
insurer’s balance sheet are technical provisions. Any variations in their
values have a great impact on the insurer’s financial strength. A large part
of the provisions is the reserve for incurred but not reported (IBNR) claims
or simply – the loss reserve, which is crucial to the insurer’s solvency. The
total loss reserve is generally determined by statistical methods based on
both deterministic techniques and stochastic models and is a sum of
outstanding claims liabilities.

A wide variety of the total loss reserve prediction tools are presented in
literature in the group of stochastic methods. For years the most popular
method with insurers was the chain-leader technique (e.g. (Mack, 1994), (Mack
and Venter, 2000), (Verrall, 2000)). Parallel methods based on the regression
model have been developed (e.g. (Kremer, 1982), (Mack, 1991)). In (Renshaw and
Verrall, 1998) it is shown that the chain-leader method is equivalent to the
well-known generalized linear model (GLM) with fixed effects assuming over-
dispersed Poisson (ODP) distribution for incremental claims. This paper puts
forward a mixed model for loss reserving as an extension of the GLM in which
data for different homogenous groups are assumed to be dependent in some way.
The model belongs to the class of HGLMs and contains both fixed effects as
well as random effects. Model’s parameters are estimated in the frequentist
approach using extended likelihood (Lee and Nelder, 1996), (Lee and Nelder,
2001), (Lee and Nelder, 2003). As the total loss reserve is actually a random
variable, its value has to be predicted. In order to measure the prediction
accuracy, the bootstrap root mean squared error of prediction (RMSEP)
estimator is proposed.

##  2 Data – single loss triangle and total loss reserve

In loss reserving, data have a specific form referred to as the loss triangle.
Consider the random variable $Y_{ij}$ with $y_{ij}$ realizations
$i,j=0,\ldots,n$ being the incremental quantity for insurance claims that
occurred in year $i$ (origin year) and were reported to the insurer after $j$
years (development year). It could be the value of claims as well as their
number. The matrix $[Y_{ij}]_{(n+1)\times(n+1)}$ represents random variables
in the loss triangle for a single LOB as shown in Table 1 below.

Table 1: Incremental loss triangle for $k$th group  i/j | 0 | … | $n-1$ | $n$ |
---|---|---|---|---|---
0 | $Y_{00}$ | … | $Y_{0n-1}$ | $Y_{0n}$ |
1 | $Y_{10}$ | … | $Y_{1n-1}$ |  |
… | … | … |  |  |
n | $Y_{n0}$ |  |  |  |

The realizations $y_{ij}$ for $i+j\leq n$ are observed data while $y_{ij}$ for
$i+j>n$ represent the future unobserved data. The goal in loss reserving is to
predict the random variable $R$ being the total loss reserve defined as the
sum of unknown future claims:

$R=\sum_{i,j;i+j>n}Y_{ij}$ (1)

In order to asses the uncertainty of predictor $\hat{R}$, usually

$RMSEP(\hat{R})=\sqrt{E[(\hat{R}-R)^{2}]}$ (2)

is adopted.

##  3 The model and the prediction

Let us assume $Y_{ij}\sim ODP(\mu_{ij},\phi)$ as independent and identically
distributed (i.i.d) with $E[Y_{ij}]=\mu_{ij}$ and $Var(Y_{ij})=\phi\mu_{ij}$
and the unknown dispersion parameter $\phi$. In the GLM the mean has a linear
form $\log(\mu_{ij})=\bm{x}_{ijk}^{\prime}\bm{\beta}$, where
$\bm{\beta}=(c,u_{1},\ldots,u_{n},\beta_{1},\ldots,\beta_{n})^{\prime}$ is the
vector of fixed effects, $\bm{x}_{ijk}^{\prime}$ is the row of the design
matrix $\bm{X}$ and $u_{0}=\beta_{0}=1$. The dispersion parameter $\phi$ is
the same for all elements of the loss triangle. It is easy to extend this
model to the general class of models assuming that incremental claims follow
the three-parametric Tweedie $T(\mu,\phi,p)$ (Jorgensen, 1987) with dispersion
parameter $\phi$ and power $p\in(0,\infty]\setminus(0,1)$. The GLM is then of
the form:

$E[Y_{ij}]=\mu_{ij}=\exp(\bm{x}_{ij}^{\prime}\bm{\beta}),Var(Y_{ij})=\phi\mu_{ij}^{p},$
(3)

In order to draw an inference about the total loss reserve $R$, MLE estimators
of parameters
$(\hat{c},\hat{u}_{1},\ldots,\hat{u}_{n},\hat{\beta}_{1},\ldots,\hat{\beta}_{n},\hat{\phi})^{\prime}$
can be plugged into Formula (1) giving the predictor:

$\hat{R}_{GLM}=\sum_{i,j;i+j>n}\hat{Y}_{ij}=\sum_{i,j;i+j>n}\exp(\hat{c}+\hat{u}_{i}+\hat{\beta}_{j}).$
(4)

The total loss reserve prediction accuracy is typically measured using the
$RMSEP$. According to (England and Verrall, 1999), the predictor defined in
Formula (4) can be expressed as:

$\displaystyle
RMSEP(\hat{R}_{GLM})=\sum_{i,j;i+j>n}\phi_{ij}\mu_{ij}^{p}+\sum_{i,j;i+j>n}\mu_{ij}^{2}Var(\bm{x}_{ij}^{\prime}\bm{\beta})+$
$\displaystyle+\sum_{i,j;i+j>n}\mu_{j_{1}j_{1}}\mu_{j_{2}j_{2}}Cov(\bm{x}_{i_{1}j_{1}}^{\prime}\bm{\beta},\bm{x}_{i_{2}j_{2}}^{\prime}\bm{\beta}).$
(5)

The fundamental assumption in the GLM is the independence between random
incremental claims $Y_{ij}$, $i,j=0,\ldots,n$ taken from a single loss
triangle. However, this assumption may be inadequate, as indicated by the fact
that origin years reflect the same process of loss development. This means
that there is a dependence between claims from the same development year $j$
but from different origin years. The dependence may be taken into account by
applying a mixed model with fixed and random effects in place of the GLM with
fixed effects only.

Let us assume that effects $u_{0},\ldots,u_{n}$ occuring in model (3) are
independent realizations of random variable $U$ following Tweedie
distributions $u_{i}\sim T(\psi_{iu},\phi_{u})$ with $E[u_{i}]=\psi_{iu}$ and
$Var(u_{i})=\phi_{u}\psi_{iu}^{p}$. This means that incremental claims from
the loss triangle are now conditionally independent, which can be described in
the following way: $Y_{ij}$ of different origin years are independent, but due
to the random effects, for any origin year $i$, the random variables $Y_{ij}$
for the different development years $j$ are dependent (cf. (Gigante et al.,
2013b), pp. 383, (a1)–(4)). This approach is examined in (Gigante et al.,
2013a), (Gigante et al., 2013b), where the hierarchical GLM (HGLM) is
described. The authors assume a conditional ODP distribution for incremental
claims $Y_{ij}|u_{i}$. Like for the GLM, this assumption is extended to any
$p$ such that $Y_{ij}|u_{i}\sim T(\mu_{ij},\phi,p)$, which gives the following
general form of the HGLM:

$\displaystyle
E[Y_{ij}|u_{i}]=\mu_{ij}(u_{i})=\exp(\bm{x}_{ij}^{\prime}\bm{\beta}+\bm{z}_{ij}^{\prime}\bm{v})$
$\displaystyle Var(Y_{ij}|u_{i})=\phi\mu_{ij}^{p}(u_{i}),$ (6)

where $\bm{\beta}=(c,\beta_{1},\ldots,\beta_{n})^{\prime}$,
$v=(\log(u_{1}),\ldots,\log(u_{n}))^{\prime}$ and $\bm{z}_{ijk}$ is the row of
the design matrix $\bm{Z}$. Under the model defined in (3) and using the
conditional independence assumption, the total loss reserve is defined as the
conditional expected value:

$R_{HGLM}=\sum_{i,j;i+j>n}E[Y_{ij}|u_{i}]=\sum_{i,j;i+j>n}\exp(\bm{x}_{ij}^{\prime}\bm{\beta}+\bm{z}_{ij}^{\prime}\bm{v})$
(7)

The estimation of the HGLM parameters is more complex than in the case of the
GLM due to the presence of random effects. Following (Gigante et al., 2013a),
the likelihood-based approach is applied using the h-likelihood function. The
idea is to treat the vector of random effects $\bm{v}$ as the vector of fixed
effects and transform the HGLM into an augmented GLM. Details concerning the
h-likelihood estimation are presented in two fundamental papers (Lee and
Nelder, 1996), (Lee and Nelder, 2001). This numerical analysis uses the
$\bm{R}$ implementation of the algorithm taken from the hglm package
(Ronnegard et al., 2010). After estimators $\hat{\bm{\beta}}$ and predictors
$\hat{\bm{v}}$ are found, the total loss reserve predictor is of the form:

$\hat{R}_{HGLM}=\sum_{i,j;i+j>n}\exp(\hat{c}+\hat{\beta}_{j})\hat{u}_{i}$ (8)

A similar mixed model-based approach is described in (Antonio et al., 2006),
but only for a lognormal distribution. The authors adopted the generalized
linear mixed model (GLMM) and the REML as the method of estimation of the
model parameters.

##  4 Bootstrap MSEP

In (Gigante et al., 2013b) the estimator of $RMSEP(\hat{R})$ is derived
analytically in the case of a single loss triangle. Unfortunately, the formula
describing the estimator has a rather complex structure. An alternative may be
to use the parametric bootstrap technique, like for the $RMSEP(\hat{R}_{GLM})$
proposed in (England and Verrall, 1999). This paper proposes another
procedure.

Let us denote the realization of $Y_{ij}$ as $y_{ij}$. As the HGLM is commonly
an extension of the GLM, Pearson residuals
$r_{ij}=\frac{\hat{y}_{ij}-\mu_{ij}}{\sqrt{\mu_{ij}}}$ are bootstrapped
instead of the data from the loss triangles. A single $b$th simulation,
$b=1,\ldots,B$ in the full bootstrap procedure has the following steps:

  * •

take a sample $r_{ij}^{*b}$ with replacement from Pearson residuals $r_{p}$
based on data from the loss triangle $i+j\leq n$

  * •

calculate new data
$y_{ij}^{*b}=r_{ij}^{*b}\sqrt{\hat{\mu}_{ij}}+\hat{\mu}_{ij}$, $i+j\leq n$

  * •

estimate the vector of parameters
$(\bm{\beta}^{*b},\bm{v}^{*b},\phi^{*b},\phi_{u}^{*b})^{\prime}$ based on new
data $y_{ijk}^{*b}$

  * •

calculate predicted values
$\hat{y}_{ij}^{*b}=\exp(\hat{c}+\hat{\beta}_{j})\hat{u}_{i}$ for $i+j>n$

  * •

generate values $y_{ij}^{*b}$ for $i+j>n$ from Tweedie distribution with the
dispersion parameter $\phi^{*b}$

After $B$ iterations of the bootstrap procedure the estimator of
$RMSEP(\hat{R}_{HGLM})$ is of the following form:

$\widehat{RMSE}(\hat{R})=\sqrt{\frac{1}{B}\sum_{b=1}^{B}(\sum_{i,j;i+j>n}\hat{y}_{ij}^{*b}-\sum_{i,j;i+j>n}{y}_{ij}^{*b})^{2}}$
(9)

Let us notice that the commonly used $RMSEP$ gives information only about the
average error of prediction. If full distribution of the prediction error is
of interest, quantiles of the absolute prediction error can be used (see
(Ża̧dło, 2013)), where the $p$th quantile is defined as
$Q_{p}(\hat{R})=Q_{p}(|\hat{R}-R|)$. Quantiles reflect the relation between
the magnitude of the error and the probability of its realization. In order to
estimate $Q_{p}(\hat{R})$, the bootstrap procedure described above can be
applied. In the last step, Formula (9) is replaced by the formula presented
below:

$\widehat{Q_{p}}(\hat{R})=Q_{p}(|\sum_{i,j;i+j>n}\hat{y}_{ijk}^{*b}-\sum_{i,j;i+j>n}{y}_{ijk}^{*b}|)$
(10)

As a result, the bootstrap estimator of the $Q_{p}(\hat{R})$ error is
obtained.

##  5 Numerical example

The subject of this analysis is the single loss triangle taken from (Wüthrich
and Merz, 2008), pp. 33, like in (Gigante et al., 2013a). In order to
demonstrate the HGLM for loss reserving, it is assumed that incremental claims
conditionally follow Tweedie distribution $Y_{ij}|u_{i}\sim
T(\mu_{ij},\phi_{u})$ with power $p=1$ and random effects follow Tweedie
distribution $u_{i}\sim T(\psi_{u},\phi_{u})$ with power $p=2$. This model is
equivalent to the ODP-Gamma HGLM. Moreover, dispersion parameters $\phi$,
$\phi_{u}$ are assumed as constant for both incremental claims and random
effects, respectively. The vector of parameters

$(c,\beta_{1},\ldots,\beta_{9},u_{0},\ldots,u_{9},\phi,\phi_{u})^{\prime}$

is estimated using the h-likelihood function. Plugging the obtained values
into Formula (8), the total loss reserve and the loss reserves for origin
years are obtained. In order to find the estimator of $RMSEP$, the bootstrap
procedure was implemented. The number of simulations was $n=1000$. The results
are presented in Table 2.

Table 2: Total loss reserve and bootsrap estimator $\widehat{RMSE}$ Origin year $i$ | $\hat{R}_{HGLM}$ | $\widehat{RMSEP}_{HGLM}$ | $\hat{R}_{GLM}$ | $\widehat{RMSEP}_{GLM}$
---|---|---|---|---
1 | 15 239 | 19 763 | 15 125 | 19 620
2 | 26 415 | 24 220 | 26 257 | 24 508
3 | 35 264 | 27 182 | 34 538 | 27 613
4 | 86 572 | 38 331 | 85 302 | 40 690
5 | 157 840 | 49 339 | 156 494 | 54 876
6 | 290 991 | 66 545 | 286 121 | 70 027
7 | 465 024 | 81 939 | 449 161 | 86 568
8 | 1 078 769 | 130 374 | 1 043 237 | 132 874
9 | 3 988 971 | 281 813 | 3 950 809 | 315 437
Total | 6 145 085 | 373 914 | 6 047 044 | 403 506
|  $\phi$=14 739 |  |  $\phi$=14 714 |
|  $\phi_{u}$=0.0054 |  |  |

It can be seen that the HGLM gives higher loss reserve values compared to the
GLM but the errors are generally lower. Naturally, errors rise for subsequent
origin years in both cases, which results directly from the fact that there
are fewer and fewer observations in the loss triangle. The two models are not
compared directly because in the case of the HGLM it is the conditional value
of reserves that is determined.

Figure 1: Measures of accuracy for origin years and the total loss reserve

The bootstrap procedure makes it possible to obtain not only the RMSEP but
also the quantiles of the absolute error of prediction $Q_{p}$. Figure 1
presents quantiles $Q_{0.5}$, $Q_{0.75}$, $Q_{0.9}$ and $Q_{0.95}$ for origin
years and for the total loss reserve. It can be seen that in earlier origin
years the RMSEP value is more or less equal to quantiles of the order of 0.75,
whereas for the total loss reserve the value is close to quantiles of the
order of 0.9.

##  6 Conclusions

GLMs are popular statistical techniques in actuarial practice, especially in
ratemaking but also in loss reserving. However, the independence assumption
needed in GLMs is generally violated in many cases. There are three basic
advantages of the mixed HGLM application. Firstly, by introducing random
effects into the model, the dependencies between development years are taken
into account. Secondly, by imposing certain constraints on random effects, it
is possible to take account of external information which does not come from
the sample directly, like in (Gigante et al., 2013a), and which has an impact
on the total loss reserve value. Thirdly, distributions can be set flexibly
within the Tweedie family, e.g. the Gamma, the inverse-Gaussian or the
compound Poisson distribution. The downside, however, is the complex form of
the error prediction. Therefore, it is proposed herein that the error should
be determined by means of the bootstrap technique. Although this solution is
not perfect, its important advantage is that full information on the absolute
error distribution can be obtained easily using quantiles.

## References

  * Antonio et al. [2006] K. Antonio, J. Beirlant, T. Hoedemakers, and R. Verlaak.  Lognormal mixed models for reported claims reserves.  _North American Actuarial Journal_ , 10(1):30–48, 2006.
  * England and Verrall [1999] P. England and R. Verrall.  Analytic and bootstrap estimates of prediction errors in claims reserving.  _Insurance: mathematics and economics_ , 25(3):281–293, 1999.
  * Gigante et al. [2013a] P. Gigante, L. Picech, and L. Sigalotti.  Claims reserving in the hierarchical generalized linear model framework.  _Insurance: Mathematics and Economics_ , 52(2):381–390, 2013a.
  * Gigante et al. [2013b] P. Gigante, L. Picech, and L. Sigalotti.  Prediction error for credible claims reserves: an h-likelihood approach.  _European Actuarial Journal_ , 3(2):453–470, 2013b.
  * Jorgensen [1987] B. Jorgensen.  Exponential dispersion models.  _Journal of the Royal Statistical Society. Series B (Methodological)_ , pages 127–162, 1987.
  * Kremer [1982] E. Kremer.  Ibnr-claims and the two-way model of anova.  _Scandinavian Actuarial Journal_ , 1982(1):47–55, 1982.
  * Lee and Nelder [2003] Y. Lee and J. Nelder.  Extended-reml estimators.  _Journal of Applied Statistics_ , 30(8):845–856, 2003.
  * Lee and Nelder [1996] Y. Lee and J. A. Nelder.  Hierarchical generalized linear models.  _Journal of the Royal Statistical Society. Series B (Methodological)_ , pages 619–678, 1996.
  * Lee and Nelder [2001] Y. Lee and J. A. Nelder.  Hierarchical generalised linear models: a synthesis of generalised linear models, random-effect models and structured dispersions.  _Biometrika_ , 88(4):987–1006, 2001.
  * Mack [1991] T. Mack.  A simple parametric model for rating automobile insurance or estimating ibnr claims reserves.  _Astin Bulletin_ , 21(1):93–109, 1991.
  * Mack [1994] T. Mack.  Which stochastic model is underlying the chain ladder method?  _Insurance: mathematics and economics_ , 15(2):133–138, 1994.
  * Mack and Venter [2000] T. Mack and G. Venter.  A comparison of stochastic models that reproduce chain ladder reserve estimates.  _Insurance: Mathematics and Economics_ , 26(1):101–107, 2000.
  * Renshaw and Verrall [1998] A. E. Renshaw and R. J. Verrall.  A stochastic model underlying the chain-ladder technique.  _British Actuarial Journal_ , 4(04):903–923, 1998\.
  * Ronnegard et al. [2010] L. Ronnegard, X. Shen, and M. Alam.  hglm: A package for fitting hierarchical generalized linear models.  _The R Journal_ , 2(2):20–28, 2010.
  * Verrall [2000] R. J. Verrall.  An investigation into stochastic claims reserving models and the chain-ladder technique.  _Insurance: Mathematics and Economics_ , 26(1):91–99, 2000.
  * Wüthrich and Merz [2008] M. V. Wüthrich and M. Merz.  _Stochastic claims reserving methods in insurance_ , volume 435.  John Wiley & Sons, 2008.
  * Ża̧dło [2013] T. Ża̧dło.  On parametric bootstrap and alternatives of mse.  In _Proceedings of the 31st International Scientific Conference Mathematical Method in Economics, (Hana Vojáčková ed.), College of Polytechnics Jihlava, Jihlava_ , pages 1081–1086, 2013.


# Large-Scale Self- and Semi-Supervised Learning for Speech Translation

###### Abstract

In this paper, we improve speech translation (ST) through effectively
leveraging large quantities of unlabeled speech and text data in different and
complementary ways. We explore both pretraining and self-training by using the
large Libri-Light speech audio corpus and language modeling with CommonCrawl.
Our experiments improve over the previous state of the art by 2.6 BLEU on
average on all four considered CoVoST 2 language pairs via a simple recipe of
combining wav2vec 2.0 pretraining, a single iteration of self-training and
decoding with a language model. Different to existing work, our approach does
not leverage any other supervision than ST data. Code and models will be
publicly released.

Index Terms: speech translation, unsupervised pretraining, self-training,
semi-supervised learning

11footnotetext: Equal contribution.

##  1 Introduction

Similar to many machine learning tasks, the amount of labeled data for speech-
to-text applications such as automatic speech recognition (ASR) and speech
translation (ST) is limited. For speech and language tasks, this problem is
exacerbated by the fact that people speak many languages around the world and
obtaining large quantities of labeled speech translation data for every
language is simply not practical. On the other hand, unlabeled speech audio or
text data is much more plentiful and various techniques to utilize it have
been explored.

Semi-supervised learning techniques for ASR such as unsupervised speech
pretraining [1, 2, 3, 4] and self-training or (iterative) pseudo-labeling [5,
6, 7, 8] have substantially improved performance on the traditional
LibriSpeech benchmark, and led to systems that can learn with little
supervision. In addition, it was recently shown that self-training and self-
supervised pretraining can be effectively combined [9, 10] for speech
recognition to achieve 4.8 WER on LibriSpeech test-other with only 10 minutes
of annotated speech data.

Similarly for speech translation, there has been work on combating data
scarcity, namely multitask learning [11, 12, 13], pretraining on ASR data [12,
14, 15, 16], data augmentation [17, 18, 19, 20], self-supervised pretraining
[21, 22], self-training [23] or multilingual speech translation [24, 25, 26,
27, 28]. However, multitask learning, pretraining, data augmentation and
multilingual ST rely on additional supervision provided by labeled ASR data,
machine translation (MT) data or ST data while self-training and self-
supervised pretraining exploit unlabeled source speech data. Unlike in ASR,
the complementarity of self-supervised learning and self-training (and other
semi-supervised techniques) has not been studied for the ST task.

On the other hand, past work in ST relied heavily on labeled data in the form
of speech-to-text translation data, ASR transcriptions [29] or MT sentence
pairs to improve performance [28]. In this work, we follow the trend of
leveraging purely unlabeled data to improve performance and show that this
obtains strong performance while leveraging only supervision from ST data.
Using wav2vec 2.0 pretraining, self-training and language model decoding, we
show that we can outperform previous work while leveraging much less
supervision.

Our contributions are as follows: we present a comprehensive study of the
impact of existing semi-supervised learning techniques on speech translation
and show that they greatly reduce the need for additional supervision in the
form of labeled ASR or translation parallel data. We show that our simple
approach obtains state-of-the-art results on all four language pairs we
evaluate on: English to German, English to Catalan, English to Arabic and
English to Turkish. We also conduct an ablation study on the impact of the
quantity of unlabeled data for self-training and self-supervised pre-training
in the context of ST.

In what follows, we describe the semi-supervised learning techniques and the
system we use for ST. Then we present our results on the CoVoST 2 benchmark
[27] on four language pairs and compare our work to the literature.

##  2 Learning from Unlabeled Data

In this section, we describe the techniques to leverage unlabeled speech or
text data which we use in this study.

###  2.1 Unsupervised Pretraining

Unsupervised pretraining has been very effective in multiple fields of machine
learning, including natural language processing [30], computer vision [31, 32]
and speech recognition [3]. In this work, we demonstrate the impact of
unsupervised pretraining for speech translation (ST) by leveraging a wav2vec
2.0 model pretrained on Libri-
Light111github.com/pytorch/fairseq/blob/master/examples/wav2vec, a unlabeled
dataset comprising 53K hours of English read audio books. The model is trained
by predicting the latent speech representations of masked time-steps using a
loss similar to SimCLR [32]. The latent speech representations are quantized
for the prediction task and there is a fixed number of latents stored in a
codebook.

For ST, we construct a sequence-to-sequence model with attention [33, 34] by
adding a randomly initialized decoder model on top of a wav2vec 2.0 encoder.
After pretraining, all parameters are fine-tuned on the CoVoST 2 ST data [27].
The decoder is also a Transformer model but smaller, with 7 layers and model
dimension 256, which we do not pretrain222Pretraining with a masked language
model did not improve performance in our setting.

###  2.2 Self-Training

Self-training is a semi-supervised learning method that first trains a teacher
model on labeled data. The teacher model is then used to synthetically
annotate unlabeled data in order to train a new student model on the
combination of labeled and pseudo-labeled data [35, 23, 36]. We further fine-
tune the student model on labeled data to alleviate the domain mismatch
between labeled and unlabeled data. Recent work showed that unsupervised
pretraining and self-training can be complementary for natural language
understanding [37] and speech recognition [9]. In this work, we adopt this
setting by considering self-training on top of a wav2vec 2.0 pretrained model
and show the complementarity of these two learning approaches for speech
translation.

###  2.3 Decoding with Language Model

With unsupervised pretraining and self-training, we leverage additional
unlabeled speech data to improve the performance of a ST system. We also make
use of monolingual text data in the target language to further improve
translation quality. Specifically, we train a language model (LM) on part of
the CommonCrawl data333http://data.statmt.org/cc-100/ that is in similar
domains as CoVoST 2 (§ 3.1). And then we combine ST model and LM scores at
every time step in beam search decoding (shallow fusion) [38].

##  3 Datasets and Training Details

In this section, we describe the datasets we use for speech translation (ST),
the setups for unsupervised pretraining, self-training as well as language
modeling. Then we give details on how we train our models.

###  3.1 Datasets

Table 1: CoVoST 2 dataset. We use four language pairs of the CoVoST 2 dataset
for our experiments, all from English due to the easy access of open-source
English unnanotated data in that source language.

Language pair | en-de | en-ca | en-ar | en-tr
---|---|---|---|---
Speech hours | 430h | 430h | 430h | 430h
Target utterances | 288K | 288K | 288K | 288K
Target words | 2.8M | 3.1M | 2.4M | 2.2M

CoVoST 2: speech translation data. CoVoST
2444https://github.com/facebookresearch/covost is a large-scale multilingual
speech translation corpus covering translations from 21 languages into English
and from English into 15 languages. This represents the largest open dataset
available to date from total volume and language coverage perspective.
Specifically, we cover four language directions from English to German (de),
Catalan (ca), Arabic (ar) and Turkish (tr), which contain 430 hours of
annotated data each. For simplicity, we choose to focus on speech translation
from English because of easy public access to unlabeled data in that language.

Libri-Light: pretraining and self-training data. We use the wav2vec 2.0
model555dl.fbaipublicfiles.com/fairseq/wav2vec/wav2vec_vox_new.pt pretrained
on the Libri-Light data [39], a dataset consisting of more than 60k hours of
unlabeled read speech data666https://github.com/facebookresearch/libri-light.
The dataset is derived from open-source audio books from the LibriVox project
and is the largest freely-available corpus of speech. Previous work in
unsupervised pretraining showed that using Libri-Light over LibriSpeech led to
better performance [3]. We use the 6000-hour subset of Libri-light (LV-6k) as
well as the full 60k hours dataset (LV-60k) for self-training. Specifically,
we train models on CoVoST 2 whose domain is read speech, and synthetically
label unannotated data from LV-60k - which is also read speech. We then fine-
tune a wav2vec 2.0 pretrained model - the student model - on both the
synthetically annotated data and the ground-truth CoVoST 2 training data.
Self-training is a good alternative to back-translation in the case of speech
translation, although the target data generated by the teacher model is not
real data. To remedy this issue, we leverage ground-truth target data through
language model decoding. Note that similar to NoisyStudent, we also inject
noise in our input during self-training through the same masking strategy as
wav2vec 2.0 [40, 3].

LibriSpeech: pretraining and self-training ablation data. In order to study
the effect of the amount of unlabeled data on both pretraining and self-
training, we use the 960 hours of LibriSpeech as a smaller-scale alternative
to Libri-light and LV-60k for wav2vec 2.0
pretraining777dl.fbaipublicfiles.com/fairseq/wav2vec/wav2vec_small.pt and for
self-training.

Table 2: BLEU on four language pairs of CoVoST-V2: English-German (en-de), English-Catalan (en-ca), English-Arabic (en-ar) and English-Turkish (en-tr). The results show that self-supervised pre-training on Librispeech (wav2vec-2.0) followed by self-training on the same data can improve performance. Using a language model during decoding improves performance further. Row | Model | en-de | en-ca | en-ar | en-tr | Avg
---|---|---|---|---|---|---
Baselines from CoVoST 2
1 |  Wang et al. [27] (w/o pre-ASR) | 13.6 | 20.2 | 8.7 | 8.9 | 12.9
2 |  Wang et al. [27] (w/ pre-ASR) | 16.3 | 21.8 | 12.1 | 10.0 | 15.1
Previous state-of-the-art results
3 |  Prev E2E SOTA [41, 27, 28] | 18.4 | 23.6 | 13.9 | 11.7 | 16.9
4 |  Cascade SOTA [28] | 19.4 | 25.0 | 14.3 | 11.7 | 17.6
5 |  Li et al. [28] (joint training) | 25.8 | 30.9 | 18.0 | 17.0 | 22.9
6 |  Li et al. [28] (+ extra MT data) | 26.6 | 30.4 | 18.6 | 16.3 | 23.0
Our results
8 | wav2vec-2.0 | 23.8 | 32.4 | 17.4 | 15.4 | 22.3
9 | wav2vec-2.0 + decoding w/ LM | 24.9 | 34.0 | 18.0 | 16.7 | 23.4
10 | wav2vec-2.0 + self-training (LV-60k) | 26.5 | 34.1 | 20.2 | 17.5 | 24.6
11 | wav2vec-2.0 + self-training (LV-60k) + decoding w/ LM | 27.2 | 35.6 | 20.8 | 18.9 | 25.6

CommonCrawl: language model data. In order to obtain language model data in
the right domain, we leverage CommonCrawl (CC) data from the CC100 corpus [42,
43] for the four target languages studied (German, Catalan, Arabic, Turkish).
First, we train 4-gram language models with the KenLM toolkit [44] on the
training set of the CoVoST 2 data and on the CC data and filter CC by
averaging the LM scores [45]. We then keep only one tenth of the original
data, and use transformer-based language models trained on both the CoVoST 2
training set and the additional CC data for language model rescoring.

###  3.2 Training Details

Unsupervised pretraining. For wav2vec 2.0 pretraining, we use the Large model,
which comprises 24 self-attention blocks with model dimension 1024, inner
dimension 4096 and 16 attention heads, resulting in a total of about 300M
parameters. The feature encoder contains seven blocks and the temporal
convolutions in each block have 512 channels with strides (5,2,2,2,2,2,2) and
kernel widths (10,3,3,3,3,2,2), resulting in a receptive field of about 25ms
and a stride of about 20ms.

Speech translation. We use a sequence-to-sequence model where the encoder is a
wav2vec 2.0 model with several layers of convolutions followed by a
Transformer network with 24 layers. The decoder is also a Transformer network
with 7 layers, an embedding size of 256, 4 attention heads and FFN dimension
of 2048. A 10K BPE vocabulary is built on the CoVoST 2 target text for each
target language. We train our model with Adam, a learning rate of 5e-5, label
smoothing with probability 0.1, an effective batch size of 6.4M tokens, layer
drop 0.05, a masking strategy similar to wav2vec 2.0 with mask length 5 and
mask probability 0.15. During the fine-tuning phase, the wav2vec 2.0 encoder
is frozen for 10K updates. Models are trained for 250K updates and the best
checkpoint is selected based on the BLEU score on the validation set. For
self-training, we use a learning rate of 3e-5 and the same setting for the
remaining hyperparameters. The pseudo-labels are generated with a beam size of
4.

Language models. For the language models, we use the Google Billion Word
Transformer architecture of [46] with 12 decoder layers and an embedding
dimension of 512. We train the model with Adam, and an inverse sqrt scheduler.
When using LM scores in decoding, we scale the scores by 0.1 and use a length
penalty of 0.7.

Table 3: Effect of the amount of unlabeled data on pretraining and self-training (cf. Table 2). Row | Model | en-de | en-ca | en-ar | en-tr | Avg
---|---|---|---|---|---|---
Increasing the amount of pretraining data
1 | wav2vec-2.0 (LS-960) | 20.5 | 27.3 | 15.4 | 14.0 | 19.3
2 | wav2vec-2.0 (LV-60k) | 23.8 | 32.4 | 17.4 | 15.4 | 22.3
Increasing the amount of self-training data
3 | wav2vec-2.0 + Self-training (LS-960) | 25.4 | 32.9 | 19.1 | 16.4 | 23.5
4 | wav2vec-2.0 + Self-training (LV-6k) | 25.9 | 33.2 | 19.9 | 17.5 | 24.1
5 | wav2vec-2.0 + Self-training (LV-60k) | 26.5 | 34.1 | 20.2 | 17.5 | 24.6

##  4 Experiments and Results

In this section, we describe our experimental results obtained in Table 2 and
3 where we combine self-supervised and semi-supervised learning techniques.

###  4.1 Improvements from Unsupervised Pretraining

We observe strong gains using wav2vec 2.0 models compared to previous
baselines which were using similar speech translation architectures without
pretraining. The Libri-light pretrained wav2vec 2.0 model (row 8) achieves
22.3 BLEU on average, which is on average 9.4 BLEU points better than the
baseline of [27] (row 1) and 7.2 BLEU higher than their model which leveraged
ASR pretraining as additional supervision (row 2). These results demonstrate
that the features learned by the wav2vec 2.0 model are very useful beyond
speech recognition, and applicable to other speech tasks such as speech
translation. The results with only wav2vec 2.0 pretraining (row 8) come very
close to the most recent state-of-the-art results (row 6) on CoVoST 2 [28]
which obtained 23 BLEU and even leads to a new state of the art on English-
Catalan (row 8), without using any other supervision than the speech
translation data.

###  4.2 Improvements from Self-Training

As previously shown for computer vision [47], natural language understanding
[37] and speech recognition [9], supervised or unsupervised pretraining can be
complementary to self-training. Combining self-supervised learning with semi-
supervised learning for speech translation in this work, we first use the
previously described Libri-Light wav2vec 2.0 models fine-tuned on the CoVoST 2
speech translation data as teacher models, and synthetically annotate the
60k-hour Libri-light dataset (LV-60k). We then leverage self-training by
finetuning a wav2vec 2.0 pretrained student model on both LV-60k synthetic
data and the CoVoST data [48, 23]. Note that we upsample the ground-truth
CoVoST data such that it has the same importance during finetuning than the
synthetically annotated LV-60k. After following this procedure, we obtain 24.6
average BLEU score (row 10) which is 2.3 BLEU better than using wav2vec 2.0
pretraining only, demonstrating the complementarity of pretraining and self-
training for speech translation. We also observe similar level of improvements
across language pairs which shows the consistency of this approach. With this
method, we reach a new state of the art on the CoVoST 2 benchmark for Catalan,
Arabic and Turkish (row 10).

###  4.3 Improvements from Decoding with Language Model

Self-training and unsupervised pretraining both leverage additional
unannotated speech data to improve performance. But self-training generates
noisy output on which the student model is fine-tuned which may lead the model
to learn incorrect patterns. To inject more prior knowledge about the target
language structure, a natural solution is to use unannotated text in the
target domain. In this work, we leverage language modeling as one way to do
that, and use it to improve generation through decoding. This improves the
wav2vec 2.0 baseline by 1.1 BLEU on average across all language pairs (row 9
vs. row 8). It also improves the stronger setting of wav2vec 2.0 + self-
training by 1 BLEU (row 11 vs. row 10). Combining wav2vec 2.0 pretraining,
self-training and language model decoding in row 11, we reach a new state of
the art on the CoVoST 2 benchmark, with an average BLEU score of 25.6 over the
four language pairs.

###  4.4 Comparison with Previous Work

The combination of pretraining, self-training and LM decoding outperforms the
prior state of the art [28] in all language directions and by 2.6 BLEU on
average. The prior state of the art uses both pretraining with wav2vec 2.0 and
mBART, as well as a minimalistic LNA (LayerNorm and Attention) finetuning. We
note that mBART was fine-tuned using additional labeled machine translation
data and therefore our approach relies on less labeled data. Moreover, this
work is conducted on a bilingual setting without additional supervision coming
from multiple language pairs. Finally, the prior state-of-the-art model uses
an adapter between the encoder and the decoder that further downsamples the
input by a factor of 8. In contrast, our model architecture was simplified by
removing the adapter module, as well as the LNA finetuning.

###  4.5 Data Ablation for Unsupervised Pretraining and Self-Training

Prior work has shown that increasing the amount of unlabeled data for
pretraining and self-training can improve ASR performance [3, 9]. To
understand whether the same holds true for ST, we compare pretraining a
wav2vec 2.0 model on the 960 hours unlabeled speech audio of the smaller
LibriSpeech corpus (LS-960) instead of the 60k hours of the LibriVox corpus
(LV-60k) used so far. Table 3 shows that pretraining the speech encoder on
more data leads to a large improvement of 3 BLEU on average across the four
language pairs (row 2 vs. row 1). Note that pretraining wav2vec 2.0 on
LibriSpeech still provides an average improvement of 4.2 BLEU over the
supervised baseline with ASR pretraining (row 1 vs. row 2 from Table 2).

Next, we examine increasing the amount of unlabeled data for self-training. We
compare pseudo-labeling LS-960 with pseudo-labeling LV-6k and LV-60k. From
Table 3 we can see that the increase of self-training data from LS-960 to
LV-6k brings an additional gain of 0.6 BLEU on average (row 4 vs. row 3).
Scaling self-training even further to LV-60k leads to an additional gain of
0.5 BLEU on average (row 5 vs. row 4). For both pretraining and self-training,
we found that gains using more unlabeled data were similar when using a
language model.

##  5 Conclusion

We pushed the limits of self-supervised and semi-supervised learning for
speech translation by leveraging pretraining with wav2vec 2.0 and self-
training. These techniques can outperform the previous state of the art by an
average of 1.3 BLEU across four language directions without using any type of
supervision other than the CoVoST 2 data. We also demonstrated the
complementarity of unsupervised pretraining, self-training and language model
decoding, outperforming previous approaches by 2.6 BLEU. Our work provides
stronger and simpler baselines for speech translation and demonstrates the
effectiveness of wav2vec 2.0 unsupervised pretraining for speech translation.

## References

  * [1] A. v. d. Oord, Y. Li, and O. Vinyals, “Representation learning with contrastive predictive coding,” _arXiv_ , vol. abs/1807.03748, 2018.
  * [2] S. Schneider, A. Baevski, R. Collobert, and M. Auli, “wav2vec: Unsupervised Pre-Training for Speech Recognition,” in _Proc. Interspeech 2019_ , 2019\.
  * [3] A. Baevski, Y. Zhou, A. Mohamed, and M. Auli, “wav2vec 2.0: A framework for self-supervised learning of speech representations,” in _Proc. NeurIPS_ , 2020.
  * [4] A. Conneau, A. Baevski, R. Collobert, A. Mohamed, and M. Auli, “Unsupervised cross-lingual representation learning for speech recognition,” _arXiv_ , vol. abs/2006.13979, 2020.
  * [5] G. S. et al., “End-to-end ASR: from supervised to semi-supervised learning with modern architectures,” _arXiv_ , vol. abs/1911.08460, 2019.
  * [6] J. Kahn, A. Lee, and A. Hannun, “Self-training for end-to-end speech recognition,” in _Proc. ICASSP_ , 2020.
  * [7] Q. Xu, T. Likhomanenko, J. Kahn, A. Hannun, G. Synnaeve, and R. Collobert, “Iterative pseudo-labeling for speech recognition,” _arXiv_ , vol. abs/2005.09267, 2020.
  * [8] T. Likhomanenko, Q. Xu, J. Kahn, G. Synnaeve, and R. Collobert, “slimipl: Language-model-free iterative pseudo-labeling,” _arXiv_ , vol. abs/2010.11524, 2020.
  * [9] Q. Xu, A. Baevski, T. Likhomanenko, P. Tomasello, A. Conneau, R. Collobert, G. Synnaeve, and M. Auli, “Self-training and pre-training are complementary for speech recognition,” in _Proc. ICASSP_ , 2020.
  * [10] Y. Zhang, J. Qin, D. S. Park, W. Han, C.-C. Chiu, R. Pang, Q. V. Le, and Y. Wu, “Pushing the limits of semi-supervised learning for automatic speech recognition,” _arXiv_ , vol. abs/2010.10504, 2020.
  * [11] R. J. Weiss, J. Chorowski, N. Jaitly, Y. Wu, and Z. Chen, “Sequence-to-sequence models can directly translate foreign speech,” _arXiv preprint arXiv:1703.08581_ , 2017.
  * [12] A. Bérard, L. Besacier, A. C. Kocabiyikoglu, and O. Pietquin, “End-to-end automatic speech translation of audiobooks,” in _Proc. ICASSP_ , 2018.
  * [13] Y. Tang, J. Pino, C. Wang, X. Ma, and D. Genzel, “A general multi-task learning framework to leverage text data for speech to text tasks,” in _Proc. ICASSP_ , 2021.
  * [14] S. Bansal, H. Kamper, K. Livescu, A. Lopez, and S. Goldwater, “Pre-training on high-resource speech recognition improves low-resource speech-to-text translation,” in _Proc. NAACL_ , 2019.
  * [15] M. C. Stoian, S. Bansal, and S. Goldwater, “Analyzing asr pretraining for low-resource speech-to-text translation,” in _Proc. ICASSP_ , 2020.
  * [16] C. Wang, Y. Wu, S. Liu, Z. Yang, and M. Zhou, “Bridging the gap between pre-training and fine-tuning for end-to-end speech translation,” in _Proc. AAAI_ , 2020.
  * [17] Y. Jia, M. Johnson, W. Macherey, R. J. Weiss, Y. Cao, C.-C. Chiu, N. Ari, S. Laurenzo, and Y. Wu, “Leveraging weakly supervised data to improve end-to-end speech-to-text translation,” in _Proc. ICASSP_ , 2019.
  * [18] J. Pino, L. Puzon, J. Gu, X. Ma, A. D. McCarthy, and D. Gopinath, “Harnessing indirect training data for end-to-end automatic speech translation: Tricks of the trade,” in _Proc. IWSLT_ , 2019.
  * [19] E. Salesky, M. Sperber, and A. W. Black, “Exploring phoneme-level speech representations for end-to-end speech translation,” in _Proc. ACL_ , 2019\.
  * [20] A. D. McCarthy, L. Puzon, and J. Pino, “Skinaugment: Auto-encoding speaker conversions for automatic speech translation,” in _Proc. ICASSP_ , 2020.
  * [21] A. Wu, C. Wang, J. Pino, and J. Gu, “Self-Supervised Representations Improve End-to-End Speech Translation,” in _Proc. Interspeech_ , 2020.
  * [22] H. Nguyen, F. Bougares, N. Tomashenko, Y. Estève, and L. Besacier, “Investigating Self-Supervised Pre-Training for End-to-End Speech Translation,” in _Proc. Interspeech_ , 2020.
  * [23] J. Pino, Q. Xu, X. Ma, M. J. Dousti, and Y. Tang, “Self-Training for End-to-End Speech Translation,” in _Proc. Interspeech_ , 2020.
  * [24] M. A. Di Gangi, M. Negri, and M. Turchi, “One-to-many multilingual end-to-end speech translation,” in _Proc. ASRU_ , 2019.
  * [25] H. Inaguma, K. Duh, T. Kawahara, and S. Watanabe, “Multilingual end-to-end speech translation,” in _Proc. ASRU_ , 2019.
  * [26] C. Wang, J. Pino, A. Wu, and J. Gu, “Covost: A diverse multilingual speech-to-text translation corpus,” in _Proc. LREC_ , 2020.
  * [27] C. Wang, A. Wu, and J. Pino, “Covost 2 and massively multilingual speech-to-text translation,” _arXiv_ , 2020.
  * [28] X. Li, C. Wang, Y. Tang, C. Tran, Y. Tang, J. Pino, A. Baevski, A. Conneau, and M. Auli, “Multilingual speech translation with efficient finetuning of pretrained models,” _arXiv_ , vol. abs/2010.12829, 2021.
  * [29] S. B. H. K. K. Livescu and A. L. S. Goldwater, “Pre-training on high-resource speech recognition improves low-resource speech-to-text translation,” in _Proc. NAACL_ , 2019.
  * [30] J. Devlin, M.-W. Chang, K. Lee, and K. Toutanova, “Bert: Pre-training of deep bidirectional transformers for language understanding,” in _Proc. NAACL_ , 2019.
  * [31] K. He, H. Fan, Y. Wu, S. Xie, and R. Girshick, “Momentum contrast for unsupervised visual representation learning,” _arXiv_ , vol. abs/1911.05722, 2019.
  * [32] T. Chen, S. Kornblith, M. Norouzi, and G. Hinton, “A simple framework for contrastive learning of visual representations,” _arXiv_ , vol. abs/2002.05709, 2020.
  * [33] I. Sutskever, O. Vinyals, and Q. V. Le, “Sequence to Sequence Learning with Neural Networks,” in _Proc. NIPS_ , 2014.
  * [34] J. K. Chorowski, D. Bahdanau, D. Serdyuk, K. Cho, and Y. Bengio, “Attention-based models for speech recognition,” in _Proc. NIPS_ , 2015\.
  * [35] D. S. Park, Y. Zhang, Y. Jia, W. Han, C.-C. Chiu, and et al., “Improved noisy student training for automatic speech recognition,” _arXiv_ , 2020.
  * [36] J. Kahn, A. Lee, and A. Hannun, “Self-training for end-to-end speech recognition,” in _Proc. ICASSP_ , 2020.
  * [37] J. Du, E. Grave, B. Gunel, V. Chaudhary, O. Celebi, M. Auli, V. Stoyanov, and A. Conneau, “Self-training improves pre-training for natural language understanding,” _arXiv_ , vol. abs/2010.02194, 2020.
  * [38] C. Gulcehre, O. Firat, K. Xu, K. Cho, L. Barrault, H.-C. Lin, F. Bougares, H. Schwenk, and Y. Bengio, “On using monolingual corpora in neural machine translation,” _arXiv preprint arXiv:1503.03535_ , 2015.
  * [39] J. Kahn, M. Rivière, W. Zheng, E. Kharitonov, Q. Xu, and et al., “Libri-light: A benchmark for asr with limited or no supervision,” _arXiv_ , vol. abs/1912.07875, 2019.
  * [40] D. S. Park, W. Chan, Y. Zhang, C.-C. Chiu, B. Zoph, E. D. Cubuk, and Q. V. Le, “Specaugment: A simple data augmentation method for automatic speech recognition,” in _Proc. Interspeech_ , 2019.
  * [41] J. Iranzo-Sánchez, J. A. Silvestre-Cerdà, and et al., “Europarl-st: A multilingual corpus for speech translation of parliamentary debates,” in _Proc. ICASSP_ , 2020.
  * [42] A. Conneau, K. Khandelwal, N. Goyal, V. Chaudhary, G. Wenzek, F. Guzmán, E. Grave, M. Ott, L. Zettlemoyer, and V. Stoyanov, “Unsupervised cross-lingual representation learning at scale,” in _Proc. ACL_ , 2020.
  * [43] G. Wenzek and et al., “CCNet: Extracting high quality monolingual datasets from web crawl data,” in _Proc. LREC_ , 2020.
  * [44] K. Heafield, I. Pouzyrevsky, and et al., “Scalable modified Kneser-Ney language model estimation,” in _Proc. ACL_ , 2013.
  * [45] R. C. Moore and W. Lewis, “Intelligent selection of language model training data,” _Proc. ACL_ , 2010.
  * [46] A. Baevski and M. Auli, “Adaptive input representations for neural language modeling,” _arXiv preprint arXiv:1809.10853_ , 2018.
  * [47] B. Zoph, G. Ghiasi, T.-Y. Lin, Y. Cui, and et al., “Rethinking pre-training and self-training,” _arXiv_ , vol. abs/2006.06882, 2020.
  * [48] J. He, J. Gu, J. Shen, and M. Ranzato, “Revisiting self-training for neural sequence generation,” in _ICLR_ , 2020.