# Data assimilation: A dynamic homotopy-based coupling approach

Sebastian Reich

###### Abstract.

Homotopy approaches to Bayesian inference have found widespread use especially
if the Kullback-Leibler divergence between the prior and the posterior
distribution is large. Here we extend one of these homotopy approach to
include an underlying stochastic diffusion process. The underlying
mathematical problem is closely related to the Schrödinger bridge problem for
given marginal distributions. We demonstrate that the proposed homotopy
approach provides a computationally tractable approximation to the underlying
bridge problem. In particular, our implementation builds upon the widely used
ensemble Kalman filter methodology and extends it to Schrödinger bridge
problems within the context of sequential data assimilation.

##  1\. Introduction

Sequential data assimilation interlaces dynamic processes with intermittent
partial state observations in order to provide reliable state estimates and
their uncertainties. A wide array of numerical methods have been proposed to
tackle this problem computationally. Popular methods include sequential Monte
Carlo, variational inference, and various ensemble Kalman filter formulations
[5, 8]. These methods can encounter difficulties whenever the predictive
distribution is incompatible with the incoming data; in other words whenever
the distance between the prior, as provided by the underlying stochastic
process, and the data informed posterior distribution is large. It has long
been realized that this challenge can be partially circumvented by altering
the underlying stochastic process through appropriate control terms or
modified proposal densities [5, 23, 16, 9]. Recently, the connection between
devising such control terms and Schrödinger bridge problems [4] has been made
explicit [16]. However, Schrödinger bridge problems are notoriously difficult
to solve numerically. The key contribution of this paper is to provide a
computationally tracktable (sub-optimal) solution via a novel extension of
established homotopy approaches [7, 14]. Similar to related homotopy
approaches for purely Bayesian inference, the solution of certain partial
differential equations (PDE) is required in order to find the desired control
terms [14, 25]. In line with standard ensemble Kalman filter (EnKF)
methodologies we approximate these PDEs via a constant gain approximation
[21]. We mention that there are also alternative approaches to sequential data
assimilation or inference which utilize ideas from optimal transportation; see
for example [15, 6, 20, 3].

The paper is structured as follows. The mathematical formulation of the data
assimilation problem, as considered in this paper, is laid out in Section 2.
The standard optimal control and Schrödinger bridge approach to data
assimilation is briefly summarized in Section 3, and the novel control
formulation based on an homotopy formulation is introduced in Section 4. A
practical implementation based on the EnKF methodology is proposed in Section
5. A series of increasingly complex data assimilation problems is considered
in Section 6 in order to demonstrate the feasibility of the proposed
methodologies. The paper concludes with some conclusions in Section 7.
Detailed mathematical derivations can be found in Appendix A and B,
respectively.

##  2\. Problem Formulation and Background

We consider drift diffusion processes given by a stochastic differential
equations (SDE)

(1) ${\rm d}X_{t}=f(X_{t}){\rm d}t+\sqrt{2\sigma}{\rm d}W_{t},$

where $X_{t}:\Omega\to\mathbb{R}^{d_{x}}$,
$f:\mathbb{R}^{d_{x}}\to\mathbb{R}^{d_{x}}$, $\sigma\in\mathbb{R}_{\geq 0}$,
and $W_{t}:\Omega\to\mathbb{R}^{d_{x}}$ denotes $d_{x}$-dimensional standard
Brownian motion [12, 13].

Assuming the law of $X_{t}$ is absolutely continuous w.r.t Lebesgue measure
with density $\pi_{t}$, this leads to the Fokker-Planck equation [13]

(2)
$\partial_{t}\pi_{t}=-\nabla\cdot\left(\pi_{t}\left(f-\sigma\nabla\log\pi_{t}\right)\right).$

The SDE (1) can be replaced by the mean field ODE

(3) $\frac{{\rm d}}{{\rm
d}t}\tilde{X}_{t}=f(\tilde{X}_{t})-\sigma\nabla\log\tilde{\pi}_{t}$

where $\tilde{\pi}_{t}$ denotes the law of $\tilde{X}_{t}$. Provided
$\tilde{\pi}_{0}=\pi_{0}$, it holds that $\tilde{\pi}_{t}=\pi_{t}$ for all
$t>0$. Note that the evolution of the random variable $\tilde{X}_{t}$ is
entirely deterministic subject to random initial conditions
$\tilde{X}_{0}\sim\pi_{0}$.

At time $t=T>0$, we have observations of the system according to

$y_{T}=h(x_{T}^{\dagger})+\nu$

from which we wish to infer the unknown state
$x_{T}^{\dagger}\in\mathbb{R}^{d_{x}}$. Here
$h:\mathbb{R}^{d_{x}}\to\mathbb{R}^{d_{y}}$ denotes the forward map and
$\nu\sim\mathcal{N}(0,R)$ is $d_{y}$-dimensional Gaussian noise with
covariance matrix $R\in\mathbb{R}^{d_{y}\times d_{y}}$.

Let $L:\mathbb{R}^{d_{x}}\to\mathbb{R}$ denote the corresponding negative log-
likelihood function. Since $\nu$ is Gaussian it is given by

(4) $L(x)=\frac{1}{2}(h(x)-y_{T})^{\top}R^{-1}(h(x)-y_{T}).$

The observations are combined with the predictive density $\pi_{T}$ at time
$t=T$ according to Bayes’ theorem,

(5) $\pi^{\rm a}_{T}=\frac{e^{-L}\pi_{T}}{\int e^{-L(x)}\pi_{T}(x){\rm d}x}.$

The process of transforming the random variable $X_{T}\sim\pi_{T}$ into a
random variable $X_{T}^{a}\sim\pi_{T}^{\rm a}$ is called data assimilation in
the context of dynamical systems and stochastic processes [10, 17, 8].

Since performing data assimilation can be difficult if the relative Kullback-
Leibler divergence

${\rm KL}(\pi_{T}|\pi_{T}^{\rm
a})=\int_{\mathbb{R}^{d_{x}}}\pi_{T}(x)(\log\pi_{T}(x)-\log\pi_{T}^{\rm
a}(x)){\rm d}x,$

also called the relative entropy [13], between the prior $\pi_{T}$ and
posterior $\pi_{T}^{a}$ is large and/or if the involved distributions are
strongly non-Gaussian [19, 1], we propose to construct a new SDE with state
process $X_{t}^{\rm h}$ such that $X_{0}^{\rm h}\sim\pi_{0}$ and $X_{T}^{\rm
h}\sim\pi_{T}^{a}$. In other words, we are looking for a stochastic process
(bridge) with initial density $\pi_{0}$ and final density $\pi^{\rm a}_{T}$.
The problem of finding the optimal process (in the sense of minimal Kullback-
Leibler divergence) is known as the Schrödinger bridge problem [4].

##  3\. Schrödinger Bridge Approach

The Bayesian adjustment (5) at final time $t=T$ leads in fact to an adjustment
over the whole solution space of the underlying diffusion process described by
(1). Let us denote the so called smoothing distribution by $\pi_{t}^{\rm a}$,
$t\in[0,T]$ [18, 5]. It is well established that these marginal distributions
can be generated from a controlled SDE

(6) ${\rm d}X_{t}^{\rm a}=f(X^{\rm a}_{t}){\rm d}t+g_{t}^{\rm a}(X^{\rm
a}_{t}){\rm d}t+\sqrt{2\sigma}{\rm d}W_{t}$

for appropriate control $g_{t}^{\rm
a}:\mathbb{R}^{d_{x}}\to\mathbb{R}^{d_{x}}$ such that $X_{0}^{\rm
a}\sim\pi_{0}^{\rm a}$ implies $X_{t}^{\rm a}\sim\pi_{t}^{\rm a}$ for all
$t>0$. It is also well known that finding a suitable $g_{t}^{\rm a}$ can be
formulated as an optimal control problem which in turn is closely related to
the backward Kolmogorov equation [13, 16]. Formulations related to (6) have
also been used in the context of sequential Monte Carlo methods [5].

As proposed in [16], an alternative perspective on sequential data
assimilation is provided by Schrödinger bridges. Given two marginal
distributions $q_{0}$ and $q_{T}$ and stochastic process $X_{t}$ (referred to
as the _reference process_), a Schrödinger bridge is another stochastic
process $\hat{X}_{t}$ such that $\hat{X}_{0}\sim q_{0}$, $\hat{X}_{T}\sim
q_{T}$ and the Kullback-Leibler divergence between the processes
$\\{\hat{X}_{t}\\}_{t\in[0,T]}$ and $\\{X_{t}\\}_{t\in[0,T]}$ is minimal. In
the language we have been using so far that means the marginals are the
initial and posterior densities, i.e. $q_{0}=\pi_{0}$ and $q_{T}=\pi^{\rm
a}_{T}$, and the reference process is the solution to (1). The solution to the
associated Schrödinger bridge problem is again of the form (6) with modified
control term denoted by $g_{t}^{\rm SB}(x)$.

A Schrödinger bridge is thus the optimal coupling as measured by the Kullback-
Leibler divergence to the underlying reference process. Unfortunately the
required control term $g_{t}^{\rm SB}$ seems rather difficult to compute in
practice. In addition to the computational complexity of solving nonlinear
Schrödinger bridge problems, the target distribution $\pi_{T}^{a}$ is
implicitly defined in the setting of data assimilation. The next section
offers a solution to both of these issues. We point to [23] for a discussion
of alternative approaches which introduce appropriate control terms into data
assimilation procedures.

##  4\. Homotopy Induced Dynamic Coupling

Since Schrödinger bridges are computationally challenging, we ask whether a
less optimal but cheaper approach might also be feasible. Indeed, in the
context of data assimilation a non-optimal coupling can be found via a
homotopy between the initial and target distribution as follows. Let

$\pi^{\rm h}_{t}(x)=Z^{-1}_{t}e^{-\frac{t}{T}L(x)}\pi_{t}(x)$

denote the homotopy in question, with $Z_{t}=\int
e^{-\frac{t}{T}L(x)}\pi_{t}(x){\rm d}x$ the time dependent normalization
constant. It clearly holds that $\pi^{\rm h}_{0}=\pi_{0}$ and $\pi^{\rm
h}_{T}=\pi_{T}^{\rm a}$. Note that the scaling $t\mapsto e^{-\frac{t}{T}L}$
was chosen for its simplicity and follows previous work on Bayesian inference
problems [7, 14]. Finding better homotopies or systematic ways of constructing
one could be an interesting direction for future research.

We can then reason backwards from the Fokker-Planck equation of $\pi^{h}_{t}$
to conclude that if it is the density of a random variable $X^{h}_{t}$ then
that random variable must satisfy the modified SDE:

(7) ${\rm d}X_{t}^{\rm h}=f(X^{\rm h}_{t}){\rm d}t-\frac{\sigma t}{T}\nabla
L(X^{\rm h}_{t}){\rm d}t+g_{t}(X^{\rm h}_{t}){\rm d}t+\sqrt{2\sigma}{\rm
d}W_{t},$

where $g_{t}$ is a solution to the PDE

(8) $\nabla\cdot(\pi^{\rm h}_{t}g_{t})=\frac{1}{T}\pi^{h}_{t}\left(L+t\nabla
L\cdot\left(f-\frac{\sigma t}{T}\nabla L-\sigma\nabla\log\pi^{\rm
h}_{t}\right)\right)+\pi^{\rm h}_{t}\frac{\dot{Z}_{t}}{Z_{t}}.$

The derivations of (8) can be found in Appendix A.

Note that (7) constitutes a mean field model since $g_{t}$ depends on the
distribution $\pi_{t}^{\rm h}$ of $X_{t}^{\rm h}$. We also wish to point out
that

(9) $\hat{g}_{t}^{\rm SB}(x):=-\frac{\sigma t}{T}\nabla L(x)+g_{t}(x)$

provides a (non-optimal) solution to the associated Schrödinger bridge
problem.

Since (8) is linear in $g_{t}$ we can decompose (8) into a set of simpler
equations $\nabla\cdot\left(\pi^{\rm h}_{t}g^{i}_{t}\right)=\pi^{\rm
h}_{t}(k^{i}-\mathbb{E}k^{i})$ such that the $k^{i}$ add up to the right hand
side of (8). In order to maintain $\int\nabla\cdot\left(\pi^{\rm
h}_{t}g_{t}\right)=0$ for the individual $g^{i}_{t}$ we can make use of the
fact that the terms in (8) are of the form $\pi^{\rm
h}_{t}\left(k-\mathbb{E}k\right)$. Separating the terms, we obtain the
following equations, the sum of whose solutions solves (8):

