# Model discovery in the sparse sampling regime

Gert-Jan Both, Georges Tod & Remy Kusters
Université de Paris, INSERM U1284
Center for Research and Interdisciplinarity (CRI)
F-75006 Paris, France
<EMAIL_ADDRESS>

###### Abstract

To improve the physical understanding and the predictions of complex dynamic
systems, such as ocean dynamics and weather predictions, it is of paramount
interest to identify interpretable models from coarsely and off-grid sampled
observations. In this work we investigate how deep learning can improve model
discovery of partial differential equations when the spacing between sensors
is large and the samples are not placed on a grid. We show how leveraging
physics informed neural network interpolation and automatic differentiation,
allow to better fit the data and its spatiotemporal derivatives, compared to
more classic spline interpolation and numerical differentiation techniques. As
a result, deep learning based model discovery allows to recover the underlying
equations, even when sensors are placed further apart than the data’s
characteristic length scale and in the presence of high noise levels. We
illustrate our claims on both synthetic and experimental data sets where
combinations of physical processes such as (non)-linear advection, reaction
and diffusion are correctly identified.

##  1 Introduction

Mathematical models are central in modelling complex dynamical processes such
as climate change, the spread of an epidemic or to design aircrafts. To derive
such models, conservation laws, physical principles and phenomenological
behaviors are key. However, some systems are too complex to model with a
purely bottom up approach Bolton & Zanna (2019); Sanchez-Pi et al. (2020).
When observational data is present, automated model discovery tools are
becoming increasingly more useful to derive partial differential equations
directly from the data. The classical method for data driven model discovery
is to apply sparse regression on a set of pre-selected features, the so-called
library. In the case of partial differential equations, this library is
constructed from a set of (higher)-order spatial derivatives. Model discovery
is thus effectively a two-step process: first construct the library, then
apply sparse regression. Numerically differentiating the data to construct the
library using finite differences is extremely sensitive to noise; in practice,
usually splines are fitted first and then differentiated. Splines model the
data as piece-wise polynomials, but this expansion breaks down when the
spacing between two sensors is large. These methods, which we refer to as
classical methods, thus fundamentally limit model discovery to densely sampled
data sets: even when no noise is present, the error incurred by the numerical
differentiation corrupts the library and renders the sparse regression
algorithm useless. The limits of classical interpolation methods have long
been known and are often cited as a reason to use neural networks instead.
Automatic differentiation can then be used to calculate the derivatives Baydin
et al. (2017), resulting in much more accurate derivatives. Previous works
Long et al. (2018); Both et al. (2021); Both & Kusters (2020) have shown that
using neural networks to create a surrogate of the data allows model discovery
in noisy and small data sets.

In this paper we systematically study how sample spacing influences model
discovery and compare neural-network based interpolation with classical
methods. Our focus is the influence of the differentiation method used to
construct the (higher-order) derivatives and its impact on model discovery, in
particular when the spacing between two sensors $\Delta x$ is larger than the
underlying equations’ characteristic length scale $l_{c}$. As NN-based model
discovery method we use DeepMoD, which is able to combine NN-based
interpolation with any sparse regression method Both et al. (2021); Both &
Kusters (2020). By using an identical sparse regression algorithm for both the
classical method and DeepMoD, we can isolate the effect of interpolation on
the library and the discovered equation. Our results show that NN-based
interpolators, in contrast to classical methods, can recover the underlying
equation when $\Delta x>l_{c}$. Furthermore, we show that NN-based
interpolation can succeed even when $\Delta x\gg l_{c}$ by either randomly
sampling or displacing the sampling grid over time. We corroborate our
findings with experimental data sets of the 2D advection-diffusion equation
and the 1D cable equation. In both cases, DeepMoD, discovers the underlying
equation in this sparsely sampled regime, contrarily to classical methods. Our
findings solidify the case for deep learning methods by showing that they
succeed in a regime where classical methods fundamentally fail.

##  2 Related works

#### Sensor placement

There exists a vast literature on determining optimal sensor placement for
control theory or signal reconstruction based on a library of features,
emphasizing the importance of sampling in the limit of sparse data Brunton et
al. (2013); Manohar et al. (2018); Wang et al. (2019). While many of these
sampling strategies have been developed to either reconstruct multi-scale
data-sets Champion et al. (2019b), flow-fields Brunton et al. (2015); Loiseau
et al. (2017) or other physical properties of a system Schaeffer et al.
(2018), research on the exact role of spatial and temporal sensor density or
distribution for model discovery has received limited attention.

#### Sparse regression-based model discovery

Using sparse regression to discover differential equations was popularized by
algorithms such as SINDY Brunton et al. (2016) and PDE-find Rudy et al.
(2017b) and has received considerable interest for both ODEs Mangan et al.
(2017); Messenger & Bortz (2020) as well as for PDEs Rudy et al. (2017a); Long
et al. (2018); Vaddireddy et al. (2020). These approaches have since been
expanded to automated hyper-parameter tuning Champion et al. (2019a); Maddu et
al. (2019); a Bayesian approach for model discovery using Sparse Bayesian
Learning Yuan et al. (2019), model discovery for parametric differential
equations Rudy et al. (2019).

#### Deep learning-based model discovery

With the advent of Physics Informed Neural Networks Raissi et al. (2017a; b),
a neural network has become one of the prime approaches to create a surrogate
of the data and perform sparse regression either on the networks prediction
Schaeffer (2017); Berg & Nyström (2019) or within the loss function of the
neural network Both et al. (2021); Both & Kusters (2020). Alternatively,
Neural ODEs were also used to discover the unknown governing equation
Rackauckas et al. (2020) from physical data-sets. Different optimisation
strategy based on the method of alternating direction is considered in Chen et
al. (2020), and graph based approaches have been developed recently Seo & Liu
(2019); Sanchez-Gonzalez et al. (2018). Finally, Cranmer et al. (2020);
Greydanus et al. (2019) directly encode symmetries in neural networks using
respectively the Hamiltonian and Lagrangian framework.

##  3 Methods

#### Sparse regression

A popular approach to discover a PDE from a spatio-temporal data set is to
apply sparse regression on a library of candidate terms $\Theta$, e.g. solve,

$u_{t}=f(1,u,u_{x},...)=\Theta\cdot\xi,$ (1)

to obtain the coefficient vector $\xi$. Here $u_{t}$ is the temporal
derivative and each column in $\Theta$ is a candidate term for the underlying
equation, typically a combination of polynomial and spatial derivative
functions (e.g. $u$, $u_{x}$, $uu_{x}$). To promote the sparsity of this
solution an $l_{1}$ regularization is added to the problem, leading to the so-
called Lasso regression:

$\xi^{*}=\min_{\xi}\left\lVert
u_{t}-\Theta\cdot\xi\right\rVert^{2}+\lambda\sum_{i}|\xi_{i}|.$ (2)

Here $\lambda$ controls the strength of the regularization, and hence the
resulting level of sparsity. In this paper we use the Lasso as a sparse
regression algorithm, with $\lambda$ determined by 5-fold cross-validation.
The resulting coefficients are normalized by the $l_{2}$ norm of the feature
vector, $\hat{\xi}_{i}=\xi_{i}\cdot||\Theta_{i}||_{2}/||u_{t}||_{2}$ and
thresholded. The exact value of the threshold can significantly influence the
resulting equation. We use exactly the same Lasso and threshold for both
DeepMoD and the classical methods so as to eliminate the influence of the
exact variable selection algorithm used. Our goal is to compare how the
feature library $\Theta$ and temporal derivative $u_{t}$ as generated by
either DeepMoD or a classical method differ, and its resulting effect on model
discovery.

#### Numerical differentiation

The features of the library $\Theta$ consists of (higher-order) derivatives,
which need to be differentiated from the observed data. Numerical
differentiation is typically performed either by finite differences or by
fitting a spline on the data and subsequently differentiating this spline.
Finite difference methods directly operate on the observed data to calculate
the derivative. In this paper, we use a standard second order accurate central
difference scheme. Finite differences is computationally cheap and easy to
scale to higher dimensions, but suffers from sensitivity to noise and requires
samples to be closely spaced for accurate results; the truncation error of the
scheme given above scales with the grid sampling, $h$, as
$\mathcal{O}\left(h^{2}\right)$. In the sparse regime where $\Delta x\to
l_{c}$, higher order schemes will not further improve this method.
Furthermore, FD requires samples on the edges of the domain to be discarded,
and for small data-sets and higher order schemes this can become a significant
fraction of the total data.

