# Fermion Sampling Made More Efficient

Haoran Sun${}^{\color{blue}\dagger}$ Cavendish Laboratory, University of
Cambridge, Cambridge, CB3 0HE, U.K.  State Key Laboratory of Surface Physics,
Institute of Nanoelectronics and Quantum Computing, and Department of Physics,
Fudan University, Shanghai 200438, China  Jie Zou${}^{\color{blue}\dagger}$
State Key Laboratory of Surface Physics, Institute of Nanoelectronics and
Quantum Computing, and Department of Physics, Fudan University, Shanghai
200438, China  Xiaopeng Li<EMAIL_ADDRESS>State Key Laboratory of
Surface Physics, Institute of Nanoelectronics and Quantum Computing, and
Department of Physics, Fudan University, Shanghai 200438, China  Shanghai Qi
Zhi Institute, Shanghai 200030, China  Shanghai Research Center for Quantum
Sciences, Shanghai 201315, China

###### Abstract

Fermion sampling is to generate probability distribution of a many-body
Slater-determinant wavefunction, which is termed “determinantal point process”
in statistical analysis. For its inherently-embedded Pauli exclusion
principle, its application reaches beyond simulating fermionic quantum many-
body physics to constructing machine learning models for diversified datasets.
Here we propose a fermion sampling algorithm, which has a polynomial time-
complexity—quadratic in the fermion number and linear in the system size. This
algorithm is about $100\%$ more efficient in computation time than the best
known algorithms. In sampling the corresponding marginal distribution, our
algorithm has a more drastic improvement, achieving a scaling advantage. We
demonstrate its power on several test applications, including sampling
fermions in a many-body system and a machine learning task of text
summarization, and confirm its improved computation efficiency over other
methods by counting floating-point operations.

The fermion sampling problem is to generate a position configuration of a
quantum system composed of $N$ fermions, ${\bf x}=(x_{1},x_{2},\ldots,x_{N})$,
according to Slater-determinants, where the position $x_{j}$ is an integer
running from $1$ to the system size $L$. The single-particle wavefunctions of
the fermions are stored in a $L\times N$ matrix $U$ as a series of column
vectors. The orthonormal condition of the wavefunctions implies
$U^{\dagger}U=\mathbb{1}$. In this sampling problem, the configuration ${\bf
x}$ is drawn from a probability distribution

$P({\bf x})=\frac{1}{N!}|{\rm det}(U_{{\bf x},{\bf n}})|^{2},$ (1)

with ${\bf n}=(1,2,\ldots,N)$. This problem arises in modeling fermionic
quantum many-body physics, whose computation complexity causes numerical
challenges in simulating many-electron interference [1] and variational Monte
Carlo calculations [2]. Such numerical simulations play important roles in
understanding strongly correlated quantum many-body physics including electron
correlation effects in molecules [3] and in quantum materials [4], and
equation-of-state of neutron stars [5]. A standard traditional approach to
perform fermion sampling in physics is through Markov chain methods, where the
computation cost of the direct calculation of Slater-determinants is avoided
by taking sequential local updates. The resultant sampling complexity is
reduced to ${\cal O}(N^{3})$. However this approach becomes inefficient in
presence of large autocorrelation length in the Markov chain, which
generically appears in the simulation of quantum criticality [6] and many-body
localization in disorder systems [7].

In statistical analysis and machine learning applications, the fermion
sampling problem has been termed “determinantal point process” (DPP) [8]. It
has deep connection with random matrix theory [9] and completely integrable
systems [10] in mathematics. Due to the built-in negative correlation by the
quantum Pauli exclusion principle, DPP has been adopted as an elegant way to
enhance diversity [11]. This has triggered tremendous interests in machine
learning, as the ingredient of diversity is fundamentally crucial in such
learning tasks as recommendation systems [12], text summarization, and image
searching [13], and is difficult to characterise with other models [11]. One
basis underlying these machine learning applications is the Hough-Krishnapur-
Peres-Virag (HKPV) algorithm [14], whose sampling complexity has been improved
from the original ${\cal O}(N^{3}L)$ to ${\cal O}(N^{2}L)$ with certain
modification [15].

Here, we propose a novel fast fermion sampling (FFS) algorithm with a time
complexity of ${\cal O}(N^{3}+N^{2}L)$. This is an explicit sampling
algorithm, free of the autocorrelation problem present in the Markov chain
sampling, and about $100\%$ more efficient than the modified HKPV algorithm in
computation time. The ${\cal O}(N^{3})$ tail in the time complexity is less
dominant in general since the system size $L$ has to be larger than the
particle number $N$ due to Pauli exclusion principle. The improvement is more
dramatic for sampling the corresponding marginal distribution—our algorithm
then has a scaling advantage in the computation time. We demonstrate our
algorithm in both quantum physics and machine learning examples. In
application to sampling free fermions in a double-well, and interacting
fermions in disorder potentials, we find the sampling error by our FFS
algorithm is much smaller (more than ten times smaller in certain cases) than
the Markov chain sampling given the same level of computation resources. Their
distinction coincides with the autocorrelation length, as FFS is an explicit
sampling algorithm. In application to text summarization, we find our FFS
algorithm outperforms the modified HKPV algorithm in computation
efficiency—the float point operations with FFS are $50\%$ smaller than the
modified HKPV for the same learning task. This $100\%$ improvement in the
computation efficiency remains with increasing text-summary length, meaning
our FFS algorithm is substantially more powerful for heavy text summarization
tasks.

Our fast fermion sampling algorithm is based on a mathematical formula
established for the Slater determinant wavefunction [1],

$\sum_{\bf m}\left[\prod_{k}P(x_{k}|x_{1},\ldots x_{k-1};{\bf m})\right]=|{\rm
det}(U_{{\bf x},{\bf n}})|^{2},$ (2)

with ${\bf m}=(m_{1},m_{2},\ldots m_{N})$ running over all permutation of
${\bf n}$, and the conditional probability distribution, $P(x_{k}\mid
x_{1},\ldots x_{k-1};{\bf m})=\frac{1}{k!}\left|{\rm det}\left[U_{(x_{1}\ldots
x_{k}),(m_{1}\ldots m_{k})}\right]\right|^{2}.$ We first generate a vector
${\bf m}$ as a random permutation of ${\bf n}$. Then the fermion position
$x_{k}$ is drawn according to the conditional probability distribution, with
the index $k$ iteratively increased from $1$ to $N$ step by step. It is
guaranteed by Eq. (2) that the sampled configuration ${\bf
x}=(x_{1},x_{2},\ldots x_{N})$ obeys the required distribution in Eq. (1).

One key observation we make here is that the sampling $x_{k}$ from the
conditional probability $P(x_{k}\mid x_{1},\ldots x_{k-1};{\bf m})$ does not
require calculation of determinants, for their geometrical interpretation. A
determinant is equal to the volume of the high dimensional parallelotope
spanned by the row vectors. Denoting the $k$-dimensional row vectors in
$U_{(x_{1}\ldots x_{k}),(m_{1}\ldots m_{k})}$ as {${\bf u}_{x_{1}}$, ${\bf
u}_{x_{2}}$, …${\bf u}_{x_{k}}$}, the geometrical interpretation of the
determinant implies an important property for the conditional probability,

$P(x_{k}\mid x_{1},\ldots x_{k-1};{\bf m})\propto|{\bf u}_{x_{k}}^{T}{\bf
h}({\bf u}_{x_{1}},{\bf u}_{x_{2}},\ldots,{\bf u}_{x_{k-1}})|^{2},$ (3)

with ${\bf h}$ a $k$-dimensional normal vector perpendicular to the $k-1$
vectors {${\bf u}_{x_{1}},{\bf u}_{x_{2}},\ldots,{\bf u}_{x_{k-1}}$}. The
normal vector can be efficiently calculated with an iterative Gaussian
elimination method, which consumes $2kN$ operations (Supplementary
Information). Calculating the inner product ${\bf u}_{x_{k}}^{T}{\bf h}$ for
all $x_{k}$ yields a computation cost of $2Lk$. We then have a sampling
complexity at $k$-th step, with $2Nk+2Lk$ operations. Running $k$ from $1$ to
$N$, the overall operations of generating one sample configuration $\bf{x}$ is
$LN^{2}+N^{3}$, to the leading order of $L$ and $N$. However, the modified
HPKV algorithm requires $2LN^{2}$ operations in total [15], which is nearly
the twice of our method in the limit of $L\gg N$. Furthermore, different
sample configurations are completely uncorrelated, in sharp contrast to the
Markov chain sampling. In particular, for a computation task of sampling the
marginal distribution of $N_{<}$ ($<N$) fermions, the computation cost of our
FFS algorithm is ${\cal O}(LN_{<}^{2}+N_{<}^{3})$, having a scaling advantage
over the modified HKPV whose computation cost goes as ${\cal O}(LNN_{<})$.

Figure 1: Stimulation of double well potential via FFS and Markov chain
sampling methods. (A) A schematic illustration of double well potential and
its wave functions in three stages: balanced, transition, and imbalanced.
Colorful lines illustrate different orbital wave functions and the relative
positions correspond to the energy levels. All the orbits are strongly
localized in one potential well except for the highest orbit during the
transition stage. (B) Particle number density difference in two potential
wells under different level of skewness. (C) Error of FFS and Markov chain
methods around the transition point, comparing to a theoretical error assuming
a Bernoulli distribution (Supplementary Information). (D) illustrates that
squared error ratio $\eta$ is strongly related to autocorrelation length
$\tau$ in the Markov chain method.

In order to benchmark our algorithm, we first consider a test case of sampling
fermions in a double-well potential, and compare our algorithm with the Markov
chain sampling. We take a standard form of a one-dimensional double-well
potential, $V(x)=a(x^{2}-b)^{2}-cx$, with $x$ the spatial coordinate. The
system Hamiltonian upon space discretization is given by

$H=-t\sum_{\langle ij\rangle}c_{i}^{\dagger}c_{j}+\sum_{i}v_{i}n_{i},$