(10a) $\displaystyle\nabla\cdot\left(\pi^{\rm h}_{t}g_{t}^{1}\right)$
$\displaystyle=\pi^{\rm h}_{t}\left(\frac{L}{T}-\mathbb{E}\frac{L}{T}\right)$
(10b) $\displaystyle\nabla\cdot\left(\pi^{\rm h}_{t}g_{t}^{2}\right)$
$\displaystyle=\pi^{\rm h}_{t}\left(\frac{t}{T}\nabla L\cdot
f-\mathbb{E}\frac{t}{T}\nabla L\cdot f\right)$ (10c)
$\displaystyle\nabla\cdot\left(\pi^{\rm h}_{t}g_{t}^{3}\right)$
$\displaystyle=-\pi^{\rm h}_{t}\left(\frac{\sigma t}{T}\nabla
L\cdot\nabla\log\pi^{\rm h}_{t}-\mathbb{E}\frac{\sigma t}{T}\nabla
L\cdot\nabla\log\pi^{\rm h}_{t}\right)$ (10d)
$\displaystyle\nabla\cdot\left(\pi^{\rm h}_{t}g_{t}^{4}\right)$
$\displaystyle=-\pi^{\rm h}_{t}\left(\frac{\sigma t^{2}}{T^{2}}\nabla
L\cdot\nabla L-\mathbb{E}\frac{\sigma t^{2}}{T^{2}}\nabla L\cdot\nabla
L\right).$

Note that

$\pi^{\rm h}_{t}\nabla L\cdot\nabla\log\pi^{\rm h}_{t}=\nabla\left(\pi^{\rm
h}_{t}\nabla L\right)-\pi^{\rm h}_{t}\Delta L,$

which can be used to avoid the computation of $\nabla\log\pi_{t}^{\rm h}$
(with $\Delta=\nabla\cdot\nabla$ the Laplacian operator). Thus the controlled
SDE (7) can be replaced by

(11) ${\rm d}X_{t}^{\rm h}=f(X^{\rm h}_{t}){\rm d}t-\frac{2\sigma t}{T}\nabla
L(X^{\rm h}_{t}){\rm d}t+\hat{g}_{t}(X^{\rm h}_{t}){\rm d}t+\sqrt{2\sigma}{\rm
d}W_{t},$

where $\hat{g}_{t}$ is a solution to

(12) $\nabla\cdot(\pi^{\rm h}_{t}\hat{g}_{t})=\frac{1}{T}\pi^{\rm
h}_{t}\left(L+t\nabla L\cdot\left(f-\frac{\sigma t}{T}\nabla L\right)+\sigma
t\Delta L\right)+\pi^{\rm h}_{t}\frac{\dot{Z}_{t}}{Z_{t}}.$

Furthermore, if $\Delta L$ is a constant (as a function of $x$) or small in
comparison to the other contributions in (12), then (12) simplifies further.
In particular, this is the case if the forward map is linear, that is,
$h(x)=Hx$.

Building upon the mean field ODE (3), one obtains the equivalent controlled
mean field ODE system

$\frac{{\rm d}}{{\rm d}t}\tilde{X}_{t}^{\rm h}=f(\tilde{X}^{\rm
h}_{t})-\frac{\sigma t}{T}\nabla L(\tilde{X}^{\rm
h}_{t})-\sigma\nabla\log\tilde{\pi}^{\rm h}_{t}(\tilde{X}_{t}^{\rm
h})+g_{t}(\tilde{X}^{\rm h}_{t}),$

with $g_{t}$ defined as before. This mean field formulation again requires
knowledge (or approximation) of $\nabla\log\tilde{\pi}_{t}^{\rm h}$. A
Gaussian approximation might be sufficient in certain circumstances giving
rise to

$\nabla\log\tilde{\pi}_{t}^{\rm h}(x)\approx-(\tilde{\Sigma}_{t}^{\rm
h})^{-1}(x-\tilde{\mu}_{t}^{\rm h})$

where $\tilde{\mu}_{t}^{\rm h}$ denotes the mean of $\tilde{X}_{t}^{\rm h}$
and $\tilde{\Sigma}_{t}^{\rm h}$ its covariance matrix.

##  5\. Numerical implementation

No analytic solution to (8) is known and we thus have to resort to
approximations. We note that a similar PDE arises in the computation of the
gain in the feedback particle filter [24] and one could use the diffusion map
based approximation [22] for the problem at hand. This method also transforms
the PDE into a Poisson equation which it then translated into an equivalent
integral equation, the semi-group form of the Poisson equation. As the name
suggests the integral equation makes use of the generator of a semi-group
which can be approximated by diffusion maps. Here we instead propose to follow
the constant gain approximation first introduced in the EnKF methodology [21].

###  5.1. Ensemble Kalman mean field approximation

Let us assume that $\Delta L\approx const$ in (12). Then we only need to deal
with the modified negative log likelihood function

$\displaystyle\tilde{L}(x)$ $\displaystyle=\frac{1}{T}\left(L(x)+t\nabla
L(x)\cdot\left(f(x)-\frac{\sigma t}{T}\nabla L(x)\right)\right)$
$\displaystyle\approx\frac{1}{T}\left(L(x)+\frac{t}{\Delta
t}\left\\{L(x)-L\left(x-\Delta tf(x)+\Delta t\frac{t\sigma}{T}\nabla
L(x)\right)\right\\}\right)$

with $\Delta t$ being the time-step also used later for time-stepping the
evolution equations (7) or (11), respectively. Since $L$ is given by (4), we
define the modified forward map

$\tilde{h}(x)=h\left(x-\Delta tf(x)+\Delta t\frac{\sigma t}{T}\nabla
L(x)\right)$

and thus

(13) $\tilde{L}(x)\approx\frac{t+\Delta t}{2\Delta
tT}(h(x)-y_{T})^{\top}R^{-1}(h(x)-y_{T})-\frac{t}{2\Delta
tT}(\tilde{h}(x)-y_{T})^{\top}R^{-1}(\tilde{h}(x)-y_{T}).$

Following the standard EnKF methodology, this suggests to approximate the
drift function $\hat{g}_{t}$ in (11) as follows:

(14) $\hat{g}_{t}^{\rm KF}(x)=-\frac{t+\Delta t}{\Delta
tT}\Sigma_{t}^{xh}R^{-1}\left(\frac{1}{2}\left(h(x)+\pi_{t}^{\rm
h}[h]\right)-y_{T}\right)+\frac{t}{\Delta
tT}\Sigma_{t}^{x\tilde{h}}R^{-1}\left(\frac{1}{2}\left(\tilde{h}(x)+\pi_{t}^{\rm
h}[\tilde{h}]\right)-y_{T}\right).$

Here we have introduced the notation $\pi_{t}^{\rm h}[l]$ to denote the
expectation value $\mathbb{E}l$ of a function $l(x)$ under the PDF
$\pi_{t}^{\rm h}$. Furthermore, $\Sigma_{t}^{xh}$ denotes the correlation
matrix between $x$ and $h(x)$ under the PDF $\pi_{t}^{\rm h}$ etc. The
derivation of (14) can be found in Appendix B.

###  5.2. Particle approximation and time-stepping

The controlled mean field equations (11) can be implemented numerically by the
standard Monte Carlo ansatz, that is, $M$ particles $X_{t}^{(i)}$ are
propagated according to

(15) ${\rm d}X_{t}^{(i)}=f(X^{(i)}_{t}){\rm d}t-\frac{2\sigma t}{T}\nabla
L(X^{(i)}_{t}){\rm d}t+\hat{g}_{t}^{\rm KF}(X^{(i)}_{t}){\rm
d}t+\sqrt{2\sigma}{\rm d}W_{t}^{(i)}$

for $i=1,\ldots,M$. The required expectation values in $\hat{g}_{t}^{\rm KF}$
are evaluated with respect to the empirical measure

$\hat{\pi}_{t}^{\rm h}(x)=\frac{1}{2}\sum_{i=1}^{M}\delta(x-X_{t}^{(i)}).$

The interacting particle system can be time-stepped using an appropriate
adaptation of (27) from Appendix B. The computation of gradients can be
avoided by applying the statistical linearisation (26).

##  6\. Examples

We now discuss a sequence of increasingly complex examples. The purpose is
both to illuminate certain aspects of the proposed control terms as well as to
indicate the computational advantages of the proposed methodology. All
examples will be based on linear forward maps $h(x)=Hx$ and, therefore $\Delta
L$ is constant and can be ignored.

###  6.1. Pure diffusion processes

We set the drift $f$ to zero in (1) and also assuming Gaussian initial
conditions. Then the homotopy approach gives rise to the controlled SDE

(16a) $\displaystyle{\rm d}X_{t}^{\rm h}$ $\displaystyle=\sqrt{2\sigma}{\rm
d}W_{t}-\frac{2\sigma t}{T}H^{\top}R^{-1}(HX_{t}^{\rm h}-y_{T}){\rm d}t$ (16b)
$\displaystyle\qquad-\,\,\Sigma_{t}^{\rm
h}H^{\top}\left\\{\frac{1}{T}R^{-1}-\frac{2\sigma
t^{2}}{T^{2}}R^{-1}HH^{\top}R^{-1}\right\\}\left(\frac{1}{2}\left(HX_{t}^{\rm
h}+H\mu_{t}^{\rm h}\right)-y_{T}\right){\rm d}t.$

We note that $X_{t}^{\rm h}\sim\pi_{t}^{\rm h}$ will remain Gaussian for all
times and we denote the mean by $\mu_{t}^{\rm h}$ and the covariance matrix by
$\Sigma_{t}^{\rm h}$. Hence, it holds that $\Sigma_{t}^{xx}=\Sigma_{t}^{\rm
h}$ and $\pi_{t}^{\rm h}[x]=\mu_{t}^{\rm h}$.

Please note that the additional drift term in (16a) is pulling $X_{t}^{\rm h}$
towards the observation $y_{T}$ regardless of the value of $\Sigma_{t}^{\rm
h}$. It should also be noted that the drift term in (16b) can be both
attractive or repulsive with regard to the observation $y_{T}$ depending on
the eigenvalues of

$\Omega(t)=\frac{1}{T}R^{-1}-\frac{2\sigma
t^{2}}{T^{2}}R^{-1}HH^{\top}R^{-1}.$

The strength of this drift term is moderated by the covariance matrix
$\Sigma_{t}^{\rm h}$.

We consider a one-dimensional problem with $R=0.01$, $\sigma=1$, $H=1$,
$y_{T}=1$ and $T=1$. The initial conditions are Gaussian with mean $\mu_{0}=0$
and variance $\Sigma_{0}=1$. It follows that $\pi_{1}$ is Gaussian with mean
$\mu_{1}=0$ and variance $\Sigma_{1}=2$ and the resulting Gaussian posterior
$\pi_{1}^{a}$ has mean and variance given by

(17) $\mu_{1}^{a}=Ky_{T}\approx 0.9524,\qquad\mu_{1}^{a}=2-2K\approx 0.0952$

with Kalman gain $K=2/(2+0.01)\approx 0.9524$.

Figure 1. Time evolution of the mean $\mu_{t}^{h}$ and the variance
$\Sigma_{t}^{h}$ under the mean field equations (16). Their values at final
agree with the posterior values provided by (17).

In Figure 1 one can find the time evolution of the mean and the variance under
the mean field equations (16). The early impact of the data driven control
term on the dynamics is perhaps surprising and quite opposite to the standard
sequential approach to data assimilation where one first propagates to final
time and only then adjusts according to the available data. It is also worth
noticing that the sign of the corresponding $\Omega(t)$ changes sign at
$t_{c}=\sqrt{2}/20$ implying that the drift term in (16) has a destabilizing
effect on the dynamics for $t>t_{c}$.

###  6.2. Purely deterministic processes

We now set $\sigma=0$ in (1). We obtain from (14) the mean field ODE system

(18a) $\displaystyle\frac{\rm d}{{\rm d}t}X_{t}^{\rm h}$
$\displaystyle=f(X_{t}^{\rm h})-\frac{t+\Delta t}{\Delta tT}\Sigma_{t}^{\rm
h}H^{\top}R^{-1}\left(\frac{1}{2}H\left(X_{t}^{\rm h}+\mu_{t}^{\rm
h}\right)-y_{T}\right)$ (18b) $\displaystyle\qquad+\,\,\frac{t}{\Delta
tT}\Sigma_{t}^{x\tilde{h}}R^{-1}\left(\frac{1}{2}\left(\tilde{h}(X_{t}^{\rm
h})+\pi_{t}^{\rm h}[\tilde{h}]\right)-y_{T}\right)$

with

$\tilde{h}(x)=Hx-\Delta tHf(x).$

These equations can be expanded giving rise to

(19a) $\displaystyle\frac{\rm d}{{\rm d}t}X_{t}^{\rm h}$
$\displaystyle=f(X_{t}^{\rm h})-\frac{1}{T}\left\\{\Sigma_{t}^{\rm
h}+t\Sigma_{t}^{xf}\right\\}H^{\top}R^{-1}\left(\frac{1}{2}H\left(X_{t}^{\rm
h}+\mu_{t}^{\rm h}\right)-y_{T}\right)$ (19b)
$\displaystyle\qquad-\,\,\frac{t}{2T}\Sigma_{t}^{\rm
h}H^{\top}R^{-1}H\left(f(X_{t}^{\rm h})+\pi_{t}^{\rm h}[f])\right)$

