– Under Review 2021 Full Paper – MIDL 2021 submission Walter Hugo Lopez Pinaya
1<EMAIL_ADDRESS>
Petru-Daniel Tudosiu 1<EMAIL_ADDRESS>
Robert Gray 2<EMAIL_ADDRESS>
Geraint Rees 3<EMAIL_ADDRESS>
Parashkev Nachev2<EMAIL_ADDRESS>
Sébastien Ourselin 1<EMAIL_ADDRESS>
M. Jorge Cardoso1<EMAIL_ADDRESS>
1 Department of Biomedical Engineering, School of Biomedical Engineering &
Imaging Sciences, King’s College London, London, UK and 2 Institute of
Neurology, University College London, London, UK and 3 Institute of Cognitive
Neuroscience, University College London, London, UK

# Unsupervised Brain Anomaly Detection and Segmentation with Transformers

###### Abstract

Pathological brain appearances may be so heterogeneous as to be intelligible
only as anomalies, defined by their deviation from normality rather than any
specific pathological characteristic. Amongst the hardest tasks in medical
imaging, detecting such anomalies requires models of the normal brain that
combine compactness with the expressivity of the complex, long-range
interactions that characterise its structural organisation. These are
requirements transformers have arguably greater potential to satisfy than
other current candidate architectures, but their application has been
inhibited by their demands on data and computational resource. Here we combine
the latent representation of vector quantised variational autoencoders with an
ensemble of autoregressive transformers to enable unsupervised anomaly
detection and segmentation defined by deviation from healthy brain imaging
data, achievable at low computational cost, within relative modest data
regimes. We compare our method to current state-of-the-art approaches across a
series of experiments involving synthetic and real pathological lesions. On
real lesions, we train our models on 15,000 radiologically normal participants
from UK Biobank, and evaluate performance on four different brain MR datasets
with small vessel disease, demyelinating lesions, and tumours. We demonstrate
superior anomaly detection performance both image-wise and pixel-wise,
achievable without post-processing. These results draw attention to the
potential of transformers in this most challenging of imaging tasks.

###### keywords:

Transformer, Unsupervised Anomaly Segmentation, Anomaly Detection,
Neuroimaging, Vector Quantized Variational Autoencoder

††editors: Under Review for MIDL 2021

##  1 Introduction

Transformers have revolutionised language modelling, becoming the de-facto
network architecture for language tasks [Radford et al. (2018), Radford et al.
(2019), Vaswani et al. (2017)]. They rely on attention mechanisms to capture
the sequential nature of an input sequence, dispensing with recurrence and
convolutions entirely. This mechanism allows the modelling of dependencies of
the inputs without regard to their distance, enabling the acquisition of
complex long-range relationships. Since the approach generalises to any
sequentially organised data, applications in other areas such as computer
vision are increasingly seen, with impressive results in image classification
[Chen et al. (2020a), Dosovitskiy et al. (2020)] and image synthesis [Child et
al. (2019), Esser et al. (2020), Jun et al. (2020)]. The power to absorb
relationships varying widely in their distance makes transformers of potential
value in the arguably the hardest of neuroimaging tasks: anomaly detection.

The detection and segmentation of lesions in neuroimaging support an array of
clinical tasks, including diagnosis, prognosis, treatment selection and
mechanistic inference. However, the fine characterisation of these lesions
requires an accurate segmentation which is generally both ill-defined and
dependent on human expertise [Kamnitsas et al. (2017)]. Manual segmentation is
expensive and time-consuming to obtain, greatly limiting clinical application,
and the scale and inclusivity of available labelled data. Qualitative,
informal descriptions or reduced measurements are often used instead in
clinical routine [Porz et al. (2014), Yuh et al. (2012)]. For this reason, the
development of accurate computer-aided automatic segmentation methods has
become a major endeavour in medical image research [Menze et al. (2014)]. Most
methods, however, depend on an explicitly defined target class, and are
sensitive to the scale and quality of available labelled data, a sensitivity
amplified by the many sources of complex variability encountered in clinical
neuroimaging. Under real-world distributional shift, such models behave
unpredictably, limiting clinical utility.

In recent years, many machine learning algorithms have been proposed for
automatic anomaly detection. To overcome the necessity of expensive labelled
data, unsupervised methods have emerged as promising tools to detect arbitrary
pathologies [Baur et al. (2018), Baur et al. (2020b), Chen et al. (2020b),
Pawlowski et al. (2018)], relying mainly on deep generative models of normal
data to derive a probability density estimate of the input data defined by the
landscape of normality. Pathological features then register as deviations from
normality, avoiding the necessity for either labels or anomalous examples in
training. The state of the art is currently held by variational autoencoder
(VAE)-based methods [Baur et al. (2020a)] which try to reconstruct a test
image as the nearest sample on the learnt normal manifold, using the
reconstruction error to quantify the degree and spatial distribution of any
anomaly. This approach’s success is limited by the fidelity of reconstructions
from most VAE architectures [Dumoulin et al. (2016)], and unwanted
reconstructions of pathological features not present in the training data,
suggesting a failure of the model to internalise complex relationships between
remote imaging features.

In an effort to address these problems, we propose a method for unsupervised
anomaly detection and segmentation using transformers, where we learn the
distribution of brain imaging data with an ensemble of Performers [Choromanski
et al. (2020)]. We create and evaluate a robust method and compare its
performance on synthetic and real datasets with recent state-of-the-art
unsupervised methods.

##  2 Background

The core of the proposed anomaly detector is a highly expressive transformer
that learns the probability density function of 2D images of healthy brains.
This requires us to express each image’s contents as a sequence of
observations on which transformers-like models can operate. Owing to the size
and complexity of brain imaging data, instead of learning the distributions on
individual pixels directly, we use the compact latent discrete representation
of a vector quantised variational autoencoder (VQ-VAE) [Razavi et al. (2019),
van den Oord et al. (2017)]. This approach allows us to compress the input
data into a spatially smaller quantised latent image, thus reducing the
computational requirements and sequence length, making transformers feasible
in neuroimaging applications.

###  2.1 Vector quantized variational autoencoder

The VQ-VAE [Razavi et al. (2019), van den Oord et al. (2017)] is a model that
learns discrete representations of images. It comprises an encoder $E$ that
maps observations $x\in\mathbb{R}^{H\times W}$ onto a latent embedding space
$\hat{z}\in\mathbb{R}^{h\times w\times n_{z}}$, where $n_{z}$ is the
dimensionality of each latent embedding vector. Then, an element-wise
quantization is performed for each spatial code
$\hat{z}_{ij}\in\mathbb{R}^{n_{z}}$ onto its nearest vector $e_{k}$ , $k\in
1,...K$ from a codebook, where $K$ denotes the vocabulary size of the
codebook. This codebook is learnt jointly with the other model parameters. A
decoder $G$ reconstructs the observations $\hat{x}\in\mathbb{R}^{H\times W}$
from the quantized latent space. We obtain the latent discrete representation
$z_{q}\in\mathbb{R}^{h\times w}$ by replacing each code by its index $k$ from
the codebook.

###  2.2 Transformers

After training the VQ-VAE, we can learn the probability density function of
the latent representation of healthy brain images using transformers. The
defining characteristic of a transformer is that it relies on attention
mechanisms to capture the interactions between inputs, regardless of their
relative position to one another. Each layer of the transformer consists of a
(self-)attention mechanism described by mapping an intermediate representation
with three vectors: query, key, and value vectors. Since the output of this
attention mechanism relies on the inner products between all elements in the
sequence, its computational costs scale quadratically with the sequence
length. Several “efficient transformers” have recently been proposed to reduce
this computational requirement [Tay et al. (2020)]. Our study uses the
Performer, a model which uses an efficient (linear) generalized attention
framework implemented by the FAVOR+ algorithm [Choromanski et al. (2020)].
This framework provides a scalable estimation of attention mechanisms
expressed by random feature map decompositions, making transformers feasible
for longer sequences, of the size needed for neuroimaging data.

To model brain images, after assuming an arbitrary ordering to transform the
latent discrete variables $z_{q}$ into a 1D sequence $s$, we can train the
transformer to maximize the training data’s log-likelihood in an
autoregressive fashion - similar to training performed on language modelling
task. This way, the transformer learns the distribution of codebook indices
for a position $i$ given all previous values $p(s_{i})=p(s_{i}|s_{<i})$.

##  3 Proposed Method

###  3.1 Anomaly Segmentation

To segment an anomaly in a previously unseen, test image, first, we obtain the
latent discrete representation $z_{q}$ from the VQ-VAE model. Next, we reshape
$z_{q}$ into a sequence $s$, and we use the autoregressive transformer to
obtain the likelihood of each latent variable value $p(s_{i})$. These
likelihood values indicate which latent variable has a low probability of
occurring in normal data. Using an arbitrary threshold (we empirically
determined a threshold of 0.005 on a holdout set), we then can select indices
with the lowest likelihood values and create a “resampling mask” indicating
which latent variables are abnormal and should be corrected to produce a
“healed” version of the image. We then replace these abnormal values with
values sampled by the transformer. This approach attenuates the influence of
the anomalies by replacing them by values that conform to the healthy
distribution in the discrete latent space. This in-painted latent space can
then reconstruct $\hat{x}^{\prime}$ without the anomalies, in “healed” form.
Finally, we obtain the pixel-wise residuals from the difference
$|x-\hat{x}^{\prime}|$. The anomaly in the new sample is finally segmented by
thresholding the highest residuals values.

###  3.2 Spatial information from the latent space

Autoencoders are known for creating blurry reconstructions [Dumoulin et al.
(2016)]. This can cause residuals with high values in areas of the image with
high frequencies. To avoid these areas being mislabelled as anomalous, we used
the spatial information present in the in-painted “resampling mask” of the
latent space.

The resampling mask indicates the spatial location of the latent values with
anomalies according to the transformer model. Since our VQ-VAE is relatively
shallow, the latent space preserves most of the spatial information of the
input data. Based on this, we used the resampling mask to avoid mislabelling
of high-frequency regions. First, we upscaled the resampling mask from the
latent space resolution to the input data resolution. Next, we multiply the
residuals with the mask. This approach cleans areas of the residuals that were
not specified as anomalous by our transformer but where the reconstructions
might be largely due to lack of VQ-VAE capacity.

###  3.3 Multiple views of the latent space through reordering