with parameter $v_{i}$ taken from the double well potential. The parameters
$a$, and $b$ in $V(x)$ controls the depth of potential well and the parameter
$c$ controls its skewness. In order to achieve strong particle localization of
the system, the depth parameter $a$ is set $10^{3}$ times of the tunneling
strength in our implementation. In small $c$ limit, fermionic particles are
evenly distributed in both potential well at their ground state. A finite
value of $c$ introduce a break of symmetry to the system, and induce a
balance-imbalance transition. At a critical $c$ value, the difference between
the particle number in two potential wells exhibits a sharp transition from
zero to two, and an increase in autocorrelation length is observed. We measure
the particle density in each potential well, and compare the performance of
the Markov chain and FFS algorithms. Given a constant total running time, a
significantly larger error is observed at the critical point for the Markov
chain sampling. In contrast, our algorithm shows no error increase other than
the physical increase of the intrinsic observable variance in the system (Fig.
1C). By measuring an error ratio
$\eta=\frac{\epsilon_{\text{Markov}}}{\epsilon_{\text{FFS}}}$, we show our
approach significantly outperforms traditional method by achieving higher
accuracy with same computational requirement (see Fig. 1(D)).

This rise of the simulation error at phase boundary due to critical slowing
down, is a key feature generically hindering the application of Markov chain
Monte Carlo sampling methods to various models at the critical point. Our
proposing fast fermion sampling algorithm is completely immune to the critical
slowing down problem.

Figure 2:  Simulation of Anderson-Hubbard model via FFS and Markov chain
methods. The error is measured according to their sampling results on the
occupation number on the first lattice site. (A) shows the errors with both
methods and the square root of autocorrelation length in the Markov chain
sampling for the Anderson model ($V=0$) with increasing disorder strength $W$.
(B), The error ratio $\eta$ between Markov chain sampling and FFS in
simulating the Anderson-Hubbard model with different interaction ($V$) and
disorder strengths ($W$). Our algorithm outperforms for all interaction $V$
and disorder strength $W$ we choose, and the comparison is more substantial
for stronger disorder and weaker interaction. All the calculations above are
implemented on a $4\times 4$ square lattice with the particle number $N=8$.

We also apply our algorithm to sample a correlated system of interacting
fermions, described by a Slater-Jastrow wave function

$|\psi_{J}\rangle=\exp(-\frac{1}{2}V\textdagger\sum_{i}n_{i,\uparrow}n_{i,\downarrow})|\phi_{0}\rangle,$
(4)

We choose the the non-interacting part $|\phi_{0}\rangle$ to be $N$-fermion
ground state of a two-dimensional Anderson model,

$H=-t\sum_{\langle
ij\rangle}\sum_{\sigma}c_{i,\sigma}^{\dagger}c_{j,\sigma}++\sum_{i}w_{i}n_{i},,$
(5)

with $w_{i}$ drawn from a uniform distribution on $[-W,W]$. This problem
arises in variational Monte Carlo calculations of Anderson-Hubbard model. The
occupation number on the site with index $(i=0,j=0)$ is calculated with both
of our algorithm and Markov chain sampling for comparison (see Fig. 2). In
Markov chain sampling, to let the dynamics equilibrate we perform $2000$ steps
of local updates first and then start to take measurements. For a fair
comparison, we let the number of sampling configurations in our algorithm and
the number of measurements in Markov chain be equal (set to be $100$ here).
The sampling error with respect to the exact values by the two methods is
shown in Fig. 2—the exact values are obtained by sampling $10^{6}$ times with
our FFS algorithm. The sampling error by Markov chain systematically increases
with the disorder strengths, whereas the error remains at the same level for
different disorder strengths. We attribute this to the localization physics or
the tendency towards localization. With larger disorder strength, the system
develops tendency towards localization, and the Markov chain sampling becomes
less efficient in exploring the configuration space for the localization
causes nonergodicity, rendering a larger autocorrelation length. This problem
is completely absent in our algorithm for it produces independent samples. As
we increase correlation effects, the distinction between our algorithm and
Markov chain becomes less dramatic, which is as expected because interactions
effects in general compete with disorder-induced localization [7]. From these
results, it is evident that our fast fermion sampling algorithm systematically
outperforms the Markov chain method, especially when autocorrelation is
significant.

Figure 3:  Application of FFS to text summarization task. (A) A diagrammatic
representation of the general approach of text summarization task using
determinantal process. Different colors mark different sentences. (B) A simple
fermion sampling task for a direct comparison of computation time cost between
FFS and DPPy algorithms. We randomly generate a $N$-particle fermionic wave
function on a $L$-site 1D chain and measure the floating point operations
(FLOPs) required for producing a sample. Here we choose $L=1024$, which is
sufficiently larger than all possible $N$. The slope of DPPy is nearly twice
that of FFS. (C) FLOPs required to generate a summary from a 100 sentences
article, in the fermion sampling step, with different required summary length.

Last, we also apply the method to text summarization, to demonstrate the
potential impact of the fast fermion sampling algorithm to machine learning
applications. Text summarization is the task of producing a concise summary of
extended texts. One approach used in machine learning is to select sentences
in the texts that reflect the most important information. A balance should be
maintained between diversity and quality of chosen information in order to
make a concise summary. For such machine learning tasks, the quantum Pauli
exclusion principle built in the fermion sampling problem, termed
determinantal point process in statistical analysis, provides an elegant
approach to impose sampling repulsion [16].

For demonstration, we adopt a standard approach of using DPP for text
summarization [11], which mainly contains the following steps. First, all
sentences are converted to unit column-vectors ${\bf b}_{1}$, ${\bf b}_{2}$,
…, ${\bf b}_{L}$—the subscript labels correspond to different sentences, with
a standard word embedding algorithm [17], where the correlation among the
sentences is encoded into inner products of those unit vectors. The magnitude
of each vector is then assigned to reflect the quality [11] of the
corresponding sentence as a summary, giving new vectors $\tilde{{\bf b}}_{1}$,
$\tilde{{\bf b}}_{2}$, …, $\tilde{{\bf b}}_{L}$. The features considered here
include sentence length, sentence position, personal pronouns, etc. The
feature to magnitude mapping is achieved by a pretrained log-linear model
(Supplementary Information). A correlation ${\bf C}$ matrix is stored with
elements ${\bf C}_{ij}=\tilde{{\bf b}}_{i}^{T}\tilde{{\bf b}}_{j}$. Its
spectral decomposition is obtained as $C=\sum_{j}\lambda_{j}{\bf u}_{j}{\bf
u}_{j}^{T}$, with ${\bf u}_{j}$ the eigenvectors of dimension $L$. Second, we
draw a sample $(x_{1},x_{2},\ldots x_{N})$ following the standard DPP taking
the matrix $C$ as a kernel. Each eigenvector ${\bf u}_{j}$ is selected with
probability $\frac{\lambda_{j}}{\lambda_{j}+1}$, representing its importance
in the correlation matrix. The selected set of vectors, {${\bf u}_{j_{1}}$,
${\bf u}_{j_{2}}$, …${\bf u}_{j_{N}}$}, form the matrix $U$ in Eq. (1),
according to which we perform the fermion sampling with fermion particle
number $N$, and system size $L$. Our fast fermion sampling algorithm is
benchmarked against the modified HKPV algorithm widely used in DPP context.
The sampled configuration $(x_{1},x_{2},\ldots x_{N})$ marks the sentences to
be kept in the summary. Third, as the DPP process occasionally produces low-
quality summary [11], a most appropriate summary is produced by taking the
minimum Bayes risk (MBR) decoding [18] (Supplementary Information).

The overall complexity of the above scheme is ${\cal
O}(L^{2}+MN^{2}L+M^{2}N)$, with $M$ the number of samples taken in the second
step. Our algorithm is expected to have a constant speedup compared to the
modified HPKV method in the sampling procedure. Before implementing the text
summarization task, we carry out a pure fermion sampling task for a clean
comparison between two sampling algorithms. The results are shown in Fig. 3.
The results of the modified HPKV algorithm are obtained by the sampling module
in the DPPy [19] open source library, which implements the modified HKPV
algorithm. As shown in Fig. 3B, the computation cost of the two algorithms at
large system size scales the same way in the floating point operations
(FLOPs). The computation cost of our FFS algorithm is about one half of the
modified HKPV by counting FLOPs. This $100\%$ improvement in the computation
efficiency remains as we increase the system size. Then we perform a test on a
public CNN/Daily Mail news data set [20]. The results are presented in Fig.
3C. We show the counting of FLOPs solely in the fermion sampling step for a
direct comparison of our FFS with the modified HKPV algorithm in the text
summarization task. The FFS algorithm outperforms the modified HKPV for all
summary lengths consistently. The improvement is now less than $100\%$, which
can be attributed to the ${\cal O}(N^{3})$ tail in the computation cost of
FFS. But still the practical improvement with our FFS algorithm is quite
substantial at large summer length. Our approach then has promising
applications in accelerating the machine learning task of multiple-text
summarization, where the summary length can be a large number to cover all the
important information of multiple texts.

#### Discussion.

We expect the developed fast fermion sampling algorithm to have wide
applications in machine learning and quantum many-body physics. In all machine
learning tasks where describing data repulsion is important, DPP is an elegant
probabilistic model, and the fast fermion sampling algorithm has potential to
reduce the computation cost substantially. In applications to simulating
strongly interacting fermions or solving quantum chemistry problems, the fast
fermion sampling algorithm may be integrated with other importance sampling
methods to reach optimal computation performance.

We remark here that in performing fermion sampling for continuous models, a
direct implementation of our fast fermion sampling algorithm could be
inefficient because the time cost grows with $L$, which is typically a large
number in discretizing a continuous model to minimize discretization error. In
Supplementary Information, we show how this algorithm adapts to sampling
continuous models with minor modification, which leads to a computation cost
of ${\cal O}(N^{3})$.

#### Acknowledgement.

We acknowledge helpful discussion with Yang Qi. This work is supported by
National Program on Key Basic Research Project of China (Grant No.
2017YFA0304204), National Natural Science Foundation of China (Grants No.
11934002, and 11774067), Shanghai Municipal Science and Technology Major
Project (Grant No. 2019SHZDZCX01). We also acknowledge the support by the
Talented Student Program on the Fundamental Disciplines from Ministry of
Education of China.

${}^{\color{blue}\dagger}$ These authors contributed equally to this work.