upon ignoring terms of order $\mathcal{O}(\Delta t)$. Unless the drift
function $f$ is linear, these mean field equations provide only an
approximation to the controlled mean field equations (11).

###  6.3. Linear Gaussian case

It is instructive to investigate the linear case

(20) $f(x)=Fx+b$

in more detail where again everything remains Gaussian provided $X_{0}^{h}$ is
Gaussian distributed, that is $\pi_{0}(x)=\mathcal{N}(x;\mu_{0},\Sigma_{0})$.
Under these conditions the densities $\pi_{t}$ and $\pi^{\rm h}_{t}$ will also
be Gaussian, we write $\pi^{\rm h}_{t}(x)=\mathcal{N}(x;\mu^{\rm
h}_{t},\Sigma^{\rm h}_{t})$. The associated mean field equations follow from
Appendix B and are given by

(21a) $\displaystyle\frac{{\rm d}}{{\rm d}t}X_{t}^{\rm h}$
$\displaystyle=FX_{t}^{\rm h}+b+\sigma(\Sigma_{t}^{\rm h})^{-1}(X_{t}^{\rm
h}-\mu_{t}^{\rm h})-\frac{2\sigma t}{T}H^{\top}R^{-1}(HX_{t}^{\rm h}-y_{T})$
(21b) $\displaystyle\quad-\,\,\left\\{\frac{1}{T}\Sigma_{t}^{\rm
h}+\frac{t}{T}\Sigma_{t}^{\rm h}F^{\top}-\frac{2\sigma
t^{2}}{T^{2}}\Sigma_{t}^{\rm
h}H^{\top}R^{-1}H\right\\}H^{\top}R^{-1}\left(\frac{1}{2}H\left(X_{t}^{\rm
h}+\mu_{t}^{\rm h}\right)-y_{T}\right)$ (21c)
$\displaystyle\quad-\,\,\frac{t}{T}\Sigma_{t}^{h}H^{\top}R^{-1}H\left(\frac{1}{2}F\left(X_{t}^{\rm
h}+\mu_{t}^{\rm h}\right)+b\right).$

A qualitative discussion can be performed in the scalar case, that is
$d_{x}=1$, $H=1$, $\sigma=1$, $b=0$, $T=1$ and $F=\lambda$. One finds that the
control terms involving $F$ stabilize the dynamics whenever $\lambda>0$. This
observation is in line with the fact that the data is crucial only if the
dynamics in $X_{t}$ is unstable, that is, $\lambda>0$.

We consider a two dimensional diffusion process with state variable
$x=(x_{1},x_{2})^{\top}$ and linear drift term (20) given by

$F=\big{(}\begin{smallmatrix}-2&1\\\ 1&-2\end{smallmatrix}\big{)},$

$b=0$, and diffusion constant $\sigma=0.1I$. The forward operator is
$H=\begin{pmatrix}1&0\end{pmatrix}$ and the variance of the noise $R=0.01$.
The initial distribution was $\pi_{0}=\mathcal{N}((1,3),0.02I)$. The observed
value at time $T=1$ is set to $y_{T}=2.5$.

The posterior mean takes values $\mu_{1}^{\rm a}\approx 2.25$ and
$\mu_{2}^{\rm a}\approx 1.50$, while the posterior covariance matrix becomes

$\Sigma^{\rm a}\approx\left(\begin{array}[]{cc}0.0086&0.0039\\\
0.0039&0.0503\end{array}\right).$

Numerical results can be found in Figure 2. The impact of the control term on
the linear diffusion process can clearly be seen and is most prominent on the
observed $x_{1}$ component of the process. The final values of the controlled
process agree well with their posterior counterparts.

Figure 2. Left panel: Time evolution of the mean in $x_{1}$ and $x_{2}$, the
two associated variances and the covariance between $x_{1}$ and $x_{2}$ under
the linear diffusion process. Right panel: Time evolution of the the same
quantities under the controlled diffusion process.

###  6.4. Nonlinear diffusion example

We consider a two-dimensional problem and denote the state variable by
$x=(x_{1},x_{2})^{\top}$. The drift term is given by

$f(x)=-\nabla V(x),\qquad V(x)=\frac{\lambda_{1}}{2}\left(x_{2}-2+\beta
x_{1}^{2}\right)^{2}+\frac{\lambda_{2}}{2}\left(\frac{x_{1}^{4}}{2}-x_{1}^{2}\right)$

with parameters $\lambda_{1}=2000$, $\lambda_{2}=5$, and $\beta=1/5$. The
diffusion constant is set to $\sigma=1$. The choice of the potential $V(x)$
has two effects: (i) there is a relative high barrier for particles to pass
from positive to negative $x_{1}$-values and vice versa; (ii) the dynamics
stay close to the parabola $x_{2}=2-\beta x_{1}^{2}$.

The initial distribution is obtained by sampling $x_{1}$ from a Gaussian with
mean $1.5$ and variance $0.0625$. The $x_{2}$ component is obtained from the
relation

$x_{2}=2-\beta x_{1}^{2}.$

We observe the first component $x_{1}$ of the state vector at time $T=1$ with
measurement error variance $R=0.01$. The observed value is set to
$y_{T}=-1.5$. Due to the tiny observation error the posterior is centered
sharply about the observed value. Furthermore, recall that the dynamics is
essentially slaved to the parabola $x_{2}=2-\beta x_{1}^{2}$ which makes the
inference problem strongly nonlinear.

All particle simulations are run with an ensemble size of $M=1000$.
Essentially identical results are obtained for $M=100$. Smaller ensemble sizes
lead to numerical instabilities.

Figure 3.  Initial (blue) and final particle positions (red) under the given
evolution process together with the posterior approximation provided by the
EnKF (yellow). The observed value is also displayed.

In Figure 3, one can find the particle distribution at time $t=1$ which
constitutes the prior distribution for the associated Bayesian inference
problem. It is obvious that a particle filter would fail to recover the
posterior distribution which is sharply centered about the observed value. We
found that increasing the ensemble size to $M=10000$ allows a particle filter
to recover the posterior distribution; but the effective ensemble size still
drops dramatically. The approximation provided by the EnKF is also displayed.
The EnKF fails to recover the posterior due to its inherent linear regression
ansatz which is inappropriate for this strongly nonlinear inference problem
even in the limit ensemble size $M\to\infty$.

Figure 4.  Left panel: Initial and final particle positions under the
controlled evolution process. Right panel: Particle positions at intermediate
times $t_{k}\in[0,1]$.

In Figure 4, the results from the controlled mean field formulation are
displayed. It can be concluded that the posterior distribution is will
approximated despite the constant gain approximation made in order to
formulate the control term $\hat{g}_{t}^{\rm KF}$ in (14).

###  6.5. Lorenz-63 example

For a more challenging example we use the standard Lorenz-63 model [11]

(22) $\frac{{\rm d}}{{\rm d}t}X_{t}=f(X_{t})$

where $X_{t}:\Omega\to\mathbb{R}^{3}$ and

$f(x,y,z)=\left(\begin{array}[]{c}a(y-x)\\\ x(b-z)-y\\\ xy-
cz\end{array}\right)$

with parameters $a=10$, $b=28$ and $c=8/3$.

In order to obtain a reference solution $X_{t}^{\dagger}$ for $t\geq 0$, the
ODE (22) is solved numerically with step-size $\Delta t=0.005$ and initial
condition

$X_{0}^{\dagger}=\left(\begin{array}[]{c}-0.587276\\\ -0.563678\\\
16.8708\end{array}\right).$

Scalar-valued observations are generated every $\Delta t_{\rm obs}>0$ units of
time using the forward model

$y_{n\Delta t_{\rm obs}}=HX_{n\Delta t_{\rm obs}}+\nu_{n},\qquad n=1,\ldots
N,$

with measurement measurement errors $\nu_{n}\sim{\rm N}(0,1)$ and forward map
$H=(1\,0\,0)\in\mathbb{R}^{1\times 3}$. We use $\Delta t_{\rm
obs}\in\\{0.05,0.1,0.12\\}$ in our experiments and perform $N=20,000$
assimilation cycles.

The initial ensemble $\\{X_{0}^{(i)}\\}_{i=1}^{M}$ is drawn from the Gaussian
distribution ${\rm N}(X_{0}^{\dagger},0.01I)$. We employ multiplicative
ensemble inflation which amounts to replacing the Lorenz-63 dynamics by

$\frac{{\rm d}}{{\rm
d}t}X_{t}^{(i)}=f(X_{t}^{(i)})+\sigma_{k}(X_{t}^{(i)}-\hat{\mu}_{t}),\qquad
i=1,\ldots,M,$

with inflation factors

$\sigma_{k}=0.025k,\qquad k=0,\ldots,9.$

Here $\hat{\mu}_{t}$ denotes the empirical mean of the ensemble
$\\{X_{t}^{(i)}\\}_{i=1}^{M}$. These equations are combined with the augmented
evolution equations (19) and solved numerically with step-size $\Delta
t=0.005$ and ensemble sizes $M\in\\{5,10,15\\}$.

We report the resulting root mean square errors

$\mbox{RMSE}=\sqrt{\frac{1}{3N}\sum_{n=1}^{N}\|\hat{\mu}_{n\Delta t_{\rm
obs}}-X_{n\Delta t_{\rm obs}}^{\dagger}\|^{2}},$

which are computed for each ensemble size $M$, observation interval $\Delta
t_{\rm obs}$ and inflation factor $\sigma_{k}$. The results are displayed in
Table 1 where the smallest RMSE over the range of inflation factors
$\\{\sigma_{k}\\}_{k=0}^{9}$ is stated for each $M$ and $\Delta t_{\rm obs}$.
We also state the corresponding RMSEs from a standard ensemble square root
filter implementation [2, 8]. We find that the proposed homotopy approach
outperforms the ensemble square root filter in terms of RMSE in all settings
considered. The improvements increase for increasing observation intervals
$\Delta t_{\rm obs}$. The homotopy approach also appears less sensitive to the
ensemble size $M$.

$M$/$\Delta t_{\rm obs}$ | 5 | 10 | 15
---|---|---|---
0.05 | 0.5712/0.5457 | 0.5620/0.5475 | 0.5659/0.5496
0.10 | 0.8466/0.7735 | 0.8171/0.7627 | 0.8229/0.7707
0.12 | 0.9606/0.8645 | 0.9515/0.8621 | 0.9375/0.8615
Table 1. Displayed are the RMSE for both a standard ensemble square root
filter implementation (top) and our homotopy approach (bottom) in terms of
ensemble sizes $M\in\\{5,10,15\\}$ and observation intervals $\Delta t_{\rm
obs}\in\\{0.05,0.1,0.12\\}$. The homotopy based data assimilation method leads
to significantly reduced RMSEs in all settings.

##  7\. Conclusions

Devising alternative proposal densities has a long history in the context of
sequential data assimilation and filtering. Here we have explored a
computationally tracktable approach which combines the concept of Schrödinger
bridges with a rather straightforward homotopy approach. A further key
ingredient is the approximate solution of the arising PDEs in terms of
constant gain approximations, which are also widely used within the EnKF
community. Numerical examples indicate that the approach is viable and can
overcome limitations of both standard sequential Monte Carlo as well as
standard EnKF methods. This has been demonstrated for single assimilation
steps as well as long-time data assimilation using the chaotic Lorenz-63 model
with only the first component observed infrequently. It remains to be seen how
the proposed methods behave for high dimensional stochastic processes.

Acknowledgment. This research has been funded by the Deutsche
Forschungsgemeinschaft (DFG)- Project-ID 318763901 - SFB1294. We thank Nikolai
Zaki for earlier work on the topic of this paper.