Recent studies have reported some limitations of likelihood models, such as
our autoregressive model, in identifying out-of-distribution samples
[Nalisnick et al. (2018)]. Inspired by Choi et al. (2018), we made the
detection of the anomalies more robust by using an ensemble of models. Using
the same VQ-VAE model, we trained an ensemble of autoregressive transformers.
However, unlike Choi et al. (2018), each one of our transformers uses a
different reordering of the 2D latent image to create a sequence. This compels
each transformer to use a different context of the latent image when
predicting the likelihood of an element.

In our study, we focused on the raster scan class ordering. We obtain
different orderings by reflecting the image horizontally, vertically, and both
ways at the same time. We also define our orderings in images rotated 90
degrees, generating 8 different orderings from a single latent representation.
Each resampled latent representation is independently reconstructed, i.e. each
model independently creates a residuals map. We use the mean residual to
segment the anomalies.

###  3.4 Image-wise Anomaly Detection

So far, the proposed methodology has been focusing on segmenting
abnormalities. However, transformers can also be used to perform image-wise
anomaly detection, i.e. detecting if an abnormality exists somewhere in the
image. To do so, we use the likelihood predicted by the transformers. Like the
segmentation approach, first, we obtain the 1D latent representation $s$.
Then, we use the transformers to obtain the likelihood $p(s)$ of each latent
variable. To obtain the log-likelihood image-wise, we compute
$logp(x)=logp(s)=\sum_{i}logp(s_{i})$. Finally, we combined the predicted log-
likelihood of each transformer (per orientation/ordering) by computing the
mean value.

##  4 Experiments and Results

###  4.1 Experiment #1 – Anomaly Segmentation on Synthetic Data

To assess anomaly segmentation performance on synthetic data, we utilized a
subsample of the MedNIST dataset111Available at https://github.com/Project-
MONAI/tutorials, where we used the images of the “HeadCT” category to train
our models. Our test set comprised 100 images contaminated with sprites
Matthey et al. (2017), thus producing ground-truth abnormality masks. We
measure the performance using the best achievable DICE-score ($\lceil
DICE\rceil$), which constitutes a theoretical upper-bound to a model’s
segmentation performance and is obtained via a greedy search for the residual
threshold which yields the highest DICE-score on the test set. We compared our
results against state-of-the-art autoencoder models based on the architectures
proposed in the Baur et al. (2020a) comparison study, assessed in the same
manner. We also performed an ablation study of the proposed method
demonstrating the values of each contribution.

tab:syntseg

Table 1: Methods on anomaly segmentation using the synthetic dataset. The
performance is measured with best achievable DICE-score on the test set.

Method | $\lceil DICE\rceil$
---|---
AE (Dense) (Baur et al., 2020a) | 0.213
AE (Spatial) (Baur et al., 2020a) | 0.165
VAE (Dense) (Baur et al., 2020a) | 0.533
VQ-VAE (van den Oord et al., 2017) | 0.457
VQ-VAE + Transformer (Ours) | 0.675
VQ-VAE + Transformer + Masked Residuals (Ours) | 0.768
VQ-VAE + Transformer + Masked Residuals + Different Orderings (Ours) | 0.895

As presented in tab:syntseg, the models based only on autoencoders had a
highest $\lceil DICE\rceil$ of 0.533 (VAE). We observed an improvement in
performance when using the transformer to in-paint the latent space, changing
the VQ-VAE only performance from 0.457 to 0.675. The spatial information in
the resampling mask also contributed by attenuating the impact of the blurry
reconstructions (Figure LABEL:fig:exp1), achieving a 0.768 score. Finally, the
variability of the generative models with different orderings gave another
boost in performance (for 8 different raster ordering models $\lceil
DICE\rceil$=0.895).

fig:exp1

Figure 1: Residual maps on the synthetic dataset from the variational
autoencoder and different steps of our approach.

###  4.2 Experiment #2 – Image-wise Anomaly Detection on Synthetic Data

Next, we evaluated our method to detect anomalous (out-of-distribution - OOD)
images, again on a synthetic setting. Using the same models trained from
Experiment #1, we obtained the mean log-likelihood for each image. We used
1,000 images from the HeadCT class as the in-distribution test set, the 100
HeadCT images contaminated by sprite anomalies as the near OOD set, and 1.000
images of each other MedNIST classes as the far OOD set (see the appendix for
details). We use the area under the receiver operating characteristic curve
(AUROC) as performance metric, with in-distribution test set and out-of-
distribution being the labels. This way, we have a threshold-independent
evaluation metric. We also measure the area under the precision recall curve
(AUCPRC), where it provides a meaningful measure for detection performance in
the presence of heavy class-imbalance. Finally, we also computed the false
positive rate of anomalous examples when the true positive rate of in-
distribution examples is at 95% (FPR95), 99% (FPR99) and 99.9% (FPR999).

tab:syntdet

Table 2: Performance of the methods on image-wise anomaly detection using the
synthetic dataset.

| AUCROC |  AUCPR In |  AUCPR Out | FPR95 | FPR99 | FPR999
---|---|---|---|---|---|---
vs. far OOD classes |  |  |  |  |  |
VAE (Dense) (Baur et al., 2020a) | 0.298 | 0.855 | 0.060 | 0.986 | 0.996 | 0.996
Our approach | 1.000 | 1.000 | 1.000 | 0.000 | 0.001 | 0.004
Our approach with general purpose VQ-VAE | 1.000 | 1.000 | 1.000 | 0.000 | 0.000 | 0.000
vs. near OOD classes |  |  |  |  |  |
VAE (Dense) (Baur et al., 2020a) | 0.111 | 0.094 | 0.672 | 1.000 | 1.000 | 1.000
Our approach | 0.921 | 0.988 | 0.707 | 0.409 | 0.885 | 0.885
Our approach with general purpose VQ-VAE | 0.932 | 0.990 | 0.721 | 0.482 | 0.882 | 0.882

tab:syntdet shows that our transformer-based method achieved an AUCROC of
0.921 and 1.000 for near OOD and far OOD, respectively. This is a improvement
compared to a method based on the error of reconstruction obtained from a VAE
model. We also evaluated our method with a VQ-VAE trained to reconstruct all
the categories from the MedNIST dataset (“general purpose VQ-VAE”) and the
ensemble of transformers trained on HeadCT images only. In this configuration,
we try to mitigate the influence of the encoder in the anomaly detection. This
approach would reduce the ability of the encoder to map an OOD image to a
familiar in-distribution latent representation, which could possibly affect
the transformer performance. This new configuration achieves a slight better
performance (AUCROC=0.932 for near OOD and AUCROC=1.000 for far OOD).

###  4.3 Experiment #3 – Anomaly Segmentation on Real Neuroimaging Data

To evaluate our method’s performance on real world lesion data, we used the
FLAIR images from the UK Biobank (UKB) dataset Sudlow et al. (2015). We
selected the 15,000 subjects, and their respective FLAIR images, with the
lowest white matter hyperintensities volume, as provided by UKB, to train our
models, as these subjects were the most radiologically normal. Then, we used
18,318 subjects from the remaining UKB dataset to evaluate our method to
detect white matter hyperintensities (WMH).

In order to test for model generalisability, we also evaluated our method on
three other datasets that also had FLAIR imaging data: the Multiple Sclerosis
dataset from the University Hospital of Ljubljana (MSLUB) dataset Lesjak et
al. (2018), which contains multiple sclerosis lesions; the White Matter
Hyperintensities Segmentation Challenge (WMH) dataset Kuijf et al. (2019); and
the Multimodal Brain Tumor Image Segmentation Benchmark (BRATS) dataset Bakas
et al. (2017, 2018); Menze et al. (2014) that contain tumours (more
information about the datasets and pre-processing on appendix A).

tab:realseg

Table 3: Results of the anomaly segmentation using real lesion data. We
compared our models against the state-of-the-art autoencoder models based on
the architecture proposed in Baur et al. (2020a). We measured the performance
using the theoretically best possible DICE-score ($\lceil DICE\rceil$) on each
dataset.

UKB Dataset | $\lceil DICE\rceil$
---|---
AE (Dense) (Baur et al., 2020a) | 0.016
AE (Spatial) (Baur et al., 2020a) | 0.054
VAE (Dense) (Baur et al., 2020a) | 0.016
VQ-VAE (van den Oord et al., 2017) | 0.028
VQ-VAE + Transformer + Masked Residuals + Different Orderings (Ours) | 0.232
MSLUB Dataset |
AE (Dense) (Baur et al., 2020a) | 0.041
AE (Spatial) (Baur et al., 2020a) | 0.061
VAE (Dense) (Baur et al., 2020a) | 0.039
VQ-VAE (van den Oord et al., 2017) | 0.040
VQ-VAE + Transformer + Masked Residuals + Different Orderings (Ours) | 0.378
BRATS Dataset |
AE (Dense) (Baur et al., 2020a) | 0.276
AE (Spatial) (Baur et al., 2020a) | 0.531
VAE (Dense) (Baur et al., 2020a) | 0.294
VQ-VAE (van den Oord et al., 2017) | 0.331
VQ-VAE + Transformer + Masked Residuals + Different Orderings (Ours) | 0.759
WMH Dataset |
AE (Dense) (Baur et al., 2020a) | 0.073
AE (Spatial) (Baur et al., 2020a) | 0.150
VAE (Dense) (Baur et al., 2020a) | 0.068
VQ-VAE (van den Oord et al., 2017) | 0.100
VQ-VAE + Transformer + Masked Residuals + Different Orderings (Ours) | 0.429

Our method showed a better performance than the autoencoder approaches from
the literature in all datasets (tab:syntseg and Figure LABEL:fig:exp3).
Compared to the numbers in Baur et al. (2020a), our autoencoder-based models
got a lower performance on the common dataset (MSLUB dataset), where they
achieved an DICE score of 0.271 with the AE (dense), 0.154 with the AE
(spatial), and 0.323 with the VAE (dense). We believe that the discrepancy
comes mostly from the significant post-processing of the Baur et al. (2020a)
work as presented in Table 8 of this reference. Differences might also arise
from the difference in resolution, as the DICE score is not invariant to
resolution.

fig:exp3

Figure 2: Residual maps on the real lesions from the variational autoencoder
and our transformer-based method.

##  5 Conclusion

Automatically determining the presence of lesion and delineating their
boundaries is essential to the introduction complex models of rich
neuroimaging features in clinical care. In this study, we propose a novel
transformer-based approach for anomaly detection and segmentation which
achieves state-of-the-art results in all tested tasks when compared to
competing methods. Transformers are making impressive gains in image analysis,
and here we show that their use to identify anomalies holds great promise. We
hope that our work will inspire further investigation of the properties of
transformers for anomaly detection in medical images, the development of new
network designs, exploration of a wider variety of conditioning information,
and the application of transformers to other medical data.