## References

  * [1] X. Li, G. Zhu, M. Han, X. Wang, Phys. Rev. A 100, 032309 (2019).
  * [2] F. Becca, S. Sorella, Variational Monte Carlo (Cambridge University Press, 2017).
  * [3] W. Kohn, Rev. Mod. Phys. 71, 1253 (1999).
  * [4] B. Keimer, J. E. Moore, Nature Physics 13, 1045 (2017).
  * [5] G. A. Baker, Phys. Rev. C 60, 054311 (1999).
  * [6] P. Coleman, A. J. Schofield, Nature 433, 226 (2005).
  * [7] D. M. Basko, I. L. Aleiner, B. L. Altshuler, Annals of Physics 321, 1126 (2006).
  * [8] O. Macchi, Advances in Applied Probability 7, 83 (1975).
  * [9] T. Tao, Topics in random matrix theory, vol. 132 (American Mathematical Soc., 2012).
  * [10] M. Jimbo, T. Miwa, Y. Môri, M. Sato, Physica D: Nonlinear Phenomena 1, 80 (1980).
  * [11] A. Kulesza, B. Taskar, et al., Foundations and Trends® in Machine Learning 5, 123 (2012).
  * [12] L. Chen, G. Zhang, E. Zhou, Advances in Neural Information Processing Systems (2018), pp. 5622–5633.
  * [13] A. Kulesza, B. Taskar, Proceedings of the 28th International Conference on International Conference on Machine Learning (2011), pp. 1193–1200.
  * [14] J. B. Hough, M. Krishnapur, Y. Peres, B. Virág, et al., Probability surveys 3, 206 (2006).
  * [15] J. A. Gillenwater, PhD thesis, University of Pennsylvania (2014).
  * [16] A. Kulesza, B. Taskar, Proceedings of the Twenty-Seventh Conference on Uncertainty in Artificial Intelligence (2011), pp. 419–427.
  * [17] R. Řehůřek, P. Sojka, Proceedings of the LREC 2010 Workshop on New Challenges for NLP Frameworks (ELRA, Valletta, Malta, 2010), pp. 45–50.
  * [18] V. Goel, W. J. Byrne, Computer Speech & Language 14, 115 (2000).
  * [19] G. Gautier, G. Polito, R. Bardenet, M. Valko, Journal of Machine Learning Research - Machine Learning Open Source Software (JMLR-MLOSS) (2019). Code at http://github.com/guilgautier/DPPy/ Documentation at http://dppy.readthedocs.io/.
  * [20] A. See, P. J. Liu, C. D. Manning, Get to the point: Summarization with pointer-generator networks (2017).

Supplementary Information

##  S-1 The geometrical approach for determinantal sampling

As described in the main text, the major computation cost of our fast fermion
sampling algorithm consists in the calculation of the determinants in the
conditional probability distribution in Eq. (2). In sampling $x_{k}$ with $k$
running from $1$ to $N$, we need to calculate $L-k+1$ number of determinants
[${\rm det}\,U_{(x_{1},\ldots x_{k}),(m_{1},\ldots m_{k})}$] in each $k$-step,
whose direct computation is too costly.

In this section, we show how these determinants are calculated efficiently.
Here we exploit the geometric interpretation of determinants, that is, the
absolute value of the determinant ${\rm det}\,U_{(x_{1},\ldots
x_{k}),(m_{1},\ldots m_{k})}$ equals to the volume spanned by its
$k$-dimensional row vectors {${\bf u}_{x_{1}}$, ${\bf u}_{x_{2}}$, …${\bf
u}_{x_{k}}$}. In the $k$-th step, since the first $k-1$ row vectors of the
determinants are fixed, we can take them as the base and the remaining task of
computing the volume is to get the height (see main text).

To this end, we can calculate the unit normal vector $\bf{h}$ first, and then
the height equals to the projection of the last row vector
$\bf{u}_{x_{k}}\equiv U_{x_{k},(m_{1},\ldots m_{k})}$ onto $\bf{h}$. Of course
this normal vector can be derived by implementing the standard Gaussian
elimination method on the $(k-1)\times k$ submatrix $U_{(x_{1},\ldots
x_{k-1}),(m_{1},\ldots m_{k})}$, which takes a computation time of $O(k^{3})$.
However, instead of a naive implementation, we can compute the normal vector
iteratively, since Gaussian elimination operations on the first $k-2$ rows of
$U_{(x_{1},\ldots x_{k-1}),(m_{1},\ldots m_{k})}$ at step $k$ are exactly the
same as those of $U_{(x_{1},\ldots x_{k-2}),(m_{1},\ldots m_{k-1})}$ at step
$k-1$. Therefore, at step $k$, only the first $k-2$ entries in the $(k-1)$-th
row vector $\bf{u}_{x_{k-1}}$ needs to be eliminated. But we have not dealt
with the $k$-th column in $U_{(x_{1},\ldots x_{k-1}),(m_{1},\ldots m_{k})}$ by
now, which is also appended at the beginning of step $k$. This problem can be
naturally solved if we always carry out Gaussian elimination procedures in the
full $N$-dimensional space, that is, entries in column
$m_{k+1},m_{k+2},...,m_{N}$ are also updated at step $k$.

With the setup above, we can give a quantitative description for the
complexity of our algorithm. It requires $\sum_{j=0}^{k-3}(N-j)\approx Nk$
multiplications and approximately the same number of additions to implement
the iterative Gaussian elimination method at the $k$-th step. Then the normal
vector $\bf{h}$ can be determined with another $k(k-1)/2$ additions and
summations. Taking into account the computation of inner products of $L-k+1$
candidate vectors for $\bf{u}_{x_{k}}$ with the normal vector $\bf{h}$, our
algorithm consumes $2Nk+2Lk$ operations (including additions and summations)
at step $k$, to the leading order of $N$ and $L$. Therefore, with $k$ running
from 1 to $N$, the total operation count of the whole algorithm is
$LN^{2}+N^{3}$.

##  S-2 Determination of error and autocorrelation length

The error is defined by the root mean square of absolute error from single
trial. More specifically, for both algorithms, we take 100 samples in each
trial and obtain an averaged observable value from these samples. The trial is
repeated for approximately $10^{4}$ times (vary for different models) to
accumulate sufficient data. The absolute error for each trial is taken to be
the difference between value from single trial and an average of all trials in
the simulation of Anderson-Hubbard model. But we replace the average of all
trials with the theoretical value in the double-well model because here we
only deal with the noninteracting case which is exactly solvable. The
presented error is a root mean square of all these errors from different
trials.

The autocorrelation length (ACL) is calculated naively from a logarithm
fitting of the autocorrelation function $\hat{c}$ (ACF). Let
$\left\\{f_{n}\right\\}_{n=1}^{N}$ be a finite chain of input data with mean
$\mu$, the ACF is defined by

$\hat{c}(i)=\frac{1}{N-i}\sum_{n=1}^{N-i}\left(f_{n}-\mu\right)\left(f_{n+i}-\mu\right)$
(S1)

The ACF in a general Monte Carlo process is modeled by a exponential decrease
$\hat{c}(i)=\exp({-\frac{i}{\tau})}$. A logarithm fitting give raise to an
unbiased estimation of the ACL. We can also compute the ACL by discretize the
exponential relation above, that is to say,

$\tau=\sum_{i=0}^{\infty}\frac{\hat{c}(i)}{\hat{c}(0)}.$ (S2)

In practice we choose a sufficiently large cutoff for $i$ instead of a
summation to infinity.

##  S-3 Simulation details on Double Well

The system is stimulated in 64-site 1D lattices, with four non-interacting
spinless fermions. The double well potential is defined by

$V(x_{i})=a(x_{i}^{2}-b)^{2}-cx_{i}$ (S3)

with $a=4096$, $b=2$ and a varying $c$. The position of each lattice, $x_{i}$,
ranges evenly from $-2$ to $2$. The system tunneling strength, $t$, was taken
to be $\frac{\pi^{2}}{6}$. These values are chosen for particles to stay
strongly localized in certain potential well. A sharp transition between
balance and imbalance states of the system was observed around $c=89.929$.

The error of both methods is measured by comparing an exact solution and an
average of 100 consecutive samples. They are consequently compared with a
theoretical estimation of the natural error aroused from inherent uncertainty.
By assuming an ideal sampling of the system be Bernoulli process, with
particles falls into either side of the potential well, the theoretical error
is expressed by

$\epsilon=\sqrt{\frac{p(1-p)}{M}}$ (S4)

where $p=\sum_{x_{i}<0}|\phi_{i}|^{2}$ be the probability of the particle with
highest energy falls into the negative side of the potential, and $M$ be the
number of samples taken.

##  S-4 Simulation details on Anderson-Hubbard model

In the main text we show the strong relevance between error $\epsilon$ and
autocorrelation length $\tau$ in variational Monte Carlo (VMC) method through
the collapse of the curves of $\epsilon$ and $\sqrt{\tau}$. This square
relation is chosen from a data fitting analysis, as shown in Figure S2A. This
relation is also reasonable from a perspective of the central limit theorem
where error is proportional to $1/\sqrt{N}$ and $N$ is the amount of
statistically independent samples. The autocorrelation length reduces the
effective sample number up to a factor $1/\tau$ and the above relation is thus
derived.

For the interacting Anderson-Hubbard model, we have calculated the error ratio
of Markov chain sampling and FFS methods within a broad parameter regime. In
in Figure S2B. we show that for interacting cases the error of FFS still
nearly remains a constant and the one with Markov chain sampling grows
linearly with the increasing disorder strength. Besides, the appearance of
interaction has nothing to do with the error of FFS but relieve the trouble
brought from the localization for Markov chain sampling. Again, we compare the
error and the square root of autocorrelation length in Figure S2C, and the
relevance between them is tested. All the facts above are consistent with the
conclusions in our main text.

Figure S1: Simulation details on Anderson-Hubbard model. (A), a linear fit
between error and the square root of autocorrelation length $\sqrt{\tau}$.
(B), shows the error of FFS and Markov chain sampling with increasing disorder
strength $W$ for both noninteracting and interacting case. (C), shows the
error and autocorrelation length in Markov chain sampling with increasing
disorder strength for different interacting strengths. All the simulations
above are implemented on a $4\times 4$ square lattice and the particle number
is 8. The error is measured according to the occupation number of the first
site.

##  S-5 Implementation of our method on the Slater-Jastrow wave function

Jastrow wavefunction is a commonly used trial wave function in VMC, which is
written as