A more accurate and widely used alternative is to fit a spline to the data and
differentiate it. When fitting using splines, the data is approximated by a
piece-wise polynomial with enforced continuity at the edges. Splines yield
more accurate results in practice, but do not scale easily to higher
dimensions, especially when using splines of higher order. This hinders model
discovery, which requires these higher orders due to the derivatives in the
feature library; by using a fifth-order spline to approximate the data, we
effectively approximate the 3rd order derivative with only a second order
polynomial, hence hindering its application to model discovery.

#### DeepMoD

Both et al. (2021); Both & Kusters (2020)111github.com/PhIMaL/DeePyMoD is a
neural network-based model discovery algorithm. It uses a neural network to
learn both a noiseless surrogate of the data $\hat{u}$ and a coefficient
vector $\xi$ by minimizing,

$\mathcal{L}=\frac{1}{N}\sum_{i=1}^{N}\left(u_{i}-\hat{u}_{i}\right)^{2}+\frac{1}{N}\sum_{i=1}^{N}\left((\hat{u}_{t})_{i}-\Theta_{i}(\xi\cdot
g)\right)^{2}.$ (3)

Here $\Theta$ and $\hat{u}_{t}$ are calculated by automatic differentiation of
the neural network output $\hat{u}$. $g$ is a mask which sets the active
terms, i.e. the terms that feature in the differential equation. The first
term learns the data mapping $(x,t)\to\hat{u}$, while the second term
constrains the network to solutions of the partial differential equation given
by $\hat{u}_{t},\Theta$ and $\xi\cdot g$. During training, the coefficients
$\xi$ are determined by solving the least squares problem corresponding to the
second part of eq. 3. The mask $g$ is updated separately by a sparse
regression algorithm. The mask $g$ thus selects which terms are in the
equation, while $\xi$ are the coefficients of these active terms. The value of
the threshold can impact the discovered equation. To remove this factor from
our comparison, we use exactly the same method to find the sparse coefficient
vector $\xi^{*}$ for DeepMoD and the classical methods. More details on
DeepMoD can be found in Both et al. (2021); Both & Kusters (2020).

We emphasize two key differences with classical methods: 1) DeepMoD uses
automatic differentiation to calculate $\Theta$, and the accuracy of the
derivatives is thus not fundamentally linked to the sample distancing as with
numerical differentiation. 2) By including the regression term within the loss
function, we regularise the solution of the neural network $\hat{u}$, with the
learned solution of the PDE. The result of fitting a spline is solely based on
the data, whereas with DeepMoD it is also influenced by the constrained of the
underlying equation (i.e. $\xi$ and $g$). We show in the next section that
these two differences allow model discovery in extremely sparse and noisy data
sets, whereas classical methods fail.

##  4 Results

###  4.1 Synthetic data - Burgers equation

We consider a synthetic data set of the Burgers equation $u_{t}=\nu
u_{xx}-uu_{x}$, with a delta peak initial condition $u(x,t=0)=A\delta(x)$ and
domain $t\in[0.1,1.1],x\in[-3,4]$. This problem can be solved analytically
(see Appendix A) to yield a solution dependent on a dimensionless coordinate
$z=x/\sqrt{4\nu t}$. We recognize the denominator as a time-dependent length
scale: a Burgers data set sampled with spacing $\Delta x$ thus has a time-
dependent dimensionless spacing $\Delta z(t)$. We are interested in the
smallest characteristic length scale, which for this data set is
$l_{c}=\sqrt{4\nu t_{0}}$, where $t_{0}=0.1$ is the initial time of the data
set. We set $A=1$ and $\nu=0.25$, giving $l_{c}=\sqrt{0.1}\approx 0.3$.

Splines do not scale effectively beyond a single dimension, making it hard to
fairly compare across both the spatial and temporal dimensions. We thus study
the effect of spacing only along the spatial axis and minimize the effect of
discretization along the temporal axis by densely sampling 100 frames, i.e.
$\Delta t=0.01$. Along the spatial axis we vary the number of samples between
4 and 40, equivalent to $0.5<\frac{\Delta x}{l_{c}}<5$. We study the relative
error $\epsilon$ as the sum of the relative errors for all the derivatives,
normalized over every frame,

$\epsilon=\sum_{i=1}^{3}\left\langle\frac{\left\lVert\partial_{x}^{i}u_{j}-\partial_{x}^{i}\hat{u}_{j}\right\rVert_{2}}{\left\lVert\partial_{x}^{i}u_{j}\right\rVert_{2}}\right\rangle_{j}$
(4)

where $i$ sums the derivatives and $j$ runs over the frames. The derivatives
are normalised every frame by the $l_{2}$ norm of the ground truth to ensure
$\epsilon$ is independent of the magnitude of $u$. $\epsilon$ does not take
into account the nature of noise (e.g. if it is correlated and non-gaussian),
nor if the correct equation is discovered. However, taken together with a
success metric (i.e if the right equation was discovered), it does serve as a
useful proxy to the quality of the interpolation.

Figure 1b) shows $\epsilon$ as a function of the relative spacing $\Delta
x/l_{c}$ and whether the correct equation was discovered. The error when using
splines (yellow) increases with $\Delta x$ and, as expected, we are unable to
discover the correct equation for $\Delta x>0.8l_{c}$ (dots indicate the
correct equation is discovered and triangles indicates it failed to do so).
Considering the NN-based DeepMoD method, sampled on a grid (green), we observe
that it is able to accurately interpolate and discover the correct equation up
to $\Delta x\approx 1.2l_{c}$. The reason for this is that NN-based
interpolation constructs a surrogate of the data, informed by both the spatial
and the temporal dynamics of the data set, while classical interpolation is
intrinsically limited to a single time frame.

In figure 1c) we consider the same graph with $20\%$ white noise on the data.
Despite smoothing, the spline is unable to construct an accurate library and
fails to discover the correct equation in every case. DeepMoD stands in stark
contrast, discovering the correct equation with comparable relative error as
in the $0\%$ noise case.

Figure 1: a) The three sampling strategies considered. b) and c) Error in the
function library (Eq. 4) as function of the distance between the senors
$\Delta x$, normalized with $l_{c}=\sqrt{4\nu t_{0}}$, for b) noise-less data
and c) 20$\%$ of additive noise. The yellow symbols correspond to a spline
interpolation and the green, blue and red correspond to the NN-based model
discovery with various sampling strategies. The circles indicate that model
discovery was successful while the triangles indicate that the incorrect model
was discovered. The horizontal dashed line indicates the smallest
characteristic length-scale in the problem: $\Delta x/l_{c}=1$.

#### Off-grid sampling

Whereas higher-order splines are constrained to interpolating along a single
dimension, DeepMoD uses a neural network to interpolate along both the spatial
and temporal axis. This releases us from the constraint of on-grid sampling,
and we exploit this by constructing an alternative sampling method. We observe
that for a given number of samples $n$, DeepMoD is able to interpolate much
more accurately if these samples are randomly drawn from the sampling domain.
We show in figure 1b and c (Red) that the relative error $\epsilon$ in the
sparse regime, can be as much as two orders of magnitude lower compared to the
grid-sampled results at the same number of samples. We hypothesize that this
is due to the spatio-temporal interpolation of the network. By interpolating
along both axes, each sample effectively covers its surrounding area, and by
randomly sampling we cover more of the spatial sampling domain. Contrarily,
sampling on a grid leaves large areas uncovered; we are effectively sampling
at a much lower resolution than when using random sampling.

To test whether or not this improvement is intrinsically linked to the
randomness of sampling, we also construct an alternative sampling method
called shifted-grid sampling. Given a sampling grid with sensor distance
$\Delta x$, shifted-grid sampling translates this grid a distance $\Delta a$
every frame, leading to an effective sample distance of $\Delta a\ll\Delta x$.
This strategy, similarly as random sampling varies the sensor position over
time, but does so in a deterministic and grid-based way. We show this sampling
strategy in figure 1a, while panels b and c confirm our hypothesis; shifted
grid sampling (Blue) performs similarly to random sampling. Shifted-grid
sampling relies on a densely sampled temporal axis ’compensating’ for the
sparsely sampled spatial axes. This makes off-grid sampling beneficial when
either the time or space axis, but not both, can be sampled with a high
resolution. In the experimental section we show that if both axes are sparsely
sampled, we do not see a strong improvement over grid sampling.