WHLP and MJC are supported by Wellcome Innovations [WT213038/Z/18/Z]. PTD is
supported by the EPSRC Research Council, part of the EPSRC DTP, grant Ref:
[EP/R513064/1]. PN is supported by Wellcome Innovations [WT213038/Z/18/Z] and
the UCLH NIHR Biomedical Research Centre. This research has been conducted
using the UK Biobank Resource (Project number: 58292).

## References

  * Alfaro-Almagro et al. (2018) Fidel Alfaro-Almagro, Mark Jenkinson, Neal K Bangerter, Jesper LR Andersson, Ludovica Griffanti, Gwenaëlle Douaud, Stamatios N Sotiropoulos, Saad Jbabdi, Moises Hernandez-Fernandez, Emmanuel Vallee, et al.  Image processing and quality control for the first 10,000 brain imaging datasets from uk biobank.  _Neuroimage_ , 166:400–424, 2018.
  * Avants et al. (2011) Brian B Avants, Nicholas J Tustison, Gang Song, Philip A Cook, Arno Klein, and James C Gee.  A reproducible evaluation of ants similarity metric performance in brain image registration.  _Neuroimage_ , 54(3):2033–2044, 2011.
  * Bakas et al. (2017) Spyridon Bakas, Hamed Akbari, Aristeidis Sotiras, Michel Bilello, Martin Rozycki, Justin S Kirby, John B Freymann, Keyvan Farahani, and Christos Davatzikos.  Advancing the cancer genome atlas glioma mri collections with expert segmentation labels and radiomic features.  _Scientific data_ , 4(1):1–13, 2017.
  * Bakas et al. (2018) Spyridon Bakas, Mauricio Reyes, Andras Jakab, Stefan Bauer, Markus Rempfler, Alessandro Crimi, Russell Takeshi Shinohara, Christoph Berger, Sung Min Ha, Martin Rozycki, et al.  Identifying the best machine learning algorithms for brain tumor segmentation, progression assessment, and overall survival prediction in the brats challenge.  _arXiv preprint arXiv:1811.02629_ , 2018.
  * Baur et al. (2018) Christoph Baur, Benedikt Wiestler, Shadi Albarqouni, and Nassir Navab.  Deep autoencoding models for unsupervised anomaly segmentation in brain mr images.  In _International MICCAI Brainlesion Workshop_ , pages 161–169. Springer, 2018.
  * Baur et al. (2020a) Christoph Baur, Stefan Denner, Benedikt Wiestler, Nassir Navab, and Shadi Albarqouni.  Autoencoders for unsupervised anomaly segmentation in brain mr images: A comparative study.  _Medical Image Analysis_ , page 101952, 2020a.
  * Baur et al. (2020b) Christoph Baur, Benedikt Wiestler, Shadi Albarqouni, and Nassir Navab.  Scale-space autoencoders for unsupervised anomaly segmentation in brain mri.  In _International Conference on Medical Image Computing and Computer-Assisted Intervention_ , pages 552–561. Springer, 2020b.
  * Chen et al. (2020a) Mark Chen, Alec Radford, Rewon Child, Jeffrey Wu, Heewoo Jun, David Luan, and Ilya Sutskever.  Generative pretraining from pixels.  In _International Conference on Machine Learning_ , pages 1691–1703. PMLR, 2020a.
  * Chen et al. (2020b) Xiaoran Chen, Suhang You, Kerem Can Tezcan, and Ender Konukoglu.  Unsupervised lesion detection via image restoration with a normative prior.  _Medical image analysis_ , 64:101713, 2020b.
  * Child et al. (2019) Rewon Child, Scott Gray, Alec Radford, and Ilya Sutskever.  Generating long sequences with sparse transformers.  _arXiv preprint arXiv:1904.10509_ , 2019.
  * Choi et al. (2018) Hyunsun Choi, Eric Jang, and Alexander A Alemi.  Waic, but why? generative ensembles for robust anomaly detection.  _arXiv preprint arXiv:1810.01392_ , 2018.
  * Choromanski et al. (2020) Krzysztof Choromanski, Valerii Likhosherstov, David Dohan, Xingyou Song, Andreea Gane, Tamas Sarlos, Peter Hawkins, Jared Davis, Afroz Mohiuddin, Lukasz Kaiser, et al.  Rethinking attention with performers.  _arXiv preprint arXiv:2009.14794_ , 2020.
  * Dosovitskiy et al. (2020) Alexey Dosovitskiy, Lucas Beyer, Alexander Kolesnikov, Dirk Weissenborn, Xiaohua Zhai, Thomas Unterthiner, Mostafa Dehghani, Matthias Minderer, Georg Heigold, Sylvain Gelly, et al.  An image is worth 16x16 words: Transformers for image recognition at scale.  _arXiv preprint arXiv:2010.11929_ , 2020.
  * Dumoulin et al. (2016) Vincent Dumoulin, Ishmael Belghazi, Ben Poole, Olivier Mastropietro, Alex Lamb, Martin Arjovsky, and Aaron Courville.  Adversarially learned inference.  _arXiv preprint arXiv:1606.00704_ , 2016.
  * Elliott and Peakman (2008) Paul Elliott and Tim C Peakman.  The uk biobank sample handling and storage protocol for the collection, processing and archiving of human blood and urine.  _International journal of epidemiology_ , 37(2):234–244, 2008.
  * Esser et al. (2020) Patrick Esser, Robin Rombach, and Björn Ommer.  Taming transformers for high-resolution image synthesis.  _arXiv preprint arXiv:2012.09841_ , 2020.
  * Griffanti et al. (2016) Ludovica Griffanti, Giovanna Zamboni, Aamira Khan, Linxin Li, Guendalina Bonifacio, Vaanathi Sundaresan, Ursula G Schulz, Wilhelm Kuker, Marco Battaglini, Peter M Rothwell, et al.  Bianca (brain intensity abnormality classification algorithm): A new tool for automated segmentation of white matter hyperintensities.  _Neuroimage_ , 141:191–205, 2016.
  * Jenkinson et al. (2012) Mark Jenkinson, Christian F Beckmann, Timothy EJ Behrens, Mark W Woolrich, and Stephen M Smith.  Fsl.  _Neuroimage_ , 62(2):782–790, 2012.
  * Jun et al. (2020) Heewoo Jun, Rewon Child, Mark Chen, John Schulman, Aditya Ramesh, Alec Radford, and Ilya Sutskever.  Distribution augmentation for generative modeling.  In _International Conference on Machine Learning_ , pages 5006–5019. PMLR, 2020.
  * Kamnitsas et al. (2017) Konstantinos Kamnitsas, Christian Ledig, Virginia FJ Newcombe, Joanna P Simpson, Andrew D Kane, David K Menon, Daniel Rueckert, and Ben Glocker.  Efficient multi-scale 3d cnn with fully connected crf for accurate brain lesion segmentation.  _Medical image analysis_ , 36:61–78, 2017.
  * Kuijf et al. (2019) Hugo J Kuijf, J Matthijs Biesbroek, Jeroen De Bresser, Rutger Heinen, Simon Andermatt, Mariana Bento, Matt Berseth, Mikhail Belyaev, M Jorge Cardoso, Adria Casamitjana, et al.  Standardized assessment of automatic segmentation of white matter hyperintensities and results of the wmh segmentation challenge.  _IEEE transactions on medical imaging_ , 38(11):2556–2568, 2019.
  * Lesjak et al. (2018) Žiga Lesjak, Alfiia Galimzianova, Aleš Koren, Matej Lukin, Franjo Pernuš, Boštjan Likar, and Žiga Špiclin.  A novel public mr image dataset of multiple sclerosis patients with lesion segmentations based on multi-rater consensus.  _Neuroinformatics_ , 16(1):51–63, 2018.
  * Matthey et al. (2017) Loic Matthey, Irina Higgins, Demis Hassabis, and Alexander Lerchner.  dsprites: Disentanglement testing sprites dataset, 2017.
  * Menze et al. (2014) Bjoern H Menze, Andras Jakab, Stefan Bauer, Jayashree Kalpathy-Cramer, Keyvan Farahani, Justin Kirby, Yuliya Burren, Nicole Porz, Johannes Slotboom, Roland Wiest, et al.  The multimodal brain tumor image segmentation benchmark (brats).  _IEEE transactions on medical imaging_ , 34(10):1993–2024, 2014.
  * Miller et al. (2016) Karla L Miller, Fidel Alfaro-Almagro, Neal K Bangerter, David L Thomas, Essa Yacoub, Junqian Xu, Andreas J Bartsch, Saad Jbabdi, Stamatios N Sotiropoulos, Jesper LR Andersson, et al.  Multimodal population brain imaging in the uk biobank prospective epidemiological study.  _Nature neuroscience_ , 19(11):1523–1536, 2016\.
  * Nalisnick et al. (2018) Eric Nalisnick, Akihiro Matsukawa, Yee Whye Teh, Dilan Gorur, and Balaji Lakshminarayanan.  Do deep generative models know what they don’t know?  _arXiv preprint arXiv:1810.09136_ , 2018.
  * Pawlowski et al. (2018) Nick Pawlowski, MC Lee, Martin Rajchl, Steven McDonagh, Enzo Ferrante, Konstantinos Kamnitsas, Sam Cooke, Susan Stevenson, Aneesh Khetani, Tom Newman, et al.  Unsupervised lesion detection in brain ct using bayesian convolutional autoencoders.  _Medical Imaging with Deep Learning_ , 2018.
  * Porz et al. (2014) Nicole Porz, Stefan Bauer, Alessia Pica, Philippe Schucht, Jürgen Beck, Rajeev Kumar Verma, Johannes Slotboom, Mauricio Reyes, and Roland Wiest.  Multi-modal glioblastoma segmentation: man versus machine.  _PloS one_ , 9(5):e96873, 2014.
  * Prados et al. (2016) Ferran Prados, Manuel Jorge Cardoso, Baris Kanber, Olga Ciccarelli, Raju Kapoor, Claudia AM Gandini Wheeler-Kingshott, and Sebastien Ourselin.  A multi-time-point modality-agnostic patch-based method for lesion filling in multiple sclerosis.  _Neuroimage_ , 139:376–384, 2016.
  * Radford et al. (2018) Alec Radford, Karthik Narasimhan, Tim Salimans, and Ilya Sutskever.  Improving language understanding by generative pre-training.  2018\.
  * Radford et al. (2019) Alec Radford, Jeffrey Wu, Rewon Child, David Luan, Dario Amodei, and Ilya Sutskever.  Language models are unsupervised multitask learners.  _OpenAI blog_ , 1(8):9, 2019.
  * Razavi et al. (2019) Ali Razavi, Aaron van den Oord, and Oriol Vinyals.  Generating diverse high-fidelity images with vq-vae-2.  _arXiv preprint arXiv:1906.00446_ , 2019.
  * Sudlow et al. (2015) Cathie Sudlow, John Gallacher, Naomi Allen, Valerie Beral, Paul Burton, John Danesh, Paul Downey, Paul Elliott, Jane Green, Martin Landray, et al.  Uk biobank: an open access resource for identifying the causes of a wide range of complex diseases of middle and old age.  _Plos med_ , 12(3):e1001779, 2015.
  * Tay et al. (2020) Yi Tay, Mostafa Dehghani, Samira Abnar, Yikang Shen, Dara Bahri, Philip Pham, Jinfeng Rao, Liu Yang, Sebastian Ruder, and Donald Metzler.  Long range arena: A benchmark for efficient transformers.  _arXiv preprint arXiv:2011.04006_ , 2020.
  * van den Oord et al. (2017) Aäron van den Oord, Oriol Vinyals, and Koray Kavukcuoglu.  Neural discrete representation learning.  _CoRR_ , abs/1711.00937, 2017.  URL http://arxiv.org/abs/1711.00937.
  * Vaswani et al. (2017) Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Lukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  _arXiv preprint arXiv:1706.03762_ , 2017.
  * Wardlaw et al. (2013) Joanna M Wardlaw, Eric E Smith, Geert J Biessels, Charlotte Cordonnier, Franz Fazekas, Richard Frayne, Richard I Lindley, John T O’Brien, Frederik Barkhof, Oscar R Benavente, et al.  Neuroimaging standards for research into small vessel disease and its contribution to ageing and neurodegeneration.  _The Lancet Neurology_ , 12(8):822–838, 2013\.
  * Yuh et al. (2012) Esther L Yuh, Shelly R Cooper, Adam R Ferguson, and Geoffrey T Manley.  Quantitative ct improves outcome prediction in acute traumatic brain injury.  _Journal of neurotrauma_ , 29(5):735–746, 2012\.