## References

  * Agapiou et al. [2017] S. Agapiou, O. Papaspiliopoulos, D. Sanz-Alonso, and A. Stuart.  Importance sampling: Intrinsic dimension and computational cost.  _Statistical Science_ , pages 405–431, 2017.
  * Asch et al. [2016] M. Asch, M. Bocquet, and M. Nodet.  _Data assimilation: Methods, algorithms, and applications_.  SIAM, 2016.
  * Bernton et al. [2019] E. Bernton, J. Heng, A. Doucet, and P. E. Jacob.  Schrödinger bridge samplers.  Technical report, arXiv:1912.13170, 2019.
  * Chen et al. [2021] Y. Chen, T. T. Georgiou, and M. Pavon.  Stochastic control liaisons: Richard Sinkhorn meets Gaspard Monge on a Schrödinger bridge.  _SIAM Review_ , 63:249–313, 2021.  doi: 10.1137/20M1339982.
  * Chopin and Papaspiliopoulos [2020] N. Chopin and O. Papaspiliopoulos.  _An Introduction to Sequential Monte Carlo_.  Springer Nature Switzerland AG, Cham, Switzerland, 2020.
  * Corenflos et al. [2021] A. Corenflos, J. Thornton, G. Deligiannidis, and A. Doucet.  Differentiable particle filtering via entropy-regularized optimal transport.  In _International Conference on Machine Learning_ , pages 2100–2111. PMLR, 2021.
  * Daum et al. [2010] F. Daum, J. Huang, and A. Noushin.  Exact particle flow for nonlinear filters.  In I. Kadar, editor, _Signal Processing, Sensor Fusion, and Target Recognition XIX_ , volume 7697, pages 92 – 110. International Society for Optics and Photonics, SPIE, 2010.  doi: 10.1117/12.839590.
  * Evensen et al. [2022] G. Evensen, F. Vossepoel, and P. van Leeuwen.  _Data Assimilation Fundamentals: A unified Formulation of the State and Parameter Estimation Problem_.  Springer Nature Switzerland AG, Cham, Switzerland, 2022.
  * Heng et al. [2020] J. Heng, A. N. Bishop, G. Deligiannidis, and A. Doucet.  Controlled sequential monte carlo.  _The Annals of Statistics_ , 48:2904–2929, 2020.  doi: 10.1214/19-AOS1914.
  * Law et al. [2015] K. Law, A. Stuart, and K. Zygalakis.  Data assimilation.  _Cham, Switzerland: Springer_ , 214, 2015.
  * Lorenz [1963] E. Lorenz.  Deterministic non-periodic flows.  _J. Atmos. Sci._ , 20:130–141, 1963.
  * Oksendal [2013] B. Oksendal.  _Stochastic differential equations: an introduction with applications_.  Springer Science & Business Media, 2013.
  * Pavliotis [2014] G. A. Pavliotis.  _Stochastic Processes and Applications: Diffusion Processes, the Fokker-Planck and Langevin Equations_ , volume 60 of _Texts in Applied Mathematics_.  Springer New York, New York, NY, 2014.  doi: 10.1007/978-1-4939-1323-7.
  * Reich [2011] S. Reich.  A dynamical systems framework for intermittent data assimilation.  _BIT Numerical Mathematics_ , 51(1):235–249, 2011\.
  * Reich [2013] S. Reich.  A nonparametric ensemble transform method for Bayesian inference.  _SIAM J. Sci. Comput._ , 35:A2013–A2024, 2013.  doi: 10.1137/130907367.
  * Reich [2019] S. Reich.  Data assimilation: The Schrödinger perspective.  _Acta Numerica_ , 28:635–711, 2019.  doi: 10/gkzmrc.
  * Reich and Cotter [2015] S. Reich and C. Cotter.  _Probabilistic forecasting and Bayesian data assimilation_.  Cambridge University Press, 2015.
  * Särkkä [2013] S. Särkkä.  _Bayesian Filtering and Smoothing_.  Cambridge University Press, Cambridge, 2013.
  * Snyder et al. [2008] C. Snyder, T. Bengtsson, P. Bickel, and J. Anderson.  Obstacles to high-dimensional particle filtering.  _Monthly Weather Review_ , 136(12):4629–4640, 2008.
  * Spantini et al. [2019] A. Spantini, R. Baptista, and Y. Marzouk.  Coupling techniques for nonlinear ensemble filtering.  _arXiv preprint arXiv:1907.00389_ , 2019.
  * Taghvaei et al. [2017] A. Taghvaei, J. de Wiljes, P. Mehta, and S. Reich.  Kalman filter and its modern extensions for the continuous-time nonlinear filtering problem.  _ASME. J. Dyn. Sys., Meas., Control._ , 140:030904–030904–11, 2017.  doi: 10.1115/1.4037780.
  * Taghvaei et al. [2020] A. Taghvaei, P. Mehta, and S. Meyn.  Diffusion map-based algorithm for gain function approximation in the feedback particle filter.  _SIAM/ASA J. Uncertain. Quantif._ , 8(3):1090–1117, 2020.  doi: 10.1137/19M124513X.
  * van Leeuwen et al. [2019] P. van Leeuwen, H. Künsch, L. Nerger, R. Potthast, and S. Reich.  Particle filters for high-dimensional geoscience applications: A review.  _Q. J. Royal Meteorol. Soc._ , 145:2335–2365, 2019.  doi: 10.1002/qj.3551.
  * Yang et al. [2013] T. Yang, P. G. Mehta, and S. P. Meyn.  Feedback particle filter.  _IEEE Trans. Automat. Control_ , 58(10):2465–2480, 2013.  ISSN 0018-9286.  doi: 10.1109/TAC.2013.2258825.
  * Yang et al. [2014] T. Yang, H. A. P. Blom, and P. G. Mehta.  The continuous-discrete time feedback particle filter.  pages 648–653. IEEE, 2014.  doi: 10.1109/ACC.2014.6859259.

##  Appendix A Derivation of Control Term Equation

For an evolution equation

(23) $dX_{t}=f(X_{t})dt+\sqrt{2\sigma}dW_{t}$

we get a Fokker-Planck equation

(24)
$\partial_{t}\pi_{t}=-\nabla\cdot\left(\pi_{t}\left(f-\sigma\nabla\log\pi_{t}\right)\right).$

Now we modify (23) by an additional drift term, i.e.

$dX^{\rm h}_{t}=f(X^{\rm h}_{t})dt+\tilde{g}_{t}(X^{\rm
h}_{t})dt+\sqrt{2\sigma}dW_{t}$

with $\tilde{g}_{t}:\mathbb{R}^{d_{x}}\to\mathbb{R}^{d_{x}}$. In that case, we
would get a Fokker-Planck equation for the new equation:

(25) $\partial_{t}\pi^{\rm h}_{t}=-\nabla\cdot\left(\pi^{\rm
h}_{t}\left(f+\tilde{g}-\sigma\nabla\log\pi^{\rm
h}_{t}\right)\right)=-\nabla\cdot\left(\pi^{\rm
h}_{t}\left(f-\sigma\nabla\log\pi^{\rm
h}_{t}\right)\right)-\nabla\cdot(\pi^{\rm h}_{t}\tilde{g}_{t})$

We can find $\tilde{g}_{t}$ in terms of known quantities as follows: we begin
by taking the derivative of $\pi^{\rm h}_{t}$ with respect to time:

$\partial_{t}\pi^{\rm h}_{t}=-\pi^{\rm
h}_{t}\left(\frac{\dot{Z}_{t}}{Z_{t}}+\frac{L}{T}\right)+\frac{1}{Z_{t}}e^{-\frac{t}{T}L}\partial_{t}\pi_{t}$

Next we substitute (24) for $\partial_{t}\pi_{t}$ and use
$\pi_{t}=Z_{t}e^{\frac{t}{T}L}\pi^{\rm h}_{t}$

$\displaystyle\partial_{t}\pi^{\rm h}_{t}=$ $\displaystyle-\pi^{\rm
h}_{t}\left(\frac{\dot{Z}_{t}}{Z_{t}}+\frac{L}{T}\right)-Z_{t}^{-1}e^{-\frac{t}{T}L}\nabla\cdot\left(\pi_{t}\left(f-\sigma\nabla\log\pi_{t}\right)\right)$
$\displaystyle=$ $\displaystyle-\pi^{\rm
h}_{t}\left(\frac{\dot{Z}_{t}}{Z_{t}}+\frac{L}{T}\right)-Z_{t}^{-1}e^{-\frac{t}{T}L}\nabla\cdot\left(Z_{t}e^{\frac{t}{T}L}\pi^{\rm
h}_{t}\left(f-\sigma\nabla\log Z_{t}e^{\frac{t}{T}L}\pi^{\rm
h}_{t}\right)\right)$ $\displaystyle=$ $\displaystyle-\pi^{\rm
h}_{t}\left(\frac{\dot{Z}_{t}}{Z_{t}}+\frac{L}{T}\right)-e^{-\frac{t}{T}L}\nabla\cdot\left(e^{\frac{t}{T}L}\pi^{\rm
h}_{t}\left(f-\sigma\frac{t}{T}\nabla L-\sigma\nabla\log\pi^{\rm
h}_{t}\right)\right)$ $\displaystyle=$
$\displaystyle-\pi^{h}_{t}\left(\frac{\dot{Z}_{t}}{Z_{t}}+\frac{L}{T}\right)-\nabla\cdot\left(\pi^{\rm
h}_{t}\left(f-\sigma\frac{t}{T}\nabla L-\sigma\nabla\log\pi_{t}^{\rm
h}\right)\right)$ $\displaystyle\qquad-\,\,\frac{t}{T}\pi^{\rm h}_{t}\nabla
L\cdot\left(f-\sigma\frac{t}{T}\nabla L-\sigma\nabla\log\pi_{t}^{\rm
h}\right).$

Comparing with (25) it follows that we require

$\nabla\cdot(\pi^{\rm h}_{t}\tilde{g}_{t})=\frac{1}{T}\pi^{\rm
h}_{t}\left(L+t\nabla L\cdot\left(f-\sigma\frac{t}{T}\nabla
L-\sigma\nabla\log\pi^{\rm h}_{t}\right)\right)+\pi^{\rm
h}_{t}\frac{\dot{Z}_{t}}{Z_{t}}-\nabla\cdot\left(\pi^{\rm
h}_{t}\sigma\frac{t}{T}\nabla L\right)$

For the $\dot{Z}_{t}$ term we have

$\displaystyle\frac{\dot{Z}_{t}}{Z_{t}}=$
$\displaystyle\frac{1}{Z_{t}}\int\partial_{t}e^{-\frac{t}{T}L(x)}\pi_{t}(x)dx=\frac{1}{Z_{t}}\int-\frac{L}{T}e^{-\frac{t}{T}L(x)}\pi_{t}(x)+e^{-\frac{t}{T}L(x)}\partial_{t}\pi_{t}(x)dx$
$\displaystyle=$ $\displaystyle-\frac{1}{T}\mathbb{E}L-\frac{1}{Z_{t}}\int
e^{-\frac{t}{T}L(x)}\nabla\cdot\left(\pi_{t}\left(f-\sigma\nabla\log\pi_{t}\right)\right)dx$
$\displaystyle=$
$\displaystyle-\frac{1}{T}\mathbb{E}L-\frac{1}{Z_{t}}\int\frac{t}{T}e^{-\frac{t}{T}L(x)}\pi_{t}\nabla
L\cdot\left(f-\sigma\nabla\log\pi_{t}\right)dx$ $\displaystyle=$
$\displaystyle-\frac{1}{T}\mathbb{E}L-\frac{t}{T}\mathbb{E}\nabla
L\cdot\left(f-\sigma\nabla\log\pi_{t}\right)$ $\displaystyle=$
$\displaystyle-\frac{1}{T}\mathbb{E}L-\frac{t}{T}\mathbb{E}\nabla
L\cdot\left(f-\sigma\nabla\log\pi^{\rm h}_{t}-\sigma\frac{t}{T}\nabla
L\right)$

where the third equality follows from integration by parts and the expected
value is with taken with respect to $\pi^{\rm h}_{t}$. We finally note that

$\tilde{g}_{t}(x)=g_{t}(x)-\sigma\frac{t}{T}\nabla L(x).$

##  Appendix B Ensemble Kalman filter approximations

We provide details on the derivation of the EnKF-like approximation (14) to
the controlled mean field equation (7) and the various simplifications that
arise from assuming a linear forward map.

We first recall that a continuous time formulation of the EnKF for a generic
likelihood function

$L(x)=\frac{1}{2}(\hat{h}(x)-y_{T})^{\top}\hat{R}^{-1}(\hat{h}(x)-y_{T})$

is provided by

$\frac{\rm d}{{\rm
d}t}X_{t}=-\Sigma_{t}^{x\hat{h}}\hat{R}^{-1}\left(\frac{1}{2}\left(\hat{h}(X_{t})+\pi_{t}[h]\right)-y_{T}\right).$

Here $\pi_{t}$ denotes the law of $X_{t}$, $\pi_{t}[g]$ the expectation value
of a function $g$ under $\pi_{t}$, and

$\Sigma_{t}^{x\hat{h}}=\pi_{t}\left[(x-\pi_{t}[x])(\hat{h}-\pi_{t}[h])^{\top}\right]$

is the covariance matrix between the state $x$ and the forward map $\hat{h}$.

Formal application of this approach to the two contributions to the likelihood
function (13) leads to (14). More precisely, the first term leads to
$\hat{h}=h$ and

$\hat{R}=\frac{t+\Delta t}{2\Delta tT}R$

while the second term results in $\hat{h}=\tilde{h}$ and

$\hat{R}=-\frac{t}{\Delta tT}R.$

The EnKF makes use of statistical linearization

$\Sigma_{t}^{xx}\pi_{t}[\nabla h]=\Sigma_{t}^{xh}$