###  4.2 Experimental data - 2D Advection-Diffusion

In an electrophoresis experiment, a charged dye is pipetted in a gel over
which a spatially uniform electric field is applied (see Figure 2a)). The dye
passively diffuses in the gel and is advected by the applied electric field,
giving rise to an advection-diffusion equation with advection in one
direction: $u_{t}=D(u_{xx}+u_{yy})+vu_{y}$. Both et al. (2021) showed that
_DeepMoD_ could discover the correct underlying equation from the full data-
set (size 120 x 150 pixels and 25 frames). Here, we study how much we can sub-
sample this data and still discover the advection-diffusion equation.

Figure 2: a) Experimental setup of the gel electrophoresis. b) Three time
frames of the density with a spatial resolution of 20x25. c) and d) Success
diagram for the experimental data indicating correct model discovery (Yellow
indicates the correct AD equation $u_{t}=D(u_{xx}+u_{yy})+vu_{y}$ is
discovered) as function of the spatial and temporal resolution for c) grid
sampling and d) random sampling. e) Obtained mask and coefficients ($D=0.025$
and and $v=(0,0.2)$) for the artificial data-set as function of the noise
level (11x11 spatial resolution). Hereby, yellow indicates the terms selected
by the algorithm and the red dashed box the terms that are expected in the
PDE. f) Success diagrams for various levels of additive noise, comparing the
result of DeepMoD with a grid and random sampling strategy and the classical
LassoCV algorithm on a Finite Difference (FD)-based library (after SVD
filtering of the different frames).

In figure 2 c) and d) we perform grid based as well as a random sub-sampling
of the data. The neural network-based method discovers the advection-diffusion
equation on as few as 6 x 8 spatial sampling points with 13 time-points, or
with 20 x 25 pixels on only 3 time-points. The minimum number of required
samples is similar for grid and random sampling, confirming that when both
axes are poorly sampled, there is no benefit to sample randomly.

The smallest characteristic length scale in the problem is the width of the
dye at $t=t_{0}$, which we estimate as $l_{c}\approx 10$ pixels. For the data
presented in figure 2c) and 2d), at a resolution below $10\times 13$ sensors
classical approaches would inherently fail, even if no noise was present in
the data set. This is indeed what we observe: using a finite difference-based
library we were unable to recover the advection-diffusion equation, even after
denoising with SVD (See Appendix A for details).

The use of a neural network and random sampling lead to non-deterministic
behaviour: the neural network training dynamics depend on its initialization
and two randomly sampled datasets of similar size might not lead to similar
results. In practice this leads to a ’soft’ decision boundary, where a
fraction of a set of runs with different initialization and datasets fail. We
discuss and study this issue in appendix B.

#### Noiseless synthetic data set

To further confirm our results from the previous section, we simulate the
experiment by solving the 2D advection-diffusion with a Gaussian initial
condition and experimentally determined parameters ($D=0.025$ and and
$v=(0,0.2)$. Figure 2e) shows the selected terms and their magnitude as
functions of the applied noise levels for a highly subsampled data-set (grid
sampling, spatial resolution of 11x11 and temporal resolution 14). The correct
AD equation is recovered up to noise levels of 100$\%$ (See figure 2e),
confirming the noise robustness of the NN-based model discovery. In panel f)
we compare the deep-learning based model discovery using grid and random
sampling with classical methods for various noise levels and sensor spacing
with a fixed temporal resolution of 81 frames (Data for the FD was pre-
processed with a SVD filter, see SI for further details). We confirm that,
similarly to the Burgers example of the previous section, the correct
underlying PDE is discovered even below the smallest characteristic length-
scale in the problem (indicated by a red dashed line in figure 2f).

This figure confirms our three main conclusions: 1) In the noiseless limit,
classical approaches are only slightly less performing than NN-based model
discovery for grid sampling. 2) Increasing the noise level dramatically
impacts classical model discovery while barely impacting NN-based model
discovery and 3) random sampling over space considerably improves performance,
allowing model discovery with roughly 4-8 times fewer sample points for this
particular data-set (depending on the noise level).

###  4.3 Experimental data - Cable equation

Applying a constant voltage to a RC-circuit with longitudinal resistance (see
figure 3 a) result in time-dependent voltage increase throughout the circuit
due to the charging of the capacitors. This rise is modeled by the cable
equation, which is essentially a reaction-diffusion equation
$u_{t}=u_{xx}/(R_{l}C)+u/(R_{m}C)$ with $C$ the capacitance, $R_{l}$ the
longitudinal resistance and $R_{m}$ the parallel resistance of the circuit.
The discrete nature of the experiment automatically gives $\Delta x=O(l_{c})$.
We consider an extreme case where we only have seven sensors throughout the
circuit (i.e. spatial axis), but take 2500 samples along the time axis. Figure
3b shows the measured voltage at these seven elements. Initially, all the
capacitors are uncharged and we observe a sharp voltage increase at the first
element. As the capacitors charge, this further propagates through the
circuit, charging the capacitors and resulting in the curves shown in the
figure. We apply both a classical approach with the library generated with
splines and DeepMoD to a varying amount of elements. Figure 3 c and d show
that the DeepMoD discovers the cable equation with as few as seven elements,
whereas classical methods are unable to find the cable equation at any number
of elements.

Figure 3: a) Schematic overview of the electronic setup to generate the cable
equation. b) The voltage drop $u$ as function of time for various positions
along the circuit for a circuit with 7 elements. The mask obtained for c) NN-
based and d) cross validated Lasso with spline based library model discovery
(Yellow indicates the term was recovered by the algorithm). The red boxes
indicate the two expected terms in the equation.

##  5 Discussion and future work

In this paper we showed how a deep learning approach allows to discover
partial differential equations from coarsely and off-grid sampled observations
in time and space. The correct equations are discovered, even when the sensor
spacing is larger than some data set’s characteristic length scale- an
inaccessible regime when using numerical differentiation procedures. We have
also shown that the presence of noise quickly deteriorates the performance of
classical methods, whereas the neural network based method is much less
affected. However, in the limit of very sparse data, model discovery can be
sensitive to the exact positioning of the sensors, hence sensitive to where
exactly on the grid the samples are drawn. Future work could investigate the
upper limit of the characteristic length scale above which the approach
consistently starts failing and how it relates to the spectrum of the data. We
will also focus in the future on including more structure in the interpolation
for better convergence and initialization robustness, for example using
Gaussian Processes.

## Acknowledgments

This work received support from the CRI Research Fellowship. We thank the
Bettencourt Schueller Foundation long term partnership and NVidia for
supplying the GPU under the Academic Grant program. We would also like to
thank the authors and contributors of Numpy (Harris et al. (2020)), Scipy
(Virtanen et al. (2020)), Scikit-learn (Pedregosa et al. (2011)), Matplotlib
(Hunter (2007)), Ipython (Perez & Granger (2007)), and Pytorch (Paszke et al.
(2019)) for making our work possible through their open-source software. The
authors declare no competing interest.