##  Appendix A Experimental Details and Further Analysis

###  A.1 Experiment #1 – Anomaly Segmentation on Synthetic Data

#### Dataset

To assess the performance on anomaly segmentation, we utilized a subsample of
the MedNIST dataset, where we used the 2D images of the HeadCT category to
train our VQ-VAE and transformer models (Figure LABEL:fig:S1). From the
original 10,000 HeadCT images, we used 8,000 images as the training set and
1,000 images for the validation set. The test set was comprised of 100 images
contaminated with sprites (i.e., synthetic anomalies) obtained from the
dsprites dataset Matthey et al. (2017). We selected the sprites images that
overlapped a significant portion of the head, and their values were set as 0
or 1.

fig:S1

Figure A1: Examples of the training set (HeadCT class) from the MedNIST
dataset.

#### Models

Our VQ-VAE models had a similar architecture from van den Oord et al. (2017).
The encoder consists of three strided convolutional layers with stride 2 and
window size 4 × 4. All these convolution layers had a ReLU activation
following them. This structure is followed by two residual 3×3 blocks
(implemented as 3×3 conv, ReLU, 1×1 conv, ReLU). The decoder similarly has two
residual 3×3 blocks, followed by three transposed convolutions with stride 2
and window size 4×4. All the convolution layers have 256 hidden units. The
inputted images have the dimension of 64x64 pixels which result in a latent
representation of 8x8 latent variables. For this experiment, we used a
codebook with 16 different codes. Our performers corresponded to the
transformer’s decoder structure with 24 layers with an embedding size of 256.
The embedding, feed-forward and attention dropout had a probability of 0.1.

#### State-of-the-art Models

We compared our models against state-of-the-art autoencoder-based methods (AE
dense, AE spatial, and VAE). In this experiment, we used an unified network
architecture adapted from a recent comparison study Baur et al. (2020a). Since
MedNIST images are smaller than those used in the comparison study, our models
did not have the first block (resolution of 64x64x32) and the last block
(resolution of 128x128x32). All the other blocks and layers are similar to the
original study. To train these models, we used the ADAM optimiser with
learning rate 5e-4, an exponential learning rate decay with gamma of 0.9999,
and we trained over 1,500 epochs with batch-size 256. Similar to Baur et al.
(2020a), we used the held-out validation set to select the models to be
assessed.

#### Training settings

To train the VQ-VAE models, we use the ADAM optimiser with learning rate 5e-4,
an exponential learning rate decay with gamma of 0.999, and we trained over
1,500 epochs with batch-size 256. We train the codebook using an exponential
moving average algorithm. To stabilise the codebook’s learning, in the first
100 epoch, we warm up the moving average decay from a gamma decay of 0.5,
gradually increasing to a gamma decay of 0.99. This allows the codes to adapt
faster to the frequent changes at the beginning of the training. To train the
performers, we used the ADAM optimiser with learning rate 5e-4, an exponential
learning rate decay with gamma of 0.9999, and we trained over 1,500 epochs
with batch-size 128. We used data augmentation to increase the number of
training images. We randomly performed affine transformations (scale,
translate, and rotate operations) and horizontally flip the images.

#### Different Ordering Classes

In our study, we also analysed 3 other classes of orderings (Figure
LABEL:fig:S2): a S-curve order that traverses rows in alternating directions,
and a Hilbert space-filling curve order that generates nearby pixels in the
image consecutively. In all these classes, and a random class, where the
sequence of latent variables is randomly sorted. Similar to the raster class,
we augmented the number of possible orderings by reflecting and transposing
the images, generating in total 8 different orderings per class.

fig:S2

Figure A2: Different orderings used to transform the 2D discrete
representation into a sequence.

In the tab:S1, we can observe the performance of each class. The orderings had
a performance varying from 0.843 to 0.895. We can observe that the random
ordering performed with the lowest performance. This was expected since we
would expect that local information of the image would help the transformer
modelling of the healthy data. Since the random ordering may not include the
local data in the context to predict a variable value autoregressively, this
might reduce its performance as anomaly detector too. Finally, we evaluated
the performance when combining all the orderings. It was observed a small gain
when using an ensemble of all four classes compared to the raster class only.
We opt to use the raster ordering in the main analysis to reduce the time of
training and processing.

tab:S1

Table A1: Performance of different classes of ordering and the ensemble with
all classes.

Method | $\lceil DICE\rceil$
---|---
8 different raster orderings | 0.895
8 different S-curve orderings | 0.883
8 different Hilbert curve orderings | 0.890
8 different random orderings | 0.843
32 different orderings | 0.899

#### Same ordering but different random seed

The ensemble of 8 transformers gave a boost in performance in the segmentation
comparing with a single transformer. To verify if it was due to the
variability of the generative models with different orderings instead of using
an ensemble with more models, we trained 8 models using the same raster
ordering but with different random seed. We can observe a drop in performance
when using an ensemble of transformers using the same ordering but different
random seeds, from 0.895 to 0.826.

#### Anomalies intensity

We also evaluated the influence of the synthetic anomalies’ intensity and
texture. For this, we varied the intensity of the sprites in the image from 0
to 1 and measured the segmentation performance (best achievable DICE-score).
We also performed this approach with the addition of an additive Gaussian
noise with standard deviation of 0.2. From Figure LABEL:fig:S3, we can observe
that our transformer-based method is more robust to the change in intensity
with a narrow drop in performance when the anomaly intensity is closer to the
tissue mean values.

fig:S3

Figure A3: Analysis of the intensity of the anomalies A) with and B) without
the additive noise.

###  A.2 Experiment #2 – Image-wise Anomaly Detection on Synthetic Data

#### Dataset

In this experiment, we used the same training set from the Experiment #1. For
evaluation, we used 1,000 images from the HeadCT class as the in-distribution
test set, the 100 HeadCT images contaminated by sprites anomalies as the near
out-of-distribution set (near OOD), and 1,000 images of each other classes
from the MedNIST dataset (AbdomenCT, BreastMRI, CXR, ChestCT, and Hand) as the
far out-of-distribution set (far OOD) (Figure LABEL:fig:S4). To train our
general purpose VQ-VAE, we added 8,000 images from each other classes to our
training set and 1,000 images to our validation set.

fig:S4

Figure A4: Examples from the near out-of-distribution set (near OOD) and far
out-of-distribution set (far OOD).

#### Models and Training Settings

To train the VQ-VAE with general purpose and its transformers, we used the
same architecture, and the same training settings from the models from
Experiment #1.

#### Performance of Experiment #1 models on Anomaly Detection

Besides the VAE, in tab:S2 we also present the performance of the other
autoencoder models from Experiment #1 to perform the anomaly detection task.

tab:S2

Table A2: Performance of the methods on image-wise anomaly detection.

| AUCROC |  AUCPR In |  AUCPR Out | FPR95 | FPR99 | FPR999
---|---|---|---|---|---|---
vs. far OOD classes |  |  |  |  |  |
AE (Dense) (Baur et al., 2020a) | 0.351 | 0.875 | 0.066 | 0.969 | 0.987 | 0.987
AE (Spatial) (Baur et al., 2020a) | 0.337 | 0.874 | 0.063 | 0.959 | 0.998 | 0.998
VAE (Dense) (Baur et al., 2020a) | 0.298 | 0.855 | 0.060 | 0.986 | 0.996 | 0.996
VQ-VAE (van den Oord et al., 2017) | 0.241 | 0.834 | 0.056 | 0.987 | 0.996 | 0.996
vs. near OOD classes |  |  |  |  |  |
AE (Dense) (Baur et al., 2020a) | 0.106 | 0.093 | 0.669 | 1.000 | 1.000 | 1.000
AE (Spatial) (Baur et al., 2020a) | 0.215 | 0.103 | 0.721 | 1.000 | 1.000 | 1.000
VAE (Dense) (Baur et al., 2020a) | 0.111 | 0.094 | 0.672 | 1.000 | 1.000 | 1.000
VQ-VAE (van den Oord et al., 2017) | 0.024 | 0.089 | 0.645 | 1.000 | 1.000 | 1.000

#### Impact of the general VQ-VAE on anomaly segmentation

The VQ-VAE with general purpose reduce the impact of the encoder to the
analysis, working mainly as a compression mechanism. In this analysis, we
evaluate the performance of this models when performing anomaly segmentation.
Using the general purpose VQ-VAE, we obtained a DICE score of 0.886, just a
small decrease compared to the models trained only on HeadCT data.