$|\psi_{J}\rangle=J_{v}|\phi_{0}\rangle=\exp(-\frac{1}{2}\sum_{i,j}v_{ij}n_{i}n_{j})|\phi_{0}\rangle,$
(S5)

where $n_{i}$ labels the particle number on site i and $|\phi_{0}\rangle$ is a
non-interacting wave function. $v_{ij}$s are variational parameters to be
optimized and the Jastrow factor $\hat{J}_{v}$ is expected to give a proper
description on the strong correlation between particles. In the main text we
choose the Gutzwiller factor, a simple version of Jastrow factor with only on-
site correlations.

Although our fermion sampling method is based on a noninteracting fermionic
system, it can still be applied in an interacting case with the reweighting
technique. For example, the expectation of a physical observerable $O$ can be
calculated as

$\langle
O\rangle=\frac{\langle\psi_{J}|O|\psi_{J}\rangle}{\langle\psi_{J}|\psi_{J}\rangle}=\frac{\sum_{x}\langle\psi_{J}|x\rangle\langle
x|O|\psi_{J}\rangle}{\sum_{x}\langle\psi_{J}|x\rangle\langle
x|\psi_{J}\rangle}.$ (S6)

Since the Jastrow factor $\hat{J}$ is diagonal in the position basis, the
R.H.S can be transformed as

$\frac{\sum_{x}J^{2}(x)|\psi_{0}(x)|^{2}(\langle x|O|\psi_{J}\rangle/\langle
x|\psi_{J}\rangle)}{\sum_{x}J^{2}(x)|\psi_{0}(x)|^{2}},$ (S7)

where $J(x)=\langle x|\hat{J}|x\rangle$ and $\psi_{0}(x)=\langle
x|\phi_{0}\rangle$ are both scalar functions. Then we implement the standard
reweighting procedure which is shown as

$\displaystyle\langle O\rangle$
$\displaystyle=\left(\frac{\sum_{x}J^{2}(x)|\psi_{0}(x)|^{2}o_{L}(x)}{\sum_{x}|\psi_{0}(x)|^{2}}\right)\bigg{/}\left(\frac{\sum_{x}J^{2}(x)|\psi_{0}(x)|^{2}}{\sum_{x}|\psi_{0}(x)|^{2}}\right)$
$\displaystyle=\langle J^{2}o_{L}\rangle_{0}/\langle J^{2}\rangle_{0}.$ (S8)

Here we denote $o_{L}(x)=\langle x|O|\psi_{J}\rangle/\langle
x|\psi_{J}\rangle$ and $\langle...\rangle_{0}$ labels the expectation value
derived from a free fermion ensemble. Therefore we can calculate the results
for an interacting system by sampling the corresponding noninteracting
ensemble and then reweighting according to S8.

##  S-6 Implementation details for text summarization

###  S-6.1 Feature to magnitude mapping

The feature to magnitude mapping is achieved by a log-linear model. For the
$i$th sentence in a given context $x$, we define a real valued feature
function $\boldsymbol{f}(x,i)\in\mathbb{R}^{m}$. This manually assigned
function concerns information of the sentence other than its content. The
magnitude, is consequently evaluated by

$r(x,i)=\exp(\boldsymbol{\theta}\cdot\boldsymbol{f}(x,i))$ (S9)

The parameter $\boldsymbol{\theta}$ is a linear weight of each terms in
feature function, and is trained with SGD method in our dataset to maximize
the overlap between machine summary and human summary. In this specific case,
SGD method is guaranteed to find the maximum as our characterization of the
overlap can be shown to be concave [13].

The features we considered are as follows:

  * •

A constant to control the summary length. The expected length of a summary
vary from different purpose, and therefore must be controlled manually.

  * •

Length of the sentence. Generally speaking, a medium length sentence makes a
good summary. Short sentences contains inadequate information and long
sentences are difficult to read at first glance. We adopted a polynomial to
characterize the quality based on sentence length, including up to the cubic
term.

  * •

Position in the article. Representative sentences usually appears at the
beginning and end of an article. We included the first three powers as well.

  * •

Dialogue and personal pronouns. Summary can usually be subjective and
narrative. We want to avoid choosing sentences with too much personal pronouns
that hinder the information of names, and also avoid objective statements. We
included a binary judgement of whether the sentence belongs to a dialogue, and
also the number of personal pronouns appears in the sentence.

###  S-6.2 MBR decoding

DPP model produce probabilistic results, so we must excluded summaries with
low quality. Here we select a summary result with minimum Bayes risk (MBR) of
producing unrelated result[18]. Let $\\{s_{i}\\}_{i=1}^{N}$ be a series of
summaries from a text, and $\\{\boldsymbol{v}_{i}\\}$ be the corresponding
unit vectors of these summaries from word embedding algorithm. The selected
summary $s_{MBR}$, according to MBR decoding method, is given by $s_{i}$ with

$\max_{i}\frac{1}{N}\sum_{j=1}^{N}\boldsymbol{v}_{i}\boldsymbol{v}_{j}$ (S10)

In the test of our program, 400 individual summaries are taken from each text
for MBR decoding. The value was picked to maintain a balance between summary
quality and running time. The program generally produce stable results (seldom
vary from different trials) at this level.

##  S-7 Fast fermion sampling for continuous models

In the main text, we provide an algorithm for a discrete lattice model having
computation complexity ${\cal O}(N^{2}L)$, with $N$ the fermion particle
number and $L$ the number of lattice sites. In application to a continuous
model where the fermions’ position take continuous values, this algorithm
becomes inefficient—we need to discretize the system for which a large number
of lattice sites $L$ is typically required to minimize the discretization
error. Here we describe how our fast fermion sampling algorithm adapts to
continuous models. In the first step, we still generate a vector ${\bf m}$ as
a random permutation of ${\bf n}$ (see main text). Then in generating the
fermion position $x_{k}$ from the conditional probability distribution
$P(x_{k}|x_{1},\ldots x_{k-1};{\bf m})$, we perform Markov chain sampling
instead of the explicit sampling as used for discrete models. A trial random
continuous update is proposed ($x_{k}\to x_{k}+\delta x$), and then accepted
with probability

$P_{\rm AC}(x_{k}\to x_{k}+\delta x)=\frac{|{\rm det}[U_{(x_{1}\ldots
x_{k}+\delta x),(m_{1}\ldots m_{k})}]|^{2}}{|{\rm det}[U_{(x_{1}\ldots
x_{k}),(m_{1}\ldots m_{k})}]|^{2}}=\frac{|{\bf u}_{x_{k}+\delta x}^{T}{\bf
h}({\bf u}_{x_{1}},{\bf u}_{x_{2}},\ldots,{\bf u}_{x_{k-1}})|^{2}}{|{\bf
u}_{x_{k}}^{T}{\bf h}({\bf u}_{x_{1}},{\bf u}_{x_{2}},\ldots,{\bf
u}_{x_{k-1}})|^{2}},$ (S11)

where {${\bf u}_{x_{1}},{\bf u}_{x_{2}},\ldots,{\bf u}_{x_{k}}$} are
$k$-dimensional row vectors in $U_{(x_{1}\ldots x_{k}),(m_{1}\ldots m_{k})}$
and ${\bf h}$ is a $k$-dimensional normal vector perpendicular to the $k-1$
vectors {${\bf u}_{x_{1}},{\bf u}_{x_{2}},\ldots,{\bf u}_{x_{k-1}}$}. The
update is iterated for $M_{\rm cond}$ times in order to let the Markov chain
equilibrate. This is the only modification we need to make for sampling
continuous models to the fast fermion sampling algorithm. The sampling
complexity is then ${\cal O}(N^{3})$, which is independent of $L$.

We emphasize that Markov chain is used here to sample a single-variable
probability distribution. The autocorrelation length in this sampling is at
most a constant, which does not scale up with the system size or the fermion
number. With $M_{\rm cond}$ chosen to be larger than the autocorrelation
length, the entire fermion sampling algorithm should still be taken as an
explicit sampling algorithm, free of any diverging autocorrelation length
problem.

It is worth mentioning here that the sampling algorithm for continuous models
described above also applies to problems with large $L$. We simply need to
replace the continuous update $x_{k}\to x_{k}+\delta x$ by a discrete update
$x_{k}\to x_{k}^{\prime}$. The sampling complexity is also ${\cal O}(N^{3})$.


Table 6:  Results of the validation procedure for the edges of the standard-configuration graph of the wind tunnel (Fig. 3c). For each edge, we show the p-value resulting from the Kolmogorov-Smirnov [52] test computed in step 5 of the procedure. The computed p-values are all below $1.8\times 10^{-5}$, except for the one corresponding to the edge $L_{\text{out}}\to\tilde{P}_{\text{int}}$ (0.082). The data to validate this edge was collected during a windy day (16.4.2024 in Zurich), and the effect on the intake barometer (a few pascals) was small compared to the large fluctuations (100 pascals) in the ambient atmospheric pressure. When using a more powerful rank-sum test [53]), we obtained a p-value of $0.029$. We will repeat the experiment with a longer running time (larger $N$) or on a day with fewer fluctuations in the atmospheric pressure. The experimental data, together with the interventions $x^{A},x^{B}$, waiting time $T$ and sample size $N$ for each edge can be found in the wt_validate_v1 dataset at causalchamber.org. Edge | p-value | Edge | p-value | Edge | p-value | Edge | p-value
---|---|---|---|---|---|---|---
$\theta_{1}\to\tilde{\text{I}}\text{m}$ | $2\times 10^{-13}\quad$ | $\theta_{2}\to\tilde{\text{I}}\text{m}$ | $1.6\times 10^{-14}\quad$ | $R\to\tilde{\text{I}}\text{m}$ | $2.3\times 10^{-14}\quad$ | $G\to\tilde{\text{I}}\text{m}$ | $1.6\times 10^{-14}\quad$
$B\to\tilde{\text{I}}\text{m}$ | $1.8\times 10^{-14}\quad$ | $T_{\text{Im}}\to\tilde{\text{I}}\text{m}$ | $2.9\times 10^{-14}\quad$ | $\text{Ap}\to\tilde{\text{I}}\text{m}$ | $2.3\times 10^{-14}\quad$ | $\text{ISO}\to\tilde{\text{I}}\text{m}$ | $4.2\times 10^{-14}\quad$
Table 7: Results of the validation procedure for the additional edges in the “camera” graph of the light tunnel (Fig. 3b). For each edge, we show the p-value resulting from the Kolmogorov-Smirnov [52] test computed in step 5 of the procedure. All p-values are below $10^{-12}$. For the image variable $\tilde{\text{I}}\text{m}$, the test is performed on the average of pixel values. The experimental data, together with the interventions $x^{A},x^{B}$, waiting time $T$ and sample size $N$ for each edge can be found in the lt_camera_validate_v1 dataset at causalchamber.org. Edge | p-value
---|---
$\tilde{P}_{\text{dw}}\to L_{\text{in}}$ | $2.3\times 10^{-14}\quad$
$\tilde{P}_{\text{dw}}\to L_{\text{out}}$ | $2.3\times 10^{-14}\quad$
Table 8: Results of the validation procedure for the additional edges of the
“pressure-control” graph of the wind tunnel (Fig. 3d). For each edge, we show
the p-value resulting from the Kolmogorov-Smirnov [52] test computed in step 5
of the procedure. The experimental data, together with the interventions
$x^{A},x^{B}$, waiting time $T$ and sample size $N$ for each edge can be found
in the wt_pc_validate_v1 dataset at causalchamber.org.