## References

  * Baydin et al. (2017) Atılım Günes Baydin, Barak A Pearlmutter, Alexey Andreyevich Radul, and Jeffrey Mark Siskind.  Automatic differentiation in machine learning: a survey.  _The Journal of Machine Learning Research_ , 18(1):5595–5637, 2017.
  * Berg & Nyström (2019) Jens Berg and Kaj Nyström.  Data-driven discovery of PDEs in complex datasets.  _Journal of Computational Physics_ , 384:239–252, May 2019\.  ISSN 00219991.  doi: 10.1016/j.jcp.2019.01.036.  URL http://arxiv.org/abs/1808.10788.  arXiv: 1808.10788.
  * Bolton & Zanna (2019) Thomas Bolton and Laure Zanna.  Applications of deep learning to ocean data inference and subgrid parameterization.  _Journal of Advances in Modeling Earth Systems_ , 11(1):376–399, 2019.
  * Both & Kusters (2020) Gert-Jan Both and Remy Kusters.  Sparsely constrained neural networks for model discovery of pdes.  _arXiv preprint arXiv:2011.04336_ , 2020.
  * Both et al. (2021) Gert-Jan Both, Subham Choudhury, Pierre Sens, and Remy Kusters.  Deepmod: Deep learning for model discovery in noisy data.  _Journal of Computational Physics_ , 428:109985, 2021.
  * Brunton et al. (2013) Bingni W Brunton, Steven L Brunton, Joshua L Proctor, and J Nathan Kutz.  Optimal sensor placement and enhanced sparsity for classification.  _arXiv preprint arXiv:1310.4217_ , 2013.
  * Brunton et al. (2015) Steven L Brunton, Joshua L Proctor, Jonathan H Tu, and J Nathan Kutz.  Compressed sensing and dynamic mode decomposition.  _Journal of computational dynamics_ , 2(2):165, 2015.
  * Brunton et al. (2016) Steven L. Brunton, Joshua L. Proctor, and J. Nathan Kutz.  Discovering governing equations from data by sparse identification of nonlinear dynamical systems.  _Proceedings of the National Academy of Sciences_ , 113(15):3932–3937, April 2016.  ISSN 0027-8424, 1091-6490.  doi: 10.1073/pnas.1517384113.  URL http://www.pnas.org/lookup/doi/10.1073/pnas.1517384113.
  * Champion et al. (2019a) Kathleen Champion, Bethany Lusch, J. Nathan Kutz, and Steven L. Brunton.  Data-driven discovery of coordinates and governing equations.  _arXiv:1904.02107 [stat]_ , March 2019a.  URL http://arxiv.org/abs/1904.02107.  arXiv: 1904.02107.
  * Champion et al. (2019b) Kathleen P Champion, Steven L Brunton, and J Nathan Kutz.  Discovery of nonlinear multiscale systems: Sampling strategies and embeddings.  _SIAM Journal on Applied Dynamical Systems_ , 18(1):312–333, 2019b.
  * Chen et al. (2020) Zhao Chen, Yang Liu, and Hao Sun.  Deep learning of physical laws from scarce data.  _arXiv:2005.03448 [physics, stat]_ , May 2020.  URL http://arxiv.org/abs/2005.03448.  arXiv: 2005.03448.
  * Cranmer et al. (2020) Miles Cranmer, Sam Greydanus, Stephan Hoyer, Peter Battaglia, David Spergel, and Shirley Ho.  Lagrangian Neural Networks.  _arXiv:2003.04630 [physics, stat]_ , March 2020.  URL http://arxiv.org/abs/2003.04630.  arXiv: 2003.04630.
  * Greydanus et al. (2019) Sam Greydanus, Misko Dzamba, and Jason Yosinski.  Hamiltonian Neural Networks.  _arXiv:1906.01563 [cs]_ , September 2019.  URL http://arxiv.org/abs/1906.01563.  arXiv: 1906.01563.
  * Harris et al. (2020) Charles R. Harris, K. Jarrod Millman, Stéfan J. van der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathaniel J. Smith, Robert Kern, Matti Picus, Stephan Hoyer, Marten H. van Kerkwijk, Matthew Brett, Allan Haldane, Jaime Fernández del Río, Mark Wiebe, Pearu Peterson, Pierre Gérard-Marchant, Kevin Sheppard, Tyler Reddy, Warren Weckesser, Hameer Abbasi, Christoph Gohlke, and Travis E. Oliphant.  Array programming with NumPy.  _Nature_ , 585(7825):357–362, September 2020\.  ISSN 0028-0836, 1476-4687.  doi: 10.1038/s41586-020-2649-2.  URL http://www.nature.com/articles/s41586-020-2649-2.
  * Hunter (2007) John D. Hunter.  Matplotlib: A 2D Graphics Environment.  _Computing in Science Engineering_ , 9(3):90–95, May 2007.  ISSN 1558-366X.  doi: 10.1109/MCSE.2007.55.  Conference Name: Computing in Science Engineering.
  * Loiseau et al. (2017) Jean-Christophe Loiseau, Bernd R Noack, and Steven L Brunton.  Sparse reduced-order modeling: sensor-based dynamics to full-state estimation.  _arXiv preprint arXiv:1706.03531_ , 2017.
  * Long et al. (2018) Zichao Long, Yiping Lu, Xianzhong Ma, and Bin Dong.  Pde-net: Learning pdes from data.  In _International Conference on Machine Learning_ , pp. 3208–3216, 2018.
  * Maddu et al. (2019) Suryanarayana Maddu, Bevan L. Cheeseman, Ivo F. Sbalzarini, and Christian L. Müller.  Stability selection enables robust learning of partial differential equations from limited noisy data.  _arXiv:1907.07810 [physics]_ , July 2019.  URL http://arxiv.org/abs/1907.07810.  arXiv: 1907.07810.
  * Mangan et al. (2017) Niall M Mangan, J Nathan Kutz, Steven L Brunton, and Joshua L Proctor.  Model selection for dynamical systems via sparse regression and information criteria.  _Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , 473(2204):20170009, 2017.
  * Manohar et al. (2018) Krithika Manohar, Bingni W Brunton, J Nathan Kutz, and Steven L Brunton.  Data-driven sparse sensor placement for reconstruction: Demonstrating the benefits of exploiting known patterns.  _IEEE Control Systems Magazine_ , 38(3):63–86, 2018.
  * Messenger & Bortz (2020) Daniel A Messenger and David M Bortz.  Weak sindy for partial differential equations.  _arXiv preprint arXiv:2007.02848_ , 2020.
  * Paszke et al. (2019) Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Köpf, Edward Yang, Zach DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala.  PyTorch: An Imperative Style, High-Performance Deep Learning Library.  _arXiv:1912.01703 [cs, stat]_ , December 2019.  URL http://arxiv.org/abs/1912.01703.  arXiv: 1912.01703.
  * Pedregosa et al. (2011) Fabian Pedregosa, Gaël Varoquaux, Alexandre Gramfort, Vincent Michel, Bertrand Thirion, Olivier Grisel, Mathieu Blondel, Peter Prettenhofer, Ron Weiss, Vincent Dubourg, et al.  Scikit-learn: Machine learning in python.  _the Journal of machine Learning research_ , 12:2825–2830, 2011.
  * Perez & Granger (2007) F. Perez and B. E. Granger.  IPython: A System for Interactive Scientific Computing.  _Computing in Science Engineering_ , 9(3):21–29, May 2007.  ISSN 1558-366X.  doi: 10.1109/MCSE.2007.53.
  * Rackauckas et al. (2020) Christopher Rackauckas, Yingbo Ma, Julius Martensen, Collin Warner, Kirill Zubov, Rohit Supekar, Dominic Skinner, and Ali Ramadhan.  Universal Differential Equations for Scientific Machine Learning.  _arXiv:2001.04385 [cs, math, q-bio, stat]_ , January 2020.  URL http://arxiv.org/abs/2001.04385.  arXiv: 2001.04385.
  * Raissi et al. (2017a) Maziar Raissi, Paris Perdikaris, and George Em Karniadakis.  Physics Informed Deep Learning (Part I): Data-driven Solutions of Nonlinear Partial Differential Equations.  _arXiv:1711.10561 [cs, math, stat]_ , November 2017a.  URL http://arxiv.org/abs/1711.10561.  arXiv: 1711.10561.
  * Raissi et al. (2017b) Maziar Raissi, Paris Perdikaris, and George Em Karniadakis.  Physics Informed Deep Learning (Part II): Data-driven Discovery of Nonlinear Partial Differential Equations.  _arXiv:1711.10566 [cs, math, stat]_ , November 2017b.  URL http://arxiv.org/abs/1711.10566.  arXiv: 1711.10566.
  * Rudy et al. (2017a) Samuel H Rudy, Steven L Brunton, Joshua L Proctor, and J Nathan Kutz.  Data-driven discovery of partial differential equations.  _Science Advances_ , 3(4):e1602614, 2017a.
  * Rudy et al. (2017b) Samuel H. Rudy, Steven L. Brunton, Joshua L. Proctor, and J. Nathan Kutz.  Data-driven discovery of partial differential equations.  _Science Advances_ , 3(4):e1602614, April 2017b.  ISSN 2375-2548.  doi: 10.1126/sciadv.1602614.  URL http://advances.sciencemag.org/lookup/doi/10.1126/sciadv.1602614.
  * Rudy et al. (2019) Samuel H. Rudy, J. Nathan Kutz, and Steven L. Brunton.  Deep learning of dynamics and signal-noise decomposition with time-stepping constraints.  _Journal of Computational Physics_ , 396:483–506, November 2019.  ISSN 00219991.  doi: 10.1016/j.jcp.2019.06.056.  URL http://arxiv.org/abs/1808.02578.  arXiv: 1808.02578.
  * Sanchez-Gonzalez et al. (2018) Alvaro Sanchez-Gonzalez, Nicolas Heess, Jost Tobias Springenberg, Josh Merel, Martin Riedmiller, Raia Hadsell, and Peter Battaglia.  Graph networks as learnable physics engines for inference and control.  _arXiv:1806.01242 [cs, stat]_ , June 2018.  URL http://arxiv.org/abs/1806.01242.  arXiv: 1806.01242.
  * Sanchez-Pi et al. (2020) Nayat Sanchez-Pi, Luis Marti, André Abreu, Olivier Bernard, Colomban de Vargas, Damien Eveillard, Alejandro Maass, Pablo A Marquet, Jacques Sainte-Marie, Julien Salomon, et al.  Artificial intelligence, machine learning and modeling for understanding the oceans and climate change.  In _NeurIPS 2020 Workshop-Tackling Climate Change with Machine Learning_ , 2020.
  * Schaeffer (2017) Hayden Schaeffer.  Learning partial differential equations via data discovery and sparse optimization.  _Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , 473(2197):20160446, January 2017.  ISSN 1364-5021, 1471-2946.  doi: 10.1098/rspa.2016.0446.  URL https://royalsocietypublishing.org/doi/10.1098/rspa.2016.0446.
  * Schaeffer et al. (2018) Hayden Schaeffer, Giang Tran, and Rachel Ward.  Extracting sparse high-dimensional dynamics from limited data.  _SIAM Journal on Applied Mathematics_ , 78(6):3279–3295, 2018.
  * Seo & Liu (2019) Sungyong Seo and Yan Liu.  Differentiable Physics-informed Graph Networks.  _arXiv:1902.02950 [cs, stat]_ , February 2019.  URL http://arxiv.org/abs/1902.02950.  arXiv: 1902.02950.
  * Vaddireddy et al. (2020) Harsha Vaddireddy, Adil Rasheed, Anne E Staples, and Omer San.  Feature engineering and symbolic regression methods for detecting hidden physics from sparse sensor observation data.  _Physics of Fluids_ , 32(1):015113, 2020.
  * Virtanen et al. (2020) Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, Stéfan J. van der Walt, Matthew Brett, Joshua Wilson, K. Jarrod Millman, Nikolay Mayorov, Andrew R. J. Nelson, Eric Jones, Robert Kern, Eric Larson, C. J. Carey, İlhan Polat, Yu Feng, Eric W. Moore, Jake VanderPlas, Denis Laxalde, Josef Perktold, Robert Cimrman, Ian Henriksen, E. A. Quintero, Charles R. Harris, Anne M. Archibald, Antônio H. Ribeiro, Fabian Pedregosa, Paul van Mulbregt, and SciPy 1 0 Contributors.  SciPy 1.0–Fundamental Algorithms for Scientific Computing in Python.  _Nature Methods_ , 17(3):261–272, March 2020\.  ISSN 1548-7091, 1548-7105.  doi: 10.1038/s41592-019-0686-2.  URL http://arxiv.org/abs/1907.10121.  arXiv: 1907.10121.
  * Wang et al. (2019) Zhi Wang, Han-Xiong Li, and Chunlin Chen.  Reinforcement learning-based optimal sensor placement for spatiotemporal modeling.  _IEEE transactions on cybernetics_ , 50(6):2861–2871, 2019.
  * Yuan et al. (2019) Ye Yuan, Junlin Li, Liang Li, Frank Jiang, Xiuchuan Tang, Fumin Zhang, Sheng Liu, Jorge Goncalves, Henning U. Voss, Xiuting Li, Jürgen Kurths, and Han Ding.  Machine Discovery of Partial Differential Equations from Spatiotemporal Data.  _arXiv:1909.06730 [physics, stat]_ , September 2019.  URL http://arxiv.org/abs/1909.06730.  arXiv: 1909.06730.