###  A.3 Experiment #3 – Anomaly Segmentation on Real Data

#### MRI Datasets

In our experiment, we used three datasets: the UK Biobank (UKB) dataset Sudlow
et al. (2015), the White Matter Hyperintensities Segmentation Challenge (WMH)
dataset Kuijf et al. (2019), the Multimodal Brain Tumor Image Segmentation
Benchmark (BRATS) dataset Bakas et al. (2017, 2018); Menze et al. (2014), and
the Multiple Sclerosis dataset from the University Hospital of Ljubljana
(MSLUB) dataset Lesjak et al. (2018) (Figure LABEL:fig:S5).

fig:S5

Figure A5: Examples from the test set of the Experiment #3.

The UKB is a study that aims to follow the health and well-being of 500,000
volunteer participants across the United Kingdom. From these participants, a
subsample was chosen to collect multimodal imaging, including structural
neuroimaging. Here, we used an early release of the project’s data comprising
of 33,318 HC participants. More details about the dataset and imaging
acquisition can be found elsewhere Alfaro-Almagro et al. (2018); Elliott and
Peakman (2008); Miller et al. (2016). The UK Biobank dataset has available a
mask for hyperintensities white matter lesions obtained using BIANCA Griffanti
et al. (2016); Jenkinson et al. (2012). We selected the 15k subjects with the
lowest lesion volume to train our VQ-VAE model.

The BRATS challenge is an initiative that aims to evaluate methods for the
segmentation of brain tumours by providing a 3D MRI dataset with ground truth
tumour segmentation annotated by expert board-certified neuroradiologists
Bakas et al. (2017, 2018); Menze et al. (2014). Our study used the 2018
version of the dataset composed by the MR scans of 420 patients with
glioblastoma or lower grade glioma. The images were acquired with different
clinical protocols and various scanners from multiple (n=19) institutions.
Note, the available images from the BRATS dataset were already skull stripped.

The WHM dataset is an initiative to directly compare automated WMH
segmentation techniques Kuijf et al. (2019). The dataset was acquired from
five different scanners from three different vendors in three different
hospitals in the Netherlands and Singapore. It is composed by 60 subjects
where the WMH were manually segmented according to the STandards for ReportIng
Vascular changes on nEuroimaging (STRIVE) Wardlaw et al. (2013). The MSLUB
dataset is a publicly available dataset for validation of lesion segmentation
methods. The dataset consists of 30 images from multiple sclerosis patients
that were acquired using conventional MR imaging sequences. For each case, a
reference lesion segmentation was created by three independent raters and
merged into a consensus. This way, we have access to a precise and reliable
target to evaluate segmentation methods. Full description regarding data
acquisition and imaging protocol can be found at Lesjak et al. (2018).

In our study, we used the FLAIR images from these three datasets to evaluate
our models. For each of these FLAIR images, the UKB and WMH datasets had a
white matter hyperintensity segmentation map, the BRATS dataset had a tumour
segmentation map, and the MSLUB dataset had a white-matter lesion segmentation
map. We also used T1w and brain mask to perform the MRI pre-processing.

#### MRI Pre-processing

We pre-process our images to be normalized in a common space. For this reason,
all scans and lesion masks were registered to MNI space using rigid + affine
transformation. This registration was performed using the Advanced
Normalisations Tools (ANTs - version 2.3.4) Avants et al. (2011). Since our
anomaly segmentation method relies on a training set composed by a population
with a low occurrence of lesions and anomalies, we tried to minimize the
occurrence of lesions on the transformers’ training set. For this reason,
after the traditional MRI pre-processing, we used the NiftySeg package
(version 1.0) Prados et al. (2016) to mitigate the influence of the lesions
present our training set. Using the seg_FillLesions function and the lesion
maps supplied by the UKB dataset, we in-painted the few white matter
hyperintensities present in the FLAIR images using a non-local lesion filling
strategy based on a patch based in-painting technique for image completion.
Since the VQ-VAE performs mainly a dimensionality reduction in our method, it
was trained using the normalized dataset without the NiftySeg in-painting. We
believe that the presence of the lesions in the VQ-VAE training set is
important to avoid that the autoencoder method avoids performing the
correction on its encoding part. If the encoder corrects the code by itself,
the transformer would not be able to detect the presence of a lesion. This
missing detection would result in a resampling mask that filters out the
encoder correction. On Experiment #2, we show that this approach does not
prejudice the performance of the segmentation. Finally, we selected 4 axials
slices (z=89, 90, 91, 92) per FLAIR image and, we center cropped these slices
to have the dimensions of 224 x 224 pixels.

#### Models

Our VQ-VAE models had a similar architecture from the Experiment #1 but with
three residual 3×3 blocks (instead 2). All the convolution layers had 256
hidden units. The inputted images had 224x224 pixels which result in a latent
representation of 28x28 latent variables. For this experiment, we used a
codebook with 32 different codes. Our performers corresponded to the
transformer’s decoder structure with 16 layers with an embedding size of 256.
The embedding, feed-forward and attention dropout had a probability of 0.3.

#### Training settings

To train the VQ-VAE models, we use the ADAM optimiser with learning rate 1e-3,
an exponential learning rate decay with gamma of 0.99995, and we trained over
500 epochs with batch-size 256. We train the codebook similar to Experiment
#1. To train the performers, we used the ADAM optimiser with learning rate
1e-3, an exponential learning rate decay with gamma of 0.99992, and we trained
over 150 epochs with batch-size 128. We used data augmentation to increase the
number of training images. We randomly performed small translation
transformations as well as random intensity shift and random adjusts in
contrast.

#### State-of-the-art Models

In this experiment, we used the network architecture from Baur et al. (2020a)
for the autoencoder only-based approaches. Even with our bigger input images,
we opted to use the same instead a version with an extra downsampling step to
avoid a model structure with a bottleneck with a much smaller resolution than
our proposed method. To train these models, we used the ADAM optimiser with
learning rate 1e-3, an exponential learning rate decay with gamma of 0.99995,
and we trained over 1,500 epochs with batch-size 256. Similar to Baur et al.
(2020a), we used the held-out validation set to select the models to be
assessed.

#### Impact of Mitigating Lesions on the Training set

In our pre-processing, we in-painted the white matter hyperintensity of the
training set using the NiftySeg package, as to simulate completely lesion-free
data. Here, we compare the performance of our method without this step. By
skipping this pre-processing step and training the transformers on the new
dataset, we can observe a drop in performance, from 0.232 to 0.051 in the UKB
dataset, from 0.378 to 0.264 in the MSLUB dataset, from 0.759 to 0.677 in the
BRATS dataset, and from 0.429 to 0.349 in the WMH dataset. We believe that the
highly expressive transformers can learn the few white matter hyperintensities
present in the dataset and associate a higher probability of occurrence,
reducing the performance in detection.

#### Performance of different stages of the method

Like Experiment #1, we evaluated the performance of each step of our method.
From tab:S3, we can observe that each step presents an incremental
improvement.

tab:S3

Table A3: Performance of each processing step of our transformer-based
approach in the Experiment #3

UKB Dataset | $\lceil DICE\rceil$
---|---
VQ-VAE (van den Oord et al., 2017) | 0.028
VQ-VAE + Transformer (Ours) | 0.079
VQ-VAE + Transformer + Masked Residuals (Ours) | 0.104
VQ-VAE + Transformer + Masked Residuals + Different Orderings (Ours) | 0.232
MSLUB Dataset |
VQ-VAE (van den Oord et al., 2017) | 0.040
VQ-VAE + Transformer (Ours) | 0.097
VQ-VAE + Transformer + Masked Residuals (Ours) | 0.234
VQ-VAE + Transformer + Masked Residuals + Different Orderings (Ours) | 0.378
BRATS Dataset |
VQ-VAE (van den Oord et al., 2017) | 0.331
VQ-VAE + Transformer (Ours) | 0.431
VQ-VAE + Transformer + Masked Residuals (Ours) | 0.476
VQ-VAE + Transformer + Masked Residuals + Different Orderings (Ours) | 0.759
WMH Dataset |
VQ-VAE (van den Oord et al., 2017) | 0.100
VQ-VAE + Transformer (Ours) | 0.205
VQ-VAE + Transformer + Masked Residuals (Ours) | 0.269
VQ-VAE + Transformer + Masked Residuals + Different Orderings (Ours) | 0.429

#### Post-processing Impact

Similar to Baur et al. (2020a), we verified the performance of the methods
using the prior knowledge that multiple sclerosis lesions would appear as
positive residuals as these lesions appear as hyper-intense in FLAIR images.
We assumed the same for the WMH lesions. By using only the positive values of
the residuals as a post-processing step, we got a gain on performance on both
the autoencoders only based methods and our approach (tab:S4).

tab:S4

Table A4: Performance of each method using post-processing step.

UKB Dataset | $\lceil DICE\rceil$
---|---
AE (Dense) (Baur et al., 2020a) | 0.079
AE (Spatial) (Baur et al., 2020a) | 0.054
VAE (Dense) (Baur et al., 2020a) | 0.071
VQ-VAE (van den Oord et al., 2017) | 0.056
VQ-VAE + Transformer + Masked Residuals + Different Orderings (Ours) | 0.297
MSLUB Dataset |
AE (Dense) (Baur et al., 2020a) | 0.106
AE (Spatial) (Baur et al., 2020a) | 0.067
VAE (Dense) (Baur et al., 2020a) | 0.106
VQ-VAE (van den Oord et al., 2017) | 0.077
VQ-VAE + Transformer + Masked Residuals + Different Orderings (Ours) | 0.465
WMH Dataset |
AE (Dense) (Baur et al., 2020a) | 0.166
AE (Spatial) (Baur et al., 2020a) | 0.151
VAE (Dense) (Baur et al., 2020a) | 0.161
VQ-VAE (van den Oord et al., 2017) | 0.143
VQ-VAE + Transformer + Masked Residuals + Different Orderings (Ours) | 0.441

##  Appendix B More Residuals maps from Experiment #1

fig:appendixb

Figure A6: More examples of residual maps on the synthetic dataset from the
variational autoencoder and different steps of our approach.

##  Appendix C More Residuals maps from Experiment #3

fig:appendixc

Figure A7: More examples of residual maps on the real lesions from the
variational autoencoder and our method.


# Insights from an adaptive variational wave function study of the Fermi-
Hubbard Model

Gaurav Gyawali  Laboratory of Atomic and Solid State Physics, Cornell
University, Ithaca, NY, 14853  Michael J Lawler  Laboratory of Atomic and
Solid State Physics, Cornell University, Ithaca, NY, 14853  Department of
Physics, Applied Physics and Astronomy, Binghamton University, Binghamton, New
York 13902