which holds provided $\pi_{t}$ is Gaussian or if $h$ is linear; a result known
as Stein’s identity. The identity can also be used to approximate derivatives
in a (weakly) non-Gaussian setting giving rise to

(26) $\nabla h(x)\approx(\Sigma_{t}^{xx})^{-1}\Sigma_{t}^{xh}.$

We also recall the robust time-stepping method

(27) $X_{t_{n+1}}-X_{t_{n}}=-\Delta t\Sigma_{t_{n}}^{x\hat{h}}\left(\Delta
t\Sigma_{t_{n}}^{\hat{h}\hat{h}}+R\right)^{-1}\left(\frac{1}{2}\left(\hat{h}(X_{t_{n}})+\pi_{t_{n}}[h]\right)-y_{T}\right)$

which again can be adjusted appropriately to (14).

We now assume a linear forward map, that is $h(x)=Hx$, and discuss the
simplifications that result in the computation of (14). Note that

$\tilde{h}(x)=Hx-\Delta tHf(x)+\Delta t\frac{\sigma t}{T}HH^{\top}R^{-1}Hx.$

Hence the covariance matrix $\Sigma_{t}^{x\tilde{h}}$ can be reformulated to

$\Sigma_{t}^{x\tilde{h}}=\Sigma_{t}^{xx}H^{\top}-\Delta
t\Sigma_{t}^{xf}H^{\top}+\Delta t\frac{\sigma
t}{T}\Sigma_{t}^{xx}H^{\top}R^{-1}HH^{\top}$

and (14) simplifies to

$\displaystyle\hat{g}_{t}^{\rm KF}(x)$
$\displaystyle=-\frac{1}{T}\Sigma_{t}^{xx}H^{\top}R^{-1}\left(\frac{1}{2}\left(Hx+H\mu_{t}^{\rm
h}\right)-y_{T}\right)-\frac{t}{T}\Sigma_{t}^{xf}H^{\top}R^{-1}\left(\frac{1}{2}\left(Hx+H\mu_{t}^{\rm
h}\right)-y_{T}\right)$ $\displaystyle\qquad+\,\,\frac{\sigma
t^{2}}{T^{2}}\Sigma_{t}^{xx}H^{\top}R^{-1}HH^{\top}R^{-1}\left(\frac{1}{2}\left(Hx+H\mu_{t}^{\rm
h}\right)-y_{T}\right)-\frac{t}{2T}\Sigma_{t}^{xx}H^{\top}R^{-1}H\left(f(x)+\pi_{t}^{\rm
h}[f]\right)$ $\displaystyle\qquad+\,\,\frac{\sigma
t^{2}}{T^{2}}\Sigma_{t}^{xx}H^{\top}R^{-1}HH^{\top}R^{-1}\left(\frac{1}{2}\left(Hx+H\mu_{t}^{\rm
h}\right)-y_{T}\right)+\mathcal{O}(\Delta t)$
$\displaystyle=-\left\\{\frac{1}{T}\Sigma_{t}^{xx}+\frac{t}{T}\Sigma_{t}^{xf}-\frac{2\sigma
t^{2}}{T^{2}}\Sigma_{t}^{xx}H^{\top}R^{-1}H\right\\}H^{\top}R^{-1}\left(\frac{1}{2}\left(Hx+H\mu_{t}^{\rm
h}\right)-y_{T}\right)$
$\displaystyle\qquad-\,\,\frac{t}{2T}\Sigma_{t}^{xx}H^{\top}R^{-1}H\left(f(x)+\pi_{t}^{\rm
h}[f]\right)+\mathcal{O}(\Delta t).$

Upon dropping terms of order $\mathcal{O}(\Delta t)$ and using
$\Sigma_{t}^{xx}=\Sigma_{t}^{\rm h}$, we obtain (16) for $f=0$ and (19) for
$\sigma=0$ as special cases. The mean field equations (21) also follow easily
from $f(x)=Fx+b$.


# Quantity vs Quality: Investigating the Trade-Off between
Sample Size and Label Reliability

Timo Bertram  Johannes Fürnkranz  Martin Müller

###### Abstract

In this paper, we study learning in probabilistic domains where the learner
may receive incorrect labels but can improve the reliability of labels by
repeatedly sampling them. In such a setting, one faces the problem of whether
the fixed budget for obtaining training examples should rather be used for
obtaining all different examples or for improving the label quality of a
smaller number of examples by re-sampling their labels. We motivate this
problem in an application to compare the strength of poker hands where the
training signal depends on the hidden community cards, and then study it in
depth in an artificial setting where we insert controlled noise levels into
the MNIST database. Our results show that with increasing levels of noise,
resampling previous examples becomes increasingly more important than
obtaining new examples, as classifier performance deteriorates when the number
of incorrect labels is too high. In addition, we propose two different
validation strategies; switching from lower to higher validations over the
course of training and using chi-square statistics to approximate the
confidence in obtained labels.

Machine Learning, Online Learning, Labeling, MNIST, Oracle

##  1 Introduction

Sample efficiency is one of the main concerns when training neural networks.
Having a better sample efficiency leads to faster training and, when using
online learning, requires less computational effort to generate new samples. A
different but related aspect of deep learning is the reliability of labels.
Noisy or incorrect labels hinder training and often lead to slower convergence
and worse performance (Song et al., 2020). Therefore, label reliability and
the efficiency of training are inherently related (Reed et al., 2015).

While the reliability of labels can not be influenced in some contexts, we
consider a problem setting where we have a fixed budget of samples that can be
used for training a classifier. At each step after receiving the label of a
sample, we can decide whether we want to sample a new example or resample the
label of the previously seen example in order to validate the correctness of
the label. This can, e.g., be used in noisy domains, where the observed label
has a certain error probability, or in probabilistic problems, where each
example is associated with a probability distribution over the labels. For
simplicity, we will in the following typically assume that we have a single
correct label which may be distorted by noise in the data, but the same
setting is also applicable to probabilistic domains. Practical applications of
this include a crowd-sourcing campaign where letting multiple participants
label a sample will increase the likelihood of the label being correct, or
probabilistic domains, where the label depends on a probabilistic event. We
use one such case, poker, as a motivating example, where the task is to
predict the winner of two hands, but the outcome depends on hidden community
cards.

In this work, we do not provide the classifier any information about the
noisiness of the data. In addition, we intentionally do not take any noise-
reducing measures apart from the validation process. For both series of
experiments, we use the neural network architecture that performed well on the
non-noisy variant for the experiments with noise. This noise-agnostic approach
was chosen such that noise-handling in the classifier and training algorithm
is decoupled from the domain we are investigating here; handling noise in the
data generation step.

After a brief discussion of related work (Section 2), we motivate the research
problem with an experiment in the domain of poker, where the question of
whether to re-sample known examples to increase label reliability or to sample
new examples for increasing the size of the training set naturally occurs
(Section 3). We then formalize the problem and our research questions in
Section 4, and empirically investigate it in a set of studies on the MNIST
data where controlled levels of noise were inserted (Section 5). The
conclusions that we draw from these experiments are summarized in Section 6.

##  2 Related work

Learning in noisy and probabilistic domains has been studied extensively in
the machine learning literature. However, such works mostly cover other
solutions than this work aims to, partly due to using different assumptions
about receiving training samples. One common approach when handling noisy
labels is to omit the questionable labels in training and later use a
classifier trained on non-omitted labels as ground truth to identify errors.
Through cross-validation, this then yields adjusted labels that can be used
for the final training of the network (Chen et al., 2019). Other works use
specific layers to handle the noisy data (Sukhbaatar & Fergus, 2015), which
however also requires knowledge about the noise distribution which has to be
learned in a previous step. Additionally, previous research aims to change
loss functions to account for noise in the data (Sukhbaatar & Fergus, 2015;
Zhang & Sabuncu, 2018).

Different related areas of research concern cost-effective and active learning
(Settles, 2012; Wang et al., 2017), where the algorithm can query an oracle to
label specific samples. There, one of the main goals is to train classifiers
as efficiently as possible with respect to required samples. Active learning
also provides an interesting parallel to this work as both use an oracle to
label samples. The main difference between active learning and what is studied
here lies in the fact that in this domain, the oracle is noisy and multiple
independent queries to the oracle can be used to reduce noise. In active
learning, the algorithm aims to identify which samples provide the highest
benefit and queries the oracle for the labels of those. While the combination
of both settings is not studied here, it provides an interesting ground for
further research. Similar to active learning, it may be especially beneficial
to validate more important samples further, while less interesting samples may
require fewer validations to efficiently progress the classifying performance.

Our work is also related to explicit noise filtering algorithms, such as the
ensemble filters introduced in (Brodley & Friedl, 1999). Like these
algorithms, we can use multiple queries to correct noisy labels, but we assume
that we can directly query the labeling oracle instead of using an ensemble of
diverse noise-tolerant learners, and we explicitly aim at studying the trade-
off between obtaining new examples and re-sampling known examples.

##  3 Motivating Example: Poker

Our work is motivated by experiments that we did in learning to evaluate hand
strength in Texas Hold ’em Poker. The task here is to predict which of two
2-card hands is more likely to win the game after all five community cards are
revealed (the ’river’), given the first three community cards (the ’flop’). In
this domain, the label quality vs. quantity problem occurs naturally: Should
we obtain multiple samples of the missing two community cards in order to get
a more reliable estimate of which of the two hands is stronger in the context
of the current flop, or should we instead increase our training set with more
samples, i.e., different combinations of hand and flop cards?

###  3.1 Experimental Setup

We investigate this domain by training neural networks. Thus, a training
sample consists of three inputs ($p1$, $p2$ and $flop$) while the output
signals whether $p1$ or $p2$ is the better hand in the context of the given
flop. To account for such a structure, we train a Siamese network (Koch et
al., 2015; Chicco, 2021). While the complete Siamese structure that was used
can not be explained fully here, Siamese networks can not only be used in one-
shot image learning but also to model which of several inputs is preferential
to another (Tesauro, 1989) in a given context (Bertram et al., 2021). Here,
the two hands and the flops are used as inputs for the Siamese network and its
output models which of the two hands is assumed to be better given the flop.

The test accuracy of the network is computed by taking the average over all
possible rivers that extend the current flop, i.e., if $P(p1\textrm{\
wins}|p1,p2,flop)>P(p2\textrm{\ wins}|p1,p2,flop)$, then the network should
predict that $p1$ is the better hand. For training, the network receives $p1$,
$p2$ and the $flop$ as an input. In this situation, the algorithm has the
option to choose how often it wants to validate the training signal, i.e.
which hand is assumed to be better. One validation consists of drawing two
random and legal community cards to form the river and see which of the two
hands wins the game. By doing this multiple times, a training signal is
constructed which, is less noisy for higher numbers of validations, i.e. as
more different rivers are seen, which increases the chance of the assumed
better hand to be the better hand in general.

To emphasize this, consider the following example. Hand
$p1=(Q\varheartsuit,J\spadesuit)$, hand $p2=(7\spadesuit,7\vardiamondsuit)$,
and $flop=(2\spadesuit,9\spadesuit,10\spadesuit)$. The algorithm then decides
how often it wants to draw two remaining cards to approximate the better hand,
e.g. for one draw the river may result in
$(2\spadesuit,9\spadesuit,10\spadesuit,Q\vardiamondsuit,K\vardiamondsuit)$,
resulting in the approximation that $P(p1\textrm{\ wins}|p1,p2,flop)=1$ and
$P(p2\textrm{\ wins}|p1,p2,flop)=0$, with the network using the information
that hand $p1$ is better than hand $p2$ given the flop as the training signal.
For two additional validations, those rivers may result in
$(2\spadesuit,9\spadesuit,10\spadesuit,A\vardiamondsuit,3\vardiamondsuit)$ and
$(2\spadesuit,9\spadesuit,10\spadesuit,10\vardiamondsuit,10\varheartsuit)$,
approximating the winning chances for p1 and p2 as 0.33 and 0.67 respectively.
In the limit of validations, this approximation will result in the true win-
probabilities (here $P(p1\textrm{\ wins}|p1,p2,flop)=0.6697$ and
$P(p2\textrm{\ wins}|p1,p2,flop)=0.3303)$.

###  3.2 Result

Figure 1: Correct comparisons between two poker hands on the test dataset over
time. Labels are sorted by highest peak-performance achieved.

Figure 1 shows the test accuracy of networks trained with different numbers of
validations over the course of training. In this experiment, labels are sorted
by their peak performance reached. We can observe that smaller validation
numbers increase the performance of the learning faster in the early training
but generally achieve lower peak performances. It seems that an appropriate
choice of validation numbers crucially decides the performance of the
learning: while too much re-sampling leads to greatly ineffective training,
training without re-sampling and relying on a larger sample size also
significantly hurts the peak performance.

Following this motivation to study the problem, we now formalize the problem
we are researching here.