##  Appendix A Reproducibility

###  A.1 Hyperparameters

#### DeepMoD:

In this paper we use the neural network based model discovery tool DeepMoD
222github.com/PhIMaL/DeePyMoD. Every experiment uses a neural network with
tanh-activation functions and 4 layers of 30 neurons with random
initialization, and an Adam optimizer with default learning rate $10^{-3}$ and
$\beta=(0.9,0.9)$. The sparsity scheduler has a patience of 500 epochs and a
periodicity of 50 epochs Both & Kusters (2020).We use a cross validated,
thresholded Lasso sparsity selection with a threshold of 0.2 and otherwise
default parameters from the Sklearn implementation Pedregosa et al. (2011)

#### Spline interpolation:

For fitting the spline interpolation in both the Burgers as well as the Cable
equation, we use a smoothing parameter of $s=0.01$ in the case of noisy data
and 5th order splines.

#### Finite difference and SVD filter:

To construct the function library of the 2D Advection diffusion equation we
use a second-order accurate central difference scheme. For the 2D advection-
diffusion data, the data was denoised using by decomposing it using the SVD
and (Harris et al. (2020)) selecting the 3 largest modes from the signal.

#### Noise on synthetic data:

We add white noise to the data with a strength relative to the standard
deviation of the data, i.e. $50\%$ noise corresponds to $0.5\cdot\sigma$.

###  A.2 Data preparation

#### Burgers equation:

Using the Cole-Hopf transform, the Burgers equation described in the main text
reduces to the heat equation and can be solved exactly for a delta peak
initial condition to give,

$u(x,t)=\sqrt{\frac{\nu}{t\pi}}\left(\frac{(e^{R}-1)e^{-z^{2}}}{1+\frac{\left(e^{R}-1\right)}{2}\text{erfc}(z)}\right).$
(5)

where $R=A/2\nu$ and $z=x/\sqrt{4\nu t}$, a dimensionless coordinate. The
characteristic length-scale is thus the smallest on in the system; for our
case study $l_{c}=\sqrt{4\nu t}|_{t=t_{0}}$. We use a function library
containing all combinations of up to third order spatial derivative and second
order polynomials in $u$, for a total of 12 terms, i.e.,

$\Theta=\left[1,u_{x},u_{xx},u_{xxx},u,uu_{x},uu_{xx},uu_{xxx},u^{2},u^{2}u_{x},u^{2}u_{xx},u^{2}u_{xxx}\right].$
(6)

#### Cable equation:

We measured the passive voltage drop across a RC-circuit coupled to a
longitudinal resistance (See Fig. 3A). This voltage drop across the circuit
typically serves to model the passive voltage transmission through a dendrite,
and is described by the so-called cable equation,

$u_{t}=\frac{1}{R_{l}C}u_{xx}-\frac{1}{R_{m}C}u.$ (7)

Here $C$ is the capacitance. $R_{l}$ the longitudinal resistance and $R_{m}$
the membrane resistance. This equation can be discretizised by an electric
circuit, consisting of a serial set of $n$ longitudinal resistors, $r_{i}$,
membrane resistors, $r_{m}$, and capacitors, $c_{m}$. Using Ohm’s ans
Kirchhoff’s law, the discretized versioan of an array of these elements read,

$\frac{du_{i}}{dt}=\frac{(u_{i-1}+2u_{i}+u_{i+1})}{c_{m}r_{l}}-\frac{u_{i}}{c_{m}r_{m}}.$
(8)

We use a breadboard using structures imitating GMMs, using only standard
electronics hardware ($r_{m}=10k\Omega$, $r_{c}=270\Omega$ and $c_{m}=680mF$).
We applied a voltage profile across the electronics structure using an
arbitrary wave generator (AWG) (mhinstek MHS-2300A) and used a dual channel
oscilloscope (Voltcraft DSO-1062D) to measure the voltage at various positions
along the circuitry. These positions along the circuitry are the the spatial
dimension of the cable equation. We varied the amount of elements between 5
and 13, mimicking various spatial discretizations. At every sensor, we
collected 2500 data-points. We trained the model discovery algorithm on a
function library up to first order polynomials and second order derivatives.

###  A.3 2D Advection Diffusion

#### Experiment:

We consider the 2D advection-diffusion process described by,

$u_{t}=-\nabla\cdot\left(-D\nabla u+\vec{v}u\right).$ (9)

Here $\vec{v}$ is the velocity vector describing the advection and $D$ is the
diffusion coefficient. We measure a time-series of images from an
electrophoresis experiment, tracking the advection-diffusion of a charged
purple loading dye under the influence of a spatially uniform electric field.
We capture a set of 25 images with a resolution of 120x150 pixels and show the
resultant 2D density field for three separate time-frames (in arbitrary units)
in Fig. 2a, by subtracting the reference image (no dye present). The dye
displays a diffusive and advective motion with constant velocity $\vec{v}$,
which is related to the strength of the applied electric field. We use this
data-set to asses the impact of temporal as well as spatial sensor density on
the model discovery task. We used a cross validated thresholded Lasso sparsity
selection with a threshold of 0.05 and a function library containing all
combinations of up to third order spatial derivative and second order
polynomials in $u$, for a total of 10 terms,

${\Theta}=\left[1,u_{x},u_{y},u_{xx},u_{yy},u_{xy},u_{xxx},u_{yyy},u_{xxy},u_{xyy}\right].$
(10)

##  Appendix B Sensitivity to the random sampling

In this Appendix we discuss the sensitivity of the deep learning based
approach, DeepMoD, w.r.t. the set of random samples selected, in particular is
the limit $\Delta x/l_{c}>1$. In order to show the impact of the random set of
samples drawn we perform 10 runs of DeepMoD with otherwise identical
parameters (1000 samples drawn and $10\%$ white noise and otherwise identical
parameters as discussed in Appendix A). In Fig. 4a we show the outcome for
$\Delta x/l_{c}=2>1$ indicated that in 7 of the 10 cases the correct equation
is discovered and in 3 of the 10 cases this is not the case. In Fig. 4b) we
repeat this as function of $\Delta x/l_{c}$ and show that the larger the
average distance between the samples becomes, the more pronounced the
discrepancy between discovered models becomes. We have also tested the impact
of the initialization of the neural network on the outcome, with identical set
of samples and parameters, but this had little impact to the obtained PDE.

Figure 4: a) Coefficients obtained for the Burgers equation with 10$\%$ white
noise for 10 separate runs with 10 sets of randomly sampled data-sets. b)
Fraction of correctly discovered equations over 10 runs (with 10$\%$ white
noise and 1000 samples per run) as function of the average distance between
the samples, $\Delta x$, relative the the smallest characteristic length-scale
$l_{c}$.


capbtabboxtable[][]

11institutetext: Technical University Munich, Munich, Germany 22institutetext:
Helmholtz AI, Neuherberg, Germany 33institutetext: Institute for Computational
Biology, HelmholtzZentrum Munich, Germany 44institutetext: Munich School of
Data Science (MuDS), Munich, Germany 55institutetext: ContextVision AB,
Stockholm, Sweden

# Structure-Preserving Multi-Domain Stain Color Augmentation using Style-
Transfer with Disentangled Representations

Sophia J. Wagner 112244 Nadieh Khalili 55 Raghav Sharma 33 Melanie Boxberg
1144 Carsten Marr 33 Walter de Back 55 Tingying Peng 112244

###### Abstract

In digital pathology, different staining procedures and scanners cause
substantial color variations in whole-slide images (WSIs), especially across
different laboratories. These color shifts result in a poor generalization of
deep learning-based methods from the training domain to external pathology
data. To increase test performance, stain normalization techniques are used to
reduce the variance between training and test domain. Alternatively, color
augmentation can be applied during training leading to a more robust model
without the extra step of color normalization at test time. We propose a novel
color augmentation technique, HistAuGAN, that can simulate a wide variety of
realistic histology stain colors, thus making neural networks stain-invariant
when applied during training. Based on a generative adversarial network (GAN)
for image-to-image translation, our model disentangles the content of the
image, i.e., the morphological tissue structure, from the stain color
attributes. It can be trained on multiple domains and, therefore, learns to
cover different stain colors as well as other domain-specific variations
introduced in the slide preparation and imaging process. We demonstrate that
HistAuGAN outperforms conventional color augmentation techniques on a
classification task on the publicly available dataset Camelyon17 and show that
it is able to mitigate present batch effects. 111Code and model weights are
available at https://github.com/sophiajw/HistAuGAN.

###### Keywords:

color augmentation style-transfer disentangled representations.

##  1 Introduction

Modern cancer diagnosis relies on the expert analysis of tumor specimens and
biopsies. To highlight its structure and morphological properties,
conventionally, the tissue is stained with hematoxylin and eosin (H&E) [5].
The path from the raw tissue to the final digitized image slide however
consists of many different processing steps that can introduce variances, such
as tissue fixation duration, the age and the composition of the H&E-staining,
or scanner settings. Therefore, histological images show a large variety of
colors, not only differing between laboratories but also within one laboratory
[3].

This variability can lead to poor generalization of algorithms that are
trained on WSIs from a single source. One strategy to account for this is
stain color normalization. Traditionally, this is either done by aligning the
color distribution of the test images to a reference tile in the training
domain [12] or by decomposing the color space of a reference tile into
hematoxylin and eosin components [10, 17]. Then, H&E components of the test
tiles can be aligned while keeping the structure intact.

Recently, the focus shifted toward the application of style-transfer methods
such as cycle-consistent generative adversarial networks, CycleGAN [19], for
stain normalization [16]. However, these models aim to match the target
distribution possibly leading to undesired changes in the morphological
structure [6]. To circumvent this, other approaches propose color space
transformations [14], structural similarity loss functions [9], or residual
learning [4].

We propose a novel histological color transfer model, HistAuGAN, based on a
GAN architecture for image-to-image translation. In contrast to previous
approaches, HistAuGAN disentangles the content of a histological image, i.e.,
the morphological tissue structure, from the stain color attributes, hence
preserving the structure while altering the color. Therefore, HistAuGAN can be
used as a stain augmentation technique during training of a task-specific
convolutional neural network (CNN). We demonstrate that this helps to render
the trained network color-invariant and makes it transferable to external
datasets without an extra normalization step at test time. Applied as an
augmentation technique, HistAuGAN significantly outperforms other color
augmentation techniques on a binary tumor-classification task. Furthermore,
clustering results suggest that HistAuGAN can capture sources of domain shifts
beyond color variations, such as noise and artifacts introduced in the
staining or digitization process, e.g., image compression or blurring.

To the best of our knowledge, HistAuGAN is the first GAN-based color
augmentation technique that generates realistic histological color variations.

##  2 Method

###  2.1 Model architecture

Figure 1: We propose HistAuGAN for structure-preserving multi-domain stain
color augmentation. (a) Histological slides from different laboratories
(domains) exhibit color variations. (b) Model architecture. Here, the domain
information flow is visualized by colored arrows. (c) At inference, HistAuGAN
can be used as an augmentation technique by sampling attribute $z_{a}$ and
domain $d$.

We build our model based on a multi-domain GAN using disentangled
representations, inspired by DRIT++ [8]. Originally designed for image-to-
image translation of natural images using a predefined style, we propose its
application on histological images to disentangle the morphological tissue
structure from the visual appearance. In contrast to previous CycleGAN-based
color normalization methods that use only a single encoder, HistAuGAN is able
to separate two essential image properties from each other as visualized in
Figure 1b: the domain-invariant content encoder $E_{c}$ encodes the
histopathological structure of the tissue, e.g., size and position of the
nuclei, whereas the domain-specific attribute encoder $E_{a}$ learns the
domain-specific color appearance. The model can be trained on data from
multiple domains and thereby captures both inter-laboratory variability
between multiple domains and intra-laboratory variability within each domain
at the same time. Finally, the generator $G$ takes as input a content vector
$z_{c}$, an attribute vector $z_{a}$, and the one-hot-encoded domain vector
$d$ and outputs a simulated histological image. The objective function is
given by

$L_{total}=w_{cc}L^{cc}+w_{c}L^{c}+w_{d}L^{d}+w_{recon}L^{recon}+w_{latent}L^{latent}+w_{KL}L^{KL},$
(1)