###### Abstract

Approximating the ground states of strongly interacting electron systems in
quantum chemistry and condensed matter physics is expected to be one of the
earliest applications of quantum computers. In this paper, we prepare highly
accurate ground states of the Fermi-Hubbard model for small grids up to 6
sites (12 qubits) by using an interpretable, adaptive variational quantum
eigensolver(VQE) called ADAPT-VQE Grimsley19 . In contrast with non-adaptive
VQE, this algorithm builds a system-specific ansatz by adding an optimal gate
built from one-body or two-body fermionic operators at each step. We show this
adaptive method outperforms the non-adaptive counterpart in terms of fewer
variational parameters, short gate depth, and scaling with the system size.
The fidelity and energy of the prepared state appear to improve asympotically
with ansatz depth.We also demonstrate the application of adaptive variational
methods by preparing excited states and Green functions using a proposed
ADAPT-SSVQE algorithm. Lower depth, asymptotic convergence, noise tolerance of
a variational approachVQE14 ; McClean2016 ; Sharma_2020  and a highly
controllable, system specific ansatz make the adaptive variational methods
particularly well-suited for NISQ devices.

##  I Introduction

Finding the ground state of a many-body system is an important problem in
condensed matter physics. There are numerous ways of doing so on a classical
computer including exact diagonalization, density functional theory(DFT)DFT
and density matrix renormalization group (DMRG) DMRG . However, all of these
methods suffer from an exponential scaling of computational resources with the
system size for strongly correlated electrons. Quantum computing has emerged
as a promising platform to mitigate the exponential scaling of classical
resources because the Hilbert space of the device also grows exponentially.
The current hype of quantum computers owes much to the discovery of sub-
exponential time quantum algorithms such as Shor’s factoring algorithm and
quantum phase estimation for classically hard problems. Shor1997 ; Mike&Ike .
It is therefore important to study the performance of new quantum algorithms
and assess their use in potential breakthrough applications.

Although significant strides have been made in the direction of simulating
quantum systems using quantum computers, as envisioned by FeynmanFeynman1982 ,
we are still in the NISQ era Preskill2018 . Noisy quantum devices with around
$100$ qubits have been available lately, and scientific efforts are already in
place to improve the noise tolerance and increase number of qubits IBM ;
Google . A hybrid quantum-classical algorithm called the variational quantum
eigensolver(VQE) emerged in the field of quantum chemistry to harness the
power of NISQ devices by approximating the ground states of interacting
HamiltoniansVQE14 ; McClean2016 . VQE requires optimization of a large number
of variational parameters but is quite robust to noise and results in a low-
depth circuit, which makes it ideal for the near term noisy devices with few
qubits.

An algorithm called Adaptive Derivative-Assembled Pseudo-Trotter ansatz VQE
(ADAPT-VQE) has been proposed as a quasi-optimal ansatz with the minimal
number of operators for a desired level of accuracy Grimsley19 ; Zhang2021 .
The ansatz is systematically grown by adding optimal fermionic gates one or
more at a time such that the maximal amount of correlation energy is recovered
at each step. This algorithm generates an ansatz with a small number of
parameters and shallow depth, and was tested for various small molecules on
quantum hardware. In this paper, we extend the application of ADAPT-VQE
algorithm to condensed matter physics by using it to prepare the ground as
well as excited states of a strongly correlated electron system called the
Fermi-Hubbard model Hubbard1963 .

The Hamiltonian of the Fermi-Hubbard model is defined as

$H=-t\sum_{<i,j>,\sigma}\left(c_{i\sigma}^{\dagger}c_{j\sigma}+c_{j\sigma}^{\dagger}c_{i\sigma}\right)+U\sum_{i}n_{i\uparrow}n_{i\downarrow},$
(1)

where $c_{i\sigma}^{\dagger}$,$c_{i\sigma}$ are fermionic creation and
annihilation operators corresponding to site $i$ and spin $\sigma$;
$n_{i,\sigma}=c_{i\sigma}^{\dagger}c_{i\sigma}$, and $\langle.\rangle$
indicates nearest neighbors. $t$ and $U$ are called the hopping and on-site
repulsion parameters respectively. Finding the ground state of a Fermi-Hubbard
model with uniform $U$ but arbitrary hoppings $(t_{ij}$) has been shown to be
in the QMA-complete complexity class BGorman2021 . Although it implies that
the worse-case scenario of the problem is difficult to solve even with a
quantum computer, there can be many interesting practical instances that can
be solved efficiently Cao2019 . We focus on a relatively simple yet practical
case in which the hopping and interaction parameters are uniform.

In addition to ground states, we can also use adaptive methods to prepare
excited states. This can be done by choosing a variational circuit and
applying subspace search VQE(SSVQE).Endo2020 ; Nakanishi2019 . The SSVQE
algorithm utilizes the conservation of orthogonality under unitary
transformation i.e. for a superposition of orthogonal input states, it outputs
another superposition of orthogonal states. The circuit is optimized in such a
way that these orthogonal output states are the desired excited states.
Combining this algorithm with adaptive variational circuits then produces
ADAPT-SSVQE, an adaptive method for generating excited states. The quality of
these excited states are then demonstrated with a computation of the
associated zero-temperature Green’s function via the Lehmann representation.

Our results show that ADAPT-VQE algorithm is highly successful at preparing
low-depth circuits to approximate ground states of the Fermi-Hubbard model for
small grids with high fidelity. We also demonstrate the application of the
adaptive variational method to prepare the excited states using ADAPT-SSVQE
and thus calculate the Green’s function. We show that ADAPT-VQE outperforms
state of the art non-adaptive ansatz in terms of performance with the system
size, number of parameters, and gate depth. We also observe asymptotic
improvement in fidelity as well as energy with each added gate.

##  II ADAPT-VQE algorithm

Figure 1:  Pictorial representation of an ADAPT-VQE circuit that prepares the
ground state of a 3 $\times$ 1 Fermi-Hubbard model with $t=1$ and $U=3$.
Creation and annihilation operators are denoted by a circle and a diamond
respectively. A set of diamond and a circle represent a one-body hopping term
($c^{\dagger}_{i,\sigma}c_{j,\sigma}$) whereas a set of four shapes represent
a two-body term
($c^{\dagger}_{i,\sigma}c_{j,\sigma}c^{\dagger}_{k,\sigma^{\prime}}c_{l,\sigma^{\prime}}$
). Two-body terms with $k=l$ are the correlated hopping terms for which a
circle and a diamond overlap. Transparency of the purple box indicates the
absolute value of the coefficient($\theta$) for that step with respect to the
maximum coefficient.

The ADAPT-VQE ansatz is built one operator at a time as opposed to the non-
adaptive ansatzs such as Hamiltonian variational ansatz and number preserving
ansatz in which the operators are fixed, and only optimization performed is
for the variational parameters HVA2015 ; Cade2019 . Each operator in this
ansatz is chosen from a pool of one-body and two-body Hamiltonians
$\hat{A}_{m}=\\{\hat{H}_{ij},\hat{H}_{ijkl}\\}$, where

$\displaystyle\hat{H}_{ij}$ $\displaystyle=c^{\dagger}_{i}c_{j}+hc$ (2)
$\displaystyle\hat{H}_{ijkl}$
$\displaystyle=c^{\dagger}_{i}c^{\dagger}_{j}c_{k}c_{l}+hc.,$ (3)

where $i,j,k,l$ are qubit labels. These Hamiltonians are chosen such that they
are spin conserving. The adaptive ansatz is applied on top of the initial-
state preparation circuit, which is chosen to be either a product state or a
Slater determinant. At the $n^{th}$ step of the ansatz, the energy gradient
for each of the operator ($\hat{A}_{m}$) i.e.

$\frac{\partial
E^{(n)}}{\partial\theta_{m}}=\left\langle\psi^{(n)}\left|\left[\hat{H},\hat{A}_{m}\right]\right|\psi^{(n)}\right\rangle$
(4)

is calculated. Here, $\ket{\psi^{(n)}}$ is the optimized wavefunction at the
$n^{th}$ step. The wavefunction at $(n+1)^{th}$ step is then given by

$\ket{\psi^{(n+1)}}=e^{i\theta_{n+1}\hat{A}_{n+1}}\ket{\psi^{(n)}},$ (5)

where $\hat{A}_{n+1}$ is the $\hat{A}_{m}$ with highest energy gradient. It is
important that we optimize the accumulated parameter set
$\vec{\theta}^{(n+1)}=\\{\vec{\theta}^{(n)},\theta_{n+1}$ } at each step
rather than using $\vec{\theta}^{(n)}$ from the previous step and just
optimizing $\theta_{n+1}$. This process is repeated until stopping criteria
are met i.e. $(E^{(n)}-E^{(n-1)})<\epsilon$ or
$\absolutevalue{\partialderivative{E^{(n)}}{\theta_{i}}}<\delta$. See Ref.
Grimsley19  for a more detailed description of this algorithm. The variational
algorithm was implemented using OpenFermion-FQE library available in Python
fqe21 . Mapping from fermionic to qubit operators was done via Jordan-Wigner
transform JW28  as described in Appendix A. We used L-BFGS-B algorithm
available in SciPy library to perform all optimizationsSciPy .

Formulation of the ADAPT-VQE algorithm in terms of fermionic operators makes
it convenient to interpret an ansatz. The $n^{th}$ step of an ansatz involving
operator $\hat{A}_{n}$ and variational parameter $\theta_{n}$ can be thought
of as the time evolution of the system under $\hat{A}_{n}$ for time
$-\theta_{n}$ i.e. $e^{-i\hat{A}_{n}(-\theta_{n})}$. At each step only the
most optimal gate is selected and so earlier steps are more important than
later ones. A pictorial representation of the ADAPT-VQE optimized ansatz for a
$3\times 1$ Fermi-Hubbard model is shown in Fig. 1. Here, the first four
operators are one-body/hopping gates between qubits corresponding to different
sites _and same spin_. They create a state close to the ground state for the
non-interacting part of the Hamiltonian with the same spin symmetry. The fifth
operator introduces two-body terms that apply the corrections due to
interaction. Of these corrections, majority are three qubit gates i.e.
correlated hopping. These gates appear more often when the interaction is
weak, and the system size is small. We have observed that the ground state of
a Hamiltonian with no interaction (not shown in the figure) can be prepared
only by using the hopping gates. Introducing the interaction therefore adds
three-qubit and four qubit gates and their relative numbers is determined by
the system size and interaction strength. Hence the ADAPT-VQE ansatz tells us
the spin symmetry, importance of interactions, and may even give us a sense of
how electrons move around in the material via correlated hopping processes.