##  4 Problem definition and terminology

The problem we want to handle can be defined as follows: The learner is
trained with examples $\mathcal{X}=\\{\mathbf{x}_{i}\\},i=1\dots n$, which are
associated with labels $\mathcal{L}=\\{\lambda_{j}\\},j=1\dots l$ according to
an unknown function $f:\mathcal{X}\rightarrow\mathcal{L}$, with
$y_{i}=f(\mathbf{x}_{i})$. We further assume that the labeling process is non-
deterministically guided by a noisy oracle $\mathcal{O}$ with a probability
vector $\mathbf{p}=(p_{1},\dots,p_{l})$, which assigns the example the correct
label $\lambda_{j}$ with a probability $p_{j}$, and one of the other labels
with a probability of $p_{k}$, where $p_{j}>p_{k}$ for all $k\neq j$. Thus,
$y_{i}=\lambda_{j}$ with $j=\arg\max_{k}p_{k}$.

We will often abbreviate this notation by using $q=p_{j}$ for the probability
of picking the correct label, and $w=1-q=\sum_{k\neq j}p_{k}$ for picking one
of the incorrect labels.

In this setting, we are given a budget $s_{\max}$ of queries to the oracle
$\mathcal{O}$, which we can use for creating a training set $\mathcal{T}$. As
explained above, each query to $\mathcal{O}$ returns the correct label for a
given example with probability $q$, and an incorrect label with probability
$w$. However, we can query the oracle $v$ times for the same sample, each time
receiving a newly generated label which may or may not be correct according to
$\mathbf{p}$. Such repeated queries for the same example are called
_validations_ , as we are trying to validate whether the label is correct.
Each validation counts as a new query to the budget $s_{\max}$ in the same way
as queries for new examples. This means that we can use the budget $s_{\max}$
to create a training set consisting of, e.g., $|\mathcal{T}|=s_{\max}$
different samples, each validated only once ($v=1$), or we can compute
$|\mathcal{T}|=\frac{s_{\max}}{3}$ samples where each label has been validated
three times. In general, the number of validations can be changed at any time.
This means that one may choose it prior to training, change it over the course
of training, or even change it depending on the current sample or its already
received labels.

Validating labels reduces the sizes of the training data that the classifier
will be able to use throughout training but increases the likelihood of
samples receiving the correct label. After $v$ validations for the same
example, in which each label has been observed with frequencies
$\mathbf{v}=(v_{1},\dots,v_{l})$, this example is assigned a label
$\lambda_{v}$ with a majority vote, i.e.,

$\lambda_{v}=\lambda_{k},\textrm{ with }k=\arg\max_{j}v_{j}.$ (1)

Ignoring ties between labels, the probability of $\lambda_{v}$ being the
correct label $\lambda_{j}$ after $v$ queries is equal to the probability of
this label being drawn more frequently than any of the other labels. The
probabilities of observing the labels with frequencies $\mathbf{v}$ after $v$
validations are multinominal distributed with the probability vector
$\mathbf{p}$, i.e.,

$P(v_{1},\dots,v_{l})=\textrm{Mult}(\mathbf{v},\mathbf{p}).$ (2)

Thus, a lower bound for the probability that the correct label $\lambda_{j}$
is selected is the probability that $v_{j}>v_{k},\forall k\neq j$:

$P(\lambda_{j}=\lambda_{v})=\sum_{\mathbf{v}:\sum_{i}v_{i}=v}\textrm{Mult}(\mathbf{v},\mathbf{p})\prod_{k\neq
j}{\mathbbm{1}({v_{j}>v_{k}})}$ (3)

where $\mathbbm{1}(.)$ is the indicator function.

Figure 2: Probability of receiving the correct label given the number of
validations for different probabilities of receiving incorrect labels. This
simulation assumes 10 classes where all incorrect classes have equal chances.

Obviously, $P(\lambda_{j}=\lambda_{v})$ depends on the number of validations
$v$, the number of different labels $l$ and the chance of $\mathcal{O}$
returning the correct label $q$. For $\lim_{q\to 1}$, $\lim_{v\to\infty}$ and
$\lim_{l\to\infty}$, the probability of receiving the correct label converges
to 1 and the function has a logarithmic shape. Figure 2 shows the probability
of receiving a correct label from the oracle for varying numbers of
validations for a dataset consisting of 10 different classes (e.g., like
MNIST). Each curve corresponds to a noise level $w$. There, we also assume
that all incorrect classes follow a uniform distribution, i.e.
$p_{k}=\frac{w}{l-1}$ for all $k\neq j$. The probabilities to receive the
correct labels in domains with unequal noise differ. From this, it is obvious
that the choice of validations crucially depends on the noise in the data, as
for lower noise few validations already result in almost perfect accuracy,
while higher noise data benefits longer from additional validations to
increase the reliability of the received labels.

#### Research Questions.

On the basis of the above formalization, we define a number of research
questions concerning the label validation that should be answered in this
work. Initially, it is not clear that any validation of labels is an efficient
strategy, as re-evaluating the label of a known example is equally as costly
as obtaining a label for a new example, which raises the question of whether
such validations, which reduce the training set size but increase the label
reliability, have any benefit at all. Intuitively, one might think that
creating a sample only once will yield the highest efficiency of training, as
the algorithm will use the maximum number of different examples for training
and not ”waste” any examples for validating samples. However, we show in later
experiments (Sections 5.2 and 5.3) that this is not the case and that
validating samples multiple times leads to faster and more efficient training
of a neural network in our experiments. In cases with very high levels of
noise, not validating labels even stopped the network from learning at all.

In summary, the research questions investigated in the remainder of the paper
are the following:

  * R1:

Can it be better to validate unreliable labels by pooling them with a majority
vote or should all examples be used for training?

  * R2:

How should the number of validations be chosen in order to most efficiently
train a neural network?

  * R3:

What is the relationship between the number of validations and the level of
noise in the data?

  * R4:

Can an adaptive number of validations result in better efficiency than a
stable number?

In the following, we will address these questions empirically in a study where
we add artificial noise to the MNIST dataset so that we can systematically
investigate the trade-off between re-sampling for increased label reliability
and additional sampling for increased data volume.

##  5 Experiments on MNIST

The domain in which we investigate the previously phrased questions is the
popular image-dataset MNIST (Deng, 2012). We first describe the experimental
setup and afterward the results that concern the research questions formulated
above.

###  5.1 Experimental Setup

In this section, the setup for the conducted experiments is outlined. The
MNIST database consists of 60,000 gray-scale training images of size 28x28,
for which we simulate a labeling process with a noisy oracle. The 10,000 test
images are not obstructed with noise and are used to compute the test accuracy
of the networks.

#### Labeling.

Each image falls into one of ten classes which represent the ten digits. For
this dataset, a noise level of $w$ is modeled by an oracle with returns the
correct label with probability $q=1-w$, while all other labels have a chance
of $\nicefrac{{w}}{{9}}$ of being returned. Note that in this setting, due to
the comparably small training set, multiple epochs of the whole training set
are required. Therefore, in contrast to the online learning setting described
in Section 3, the network will see the same training examples multiple times.
Due to the noisy setting, the network may also receive different labels for
the same examples in each epoch. However, in order not to disturb the
validation parameters, we do not save labels from previous epochs. Validations
are independent each time a sample is generated, even though it may have been
generated at a previous epoch.

#### Training procedure.

To provide a baseline performance, as well as providing a structure for the
experiments with noise, a neural network is constructed by optimizing its
architecture on the basic MNIST-dataset without noise. In the end, the
architecture of the network consists of 6 convolutional layers organized in 3
blocks. For the full architecture, refer to Appendix D. With this, a final
testing accuracy of 99.4% is achieved, which is deemed high enough compared to
the current state-of-the-art performance of 99.9% (An et al., 2020). The
specific choices of hyperparameters, while of no great importance to the
questions and results of this work, are outlined in Table 1.

Table 1: Chosen hyperparameters for the MNIST network

Hyperparameter | Choice
---|---
Learning rate | 0.01
Optimizer | AdamW
Loss function | Cross Entropy Loss
Dropout | 0.25
Batch normalization | Yes
Number of convolutional layers | 6
Activation function | ELU
Batch size | 256

In all experiments with noisy MNIST, the main performance metric is the test
accuracy of the network in relation to the number of seen training samples.
Note that, as outlined in Section 4, multiple validations of the same sample
count the same as seeing a new sample. For example, if $v=3$ validations are
performed on each example, the algorithm needs $3\times 256=768$ examples in
order to produce a single batch update, whereas when no validations are
performed ($v=1$), three batch updates can be performed from the same number
of examples. In our plots, all versions will consume the same amount of
examples (shown on the $x$-axis). It is important that this is a conservative
assumption and overstates the computational effort of validation. In reality,
validating a label is much faster than generating a new sample, passing it
through the network, computing its loss, and back-propagating it. Therefore,
the advantages of validation are even higher than the results show when the
true wall-clock time required for training is used as a metric.

With the assumption that validating samples with a majority vote is the same
as viewing a new sample, some considerations have to be made to allow for a
fair comparison between different validation numbers. As all performances here
are related to the number of seen samples, and as the batch size was held
consisted in the experiments, a network trained with multiple validations will
have made fewer gradient descent steps in the same ”time” as one with one
validation by a factor of $v$. Therefore, we multiply the base learning rate
of 0.01 by the number of validations that are used for the respective training
run. This again can be regarded as a disadvantage for the networks with
validations, as a higher learning rate will make it harder to efficiently
descend in the loss landscape. However, this adaptation was made to normalize
the potential decrease in loss per number of seen samples.

###  5.2 Utility of validation

We first investigate whether validations have any advantage over directly
training on all of the sampled data. To this end, we compare two learners that
see exactly the same data, which consists of examples each of which have been
sampled $v$ times. The first version of the network uses all of the samples
directly for training the network, whereas the second version first processes
the samples in an internal validation process. Thus, the first network is
trained on $N$ examples, whereas the second network is trained on
$\nicefrac{{N}}{{v}}$ examples, which have been derived from the original $N$
examples via majority voting. In the plots, however, the performances of both
versions are aligned on the $N$ examples that they both have received.

Figure 3: Comparison of processing samples with a majority vote (solid) and
using all samples individually (dashed). In all cases, voting leads to faster
training with better or similar peak performance.

Figure 3 shows this comparison of using all samples against first validating
them with a majority vote for different noise values and different validation
rates $v$. There, we see that in all settings, voting for the final label
leads to better results than directly using all examples for training. This
result makes intuitive sense, as the learner receives numerous incorrect
labels in the former case and has to converge to a good classification
performance by learning to ignore them. As such, it learns slower at the
beginning of training. However, one may also make the argument that validating
a single sample with a majority vote leads to a loss of information, as the
classifier will only see the most often obtained label without any information
about how often others were returned in queries. However, this experiment
indicates that the former argument weighs stronger than the second one in
practice.

One may argue that the learner would perform better if it could use $N$
different examples instead of repeated samples from validating. This is the
setting that will be investigated in the next section, where we also try to
investigate the relationship between validation amount and noise.

###  5.3 Fixed validation choices

First, the performances of different constant choices for the validation
parameter $v$ are compared for different noise levels in the training data. As
learning curves for closely related performances can be hard to read, we use a
different metric to evaluate how well a specific choice performed. Before
training, several levels of test accuracy are chosen and performances are
compared by the speed with which those levels were achieved. Levels ranged
from very easy to reach (80% test accuracy) to levels similar to the
performance on non-noisy MNIST (99.3% test accuracy). This metric of
comparison also makes it possible to disregard overfitting as only the first
time of reaching a level was important.

Table 2: Seen samples required to first reach the specified test accuracies
(lower is better). Only the five best fixed values of $v$ are shown for
conciseness. Full table in Appendix A