##  Appendix VI Component Datasheets

Component | Chamber |  Description and datasheets
---|---|---
light_source | wt/lt |  The light source of the light tunnel.  - LED array: datasheets/light_source.pdf - Individual LED (e.g. for wavelengths): datasheets/light_source_led.pdf
current_sensor | wt/lt |  Current sensor for the light source and fans. datasheets/current_sensor.pdf
motor | wt/lt |  The stepper motors to control the polarizers and the hatch.  - Motor: Model 17HS8401 in datasheet datasheets/motor.pdf - A4988 motor driver: datasheets/motor_driver.pdf
angle_sensor | lt |  The rotary potentiometer used to measure the polarizer angle. Model 3590S-6-502L in datasheet datasheets/angle_sensor.pdf
light_sensor | lt |  The light-intensity sensor of the light tunnel. Model Si1151-AB00-GMR in datasheet datasheets/light_sensor.pdf
led | lt |  The LEDs placed by each light sensor in the light tunnel.  - LEDs: datasheets/led.pdf - Digital rheostats: model MCP4151 in datasheets/potentiometer.pdf
camera | lt |  The light-tunnel camera, a Sony $\alpha$6100: datasheets/camera.pdf
arduino | wt/lt |  The control computer used in the chambers, i.e., an Arduino Mega Rev3: datasheets/arduino_mega.pdf
fan | wt |  The wind tunnel fans:  - Fan specifications: datasheets/fan.pdf - Speed control white paper: datasheets/fan_pwm.pdf - Specifications of a similar motor: datasheets/fan_motor.pdf
barometer | wt |  The high precision barometers used in the wind tunnel: datasheets/barometer.pdf
potentiometer | wt |  Digital potentiometer used to regulate amplitude in the amplification circuit of the wind tunnel. It is the same component as the rheostats of the light tunnel, i.e., model MCP4151 in datasheets/potentiometer.pdf
microphone | wt |  The wind tunnel microphone consists of an electret microphone with the LM385-DIP8 amplifier in datasheets/microphone_amplifier.pdf
speaker | wt |  The wind-tunnel speaker.  - Amplifier: model LM386-D08-T in datasheets/speaker_amplifier.pdf - Speaker: datasheets/speaker.pdf
Table 9: Physical components of the chambers and their related datasheets. The
datasheets can be found at github.com/juangamella/causal-chamber/.


11institutetext:  Data Science Institute, Imperial College London, London, UK
22institutetext: Department of Medicine, Imperial College London, London, UK
33institutetext: BioMedIA, Department of Computing, Imperial College London,
London, UK 44institutetext: NIHR Barts BRC, Queen Mary University of London,
London, UK 55institutetext: UK Dementia Research Institute, Imperial College
London, London, UK

# Self-Supervised Learning for Cardiac MR Image Segmentation by Anatomical
Position Prediction

Wenjia Bai 1122 Chen Chen 33 Giacomo Tarroni 33
Jinming Duan 33 Florian Guitton 11 Steffen E. Petersen 44
Yike Guo 11 Paul M. Matthews 2255 Daniel Rueckert 33

###### Abstract

In the recent years, convolutional neural networks have transformed the field
of medical image analysis due to their capacity to learn discriminative image
features for a variety of classification and regression tasks. However,
successfully learning these features requires a large amount of manually
annotated data, which is expensive to acquire and limited by the available
resources of expert image analysts. Therefore, unsupervised, weakly-supervised
and self-supervised feature learning techniques receive a lot of attention,
which aim to utilise the vast amount of available data, while at the same time
avoid or substantially reduce the effort of manual annotation. In this paper,
we propose a novel way for training a cardiac MR image segmentation network,
in which features are learnt in a self-supervised manner by predicting
anatomical positions. The anatomical positions serve as a supervisory signal
and do not require extra manual annotation. We demonstrate that this seemingly
simple task provides a strong signal for feature learning and with self-
supervised learning, we achieve a high segmentation accuracy that is better
than or comparable to a U-net trained from scratch, especially at a small data
setting. When only five annotated subjects are available, the proposed method
improves the mean Dice metric from 0.811 to 0.852 for short-axis image
segmentation, compared to the baseline U-net.

##  1 Introduction

Cardiac MR image segmentation plays a central role in characterising the
structure and function of the heart. Quantitative phenotypes derived from the
segmentations provide important biomarkers for diagnosing and managing
cardiovascular diseases. In the recent years, convolutional neural networks
have greatly advanced the performance of cardiac MR image segmentation due to
their capacity in learning discriminative image features for the segmentation
task [1, 2, 3]. Most successful methods are fully supervised and rely on a
large amount of annotated data to learn the features. However, annotated
medical imaging data may not always be available. The annotations are
expensive to acquire and often limited by the available resource of expert
image analysts. To address this challenge, we propose a novel way for training
a cardiac MR image segmentation network, which formulates a self-supervised
task for feature learning and alleviates the cost of data annotation.

There are two major contributions of this work. First, the proposed method
learns image features from anatomical positions automatically defined by
cardiac chamber view planes, which is a novel pretext task for self-supervised
learning and provides a strong supervisory signal. Importantly, the chamber
view plane information is freely available from standard cardiac MR scans,
which means the method has the potential to be extended to a clinical setting,
where a lot of unannotated MR scans stored on the PACS in hospitals can be
utilised for feature learning. Second, we demonstrate the learning performance
on two tasks, namely short-axis image and long-axis image segmentations. For
both tasks, self-supervised learning demonstrates a strong boost to
segmentation accuracy, especially at a small data setting.

###  1.0.1 Related Works:

To address the challenge of limited data annotations, there is increased
interest in developing methods that do not require a large amount of
annotations for feature learning. Directions of research include transfer
learning, domain adaptation, semi-supervised, weakly-supervised, unsupervised
and self-supervised learning [4]. Here, we focus on self-supervised learning,
which formulates a pretext task based on unannotated data for feature
learning.

For natural image and video analysis problems, a number of pretext tasks have
been explored, including prediction of image rotation [5], relative position
[6], colorisation [7] and image impainting [8] etc. In medical imaging domain,
self-supervised learning has also been explored but to a less extent.
Jamaludin et al. proposed a pretext task for subject identification [9]. A
Siamese network was trained to classify whether two spinal MR images came from
the same subject or not. The pretrained features were used to initialise a
disease grade classification network. Ross et al defined re-colourisation of
surgical videos as a pretext task and used the pretrained features to
initialise a surgical instrument segmentation network [10]. Tajbakhsh et al.
used rotation prediction as a pretext task and the self-learnt features were
transferred to lung lobe segmentation and nodule detection tasks [11].
Different from previous works in the medical imaging domain, we propose a
novel pretext task, which is to predict anatomical positions. In particular,
we leverage the rich information encoded in the cardiac MR scan view planes
and DICOM headers to define the anatomical positions for the task.

##  2 Methods

Here we describe the cardiac MR view planes, the pretext task for self-
supervised learning and architectures for transferring a self-trained network
to a new task.

(a) Cardiac MR view planes (b) Short-axis image (c) Long-axis 4Ch image

Figure 1: Cardiac MR view planes and anatomical positions. (a) Short-axis and
long-axis view planes with regard to the heart. (b) Short-axis image with
overlaid 2Ch and 4Ch view planes (yellow lines) and anatomical positions
defined by view planes (coloured boxes). (c) Long-axis 4Ch view image with
overlaid 2Ch and mid short-axis view planes and anatomical positions.

###  2.0.1 Cardiac MR view planes:

A standard cardiac MR scan consists of images acquired at different angulated
planes with regard to the heart, including short-axis, long-axis 2 chamber
(2Ch, vertical long-axis), 4 chamber (4Ch, horizontal long-axis) and 3 chamber
(3Ch) views. They are used for evaluating different anatomical regions of the
heart. For example, the short-axis view shows the cross-sections of the left
ventricle (LV) and right ventricle (RV). The long-axis views show the septal
and lateral walls of the ventricles, as well as the atrial chambers, including
the left atrium (LA) and right atrium (RA). Figure 1(a) illustrates how the
short-axis and long-axis 2Ch, 4Ch views are oriented with regard to the heart.

Most previous works on cardiac MR image segmentation [1, 2, 3] consider short-
axis and long-axis views separately and disregard the relative orientation of
different views. Typically, images from a specific view and the corresponding
label maps are used to train a segmentation network from scratch. In this
work, we propose that the relative orientation of short-axis and long-axis
views and the anatomical positions, defined by the view planes, can be used to
formulate a pretext task for training the network in a self-supervised manner
and increasing data efficiency.

###  2.0.2 Self-Supervised Learning (SSL):

Figure 1(b) shows a short-axis image, with the overlaid 2Ch view and 4Ch view
(yellow lines). As it shows, the 2Ch view bisects both the LV and RV, whereas
the 4Ch view bisects the LV. They intersect at the LV. Along the chamber view
lines, we define nine anatomical positions, represented by bounding boxes,
including the intersection, two boxes on the left, two on the right, two at
the anterior and two at the posterior. The orientations from left to right and
from posterior to anterior are available from the DICOM headers. The pretext
task is to predict the anatomical positions defined by these nine bounding
boxes. The intuition here is that for the network to recognise these
anatomical positions, it has to learn features for understanding not only
where the left and right ventricles roughly are but also what their
neighbouring regions look like. The learnt features can be transferred to a
related but more demanding task, which is accurate segmentation of the
ventricles.