##  III Green’s function

A fundamental object of interest in many-body physics is the Green’s function.
It contains many important dynamic and static properties of a many-body system
including the quasi-particle dispersionColeman2015 . Throughout this paper,
we’ll focus on Green’s function at zero temperature. Consider a fermionic
system described by Hamiltonian $H$. The retarded-time Green’s function for
this system is defined as follows Mahan1990 ; Coleman2015 :

$G(t)=-i\Theta(t)\left\langle
c_{a}(t)c_{b}^{\dagger}(0)+c_{b}^{\dagger}(0)c_{a}(t)\right\rangle,$ (6)

where $c^{\dagger}_{a}$, $c_{a}$ are creation and annihilation operators for a
fermionic mode $a$; for instance $a=(k,\sigma)$, with $k$ being the momentum,
and $\sigma$ the spin. $\Theta(t)$ is the Heaviside step function, and
$c_{a}(t)$ is a Heisenberg operator defined as
$c_{a}(t)=e^{iHt}c_{a}(0)e^{-iHt}$. The expectation value is calculated with
respect to the ground state $\ket{G}$ of $H$. Taking the Fourier transform of
Eq. 6 gives us the propagator $G(k,\omega)$. The imaginary part of the
propagator gives us information about the spectrum of the many-body system,
and thus known as the spectral function:

$A(k,\omega)=-\frac{1}{\pi}\operatorname{Im}G(k,\omega).$ (7)

We can write the retarded Greens function in Lehmann representation, which
relates the definition in Eq. 6 to the eigenspectrum of $H$, as follows:

$\displaystyle G(k,\omega)=$
$\displaystyle\sum_{n}\left(\frac{\left|\left\langle
E_{n}\left|c^{\dagger}_{k}\right|\mathrm{G}\right\rangle\right|^{2}}{\omega+E_{\mathrm{G}}-E_{n}+i\nu}+\frac{\left|\left\langle
E_{n}\left|c_{k}\right|\mathrm{G}\right\rangle\right|^{2}}{\omega-
E_{\mathrm{G}}+E_{n}+i\nu}\right),$ (8)

where $E_{n}$ is the $n^{th}$ excited energy of $H$ corresponding to the
eigenstate $\ket{E_{n}}$, and $E_{G}$ is the ground state energy.

### Green’s function by ADAPT-SSVQE

The SSVQE algorithm can prepare a desired number of excited states on a
quantum computer with a single variational circuit by using the orthogonality
conservation property of a unitary transformationNakanishi2019 ; Endo2020 .
The algorithm results in a variational quantum circuit $U(\vec{\theta})$ which
acts on a set of orthonormal states $\\{\ket{\psi_{j}}_{j=0}^{K-1}\\}$ to
approximate the lowest $K$ eigenstates of a Hamiltonian $H$ i.e.
$\ket{E_{j}}=U(\vec{\theta})\ket{\psi_{j}},0<j<K$. In this paper, we focus on
an adaptive construction of $U(\vec{\theta})$ by using ADAPT-VQE with the
following cost function:

$\mathcal{C}_{0}(\vec{\theta})=\sum_{j=0}^{K-1}w_{j}\left\langle\psi_{j}\left|U(\vec{\theta})^{\dagger}HU(\vec{\theta})\right|\psi_{j}\right\rangle=\sum_{j=0}^{K-1}w_{j}E_{j},$
(9)

where $\vec{\theta}$ represents the variational parameters, and
$\mathrm{w}_{0}>\cdots\mathrm{w}_{K-1}>0$ are the weights which are chosen to
ensure that the eigenenergies are in ascending order i.e. $E_{0}<E_{1}\cdots
E_{K-1}$. Once we have the approximate eigenstates $\ket{E_{n}}$ and
eigenenergies $E_{n}$ for both particle and hole excitations, we can obtain
the transition amplitudes of creation and annihilation operators, and thus use
Lehmann representation to compute the propagator. A detailed description on
how to calculate the transition amplitudes is given in Ref. Endo2020 .

##  IV Results and discussion

###  IV.1 Ground state results

Figure 2: Scaling of the number of parameters required to achieve $0.99$ fidelity with the system size for 1-dimensional Fermi-Hubbard model. Solid lines indicate ADAPT-VQE algorithm whereas dashed line is for non-adaptive algorithm using a number preserving ansatz, taken from Ref.Cade2019 . Figure 3: Scaling of infidelity(1-fidelity) and error in energy ($\Delta E=E-E_{exact}$) with the number of steps for a 3 $\times$ 2 grid. Dashed black line indicates 0.99 fidelity. Notice that the accuracy increases asymptotically with the number of gates.  Grid | U | $\text{E}_{\text{exact}}$ | Depth | Fidelity | $\Delta E$
---|---|---|---|---|---
$2\times 1$ | 1 | -2.56 | 4 | 1.00 | $3.80\times 10^{-12}$
3 | -4.00 | 4 | 1.00 | $1.03\times 10^{-11}$
6 | -6.61 | 4 | 1.00 | $1.15\times 10^{-11}$
$3\times 1$ | 1 | -3.51 | 9 | 1.00 | $1.17\times 10^{-09}$
3 | -5.12 | 9 | 1.00 | $2.81\times 10^{-09}$
6 | -7.85 | 10 | 1.00 | $1.52\times 10^{-10}$
$4\times 1$ | 1 | -5.58 | 19 | 1.00 | $2.90\times 10^{-04}$
3 | -8.35 | 24 | 1.00 | $1.17\times 10^{-04}$
6 | -13.43 | 29 | 1.00 | $1.42\times 10^{-03}$
$2\times 2$ | 1 | -5.34 | 35 | 1.00 | $2.87\times 10^{-08}$
3 | -8.42 | 32 | 1.00 | $6.52\times 10^{-05}$
6 | -13.63 | 32 | 1.00 | $2.19\times 10^{-05}$
$5\times 1$ | 1 | -6.73 | 32 | 1.00 | $7.88\times 10^{-04}$
3 | -9.74 | 53 | 1.00 | $1.50\times 10^{-03}$
6 | -14.99 | 75 | 1.00 | $7.15\times 10^{-04}$
$6\times 1$ | 1 | -8.63 | 51 | 1.00 | $1.56\times 10^{-03}$
3 | -12.72 | 84 | 1.00 | $6.09\times 10^{-03}$
6 | -20.27 | 79 | 0.98 | $2.19\times 10^{-02}$
$3\times 2$ | 1 | -9.28 | 52 | 1.00 | $1.12\times 10^{-03}$
3 | -13.28 | 80 | 0.99 | $2.93\times 10^{-02}$
6 | -20.73 | 71 | 0.88 | $6.74\times 10^{-02}$
Table 1: Final depth, fidelity and error in energy ($\Delta E=E-E_{exact}$)
reached for the adaptive ansatz optimized variationally using the L-BFGS-B
method available in SciPy. Stopping criteria were
$(E^{(n)}-E^{(n-1)})<10^{-3}$ and
$\absolutevalue{\partialderivative{E^{(n)}}{\theta_{i}}}<10^{-4}$, where the
partial derivative is taken with respect to the $i^{th}$ parameter at depth
$n$.

We successfully used the ADAPT-VQE algorithm to create quantum circuits that
produce high fidelity ground states of the Fermi-Hubbard model with various
system sizes up to 6 sites (12 qubits). We explored these systems for
interaction strength($U$) up to 6, and the method works well for the range
tested. Table 1 summarizes the final depth of the variational circuit as well
as achieved fidelity and the error in energy ($\Delta E=E-E_{exact}$) for the
following stopping criteria: $(E^{(n)}-E^{(n-1)})<10^{-3}$ and
$\absolutevalue{\partialderivative{E^{(n)}}{\theta_{i}}}<10^{-4}$. Depth here
refers to the number of fermionic gates constituting the ansatz, which is also
equal to the number of variational parameters. The scaling of the ansatz depth
required to achieve a fidelity of $0.99$ with the system size for a
1-dimensional Fermi-Hubbard model is shown in Fig. 2.

In Ref. Cade2019 , the authors constructed ground states for the Fermi-Hubbard
model by using non-adaptive methods such as the Hamiltonian variational ansatz
and the number-preserving ansatz to reach 0.99 fidelity. Our circuits use
significantly less parameters to reach the same fidelity. For instance, they
report that 160 parameters are required to reach 0.99 fidelity for 3 $\times$
2 grid, $U=2$, and we obtain a better fidelity with just 62 parameters. In
Fig. 2, we compare the scaling of the number of parameters between number
preserving ansatz from Ref.Cade2019  and ADAPT-VQE ansatz for $U=2$, and it
can be seen that ADAPT-VQE scales better. Although for the system sizes
studied, the scaling appears to be exponential, the number preserving ansatz
scales sub-exponentially for larger system sizesCade2019 , so the adaptive
ansatz can be expected to exhibit a similar behavior since it outperforms
number-preserving ansatz in every single study we performed.

We also found that fidelity as well as error in energy improve consistently
with the number of gates, and the wavefunction becomes exact in the limit of
infinite depth. Such scaling for the largest grid that we studied is shown in
Fig. 3. As expected, the number of steps required to produce the ground state
for a given fidelity increases with the strength of the interaction ($U$).
Convergence was not an issue for any of the systems that we studied although
the runtime tended to increase with increased system size, interaction, and
desired accuracy. For those problems in which estimating the ground state
energy of a strongly interacting system is the main goal, reasonable estimates
can be obtained with just a few steps. The first few ($\sim N$) gates in our
optimized ansatz are predominantly hopping terms involving only two qubits,
which further adds to the shallowness of the circuit. In total, it appears
from these results that arbitrary accuracy can be obtained at the cost of
increased gate depth and number of parameters.

The ADAPT-VQE method is particularly well suited for NISQ devices compared to
the non-adaptive Hamiltonian inspired variational ansatz Cade2019 . In our
adaptive method, a step involves adding one (or a specified number) operator
to the ansatz from a given pool. Non-adaptive methods usually require addition
of all the interaction terms present in the Hamiltonian at each step. The
order and type of gates used in the non-adaptive methods is not optimized at
each step, so they are not as flexible as the adaptive methods in terms of
when we want to terminate circuit. Since the NISQ devices are noisy, having a
deeper circuit will add more noise than add to the accuracy of ground state.
The stopping criteria $\epsilon$ makes ADAPT-VQE aware of when exactly adding
a gate stops being advantageous. Since there is always going to be some number
of layers before the stopping criteria are met, progress is always made
towards some understanding of the ground state due to the interpretability of
our ansatz.