where $L^{cc}$ is the cycle-consistency loss, $L^{c}$ and $L^{d}$ are
adversarial losses for the content and the attribute encoder, $L^{recon}$ is
an $L_{1}$-loss for image reconstruction, $L^{latent}$ is an $L_{1}$-loss for
latent space reconstruction, and $L^{KL}$ enforces the latent attribute space
to be distributed according to the standard normal distribution. Please refer
to [8] for a detailed explanation of each loss and the precise hyperparameter
setting.

Figure 2: Overview of the color variation in the dataset and the augmentation
techniques used in this paper using the framed image as example tile.

At inference, using the fixed content encoding of the input image $z_{c}$, we
can sample the attribute vector $z_{a}$ and the one-hot encoded domain vector
$d$ as visualized in Figure 1c. Hence, we can map one image to many different
structure-preserving augmentations. More specifically, we sample a random
color attribute $z_{a}$ from a normal distribution that parametrizes the stain
color variabilities in one domain. Figure 2b shows randomly sampled outcomes
of intra-domain augmentations. Additionally, we can change the one-hot-encoded
domain vector $d$ to project the input image into multiple target domains as
visualized in Figure 2c. In addition to sampling from the training domains, we
can also interpolate between these domains to obtain an even broader variety
of realistic color appearances for histopathological images. Figure 2d
demonstrates this by linearly interpolating the domain from domain RUMC to
domain UMCU according to

$d=(1-t)\cdot d_{\mathrm{RUMC}}+t\cdot d_{\mathrm{UMCU}},\quad\mathrm{for}\
t\in[0,1].$ (2)

###  2.2 Competing methods for stain color augmentation

Most existing stain color transfer methods are used for stain normalization,
i.e., to transfer the stain color of the test domain to that of the training
domain. Recently, it has been shown that simple stain color augmentations,
such as perturbing the HSV color space of the histological images, perform
better and lead to more robust models than traditional and network-based
normalization techniques [15]. Therefore, we compare our HistAuGAN to the HSV
augmentations used in [15]. Besides HSV augmentation, there is a more
complicated augmentation technique based on the Wasserstein distance of
different domains [11]. But the method is much slower than HSV and HistAuGAN,
thus difficult to be used as an on-the-fly augmentation technique.

For a quantitative evaluation of our augmentation technique, we consider the
following augmentation methods:

  * •

Geometric augmentations: vertical and horizontal flipping, as well as
$90^{\circ}$, $180^{\circ}$, and $270^{\circ}$ rotations.

  * •

HSV color augmentations: geometric augmentations with Gaussian blur and
contrast and brightness perturbations applied with probability 0.25 and 0.5,
respectively. We tried both light and strong color augmentations, as suggested
in [15]. Strong color augmentations can generate unrealistic color
appearances. However, applying hue and saturation jittering with factor 0.5
and probability 0.5, which results in relatively strong color perturbance as
shown in Figure 2e, performed best for us.

  * •

HistAuGAN: geometric augmentations combined with our augmentation technique
applied to half of the images during training. For each image, we randomly
pick a target domain from the training domains and sample a color attribute
vector $z_{a}\in\mathbb{R}^{8}$ from the standard normal distribution.

###  2.3 Evaluation

We evaluate HistAuGAN on three different aspects, in particular, i) whether it
can remove batch effects present in histological images collected from
multiple medical laboratories, ii) how it affects the out-of-domain
generalization of a deep learning model trained for a specific down-stream
task, and iii) how HistAuGAN preserves morphological structure during
augmentation. For ii), we choose a binary classification task of classifying
WSI tiles into the classes tumor versus non-tumor as described in more detail
in Section 3.3. Question iii) is evaluated by asking a pathology expert to
check image similarity before and after augmentation. To explore how
generalizable our model is, we extend the HistAuGAN training data (lymph
nodes) by tiles from unseen tissue and tumor types, in particular, breast
tissue [13].

##  3 Results and Discussion

###  3.1 Dataset

For the quantitative evaluation of HistAuGAN, we choose the publicly available
Camelyon17 dataset [1] that provides WSIs from five different medical centers
(denoted by RUMC, CWZ, UMCU, RST, and LPON) with different scanning properties
and stain colors as shown in Figure 2a. Pixel-wise annotations are given for
50 WSIs in total, 10 from each medical center. To create the training patches,
we first threshold the images with naive RGB thresholding combined with Otsu
thresholding and then patch the tissue regions of each WSI at the highest
resolution based on a grid into tiles of size $512\times 512$ pixels. Each
tile is labeled as tumor if the ratio of pixels annotated as tumor pixels is
larger than 1%, otherwise, it is labeled as non-tumor. The tiled dataset has
an imbalanced class distribution, i.e., overall, 7% of the tiles are labeled
as tumor and the ratio of tumor tiles is in the same order of magnitude across
all medical centers.

###  3.2 Evaluation of batch-effect removal

Figure 3: Effect of color augmentation on batch effects in color statistics.
(a-d) UMAP embeddings of color statistics of training data, color-coded by
source domains. (e) The quantification of mixing based on mean local diversity
(mLD, higher is better) suggests HistAuGAN effectively mitigates batch
effects.

To evaluate how color augmentation mitigates batch effects, we quantify the
mixing of images from different medical centers with respect to their color
statistics. A random set of 1,000 image tiles were extracted from the WSIs
from each center and analyzed in terms of the average values of each component
after transformation to various color spaces (RGB, HSV, LAB, HED, grayscale).
To visually observe batch effects, we reduced the dimensionality to 2D using
UMAP [2] and labeled points according to their domain as shown in Figure 3a-d.
To quantify the mixing of different domains, we measured the mean over the
local diversity (mLD) for all $k$-nearest neighborhoods ($k=10$) in the 2D
projection using Shannon’s equitability which varies between 0 for non-mixed
and 1 for perfectly mixed populations (cf. Figure 3e).

Without color augmentation, we observe a clear batch effect: tiles from
different domains form distinct clusters ($\mathrm{mLD}=0.2$, Figure 3a). HSV
augmentations improve data mixing, but domain-correlated clusters are still
visible ($\mathrm{mLD}=0.48$, Figure 3b) and single domains, e.g. LPON, are
not mixed with other domains. In contrast, HistAuGAN mixes data from multiple
domains (Figure 3c,d) with a high local diversity ($\mathrm{mLD}=0.85$). If
HistAuGAN is used to transfer colors to discrete domains, the distinct domain
clusters are retained, but each cluster contains well-mixed image samples
transferred from all domains (Figure 3c). When HistAuGAN is used to randomly
interpolate between domains, a continuous well-mixed color subspace is
obtained without any clustering structure (Figure 3d).

These results show that HistAuGAN is highly effective in removing batch
effects present in color statistics of images sampled from different medical
centers.

###  3.3 Evaluation on a down-stream classification task

Figure 4: Precision-recall AUC (left) and F1-score (right) of our binary
classification task. The bold bars depict the results on the out-of-domain
centers averaged across all runs. The most-right, pale bars denote the in-
domain test performance of the classifiers trained with geometric
augmentations.

To evaluate the effect of our proposed augmentation method, we train a CNN on
a binary tumor classification task and compare the performance on different
out-of-domain test sets based on the Camelyon17 dataset. Due to the relatively
small size of our dataset, in particular the small number of tumor tiles, we
choose a small CNN, namely, a pre-trained ResNet18 [7], and fine-tune the last
two ResNet-blocks together with the fully-connected layer on our dataset. For
training, we use weighted cross-entropy-loss to rebalance the contribution of
each class, with a learning rate of 1e-5 and an $L_{2}$-regularization of 1e-5
across all runs and for all augmentation techniques. Furthermore, we used
random erasing as regularization on all augmentation techniques [18]. Since
our dataset is highly imbalanced, we report the F1-score of the tumor class in
addition to the area under the precision-recall curve (PR-AUC).

Figure 4 shows the results of the quantitative evaluation of different
augmentation techniques on the binary tumor-classification task. For each
medical center, we trained three classifiers, one for each augmentation type,
and aggregated the results evaluated on the test domains. All experiments were
repeated three times. On both metrics, HistAuGAN shows better performance on
all of the out-of-domain test sets. As visualized in Figure 2, the appearance
of images from medical center UMCU and LPON deviates strongly from the other
centers, explaining their lower scores. In comparison to HSV color
augmentation, HistAuGAN performs better in handling the stain color
discrepancy between training and test domain and is therefore able to generate
a more robust classification model that generalizes better to out-of-domain
test sets. This can also be measured in the standard deviation of the results
across the out-of-domain test sets centers. For our model, the standard
deviation of the PR-AUC for the tumor class amounts to 0.08, whereas it higher
for geometric (0.22) and color (0.14) augmentations, respectively, which
demonstrates that our model is more robust to underlying stain color
variations. The right-most group shows the in-domain test results for
geometric augmentations. It can be seen as an upper bound for any stain
normalization technique, and thus shows that HistAuGAN can even outperform
stain normalization techniques on some of the five domains.