Similarly, for a long-axis 4Ch view image, we can overlay the 2Ch view and mid
short-axis view on it, shown by Figure 1(c). Along the 2Ch view and mid short-
axis view lines, we define nine anatomical positions, including the
intersection, two boxes on the left, two on the right, two at the superior and
two at the inferior. The pretext task for long-axis image analysis is to
predict these anatomical positions. To learn from the pretext task, we train a
10-way segmentation network, which segments the nine bounding boxes and the
background. A standard U-net architecture [12] is used, which consists of the
encoder part, decoder part, skip connections between them and the task head
(the last convolutional layer), depicted by Figure 2(a). Cross-entropy is used
as the loss function.

(a) SSL (b) SSL+Decoder (c) SSL+All (d) SSL+MultiTask

Figure 2: Network architectures for self-supervised learning (SSL) and three
different ways for transfer learning. The gray area in (b) denotes the freezed
encoder.

###  2.0.3 Transfer Learning:

After the network is self-trained on the pretext task (task 1), it is
transferred to a new task (task 2), which is accurate segmentation of the
anatomical structures, e.g. the LV cavity, myocardium and RV cavity. To
achieve this, we can simply replace the head of task 1 by a new head for task
2. The task head refers to the last convolutional layer of the U-net, which is
1$\times$1 convolution with $K$-channel output, $K$ denoting the number of
classes.

We investigate three different ways for transfer learning. The first way is to
freeze the weights learnt at the encoder and only finetunes the decoder and
task head, using the annotations for task 2. This method is named as
“SSL+Decoder” and illustrated by Figure 2(b). The second way is to finetune
all the weights [4], including the encoder, decoder and task head. This method
is named as “SSL+All” and illustrated by Figure 2(c). The third way is perform
multi-task learning for finetuning. Two task heads are used, one for the
pretext task and the other for the new task. This is to avoid forgetting about
the pretext task while learning for the new task. This method is named as
“SSL+MultiTask” and illustrated by Figure 2(d). The loss function for multi-
task learning is formulated as,
$L(\theta)=L_{task_{1}}(x,y_{1}|\theta)+\beta\cdot
L_{task_{2}}(x,y_{2}|\theta)$, where $\theta$ denotes the network parameters,
$x$ denotes the image, $y_{1}$ denotes the label map of nine anatomical
positions for task 1, $y_{2}$ denotes the label map of anatomical structures
manually annotated by human experts for task 2, $\beta$ denotes the weight.

##  3 Experiments and Results

###  3.0.1 Data:

For self-supervised learning, short-axis and long-axis images of 3,825
subjects were used, acquired from the UK Biobank. The typical image dimension
is 208$\times$180, with 1.82$\times$1.82 mm2 in-plane resolution. The short-
axis image stack consists of $\sim$10 slices. There is 1 slice for long-axis
4Ch view and 1 slice for 2Ch view. Bounding boxes were automatically placed at
nine anatomical positions on the short-axis and long-axis 4Ch view images at
end-diastole (ED) and end-systole (ES). Each box was empirically set to
11$\times$11 pixels and adjacent boxes were 30 pixels apart. For transfer
learning, 200 subjects with manual annotations were used, which were randomly
split into 100 subjects for training and 100 subjects for test. For short-axis
images, LV, myocardium and RV were manually annotated by experienced image
analysts at ED and ES frames. For long-axis 4Ch view images, LV, myocardium,
RV, LA and RA were manually annotated.

###  3.0.2 Implementation:

The method was implemented using Tensorflow. For self-supervised learning, the
Adam optimiser was used, with a learning rate of 0.001, a batch size of 20
image slices and 50,000 iterations. For transfer learning, the same setting
was used. For both cases, data augmentation was performed online, including
random rotation and scaling. For multi-task transfer learning, the weight
$\beta$ was empirically set to 10 to emphasise the new task. As the input
datasets were not consistent for the two tasks (task 1 with 3,825 subjects,
task 2 with much fewer and different subjects), training was implemented as
that at each iteration, task 1 was optimised for one sub-iteration, followed
by optimising task 2 for $\beta$ sub-iterations. Since stochastic optimisation
was performed, this was approximately equivalent to assigning a weight to task
2. Training multiple tasks alternately is a commonly adopted practice for
inconsistent data input [4]. We made sure that task 2 was trained for 50,000
sub-iterations to enable a fair comparison. Finally, a U-net was also trained
with exactly the same setting, but initialised with random weights. This is
our baseline method, “U-net-scratch”.

Table 1: Comparison of the Dice metrics for short-axis image segmentation. Column 1 lists the number of training subjects and manually annotated image slices. Columns 2 to 5 report the performance of the baseline method and different self-supervised learning methods. Values are mean (standard deviation).  #subjects (#slices) | U-net-scratch | SSL+Decoder | SSL+All | SSL+MultiTask
---|---|---|---|---
1 (18) | 0.361 (0.047) | 0.515 (0.099) | 0.618 (0.068) |  0.704 (0.065)
5 (102) | 0.811 (0.037) | 0.837 (0.048) | 0.844 (0.046) |  0.852 (0.046)
10 (208) | 0.859 (0.037) | 0.860 (0.039) | 0.873 (0.036) |  0.875 (0.036)
50 (980) | 0.876 (0.035) | 0.871 (0.038) |  0.884 (0.033) | 0.873 (0.037)
100 (1,936) | 0.886 (0.030) | 0.867 (0.037) | 0.883 (0.035) |  0.887 (0.031)

###  3.0.3 Short-axis image segmentation:

We evaluated the performance on two transfer learning tasks, which are
segmentations for short-axis and long-axis images. Table 1 compares the Dice
metrics (averaged across LV, myocardium and RV) for short-axis image
segmentation between the U-net trained from scratch and self-supervised
learning methods. As it shows, even if we freeze the encoder and only tune the
decoder (SSL+Decoder), we can achieve a high accuracy comparable to training a
U-net from scratch. This indicates that SSL is able to learn good features at
the encoder which are transferrable for the segmentation task. The table also
shows when we tune all the weights (SSL+All and SSL+MultiTask), the
segmentation accuracy is generally better than U-net-scratch, especially when
the number of training subjects is small. On average, SSL+MultiTask performs
the best.

Figure 3: Short-axis image segmentations for U-net-scratch and SSL+MultiTask
with an increasing number of training subjects (slices), as well as manual
segmentations. The yellow arrows indicate segmentation errors.

Figure 3 visualises exemplar segmentations for U-net-scratch and SSL+MultiTask
with an increasing number of training subjects. It shows that when $n=1$, due
to the extremely small training set, U-net-scratch completely fails to segment
the image. On the contrary, SSL+MultiTask is still able to segment some part
of the LV and myocardium. When $n$ increases to 5 or 10, SSL+MultiTask
outperforms U-net-scratch at details, for example, without RV under-
segmentation errors. However, when $n$ increases to 50 or 100, the two methods
perform similar to each other. Figure 4 plots quantitative metrics including
the Dice metric and mean contour distance error for each anatomical structure
for the two methods. It shows a similar trend that at a small data setting
($n\leq 10$), SSL+MultiTask outperforms U-net-scratch for all the structures.
When there are more training data ($n\geq 50$), their performances become
close to each other.

Figure 4: Comparison of the Dice metrics and mean contour distance errors on short-axis image segmentation for U-Net-scratch and SSL+MultiTask.  Figure 5: Long-axis image segmentations for U-net-scratch and SSL+MultiTask with an increasing number of training subjects (slices).  Table 2: Comparison of Dice overlap metrics for long-axis image segmentation. Values are mean (standard deviation).  #subjects (#slices) | U-net-scratch | SSL+Decoder | SSL+All | SSL+MultiTask
---|---|---|---|---
1 (2) | 0.678 (0.132) | 0.699 (0.069) |  0.733 (0.081) | 0.600 (0.122)
5 (10) | 0.848 (0.080) | 0.850 (0.068) |  0.875 (0.051) | 0.861 (0.049)
10 (20) | 0.875 (0.044) | 0.888 (0.039) |  0.905 (0.039) | 0.889 (0.031)
50 (100) | 0.922 (0.031) | 0.914 (0.035) | 0.925 (0.028) |  0.926 (0.029)
100 (200) | 0.930 (0.032) | 0.924 (0.037) | 0.933 (0.031) |  0.934 (0.029)

###  3.0.4 Long-axis image segmentation:

We performed similar experiments for long-axis image segmentation. Table 2
reports the mean Dice overlap metrics. It shows that with SSL, for most of the
cases, the segmentation accuracy is increased compared to U-net-scratch.
Figure 5 visualises exemplar segmentation results. It demonstrates that with
limited training data, SSL+MultiTask generally produces better segmentations
compared to U-net-scratch.

##  4 Conclusions

In this paper, we propose a novel method that leverages self-supervised
learning for cardiac MR image segmentation. We formulate anatomical position
prediction as the pretext task. Experiments on short-axis and long-axis image
segmentation tasks demonstrate that with self-supervised learning, the
proposed method outperforms a standard U-net trained from scratch at the small
data setting and is of comparable performance at the large data setting. For
future work, we will explore other anatomically meaningful pretext tasks to
increase data efficiency in medical imaging applications.

### Acknowledgements:

This research has been conducted using the UK Biobank Resource under
Application Numbers 2964 and 18545 and supported by the SmartHeart EPSRC
Programme Grant (EP/P001009/1). We would like to thank NVIDIA Corporation for
donating a Titan Xp for this research.