It would seem that because we treat the one-body and two-body fermionic
operators as our fundamental gates rather than focus on their implementation
in terms of single-qubit and entangling gates, our approach would be
inefficient on a real NISQ device. For example, the number of such gates
required to perform these operations on quantum hardware would greatly depend
on the specific architecture i.e. connectivity between the qubits and natural
gates available. But for real-device implementation, the ADAPT-VQE algorithm
allows the cost of implementing these fermionic operators to be incorporated
into the cost function. So we can optimize the ansatz not only for the
accuracy of the result but also for the depth of the resulting circuit.

#### Choice of initial state

Figure 4: Comparision of fidelity vs. starting configuration for 2 $\times$
2\. grid with $U=3$. Red line corresponds to configurations 1 and 2 for which
two of the sites are doubly occupied and remaining two are unoccupied whereas
green line corresponds to configurations 3 and 4 for which all sites are
singly occupied.

We performed variational optimization of adaptive ansatz using both non-
interacting ground state (i.e. Slater determinants) and product state as the
initial wavefunction, and found that the product state always produces the
ground state with desired accuracy with less number of total gates. Since the
operator pool is formed by number-preserving operators, it is important to
start in the subspace with desired occupancy, which can be done by applying
X-gates on the number of qubits equal to the occupation number. For all the
grids that we studied, the best choice of the product state was to evenly
spread out the fermions in the grid because of being close to the half
filling. Although with enough steps, all starting distributions can result in
a ground state with arbitrary accuracy, spread-out initial states reach the
desired accuracy faster. One such example for 2 $\times$ 2 grid with $U=3$ is
shown in Fig. 4. We note that the fidelity-scaling is dependent only on the
number of singly-occupied vs. doubly-occupied sites and not on their
permutations. We thus note any prior understanding of the ground state
electron density is rewarded in the form of lower depth.

###  IV.2 Excited states and Green’s function

Figure 5: Spectral functions calculated for Fermi-Hubbard model on a 2
$\times$ 1 grid at half filling for $U=3$ (top) and $U=6$ (bottom) by using
ADAPT-SSVQE method.

We also show that ADAPT-SSVQE (i.e. ADAPT-VQE + SSVQE) is effective at
preparing the excited states of the Fermi-Hubbard model, and hence helps us
predict the spectral properties. We used the SSVQE method with different
weights to compute the ground state for the simple case of 2 $\times$ 1 grid
with $U=3$ and $U=6$ Endo2020 . The ADAPT-SSVQE-computed spectral functions
for these specific cases, along with the ones obtained by exact
diagonalization, are shown in Fig. 5. Ref. Endo2020  reports using 48
parameters to calculate the excited states for calculating the Green’s
function using a non-adaptive ansatz whereas the adaptive algorithm was able
to do so with just 5 parameters. The number of parameters caps the number of
excited states that can be obtained using an ansatz. The effectiveness of the
adaptive variational algorithm to calculate Green’s function for this rather
small system points towards a promising future in terms of simulating
condensed matter systems using quantum computers.

We note that the excitation spectrum of the Hamiltonian is related to the
dynamics, hence the ADAPT-VQE algorithm also has potential application to
variational quantum simulation(VQS)Li_2017 ; Endo2020 .

##  V Conclusion

In this paper, we demonstrate the application of adaptive variational methods
(ADAPT-VQE and proposed ADAPT-SSVQE) to prepare the ground as well as excited
states of the Fermi-Hubbard model of strongly correlated electron systems on
small grids. In contrast with the non-adaptive methods, a system-specific
ansatz is built by adding one or more fermionic operators at a time, resulting
in fewer variational parameters and gate depth. The energy and fidelity of the
prepared states improve asymptotically with the addition of each operator.
Lower depth, consistent convergence, and highly controllable system-specific
ansatz make this method suitable to studying strongly interacting problems in
condensed matter and quantum chemistry on NISQ devices.

## Acknowledgement

The authors would like to thank Nicholas Mayhall for his helpful suggestions
on the manuscript.

## References

  * [1] Harper R. Grimsley, Sophia E. Economou, Edwin Barnes, and Nicholas J. Mayhall.  An adaptive variational algorithm for exact molecular simulations on a quantum computer.  Nature Communications, 10(1):3007, 2019.
  * [2] Alberto Peruzzo, Jarrod McClean, Peter Shadbolt, Man-Hong Yung, Xiao-Qi Zhou, Peter J. Love, Alán Aspuru-Guzik, and Jeremy L. O’Brien.  A variational eigenvalue solver on a photonic quantum processor.  Nature Communications, 5(1):4213, 2014.
  * [3] Jarrod R McClean, Jonathan Romero, Ryan Babbush, and Alán Aspuru-Guzik.  The theory of variational hybrid quantum-classical algorithms.  New Journal of Physics, 18(2):023023, feb 2016.
  * [4] Kunal Sharma, Sumeet Khatri, M Cerezo, and Patrick J Coles.  Noise resilience of variational quantum compiling.  New Journal of Physics, 22(4):043006, apr 2020.
  * [5] L.J. Sham W. Kohn.  Self-consistent equations including exchange and correlation effects.  American Physical Society, 140, 1965.
  * [6] Steven R. White.  Density matrix formulation for quantum renormalization groups.  Phys. Rev. Lett., 69:2863–2866, Nov 1992.
  * [7] Peter W. Shor.  Polynomial-time algorithms for prime factorization and discrete logarithms on a quantum computer.  SIAM Journal on Computing, 26(5):1484–1509, Oct 1997.
  * [8] Michael A. Nielsen and Isaac L. Chuang.  Quantum Computation and Quantum Information: 10th Anniversary Edition.  Cambridge University Press, 2010.
  * [9] Richard P. Feynman.  Simulating physics with computers.  International Journal of Theoretical Physics, 21(6):467–488, 1982\.
  * [10] John Preskill.  Quantum Computing in the NISQ era and beyond.  Quantum, 2:79, August 2018.
  * [11] Jay Gambetta.  Ibm’s roadmap for scaling quantum technology.  https://research.ibm.com/blog/ibm-quantum-roadmap, Feb.
  * [12] Our quantum computing journey.  https://quantumai.google/learn/map.
  * [13] Zi-Jian Zhang, Thi Ha Kyaw, Jakob S Kottmann, Matthias Degroote, and Alán Aspuru-Guzik.  Mutual information-assisted adaptive variational quantum eigensolver.  Quantum Science and Technology, 6(3):035001, Jul 2021.
  * [14] J. Hubbard.  Electron correlations in narrow energy bands.  Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences, 276(1365):238–257, 1963.
  * [15] Bryan O’Gorman, Sandy Irani, James Whitfield, and Bill Fefferman.  Electronic structure in a fixed basis is qma-complete.  2021\.
  * [16] Yudong Cao, Jonathan Romero, Jonathan P. Olson, Matthias Degroote, Peter D. Johnson, Mária Kieferová, Ian D. Kivlichan, Tim Menke, Borja Peropadre, Nicolas P. D. Sawaya, Sukin Sim, Libor Veis, and Alán Aspuru-Guzik.  Quantum chemistry in the age of quantum computing.  Chemical Reviews, 119(19):10856–10915, 10 2019.
  * [17] Suguru Endo, Iori Kurata, and Yuya O. Nakagawa.  Calculation of the green’s function on near-term quantum computers.  Phys. Rev. Research, 2:033281, Aug 2020.
  * [18] Ken M. Nakanishi, Kosuke Mitarai, and Keisuke Fujii.  Subspace-search variational quantum eigensolver for excited states.  Physical Review Research, 1(3), Oct 2019.
  * [19] Dave Wecker, Matthew B. Hastings, and Matthias Troyer.  Progress towards practical quantum variational algorithms.  Phys. Rev. A, 92:042303, Oct 2015.
  * [20] Chris Cade, Lana Mineh, Ashley Montanaro, and Stasja Stanisic.  Strategies for solving the fermi-hubbard model on near-term quantum computers.  Phys. Rev. B, 102:235122, Dec 2020.
  * [21] Nicholas C Rubin, Toru Shiozaki, Kyle Throssell, Garnet Kin-Lic Chan, and Ryan Babbush.  The fermionic quantum emulator.  arXiv preprint arXiv:2104.13944, 2021.
  * [22] P. Jordan and E. Wigner.  Über das paulische äquivalenzverbot.  Zeitschrift für Physik, 47(9):631–651, 1928.
  * [23] Pauli Virtanen et al.  SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python, 2020.
  * [24] P. Coleman.  Introduction to Many-Body Physics.  Cambridge University Press, 2015.
  * [25] G.D. Mahan.  Many-Particle Physics.  Physics of Solids and Liquids. Springer US, 1990.
  * [26] Ying Li and Simon C. Benjamin.  Efficient variational quantum simulator incorporating active error minimization.  Phys. Rev. X, 7:021050, Jun 2017.

##  Appendix A Jordan-Wigner transform

We used the Jordan-Wigner transform[22] to map the fermionic creation and
annihilation operators into Pauli operators (i.e. tensor product of Pauli
matrices). The Jordan-Wigner transform is used to map a $n-$fermionic system
to a $n-$qubit system as follows:

$\displaystyle a_{k}$
$\displaystyle\mapsto\left(\prod_{j=0}^{k-1}Z_{j}\right)\sigma_{k}^{-}$ (10)
$\displaystyle=\frac{1}{2}\left(Z_{0}\otimes\cdots\otimes
Z_{k-1}\otimes\left(X_{k}-iY_{k}\right)\right),$

where $X_{j},Y_{j},Z_{j}$ are the Pauli spin operators, and $\sigma_{j}^{+/-}$
are the ladder operators acting on qubit $j$. By virtue of this transform, the
qubit operators, which follow commutation relations, are mapped to
creation/annihilation operators, which follow anticommutation relations.

##  Appendix B Choice of Fillings

Number of electrons | Grid
---|---
2 | 2$\times$ 1, 3 $\times$ 1
4 | 4 $\times$ 1, 2 $\times$ 2, 5 $\times$ 1
6 | 6 $\times$ 1, 3 $\times$ 2
Table 2: Number of electrons/fillings and the corresponding grid sizes chosen
such that the Fermi-Hubbard Hamiltonian has non-degenerate ground states

The fillings were chosen such that the Fermi-Hubbard Hamiltonian has a non-
degenerate ground state. Degeneracy gives rise to ambiguity in the choice of
ground state to use for calculating the fidelity.