| 80% | 90 % | 95% | 98% | 99% | 99.2% | 99.3%
---|---|---|---|---|---|---|---
20% Noise |
Fixed ($v=5$) | 250k | 550k | 950k | 1651k | 8150k | 20750k | 62950k
Fixed ($v=3$) | 50k | 150k | 350k | 1050k | 5450k | 13150k | 63750k
Fixed ($v=11$) | 652k | 1151k | 2051k | 3851k | 21051k | 28751k | NaN
Fixed ($v=7)$ | 451k | 851k | 1351k | 3250k | 12951k | 40951k | NaN
Fixed ($v=15)$ | 1053k | 1552k | 2651k | 5351k | 22352k | 48252k | NaN
40% Noise |
Fixed ($v=5$) | 350k | 650k | 950k | 2350k | 15350k | 28250k | 52450k
Fixed ($v=11$) | 551k | 851k | 1450k | 3250k | 13452k | 26851k | 60850k
Fixed ($v=7$) | 351k | 550k | 951k | 2350k | 12250k | 25250k | 74650k
Fixed ($v=25$) | 2953k | 3853k | 4352k | 6553k | 29152k | 48752k | NaN
Fixed ($v=15)$ | 2953k | 3853k | 4352k | 6952k | 19552k | 60251k | NaN
60% Noise |
Fixed ($v=11$) | 352k | 851k | 1351k | 3151k | 11750k | 29451k | 45751k
Fixed ($v=15$) | 951k | 1353k | 2153k | 4053k | 28153k | 40953k | NaN
Fixed ($v=7$) | 1450k | 1750k | 2150k | 4050k | 22450k | NaN | NaN
Fixed ($v=25$) | 1551k | 2754k | 4152k | 7551k | 27752k | NaN | NaN
Fixed ($v=5$) | 1750k | 1850k | 2251k | 7550k | 45950k | NaN | NaN
80% Noise |
Fixed ($v=15$) | 1552k | 2552k | 3953k | 8353k | 42952k | NaN | NaN
Fixed ($v=25$) | 2053k | 3553k | 4551k | 8754k | 45352k | NaN | NaN
Fixed ($v=11$) | 1351k | 2351k | 3952k | 13651k | NaN | NaN | NaN
Fixed ($v=51$) | 1958k | 3458k | 5553k | 13858k | NaN | NaN | NaN
Fixed ($v=99$) | 3852k | 6852k | 11160k | 28752k | NaN | NaN | NaN

Table 2 shows an overview of this performance for four noise levels, 20%, 40%,
60% and 80%. Note that all test accuracies are computed on the full MNIST test
dataset without any noise applied to it. The choices for the constant
validation parameter were pre-defined as 1, 3, 5, 7, 11, 15, 25, 51, and 99
validations per sample.

For this comparison, the five best-performing choices of $v$ are outlined.
While some parameters performed obviously better than others, it is possible
to interpret what makes one parameter better than another in several different
ways, e.g. averaging over all obtained values, speed by which a specific level
was reached, highest peak performance obtained, and others. We decided to
order the performances by first using the speed by which the highest level,
99.3% test accuracy, was reached. Then, all choices which did not reach this
were ordered by the speed by which they achieved the next-lower level and
continuing this process until a complete ordering is obtained. As such, we
value peak performance the most, followed by the speed with which that
performance was achieved. For the full overview of all tested parameters,
refer to the Appendix A.

In the 20% noise setting, 5 and 3 validations performed best. Comparing the
two, we see that $v=3$ was able to achieve the early goals more efficiently
than $v=5$. However, the difference between both decreases throughout the
different checkpoints until the 5-validation setting reaches the final
performance slightly faster. We also see that the setting without validations,
i.e. $v=1$, performed comparably poorly and was not able to reach the top 5.
This confirms that investing some examples to obtain an increase in label
quality at the expense of data volume provided a clear benefit here. In
addition, we also see that the high-validation settings, i.e. 25, 51, and 99,
were too inefficient for this environment. This observation can be explained
by Figure 2, which shows how quickly label reliability increases with more
validations in the 20% noise setting. Therefore, querying the oracle for too
many additional validations in a comparably low-noise environment wastes
samples that should be used for new samples instead.

In the experiments with 40% noise, using 5 validations per sample still
performed best, but the order of efficiency for other choices differs. By
increasing the noise of the oracle, $v=3$, the previous second-best choice,
did not reach the top 5 and is replaced with 25. This observation underlines
the hypothesis that with increasing noise of the oracle, the effectiveness of
a higher number of validations for examples increases. We again see that in
the earlier stages of training, lower validation numbers obtain vastly higher
efficiencies. This advantage however is lost when training gets more
difficult.

This trend continues when regarding the oracles with 60% and 80% noise. With
60%, the previous top performer 5 falls drastically in performance. This is a
result of the fact that the highest levels, 99.2%, and 99.3% test accuracy
could not be achieved with lower than 11 validations per sample. In the
highest noise setting, all of the five highest choices for $v$ performed best,
again strengthening the previous observation that appropriate validation
choice increases as noise increases. In addition, we find that in this
extremely noisy environment, the low validation choices 1 and 3 were not able
to learn at all and had the same test accuracy as classifying randomly. This
shows that choosing an appropriate number of validation queries can make the
difference between reaching test accuracies higher than 99% and completely
random classification of the learner.

#### Observations.

This series of experiments provides intriguing insights into the influence of
the validation choice on the classifying performance. We observe that there
seems to be an approximate range of appropriate parameter choice which differs
depending on the noisiness of the oracle. In addition, we observe that these
good choices seem to follow a normal distribution with decreasing performance
when choosing too few or too many validations.

Secondly, we observe that too low validations do not reach the same peak
performances but generally reach earlier performance levels more efficiently.
Thus, we speculate that beginning training with lower validations and
increasing them as classifying accuracy increases may be beneficial. We
investigate this further in Section 5.4.

Lastly, we want to emphasize the correlation between noise and validations. In
contrast to the artificial environment studied here, we generally can not
assume that the noisiness of the data is known in advance. Thus, approaches
that are able to approximate noise and adapt validations based on that may
perform better than fixed choices.

This is also the key difference between the experiments on Poker and MNIST.
With MNIST, an even noise was applied to each training sample. In the Poker
domain, noise naturally varies between training samples. While some hands have
a very high chance of winning upward of 90%, which results in many validations
showing their win, other hands have equal winning chances. As such, the
effectiveness of validations is different for different hands, which may make
an appropriate choice much harder, again resulting in a need for adaptive
solutions.

###  5.4 Scheduled validation policies

Based on the results shown in Table 2, we hypothesized that the importance of
label reliability increases with improving network performance and training
duration. In the early stages of training, the classifiers trained on the
nosier samples increase in performance faster but are not able to reach higher
performances in the later stages. Therefore, we speculate that starting
training with a higher amount of lower-validated labels and switching to fewer
but more reliable samples later can outperform the previously shown fixed
approaches.

To provide a first investigation of this, we define several experiments in
which the choice of label validation was switched at pre-defined points in
training. These points do not depend on the returned label variance or network
accuracy, as both of those measures can be highly unreliable, but we rather
use fixed numbers of samples after which the choice of validations is changed
to obtain more reliable labels. We call such variants, which have pre-defined
switching points, ”scheduled” policies.

For such scheduled approaches, we defined a starting number of validations to
perform per sample, a range over which the validation number increases over
time, and a maximum number of validations to perform. After the first 10% of
training on the initial parameter, the policy switches to the next-higher
choice. There, again 10% of training time is spent and the process repeats
until it reaches the pre-defined maximum number where the remaining training
time is used. As such, a classifier that is labeled ”v = (1,3,5,7)” or short
”v=($1\dots 7$)” will have used labels that were validated with one, three,
and five queries for the first, second, and third 10% of training respectively
while using labels with seven validations for the remaining 70% of training
time. We compare such scheduled approaches, and the statistical approach
outlined in the next section, with the previously shown fixed variants in
Table 3.

Table 3: Seen samples required to first reach the specified test accuracies
(lower is better). Only the five best overall results are shown. The schedules
progress in the specified subranges of $v=(1,3,5,7,11,15,25,51,99)$. Full
table in Appendix A

| 80% | 90 % | 95% | 98% | 99% | 99.2% | 99.3%
---|---|---|---|---|---|---|---
20% Noise |
Chi ($p<0.05$) | 150k | 250k | 450k | 950k | 3450k | 12550k | 14950k
Scheduled $(1\dots 5)$ | 50k | 150k | 250k | 850k | 5850k | 17350k | 27250k
Scheduled ($1\dots 3$) | 50k | 150k | 250k | 850k | 5550k | 12250k | 41950k
Fixed ($v=5$) | 250k | 550k | 950k | 1651k | 8150k | 20750k | 62950k
Fixed ($v=3$) | 50k | 150k | 350k | 1050k | 5450k | 13150k | 63750k
40% Noise |
Chi ($p<0.05$) | 150k | 250k | 450k | 1251k | 6250k | 18051k | 34950k
Scheduled ($1\dots 11$) | 50k | 650k | 850k | 2550k | 20450k | 26450k | 38350k
Scheduled ($1\dots 15$) | 50k | 250k | 350k | 2350k | 20350k | 27850k | 40951k
Fixed ($v=5$) | 350k | 650k | 950k | 2350k | 15350k | 28250k | 52450k
Scheduled ($3\dots 11$) | 50k | 450k | 650k | 2950k | 10650k | 22350k | 52450k
60% Noise |
Fixed ($v=11$) | 352k | 851k | 1351k | 3151k | 11750k | 29451k | 45751k
Scheduled ($1\dots 25$) | 50k | 1650k | 2150k | 11250k | 20050k | 21551k | NaN
Scheduled ($1\dots 15$) | 50k | 650k | 1250k | 10350k | 20252k | 39750k | NaN
Fixed ($v=15$) | 951k | 1353k | 2153k | 4053k | 28153k | 40953k | NaN
Scheduled ($1\dots 7$) | 50k | 1350k | 1550k | 10250k | 21951k | Nan | NaN
80% Noise |
Scheduled ($11\dots 51$) | 1650k | 2351k | 4051k | 10051k | 22554k | NaN | NaN
Scheduled ($7\dots 25$) | 3250k | 8751k | 10051k | 12351k | 29052k | NaN | NaN
Scheduled ($5\dots 25$) | 7451k | 10050k | 10150k | 15052k | 29353k | NaN | NaN
Scheduled ($3\dots 25$) | 11250k | 13150k | 15251k | 19850k | 34654k | NaN | NaN
Fixed ($v=15$) | 1552k | 2552k | 3953k | 8353k | 42952k | NaN | NaN

These results seem to underline the hypothesis that starting with noisier
samples and saving more validations for later improves the efficiency of
training. In all noise categories, apart from a single outlier in the 60%
noise setting, using the scheduled switches improved the efficiency with which
the accuracy levels were reached. Additionally, in almost all cases, the
scheduled variant ending on $X$ validations performed better than the one with
$X$ fixed validations. We also observe that in general, choosing higher
parameters for the last validation value performed better than lower ones, but
this may be an artifact due to the specific choices made. Choosing too extreme
values would likely break this observation.

###  5.5 Dynamic validation policies

As the last variant, we investigate how we can use the already obtained
validation results at a given point to decide whether additional validations
are required. The goal here is to query for additional labels until enough
confidence in the obtained samples is achieved, such that further validations
are not required. As is assumed in Section 4, $\forall_{k\neq j}:p_{j}>p_{k}$.
Therefore, one policy is to stop querying the oracle for labels if the
obtained labels v give sufficient reason to believe that
$\arg\max_{i}v_{i}=j$, which also implies that the distribution from which v
was samples has the highest probability at $p_{j}$.

As this problem is not trivial to solve, we approximate it with the problem of
deciding whether v’s distribution is sufficiently different from a uniform
distribution, which can be achieved with the chi-square goodness of fit test.
Thus, we pre-define a confidence threshold and then query the oracle for
labels until the p-value of the chi-square goodness of fit test on the
returned labels falls below the threshold, thereby modeling that we are
sufficiently certain that the underlying distribution of query results is not
uniform. See Algorithm 1 for a pseudo-code of this process.

Algorithm 1 Chi-square validation

1: $s\leftarrow\text{sample to validate}$

2: $c\leftarrow\text{confidence threshold}$

3: $l\leftarrow\text{number of classes}$

4: $counts\leftarrow\text{array with zeros of length l}$

5: $p\leftarrow 1$

6: while $p>\text{confidence threshold}$ do

7: $\lambda\leftarrow\mathcal{O}(s)$

8: $\text{counts}[\lambda]\text{+=}1$

9: $p\leftarrow\text{chi square(counts)}$

10: end while

11: return argmax(counts)

This algorithm has two key disadvantages; sample sizes can be extremely small
and the algorithm will conclude querying if the label distribution is assumed
to be non-uniform, but does not have a clear peak. As appropriate validation
numbers in all experiments ranged between 3 and 25, the former problem can not
be solved with any statistical test, as they will generally not work reliably
with such small sample sizes. However, we do not require very high accuracies
of when to stop querying, and the network can learn while receiving numerous
incorrect labels. The task of this algorithm is not to be as accurate as
possible, but rather to provide an approximation of how many validations are
required for the given setting. To illustrate the second problem, consider the
following two counts of labels received by validating a sample of the MNIST-
dataset ten times: [0,0,7,0,1,0,2,0,0,0] and [0,0,5,0,0,0,0,0,5,0].