###  3.4 Qualitative evaluation by an expert pathologist

We further check the quality of HistAuGAN by an expert pathologist on the
structural similarity of original and augmented WSI tiles from the training
set, i.e., the Camelyon17 dataset, and an unseen dataset of breast tissue
[13]. We define three levels of similarity: a) “High similarity”: a
pathologist would find it difficult to distinguish the original tile from the
augmented tile. b) “Moderate similarity”: some structural variations are
observed, but do not affect pathological diagnosis. c) “Low similarity”: the
augmentated tiles can not be used for diagnostic purposes. As shown in Table
6, most of the augmented images do not have a structural modification that
affects diagnosis and over half of them can even fool an expert pathologist.
It is worth mentioning that HistAuGAN is not trained on any of the breast
cancer images but is still able to transfer its color in a structure-
preserving manner as shown in Figure 6 on a sample tile.

Figure 5: Expert evaluation.

Tissue type | High | Moderate | Low | Total
---|---|---|---|---
Lymph nodes | 10 | 7 | 3 | 20
Breast | 14 | 4 | 2 | 20

Figure 6: HistAuGAN on unseen tissue.

##  4 Conclusion

In summary, we propose a novel GAN-based technique, HistAuGAN, for color
augmentation of histopathological images. Based on the disentangled
representations of content and style, HistAuGAN is able to change the color
appearance of an histological image while preserving its morphological
structure. Moreover, HistAuGAN captures both intra-domain and inter-domain
color variations. It is able to interpolate between domains and can therefore
span a continuous color space covering a large variety of realistic stain
colors. When applied as an augmentation technique, HistAuGAN yields a robust
down-stream classifier that generalizes better to out-of-domain test sets than
other color augmentations techniques and, therefore, renders additional stain
normalization steps unnecessary. Finally, HistAuGAN can mitigate batch effects
present in histopathological data which suggests that it is also able to cover
domain shifts beyond color variations, such as noise and artifacts introduced
in image compression. The code is publicly available at
https://github.com/sophiajw/HistAuGAN together with a model of HistAuGAN
trained on the five medical centers of the Camelyon17 dataset.

## References

  * [1] Bandi, P., Geessink, O., Manson, Q., Van Dijk, M., Balkenhol, M., Hermsen, M., Ehteshami Bejnordi, B., Lee, B., Paeng, K., Zhong, A., Li, Q., Zanjani, F.G., Zinger, S., Fukuta, K., Komura, D., Ovtcharov, V., Cheng, S., Zeng, S., Thagaard, J., Dahl, A.B., Lin, H., Chen, H., Jacobsson, L., Hedlund, M., Cetin, M., Halici, E., Jackson, H., Chen, R., Both, F., Franke, J., Kusters-Vandevelde, H., Vreuls, W., Bult, P., van Ginneken, B., van der Laak, J., Litjens, G.: From detection of individual metastases to classification of lymph node status at the patient level: The CAMELYON17 challenge. IEEE Trans. Med. Imaging 38(2), 550–560 (Feb 2019)
  * [2] Becht, E., McInnes, L., Healy, J., Dutertre, C.A., Kwok, I.W.H., Ng, L.G., Ginhoux, F., Newell, E.W.: Dimensionality reduction for visualizing single-cell data using UMAP. Nat. Biotechnol. (Dec 2018)
  * [3] Bejnordi, B.E., Litjens, G., Timofeeva, N., Otte-Holler, I., Homeyer, A., Karssemeijer, N., van der Laak, J.A.: Stain specific standardization of Whole-Slide histopathological images (2016)
  * [4] de Bel, T., Bokhorst, J.M., van der Laak, J., Litjens, G.: Residual cyclegan for robust domain transformation of histopathological tissue slides. Med. Image Anal. 70, 102004 (May 2021)
  * [5] Chan, J.K.C.: The wonderful colors of the Hematoxylin–Eosin stain in diagnostic surgical pathology. Int. J. Surg. Pathol. 22(1), 12–32 (Feb 2014)
  * [6] Cohen, J.P., Luck, M., Honari, S.: Distribution matching losses can hallucinate features in medical image translation. In: Frangi, A.F., Schnabel, J.A., Davatzikos, C., Alberola-López, C., Fichtinger, G. (eds.) Medical Image Computing and Computer Assisted Intervention - MICCAI 2018 - 21st International Conference, Granada, Spain, September 16-20, 2018, Proceedings, Part I. Lecture Notes in Computer Science, vol. 11070, pp. 529–536. Springer (2018). https://doi.org/10.1007/978-3-030-00928-1_60, https://doi.org/10.1007/978-3-030-00928-1_60
  * [7] He, K., Zhang, X., Ren, S., Sun, J.: Deep residual learning for image recognition. In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR) (June 2016)
  * [8] Lee, H.Y., Tseng, H.Y., Mao, Q., Huang, J.B., Lu, Y.D., Singh, M., Yang, M.H.: DRIT : Diverse Image-to-Image translation via disentangled representations (2020)
  * [9] Liang, H., Plataniotis, K.N., Li, X.: Stain style transfer of histopathology images via Structure-Preserved generative learning. In: Machine Learning for Medical Image Reconstruction. pp. 153–162. Springer International Publishing (2020)
  * [10] Macenko, M., Niethammer, M., Marron, J.S., Borland, D., Woosley, J.T., Guan, X., Schmitt, C., Thomas, N.E.: A method for normalizing histology slides for quantitative analysis. In: 2009 IEEE International Symposium on Biomedical Imaging: From Nano to Macro. pp. 1107–1110. IEEE (2009)
  * [11] Nadeem, S., Hollmann, T., Tannenbaum, A.: Multimarginal wasserstein barycenter for stain normalization and augmentation. Med. Image Comput. Comput. Assist. Interv. 12265, 362–371 (Oct 2020)
  * [12] Reinhard, E., Adhikhmin, M., Gooch, B., Shirley, P.: Color transfer between images. IEEE Comput. Graph. Appl. 21(5), 34–41 (Jul 2001)
  * [13] Roux, L.: Mitos-atypia-14 grand challenge. https://mitos-atypia-14.grand-challenge.org/, accessed: 2021-03-03
  * [14] Shaban, M.T., Tarek Shaban, M., Baur, C., Navab, N., Albarqouni, S.: Staingan: Stain style transfer for digital histological images (2019)
  * [15] Tellez, D., Litjens, G., Bándi, P., Bulten, W., Bokhorst, J.M., Ciompi, F., van der Laak, J.: Quantifying the effects of data augmentation and stain color normalization in convolutional neural networks for computational pathology. Med. Image Anal. 58, 101544 (Dec 2019)
  * [16] Tschuchnig, M.E., Oostingh, G.J., Gadermayr, M.: Generative adversarial networks in digital pathology: A survey on trends and future potential. Patterns (N Y) 1(6), 100089 (Sep 2020)
  * [17] Vahadane, A., Peng, T., Sethi, A., Albarqouni, S., Wang, L., Baust, M., Steiger, K., Schlitter, A.M., Esposito, I., Navab, N.: Structure-Preserving color normalization and sparse stain separation for histological images. IEEE Trans. Med. Imaging 35(8), 1962–1971 (Aug 2016)
  * [18] Zhong, Z., Zheng, L., Kang, G., Li, S., Yang, Y.: Random erasing data augmentation. Proceedings of the AAAI Conference on Artificial Intelligence 34(07), 13001–13008 (Apr 2020). https://doi.org/10.1609/aaai.v34i07.7000, https://ojs.aaai.org/index.php/AAAI/article/view/7000
  * [19] Zhu, J.Y., Park, T., Isola, P., Efros, A.A.: Unpaired image-to-image translation using cycle-consistent adversarial networks. In: Proceedings of the IEEE international conference on computer vision. pp. 2223–2232 (2017)