## References

  * [1] O. Bernard et al.  Deep learning techniques for automatic MRI cardiac multi-structures segmentation and diagnosis.  IEEE Trans Med Imaging, 37(11):2514–2525, 2018.
  * [2] W. Bai et al.  Automated cardiovascular magnetic resonance image analysis with fully convolutional networks.  J Cardiovasc Magn Reson, 20(1):65, 2018.
  * [3] Q. Tao et al.  Deep learning-based method for fully automatic quantification of left ventricle function from cine MR images.  Radiology, 290(1):81–88, 2019.
  * [4] C. Doersch et al.  Multi-task self-supervised visual learning.  In ICCV, 2017.
  * [5] S. Gidaris et al.  Unsupervised representation learning by predicting image rotations.  In ICLR, 2018.
  * [6] C. Doersch et al.  Unsupervised visual representation learning by context prediction.  In ICCV, 2015.
  * [7] R. Zhang et al.  Colorful image colorization.  In ECCV, 2016.
  * [8] D. Pathak et al.  Context encoders: Feature learning by inpainting.  In CVPR, 2016.
  * [9] A. Jamaludin et al.  Self-supervised learning for spinal MRIs.  In MICCAI DLMIA Workshop, 2017.
  * [10] T. Ross et al.  Exploiting the potential of unlabeled endoscopic video data with self-supervised learning.  Int J Comput Assist Radiol Surg, 13(6):925–933, 2018.
  * [11] N. Tajbakhsh et al.  Surrogate supervision for medical image analysis: Effective deep learning from limited quantities of labeled data.  In ISBI, 2019.
  * [12] O. Ronneberger et al.  U-Net: convolutional networks for biomedical image segmentation.  In MICCAI, 2015.


# The Mind’s Eye: Visualizing Class-Agnostic Features of CNNs

###### Abstract

Visual interpretability of Convolutional Neural Networks (CNNs) has gained
significant popularity because of the great challenges that CNN complexity
imposes to understanding their inner workings. Although many techniques have
been proposed to visualize class features of CNNs, most of them do not provide
a correspondence between inputs and the extracted features in specific layers.
This prevents the discovery of stimuli that each layer responds better to. We
propose an approach to visually interpret CNN features given a set of images
by creating corresponding images that depict the most informative features of
a specific layer. Exploring features in this class-agnostic manner allows for
a greater focus on the feature extractor of CNNs. Our method uses a dual-
objective activation maximization and distance minimization loss, without
requiring a generator network nor modifications to the original model. This
limits the number of FLOPs to that of the original network. We demonstrate the
visualization quality on widely-used architectures.111Code is available at
https://git.io/JL9Wg
and our demo video: https://youtu.be/Au3jaUdnPKM

Index Terms—  Feature visualization, CNN explainability, convolutional
features

##  1 Introduction

Fig. 1: Top 50 extracted features. ResNet-50 [1] was used with features
visualized from ‘layer6.5.conv1’.

Deep learning architectures have achieved substantial breakthroughs in
comparison to hand-coded feature extractors, for a wide variety of image and
video tasks. While their performance is high, their interpretability remains
limited.

For this reason, methods on visualizing features of CNNs have received
significant research interest over the last few years. We identify two main
approaches. The first is to consider the image regions that networks find
informative [2, 3, 4]. This approach allows for the selection of salient
regions, but it does not provide a description for the feature’s appearance.
The second set of methods addresses this shortcoming by explicitly creating
visualizations that activate features of specific classes [5, 6].

While both classes of approach have shown great promise to establish robust
visual explanations for CNN features, one key aspect of deep learning method
has not been explicitly addressed. Although early extracted features can be
easily interpreted based on edges and textural patterns, features of deeper
layers are significantly more complicated while they mostly do not correspond
to singular concepts. We refer to features with polysemantic interpretations
as entangled [7].

To visualize neurons that encapsulate entangled features, we propose a multi-
objective method that creates class-agnostic visual representations of image
features. As a result, we can show the top $k$ descriptive features in sets of
images. Our approach further addresses the common problem of different
interpretations through perturbations [8], as the method is not constrained by
visualizing features in a per-class fashion. An example appears in Figure 1.

Our contributions are as follows:

  * •

We propose a class-agnostic method for visual explanations of CNN features.
Our method uses a multi-objective loss based on activation maximization that
optimizes an input image through the excitation of a user-defined number of
layer neurons.

  * •

We design an axiomatic distance objective to address entangled features by
minimizing the distance between produced image activations and the averaged
target activations of real images.

  * •

We test our approach on different layers and neurons of popular CNNs and show
that our visualizations can uncover interesting features in sets of images.

The remainder of the paper is structured as follows. We first summarize
related work on visual interpretations. We detail our method in Section 3. We
demonstrate the produced feature visualizations in Section 4 and conclude in
Section 5.

##  2 Related Work

Recent works have argued the importance of interpretability in CNNs [9, 10]
and how it can further lead to CNN improvements [11]. However, creating
methods that capture the inner workings of CNNs has been proven challenging.

One widely used method to visualize CNN features is to maximize neurons that
correspond to a specific class [5]. To maximize a class neuron, the input is
composed of trainable parameters that are updated based on the gradients. As
the sole consideration of class activations does not give an intuitive
representation of the features that correspond to classes, Zeiler and Fergus
[6] proposed a de-convolution approach that aims at approximating layer
features. Later works of Simonyan et al. [12] and Nguyen et al. [13] have
shown how the exclusive use of an activation maximization objective can lead
to the creation of unrealistic images, since the space of possible images and
patterns that can be produced and that are close to extracted patterns, is
extremely vast. This motivated the exploration of regularization techniques
aimed at constraining the space of possible visual representations. Some of
these techniques include the use of Gaussian filters [14, 15] during the image
optimization process, jitter effect [16], and creating center-biased gradient
masks [17].

Other approaches to improve the realism in images consider using a separate
network that is capable of synthesizing feature visualizations [18, 19] based
on adversarial training similar to that of generative networks. Although the
visual quality in generative models is higher, generators lack in terms of
representing the causality of learned features [20], as they include an
additional factor of ambiguity through the generator sub-network.

To address the problems associated with activation-maximization, we propose a
method that is inspired by recent distance-minimization-based generative
networks. Our method optimizes an image-based dot product of activations from
generated and real-world images while additionally decreasing their activation
distances within the feature space.

Fig. 2: Feature space distances. Point corresponds to top-10 embeddings
($out_{n},\;n\in N$) with center $out_{c}$.

##  3 Visualization methodology

We use a dot-product activation maximization with an additional distance
minimization regression objective to optimize a trainable input image to
represent the most informative features for a specific layer in the model.

###  3.1 Multi-faced clustered neuron selection

We use a multi-faceted technique similar to the one proposed by Nguyen et al.
[17] to optimize the creation of the initial images ($\overline{img}$) as well
as the target activations ($t_{i}$) of layer $i$ used for regressions. We
define $C$ target facets. We use real images and perform a normal forward-pass
in the network until layer $i$. We then reduce the original channel
dimensionality of activations $a_{i}$ in layer $i$ through PCA [21] and t-SNE
[22] to create 2-dimensional embeddings $out$ that are then clustered into $C$
clusters with k-Means [23]. Instead of using the average within each of the
$C$ clusters $c$ ($1\leq c\leq C$), we consider the $N=10$ 2D embeddings
closest to each cluster center ($out_{c}$) to allow for a better
correspondence to feature activations that are characteristic for cluster $c$.
Based on the euclidean distance between the 2D embeddings and the cluster
center ($edist(out_{c},out_{n})\forall n\in N$) we create a weight penalty
$w_{c,n}$. The weight corresponds to the effect of each activation ($a_{i,n}$)
that is exponentially counter-equal to the euclidean distance between its 2D
embedding $out_{n}$ and the cluster center $out_{c}$. This is summarized in
Eq. 1 with the image ($\overline{img}_{c}$) initialization and the discovered
target activations ($t_{i,c}$). We include a constant ($\gamma=1e^{-5}$) for
numeric stability.

$\begin{split}{w}_{c,n}=\frac{e^{edist(out_{c},out_{n})^{-1}+\gamma}}{\sum\limits_{n\in
N}e^{edist(out_{c},out_{n})^{-1}+\gamma}}\qquad\qquad\\\
\overline{img}_{c}=\sum\limits_{n\in
N}img_{n}*w_{c,n}\;\;\&\;\;t_{i,c}=\sum\limits_{n\in
N}a_{i,n}*w_{c,n}\end{split}\vspace{-1.5mm}$ (1)

The effect is visualized in Figure 2, where the distance between each 2D
embedding $out_{n}$ and the corresponding cluster center $out_{c}$ determines
the contribution of each activation $a_{i,n}$ to the final target activation
$t_{i,c}$.

Fig. 3: Image optimization iterations. The target image from ImageNet’s [24]
padlock class is the image closest to the cluster center. The top 30 features
of wide-ResNet101’s [25] ‘layer6.12.conv3’ are used.

###  3.2 Objective formalization

We define our loss function based on two additional auxiliary objectives to
improve the feature clarity while simultaneously reduce the effects of feature
entanglement. To visualize specific features, we select the top $k$ features
of each target activation $t_{i,c}$ of each cluster $c$. This creates an
averaged overview over the most informative top $k$ features that should be
visualized.

We define a dot-product activation maximization loss ($DM$) as the channel-
wise dot product of the produced activation maps $a_{i}$ and the discovered
target layer activations $t_{i,c}$. This calculation is performed for all top
$k$ channels:

$DM=\sum\limits_{j\in k}a_{i,j}*t_{i,c,j}$ (2)

Through the maximization of the dot-product of the produced activations and
the target layer activations, we create a path towards meaningful features for
the gradients. However, this also corresponds to larger feature entanglement
as the span of possible gradient directions to provide positive improvements
can be extensive.

To address this issue, we include a second objective: a multi-dimensional
distance minimization between the produced activation maps $a_{i}$ and the
target activations $t_{i}$. For the distance loss function, we use a super-set
of distance methods as proposed by Barron [26] (denoted as $AD$). Because of
the heteroscedastic nature of the produced distances, using trainable
parameters ($r,b$), shown in Eq. 3, can better fit the produced multivariate
distances.

$\begin{split}AD=\begin{cases}\frac{1}{2}(\frac{mdist}{r})^{2},\;for\;b=2\\\
log(\frac{1}{2}(\frac{mdist}{r})^{2}+1)\;for\;b=0\\\
1-exp(-\frac{1}{2}(\frac{mdist}{r})^{2})\;for\;b\rightarrow-\infty\\\
\frac{|b-2|}{b}\left(\left(\frac{(\frac{mdist}{r})^{2}}{|b-2|}+1\right)^{b/2}+1\right)\;otherwise\end{cases}\\\
where,\>\>mdist=\left(\sum\limits_{j\in
k}||a_{i,j}-t_{i,c,j}||\right)\qquad\quad\end{split}$ (3)