Both of those provide very high confidences, i.e., low p-values of 0.000001411
and 0.000007599 respectively, in rejecting a uniform distribution. However,
while we can reasonably predict label 2 for the first case, it is not clear
which label should be used for the latter one. While making incorrect
decisions in such cases will likely influence the performance of the learning,
much more computationally intensive practices would be required to avoid this
problem. Table 3 includes the performances of this approach with a confidence
threshold of 0.05. In the settings with lower noise, this worked extremely
well, which makes intuitive sense. In the 20% noise setting, labels are very
often correct. If the oracle returns two of the same labels at the start, it
is already unlikely that it was not the true label of that sample
($9\cdot(\frac{0.2}{9})^{2}\approx 0.4\%$) and only two validations are
required to identify the likely true label. Additional queries are therefore
only needed if at least one incorrect label was retrieved among the first two.
We find that, while not perfect, mean validations (2.99, 4.93, 10.59, and
58.30 for the four levels of noise) are reasonable validation choices for each
setting. However, the chi-square method was only able to achieve the best
performances in the two lower-noise experiments, which leads us to speculate
that the spread of obtained validations per sample was too large in the other
two, as standard deviations were much higher there (9.2 and 64.36 compared to
1.55 and 3.49). A deeper investigation of this problem and other potential
solutions is required. ´

##  6 Conclusion

In this paper, we have investigated the trade-off between label reliability
and sample size. Our results have shown that in noisy or probabilistic
domains, it is beneficial to re-sample the labels of known examples in order
to increase their reliability, rather than obtaining new examples. Conversely,
we have also seen that this advantage is lost if too many replicas of the same
examples are re-sampled and that the optimal validation rate increases with
increasing noise levels. Moreover, we have demonstrated the effectiveness of a
scheduled switching of the validation policy, which increases the number of
validations over time. This variant combines the advantages of low validation
rates, which converge much faster in the initial training phases, and high
validation rates, which are most effective for the final tuning of the
network. Finally, we have also obtained first encouraging results on a
dynamic, sample-dependent validation policy, which worked especially well in
low and medium noise settings, but more work in this area is required. A more
robust method of computing the reliability of obtained validations could
increase the performance of such a dynamic variant further and reduce
currently seen problems when handling very high-noise settings.

## References

  * An et al. (2020) An, S., Lee, M. J., Park, S., Yang, H., and So, J.  An ensemble of simple convolutional neural network models for MNIST digit recognition.  _arXiv_ , 2008.10400, 2020.
  * Bertram et al. (2021) Bertram, T., Fürnkranz, J., and Müller, M.  A comparison of contextual and non-contextual preference ranking for set addition problems.  In _Proceedings of ICML SubsetML workshop_ , 2021.  arXiv:2107.04438.
  * Brodley & Friedl (1999) Brodley, C. E. and Friedl, M. A.  Identifying mislabeled training data.  _Journal of Artificial Intelligence Research_ , 11:131–167, 1999.  doi: 10.1613/jair.606.
  * Chen et al. (2019) Chen, P., Liao, B., Chen, G., and Zhang, S.  Understanding and utilizing deep neural networks trained with noisy labels.  In Chaudhuri, K. and Salakhutdinov, R. (eds.), _Proceedings of the 36th International Conference on Machine Learning (ICML)_ , pp. 1062–1070, Long Beach, California, USA, 2019. PMLR.
  * Chicco (2021) Chicco, D.  Siamese neural networks: An overview.  In Cartwright, H. M. (ed.), _Artificial Neural Networks_ , pp. 73–94. Springer, 3rd edition, 2021.
  * Deng (2012) Deng, L.  The MNIST database of handwritten digit images for machine learning research.  _IEEE Signal Processing Magazine_ , 29(6):141–142, 2012.
  * Koch et al. (2015) Koch, G., Zemel, R., and Salakhutdinov, R.  Siamese neural networks for one-shot image recognition.  In _Proceedings of the ICML Deep Learning Workshop_. Lille, 2015\.
  * Reed et al. (2015) Reed, S. E., Lee, H., Anguelov, D., Szegedy, C., Erhan, D., and Rabinovich, A.  Training deep neural networks on noisy labels with bootstrapping.  In Bengio, Y. and LeCun, Y. (eds.), _Workshop Track Proceedings of the 3rd International Conference on Learning Representations (ICLR)_ , San Diego, CA, USA, 2015.  arXiv:1412.6596.
  * Settles (2012) Settles, B.  _Active Learning_.  Synthesis Lectures on Artificial Intelligence and Machine Learning. Morgan & Claypool Publishers, 2012.  doi: 10.2200/S00429ED1V01Y201207AIM018.
  * Song et al. (2020) Song, H., Kim, M., Park, D., and Lee, J.  Learning from noisy labels with deep neural networks: A survey.  _arXiv_ , 2007.08199, 2020.
  * Sukhbaatar & Fergus (2015) Sukhbaatar, S. and Fergus, R.  Learning from noisy labels with deep neural networks.  In Bengio, Y. and LeCun, Y. (eds.), _Workshop Track Proceedings of the 3rd International Conference on Learning Representations (ICLR)_ , San Diego, CA, USA, 2015.  arXiv:1406.2080.
  * Tesauro (1989) Tesauro, G.  Connectionist learning of expert preferences by comparison training.  In Touretzky, D. S. (ed.), _Advances in Neural Information Processing Systems 1_ , pp. 99–106, Denver, Colorado, 1989. Morgan Kaufmann.
  * Wang et al. (2017) Wang, K., Zhang, D., Li, Y., Zhang, R., and Lin, L.  Cost-effective active learning for deep image classification.  _IEEE Transactions on Circuits and Systems for Video Technology_ , 27(12):2591–2600, 2017.
  * Zhang & Sabuncu (2018) Zhang, Z. and Sabuncu, M. R.  Generalized cross entropy loss for training deep neural networks with noisy labels.  In Bengio, S., Wallach, H. M., Larochelle, H., Grauman, K., Cesa-Bianchi, N., and Garnett, R. (eds.), _Advances in Neural Information Processing Systems 31_ , pp. 8792–8802, Montréal, Canada, 2018\.

##  Appendix A Full MNIST results

Table 4: Seen samples required to first reach the specified test accuracies
(lower is better).

| 80% | 90 % | 95% | 98% | 99% | 99.2% | 99.3%
---|---|---|---|---|---|---|---
20% Noise |
Chi ($p<0.05$) | 150k | 250k | 450k | 950k | 3450k | 12550k | 14950k
Scheduled ($v=(1,3,5)$) | 50k | 150k | 250k | 850k | 5850k | 17350k | 27250k
Scheduled ($v=(1,3)$) | 50k | 150k | 250k | 850k | 5550k | 12250k | 41950k
Fixed ($v=5$) | 250k | 550k | 950k | 1651k | 8150k | 20750k | 62950k
Fixed ($v=3$) | 50k | 150k | 350k | 1050k | 5450k | 13150k | 63750k
Scheduled ($v=(1,3,5,7,11)$) | 50k | 150k | 250k | 1250K | 5850k | 24351k | NaN
Fixed ($v=11$) | 652k | 1151k | 2051k | 3851k | 21051k | 28751k | NaN
Fixed ($v=7)$ | 451k | 851k | 1351k | 3250k | 12951k | 40951k | NaN
Fixed ($v=15)$ | 1053k | 1552k | 2651k | 5351k | 22352k | 48252k | NaN
Fixed ($v=25)$ | 1450k | 2453k | 3652k | 7353k | 27953k | 71053k | NaN
Fixed ($v=1)$ | 50k | 250k | 350k | 1050k | 28350k | NaN | NaN
Fixed ($v=51)$ | 2258k | 3360k | 4953k | 10759k | 41855k | NaN | NaN
Fixed ($v=99)$ | 2255k | 5170k | 8751k | 20759k | 98765k | NaN | NaN
40% Noise |
Chi ($p<0.05$) | 150k | 250k | 450k | 1251k | 6250k | 18051k | 34950k
Scheduled ($v=(1,3,5,7,11)$) | 50k | 650k | 850k | 2550k | 20450k | 26450k | 38350k
Scheduled ($v=(1,3,5,7,11,15)$) | 50k | 250k | 350k | 2350k | 20350k | 27850k | 40951k
Fixed ($v=5$) | 350k | 650k | 950k | 2350k | 15350k | 28250k | 52450k
Scheduled ($v=(3,5,7,11)$) | 50k | 450k | 650k | 2950k | 10650k | 22350k | 52450k
Fixed ($v=11$) | 551k | 851k | 1450k | 3250k | 13452k | 26851k | 60850k
Fixed ($v=7$) | 351k | 550k | 951k | 2350k | 12250k | 25250k | 74650k
Fixed ($v=25$) | 2953k | 3853k | 4352k | 6553k | 29152k | 48752k | NaN
Fixed ($v=15)$ | 2953k | 3853k | 4352k | 6952k | 19552k | 60251k | NaN
Fixed ($v=3)$ | 50k | 350k | 650k | 2650k | 3055k | NaN | NaN
Fixed ($v=51)$ | 1958k | 3458k | 4951k | 13454k | 48556k | NaN | NaN
Fixed ($v=1)$ | 50k | 250k | 350k | 2450k | NaN | NaN | NaN
Fixed ($v=99)$ | 3472k | 8170k | 13273k | 29852k | NaN | NaN | NaN

| 80% | 90 % | 95% | 98% | 99% | 99.2% | 99.3%
---|---|---|---|---|---|---|---
60% Noise |
Fixed ($v=11$) | 352k | 851k | 1351k | 3151k | 11750k | 29451k | 45751k
Scheduled ($v=(1,3,5,7,11,15,25)$) | 50k | 1650k | 2150k | 11250k | 20050k | 21551k | NaN
Scheduled ($v=(1,3,5,7,11,15)$) | 50k | 650k | 1250k | 10350k | 20252k | 39750k | NaN
Fixed ($v=15$) | 951k | 1353k | 2153k | 4053k | 28153k | 40953k | NaN
Scheduled ($v=(1,3,5,7)$) | 50k | 1350k | 1550k | 10250k | 21951k | Nan | NaN
Fixed ($v=7$) | 1450k | 1750k | 2150k | 4050k | 22450k | NaN | NaN
Fixed ($v=25$) | 1551k | 2754k | 4152k | 7551k | 27752k | NaN | NaN
Chi ($p<0.05$) | 1751k | 1951k | 2451k | 4551k | 37050k | NaN | NaN
Fixed ($v=5$) | 1750k | 1850k | 2251k | 7550k | 45950k | NaN | NaN
Scheduled ($v=(1,3,5)$) | 50k | 1950k | 2650k | 10350k | NaN | NaN | NaN
Fixed ($v=51$) | 2754k | 4254K | 7157k | 13454k | NaN | NaN | NaN
Fixed ($v=99$) | 41564 | 8170k | 14160k | 30054k | NaN | NaN | NaN
Fixed ($v=1$) | 50k | 850k | 1850k | NaN | NaN | NaN | NaN
Fixed ($v=3$) | 1150k | 1450k | 3150k | NaN | NaN | NaN | NaN
80% Noise |
Scheduled ($v=(11,15,25,51)$) | 1650k | 2351k | 4051k | 10051k | 22554k | NaN | NaN
Scheduled ($v=(7,11,15,25)$) | 3250k | 8751k | 10051k | 12351k | 29052k | NaN | NaN
Scheduled ($v=(5,7,11,15,25)$) | 7451k | 10050k | 10150k | 15052k | 29353k | NaN | NaN
Scheduled ($v=(3,5,7,11,15,25)$) | 11250k | 13150k | 15251k | 19850k | 34654k | NaN | NaN
Fixed ($v=15$) | 1552k | 2552k | 3953k | 8353k | 42952k | NaN | NaN
Fixed ($v=25$) | 2053k | 3553k | 4551k | 8754k | 45352k | NaN | NaN
Fixed ($v=11$) | 1351k | 2351k | 3952k | 13651k | NaN | NaN | NaN
Fixed ($v=51$) | 1958k | 3458k | 5553k | 13858k | NaN | NaN | NaN
Fixed ($v=99$) | 3852k | 6852k | 11160k | 28752k | NaN | NaN | NaN
Chi ($p<0.05$) | 3655k | 7759k | 13056k | 43557k | NaN | NaN | NaN
Fixed ($v=7$) | 3350k | 4450k | 28450k | NaN | NaN | NaN | NaN
Fixed ($v=5$) | 7850k | 16750k | NaN | NaN | NaN | NaN | NaN
Fixed ($v=1$) | NaN | NaN | NaN | NaN | NaN | NaN | NaN
Fixed ($v=3$) | NaN | NaN | NaN | NaN | NaN | NaN | NaN

##  Appendix B MNIST level results as graphs

(a) Accuracy checkpoints for 20% noise

(b) Accuracy checkpoints for 40% noise

(a) Accuracy checkpoints for 60% noise

(b) Accuracy checkpoints for 80% noise

Figure 5: Comparing the different approaches by the speed with which they
achieved pre-defined test accuracy checkpoint. This is a different
representation of the data seen in Appendix A

.

##  Appendix C MNIST training curves

(a) Training curves for 20% noise

(b) Training curves for 40% noise

(a) Training curves for 60% noise

(b) Training curves for 80% noise

Figure 7: Test accuracy of trained networks over the course of training

##  Appendix D MNIST network architecture