Finally, we synthesize our loss function $\mathcal{L}$ from Eqs. 2 and 3 by
including a penalty for the activations of the previous layers ($a_{i-1}$)
with channel size $D$, through $L_{1}$ regularization. The scaling value
$\lambda$ has an initial value of $1e^{-3}$ which is linearly decreased to
$1e^{-4}$ during training. The final loss is:

$\mathcal{L}=AD-DM+\lambda*\sum\limits_{j\in D}||a_{i-1,j}||\vspace{-1.5mm}$
(4)

###  3.3 Parameterization setup

We include parameterized noise functions within our workflow as constraints
for the high-frequency gradients during back-propagation. This allows for the
minimization of possible feature imbalances [15] and improves the final visual
representation quality of CNN features. We include popular techniques used in
visualization methods such as image blurring with low-variance Gaussian
kernels [14, 15] in combination with a denoising split Bregman algorithm [27].
We additionally include a center-based gradient mask [17] to limit feature
cluttering as well as to limit feature duplication during training. Learning
rate $lr$, standard deviations of the Gaussian blur, denoising, and center-
based gradient mask ($\sigma$) are decreased linearly over each iteration $t$
($1\leq t\leq T$), based on user-set starting ($start$) and final ($end$)
values:

$lr,\sigma=\frac{((end-start)*t)+(start*T)}{T-1}\vspace{-1.5mm}$ (5)

We present iterative changes in Figure 3 where an image is optimized to
visualize features that correspond to a target.

Fig. 4: Cross-layer visualizations of top-100 features. Images were sampled
from ImageNet’s ‘screw’ class. The targets include variations in orientation,
number of objects and shape. We use ResNeXt-101 [28] with the corresponding
optimization layers being displayed at the top of each column.

##  4 Visualizations

We demonstrate in Figure 1 three cases of the same object where the produced
feature visualizations present some degree of variation. Although the target
images show the same object (from class ‘banana’ in ImageNet), the number of
objects present may correspond to differences in terms of the most dominant
feature activations. We note however, that differences in the image targets do
not present detrimental effects in the overall feature combination.

In Figure 3, we present how features are visualized over time. This shows the
variations in regions of the image being optimized and how feature-inclusive
regions change over time. This demonstrates the ability of a specific
architecture to combine only a certain number of features in order to
encapsulate the general appearance of an object. For example, for images from
the ‘padlock’ class, the 30 most highly activated features can be seen as
sufficient for describing the overall appearance of the object. This can aid
in determining the images and classes that are easier or harder for the
network to extract meaningful features from.

Lastly, we show in Figure 4 the features that are extracted among different
layers of the same network. Using images from ImageNet’s ‘screw’ class, we
sample target sets of images varying in terms of their general orientation
with objects being perpendicular or at an angle, as well as differences in
screw head types and number of objects. Visualizations of later layers (7.2
and 6.20) include the metallic look of the object as well as screw body
patterns and detail in terms of the screw head type (circular and hexagonal).
Such specific features however are absent in earlier layers as their feature
extractors seem to correspond to basic features in terms of the general object
appearance with the number of information-rich entangled features being lower
from that of later layers.

##  5 Conclusions

We have proposed a novel feature visualization method aimed at providing
visual explanations for the top features extracted by CNN layers. Images are
optimized based on a dual-objective loss which includes the dot-product
activation maximization and the distance minimization between the produced and
target layer activations. To reduce feature overlap and to improve the overall
visualization clarity, we apply blurring and de-blurring filters as well as a
gradient mask to the generated images during optimization.

We present corresponding CNN features that are associated with specific
classes and images. Based on this, we believe that the produced visual
explanations can improve the in-depth understanding of trained CNNs and the
features with polysemantic interpretations that are associated with a specific
image or set of images.

## References

  * [1] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun,  “Deep residual learning for image recognition,”  in Computer Vision and Pattern Recognition (CVPR), 2016, pp. 770–778.
  * [2] Ruth C Fong and Andrea Vedaldi,  “Interpretable explanations of black boxes by meaningful perturbation,”  in International Conference on Computer Vision (ICCV), 2017, pp. 3429–3437.
  * [3] Grégoire Montavon, Sebastian Lapuschkin, Alexander Binder, Wojciech Samek, and Klaus-Robert Müller,  “Explaining nonlinear classification decisions with deep Taylor decomposition,”  Pattern Recognition, vol. 65, pp. 211–222, 2017.
  * [4] Ramprasaath R Selvaraju, Michael Cogswell, Abhishek Das, Ramakrishna Vedantam, Devi Parikh, Dhruv Batra, et al.,  “Grad-CAM: Visual explanations from deep networks via gradient-based localization.,”  in Internation Conference on Computer Vision (ICCV), 2017, pp. 618–626.
  * [5] Dumitru Erhan, Yoshua Bengio, Aaron Courville, and Pascal Vincent,  “Visualizing higher-layer features of a deep network,”  Tech. Rep. 1341-3, University of Montreal, 2009.
  * [6] Matthew D Zeiler and Rob Fergus,  “Visualizing and understanding convolutional networks,”  in European conference on computer vision (ECCV). Springer, 2014, pp. 818–833.
  * [7] Jesse Mu and Jacob Andreas,  “Compositional explanations of neurons,”  Advances in Neural Information Processing Systems (NeurIPS), vol. 33, 2020.
  * [8] Amirata Ghorbani, Abubakar Abid, and James Zou,  “Interpretation of neural networks is fragile,”  in Conference on Artificial Intelligence (AAAI), 2019, vol. 33, pp. 3681–3688.
  * [9] David Bau, Bolei Zhou, Aditya Khosla, Aude Oliva, and Antonio Torralba,  “Network dissection: Quantifying interpretability of deep visual representations,”  in Conference on Computer Vision and Pattern Recognition (CVPR), 2017, pp. 6541–6549.
  * [10] Dong Huk Park, Lisa Anne Hendricks, Zeynep Akata, Anna Rohrbach, Bernt Schiele, Trevor Darrell, and Marcus Rohrbach,  “Multimodal explanations: Justifying decisions and pointing to the evidence,”  in Conference on Computer Vision and Pattern Recognition (CVPR), 2018, pp. 8779–8788.
  * [11] Lisa Anne Hendricks, Ronghang Hu, Trevor Darrell, and Zeynep Akata,  “Grounding visual explanations,”  in European Conference on Computer Vision (ECCV), 2018, pp. 269–286.
  * [12] Karen Simonyan, Andrea Vedaldi, and Andrew Zisserman,  “Deep inside convolutional networks: Visualising image classification models and saliency maps,”  arXiv preprint arXiv:1312.6034, 2013.
  * [13] Anh Nguyen, Jason Yosinski, and Jeff Clune,  “Deep neural networks are easily fooled: High confidence predictions for unrecognizable images,”  in Conference on Computer Vision and Pattern Recognition (CVPR), 2015, pp. 427–436.
  * [14] Jason Yosinski, Jeff Clune, Anh Nguyen, Thomas Fuchs, and Hod Lipson,  “Understanding neural networks through deep visualization,”  in International Conference of Machine Learning Workshops (ICML), 2015.
  * [15] Feng Wang, Haijun Liu, and Jian Cheng,  “Visualizing deep neural network by alternately image blurring and deblurring,”  Neural Networks, vol. 97, pp. 162–172, 2018.
  * [16] Alexander Mordvintsev, Christopher Olah, and Mike Tyka,  “Inceptionism: Going deeper into neural networks,” 2015.
  * [17] Anh Nguyen, Jason Yosinski, and Jeff Clune,  “Multifaceted feature visualization: Uncovering the different types of features learned by each neuron in deep neural networks,”  in International Conference of Machine Learning Workshops (ICML), 2016.
  * [18] Christian F Baumgartner, Lisa M Koch, Kerem Can Tezcan, Jia Xi Ang, and Ender Konukoglu,  “Visual feature attribution using Wasserstein GANs,”  in Conference on Computer Vision and Pattern Recognition (CVPR), 2018, pp. 8309–8319.
  * [19] Anh Nguyen, Alexey Dosovitskiy, Jason Yosinski, Thomas Brox, and Jeff Clune,  “Synthesizing the preferred inputs for neurons in neural networks via deep generator networks,”  in Advances in Neural Information Processing Systems (NeuIPS), 2016, pp. 3387–3395.
  * [20] Andreas Holzinger, Georg Langs, Helmut Denk, Kurt Zatloukal, and Heimo Müller,  “Causability and explainability of artificial intelligence in medicine,”  Wiley Interdisciplinary Reviews: Data Mining and Knowledge Discovery, vol. 9, no. 4, pp. e1312, 2019.
  * [21] IT Jolliffe,  “Principal component analysis,”  Technometrics, vol. 45, no. 3, pp. 276, 2003.
  * [22] Laurens van der Maaten and Geoffrey Hinton,  “Visualizing data using t-SNE,”  Journal of machine learning research, vol. 9, no. 86, pp. 2579–2605, 2008.
  * [23] Stuart Lloyd,  “Least squares quantization in PCM,”  IEEE Transactions on Information Theory, vol. 28, no. 2, pp. 129–137, 1982.
  * [24] Olga Russakovsky, Jia Deng, Hao Su, Jonathan Krause, Sanjeev Satheesh, Sean Ma, Zhiheng Huang, Andrej Karpathy, Aditya Khosla, Michael Bernstein, Alexander C. Berg, and Li Fei-Fei,  “ImageNet Large Scale Visual Recognition Challenge,”  International Journal of Computer Vision (IJCV), vol. 115, no. 3, pp. 211–252, 2015.
  * [25] Sergey Zagoruyko and Nikos Komodakis,  “Wide residual networks,”  in British Machine Vision Conference (BMVC), 2016, pp. 87.1–87.12.
  * [26] Jonathan T Barron,  “A general and adaptive robust loss function,”  in Conference on Computer Vision and Pattern Recognition (CVPR), 2019, pp. 4331–4339.
  * [27] Pascal Getreuer,  “Rudin-Osher-Fatemi total variation denoising using split Bregman,”  Image Processing On Line, vol. 2, pp. 74–95, 2012.
  * [28] Saining Xie, Ross Girshick, Piotr Dollár, Zhuowen Tu, and Kaiming He,  “Aggregated residual transformations for deep neural networks,”  in Computer Vision and Pattern Recognition (CVPR), 2017, pp. 5987–5995.