11institutetext: Department of Pathology, Brigham and Women’s Hospital,
Harvard Medical School, Boston, MA
Department of Pathology, Massachusetts General Hospital, Harvard Medical
School, Boston, MA
11email: {icarmichael, asong, rchen14, dwilliamson, tchen25,
<EMAIL_ADDRESS>
# Incorporating intratumoral heterogeneity into weakly-supervised deep
learning models via variance pooling

Iain Carmichael These authors contributed equally Andrew H. Song* Richard J.
Chen  Drew F.K. Williamson  Tiffany Y. Chen  Faisal Mahmood

###### Abstract

Supervised learning tasks such as cancer survival prediction from gigapixel
whole slide images (WSIs) are a critical challenge in computational pathology
that requires modeling complex features of the tumor microenvironment. These
learning tasks are often solved with deep multi-instance learning (MIL) models
that do not explicitly capture intratumoral heterogeneity. We develop a novel
variance pooling architecture that enables a MIL model to incorporate
intratumoral heterogeneity into its predictions. Two interpretability tools
based on “representative patches” are illustrated to probe the biological
signals captured by these models. An empirical study with 4,479 gigapixel WSIs
from the Cancer Genome Atlas shows that adding variance pooling onto MIL
frameworks improves survival prediction performance for five cancer types.

###### Keywords:

Computer Vision, Computational Pathology, Multiple Instance Learning,
Intratumoral Heterogeneity, Interpretability

##  1 Introduction

Cancer diagnosis and prognosis using hematoxylin and eosin (H&E) stained whole
slide images (WSIs) are central tasks in computational pathology. These
supervised learning problems are made difficult by the size and complexity of
WSIs; these images can be gigapixels in scale and involve localizing subtle
visual signals that expert pathologists are trained for years to identify. In
recent years, weakly supervised deep learning algorithms have made significant
progress in addressing some of these challenges [Coudray et al., 2018, Couture
et al., 2018, Campanella et al., 2019, Courtiol et al., 2019, Wulczyn et al.,
2021, Lu et al., 2021].

Cancer cell morphology can vary greatly within a single tumor and it is
increasingly recognized that intratumoral heterogeneity (ITH) plays an
important role in prognostication for certain cancer types [Andor et al.,
2016, Marusyk et al., 2020], though its exact role is still being elucidated.
Furthermore, ITH complicates prognostication based on WSIs by increasing the
range of morphologies that correspond to a particular prognosis. For example,
pathologists often encounter clear cell renal cell carcinoma displaying
significant heterogeneity within the same tumor, with some areas demonstrating
classic clear cell morphology and others displaying morphologies bordering on
diagnoses such as chromophobe renal cell carcinoma or clear cell
tubulopapillary renal cell carcinoma.

Existing supervised algorithms that scale to WSIs do not generally incorporate
histological ITH. These algorithms typically approach a supervised learning
task such as survival prediction through a multiple instance learning (MIL)
framework [Ilse et al., 2018]. First, a WSI is represented as a bag of many
smaller image patches. Next, neural network features are extracted from each
patch. Finally, these patch features are aggregated for prediction. The
standard MIL aggregation approaches capture only simple statistics (e.g. mean
or max) and neglect higher-order information (e.g. variance), which can
capture heterogeneity.

To this end, we develop a novel MIL aggregation component that scales to
gigapixel WSIs and explicitly incorporates ITH into its predictions. We
accomplish this through a variance pooling operation that quantifies ITH as
the variance along a collection of low rank projections of patch features
(Section 2.2). Moving beyond pre-specified measures of heterogeneity (e.g. the
variance of tumor cell morphology), this architecture uses the data to
determine how to measure heterogeneity (e.g. the projection vectors, along
which the variance of the features are measured, are learned.) This operation
can be seamlessly incorporated into existing MIL architectures, as we now
demonstrate.

We demonstrate empirically in Section 3 that including this operation improves
baseline MIL models’ performance on survival prediction tasks across five
cancer types from The Cancer Genome Atlas (TCGA). Through attention and a
novel “variance projection contrast” visualization, we provide interpretable
insights into what biological signals these architectures utilize to make the
predictions in Section 4. Code: https://github.com/mahmoodlab/VARPOOL.

###  1.1 Related Work

Survival Analysis using WSIs: Survival analysis in histopathology has been a
longstanding problem in medical imaging. Early work required manual
identification of regions of interest [Gurcan et al., 2009]. Other methods
make use of hand-crafted tumor cell features (e.g. shape) to predict patient
survival [Lu et al., 2018]. Most recent survival methods use MIL algorithms
that are directly trained on WSIs without additional annotations [Zhu et al.,
2017, Mobadersany et al., 2018, Wulczyn et al., 2021, Zhao et al., 2020]. Some
of these methods use graph convolutional networks (GCNs), which learn more
context-aware features [Zhao et al., 2020, Chen et al., 2021].

Intratumoral Heterogeneity: Much of the previous work on ITH uses genomic data
to identify and quantify subclones of tumor cells that share mutations
[Raynaud et al., 2018]. There is growing interest in using machine learning to
quantify the histologic ITH that pathologists observe on a daily basis [Faust
et al., 2020].

##  2 Incorporating heterogeneity with variance pooling

This section presents the novel variance pooling operation for weakly-
supervised MIL tasks. In this setting, we observe a set of bags – an unordered
collection of vectors (instances) – that we use to predict some bag level
response. In computational pathology, each patient has a WSI which is used to
predict the response (e.g. class label or survival outcome). Each patent’s
WSI111When a patient has multiple WSIs, the patches are unioned across the
WSIs. (the bag) is broken up into many small image patches and the instances
are either the raw image patches or patch features extracted by a neural
network.

We first define some notation. Suppose we have $N$ patients. For the $i$th
patient, we observe a response $y^{(i)}$ and a bag of instances (image patch
features) $\mathcal{X}^{(i)}=\\{x^{(i)}_{1},\ldots,x^{(i)}_{n^{(i)}}\\}$,
where each $x^{(i)}_{j}\in\mathbb{R}^{d}$ and $n^{(i)}$ is the size of the
bag. For classification tasks $y^{(i)}$ is a patient-level class label; for
survival tasks $y^{(i)}=(t^{(i)},c^{(i)})$ where $t^{(i)}\in\mathbb{R}_{+}$ is
the observed survival time and $c^{(i)}\in\\{0,1\\}$ indicates whether or not
this time was censored. We want to learn a mapping
$f(\mathcal{X}^{(i)})=z^{(i)}$ that outputs a patient-level prediction
minimizing a supervised loss function.

In the following two sections, we present the variance pool architecture in
the context of the attention mean pooling model of Ilse et al. [2018]. We
emphasize variance pooling can be easily incorporated into other MIL
architectures. For notational simplicity we generally drop the patient
superscript $i$.

Figure 1: Variance pooling framework. Zoomed-in WSI patches show examples of
intratumoral heterogeneity, with the colonic adenocarcinoma displaying
multiple moderately- to poorly-differentiated morphologies. The mean and
variance pooling branches capture different statistical characteristics of the
data and are aggregated to yield the final prediction.

###  2.1 Attention mean pooling architecture

The attention mean pool architecture of Ilse et al. [2018] reduces a bag of
instances, $\\{x_{j}\\}_{j=1}^{n}$, to a vector through a weighted mean
operation. Specifically, each instance is first passed through an (optional)
encoding network, as $h_{j}=\phi(x_{j})\in\mathbb{R}^{d_{\text{enc}}}$. The
attention neural network $\alpha$, takes in $\\{h_{j}\\}_{j=1}^{n}$ and
computes the normalized attention weights
$a=\alpha\left(\\{h_{j}\\}_{j=1}^{n}\right)\in\mathbb{R}^{n}_{+}$, which sum
to 1. The attention mean pool vector $p_{\text{mean}}$ is obtained as

$p_{\text{mean}}=\sum_{j=1}^{n}a_{j}h_{j}.$ (1)

Finally, the bag prediction (e.g. predicted survival risk) is obtained by
passing $p_{\text{mean}}$ through an additional multi-layer perceptron (MLP),
$z=\rho(p_{\text{mean}})$.

###  2.2 Attention variance pooling

We incorporate within-bag heterogeneity (intratumoral heterogeneity) by
computing the variance of learned projections of the instances. First, let us
define the attention variance function,
$\text{Var}\left(\cdot;a\right):\mathbb{R}^{n}\to\mathbb{R}$, of $n$ numbers,
$u\in\mathbb{R}^{n}$, as

$\text{Var}\left(u;a\right):=\sum_{j=1}^{n}a_{j}\left(u_{j}-\sum_{\ell=1}^{n}a_{\ell}u_{\ell}\right)^{2}$
(2)

with attention weights $a\in\mathbb{R}^{n}_{+}$ summing to one. Note, if the
attention weights are all equal this is the empirical variance of the entries
of $u$.

We then define a $K$-dimensional variance pool of $\\{h_{j}\\}_{j=1}^{n}$ as
follows. Let $v_{k}\in\mathbb{R}^{d_{\text{enc}}},k\in[K]:=\\{1,\dots,K\\},$
be variance projection vectors, which are learnable parameters in the network.
Denote the matrix of encoded instance features as $H\in\mathbb{R}^{n\times
d_{\text{enc}}}$, i.e. the $j$th row of $H$ is given by $h_{j}$. The $k$th
entry222Without the attention weights this would be the quadratic form
$v_{k}^{T}\text{cov}(H)v_{k}$. of the variance pool vector
$p_{\text{var}}\in\mathbb{R}^{K}_{+}$ is obtained as

$p_{\text{var},k}=\text{Var}\left(Hv_{k};\alpha\right)=\sum_{j=1}^{n}a_{j}\left(h_{j}^{T}v_{k}-\sum_{\ell=1}^{n}a_{\ell}h_{\ell}^{T}v_{k}\right)^{2}.$
(3)

This can be viewed as dimensionality reduction on the (attention weighted)
covariance matrix of $H$, i.e., from $O(d_{\text{enc}}^{2})$ to $K$.

We pass $p_{\text{var}}$ through an entrywise non-linearity, $\eta(\cdot)$,
such as $\sqrt{\cdot}$ or $\log(\epsilon+\cdot)$. This non-linearity operation
is required since linear combinations of several variance pools is equivalent
to just a single variance pool. Finally, the first and second moment
information of the attended $\\{h_{j}\\}_{j=1}^{n}$ are combined and the
prediction is obtained by passing the concatenated mean and variance pool
vector,

$p_{\text{cat}}=\left[p_{\text{mean}};\eta(p_{\text{var}})\right]\in\mathbb{R}^{d_{\text{enc}}+K}$
(4)

through the MLP, $z=\rho(p_{\text{cat}})$333A similar approach, which we were
unaware of until after publication, has been successfully used for
microsatellite instability (MSI) prediction Schirris et al. [2022]. This
approach is most related to variance pooling with K=d and the variance
projection coefficients fixed to be the standard basis vectors. .

Figure 1 shows the entire architecture. The learnable parameters are the
networks $\\{\phi,\alpha,\rho\\}$ and $\\{v_{k}\\}_{k=1}^{K}$. Several design
choices can be modified including: no attention mechanism, separate attention
mechanisms for the mean and variance pools, and different combination schemes
for the mean/variance pools.

##  3 Experiments with survival prediction

Data: From the public TCGA dataset, we use: Breast Invasive Carcinoma (BRCA,
$N=1,061$), Glioblastoma & Lower Grade Glioma (GBMLGG, $N=872$), Uterine
Corpus Endometrial Carcinoma (UCEC, $N=504$), Colon Adenocarcinoma & Rectum
Adenocarcinoma (COADREAD, $N=612$), and Bladder Urothelial Carcinoma (BLCA,
$N=386$). The choice of these types was based on the number of patients. We
use progression-free interval (PFI) [Liu et al., 2018].

We use the concordance index (c-Index) to evaluate survival prediction models.
The c-Index, which lies between 0 and 1 (larger is better), measures the
agreement between predicted risk scores and actual survival outcomes using
only order information (see discussion around (5) below.) We split each
dataset into a $70\%/30\%$ train/test split. To ensure that our results are
not sensitive to a particular data split, we randomly split the data into
train/test folds 10 times.

Implementation Details: We extract $d=1,024$ features from each WSI image
patch ($256\times 256$ pixels at $20\times$ magnification) using an ImageNet-
pretrained ResNet-50 CNN encoder, truncated after the third residual block and
spatially averaged. We use a batch size of $B=32$ patients444Note the patient
batch size is not equal to the number of summands in (5); there are between 0
and ${B\choose 2}$ summands depending on the number of comparable pairs.. To
fit each mini-batch into GPU memory during training we cap–by random
subsampling– the number of instances in each bag. For each cancer type, the
maximum number of instances per bag is set to be the $75\%$ quantile of the
overall bag sizes (WSIs typically have $15-20,000$ patches). In the
Supplemental information, we report the summary statistics for number of
instances per WSI slide in the TCGA dataset. All patches in the bag are used
at test time. We use the Adam optimizer with a learning rate of $2\times
10^{-4}$, weight decay of $1\times 10^{-5}$, and train for 30 epochs.

We use three MIL architectures: Deep Sets [Zaheer et al., 2017], Attention
MeanPool [Ilse et al., 2018], and DeepGraphConv [Zhao et al., 2020]. For the
Deep Sets framework, the attention weights are not learned and instead set to
identical values of $1/n$. For the DeepGraphConv framework, the bag of
features is first passed through two layers of graph convolutional neural
network layers, the output of which are passed through the attention
mean/variance pooling modules. These are some of the most popular MIL
architectures used in computational pathology and thus provide great testbeds
for variance pooling module extension.

We use $K=10$ variance pooling projections and an
$\eta(\cdot)=\log(0.01+\cdot)$ nonlinearity. We performed an internal ablation
study for the nonlinearity comparing the logarithm, square root, and sigmoid
functions, but did not find any significant difference in performance. We
initialize the projection vectors $v_{k}\sim\mathcal{N}(0,1/d_{\text{enc}})$
to make them approximately orthonormal with high probability.

Loss function: Suppose the bag prediction network, $f(\cdot)$, outputs a risk
score for each patient. We train the network with the survival ranking loss
[Steck et al., 2007, Luck et al., 2018], which is a continuous approximation
of the negative c-Index given by

$\mathcal{L}\left(\\{\mathcal{X}^{(i)},t^{(i)},c^{(i)}\\}_{i=1}^{N}\right)=-\frac{1}{|\mathcal{C}|}\sum_{(w,b)\in\mathcal{C}}\psi\left(f(\mathcal{X}^{(w)})-f(\mathcal{X}^{(b)})\right),$
(5)

where $\mathcal{C}=\left\\{(w,b)|c^{(w)}=0\text{ and }t^{(w)}<t^{(b)},\text{
for }w,b\in[N]\right\\}$ is the set of comparable pairs i.e. pairs of patients
where we know one patient has worse survival than the other. If
$\psi(\cdot)=\mathbf{1}(\cdot>0)$ then (5) would be exactly the negative
c-Index. We set $\psi(\cdot)$ to be a sigmoid, which is a differentiable
approximation of this indicator.

We chose the ranking loss because it is separable, which allows the use of
proper mini-batches unlike the Cox loss. Additionally, the c-Index is a
natural objective for survival analysis. An alternative is the negative log-
likelihood of a discrete time-to-event model [Zadeh and Schmid, 2020]. This
loss function coarsely bins the survival times and does not fully utilize the
available survival ordering information.

Table 1: Test set c-Index performance (mean $\pm$ std) across 10 re-sampled data splits on the 100 point scale. Frameworks |  BRCA |  GBMLGG |  COADREAD
---|---|---|---
Deep Sets |  $59.5\pm 3.1$ |  $72.8\pm 1.0$ |  $50.8\pm 4.1$
Deep Sets w/ VarPool |  $\mathbf{60.0\pm 3.8}$ |  $\mathbf{73.8\pm 0.9}$ |  $\mathbf{56.1\pm 3.0}$
Attn MeanPool |  $59.2\pm 3.7$ |  $73.0\pm 1.2$ |  $52.0\pm 2.3$
Attn MeanPool w/ VarPool |  $\mathbf{60.4\pm 3.9}$ |  $\mathbf{73.6\pm 1.0}$ |  $\mathbf{56.3\pm 2.3}$
DeepGraphConv |  $57.9\pm 4.2$ |  $70.7\pm 2.3$ |  $51.4\pm 3.0$
DeepGraphConv w/ VarPool |  $\mathbf{58.0\pm 5.2}$ |  $\mathbf{70.9\pm 2.3}$ |  $\mathbf{51.6\pm 3.2}$
|  BLCA |  UCEC |  Overall
Deep Sets |  $56.6\pm 5.3$ |  $51.3\pm 4.3$ |  $58.2$
Deep Sets w/ VarPool |  $\mathbf{57.5\pm 5.1}$ |  $\mathbf{54.1\pm 3.0}$ |  $\mathbf{60.3}$
Attn MeanPool |  $57.2\pm 5.3$ |  $50.7\pm 3.9$ |  $58.4$
Attn MeanPool w/ VarPool |  $\mathbf{57.8\pm 4.6}$ |  $\mathbf{52.5\pm 3.7}$ |  $\mathbf{60.1}$
DeepGraphConv |  $52.1\pm 7.7$ |  $54.2\pm 2.1$ |  $57.3$
DeepGraphConv w/ VarPool |  $\mathbf{53.1\pm 7.6}$ |  $\mathbf{55.0\pm 2.0}$ |  $\mathbf{57.7}$

Table 1 shows the results. We observe that for every cancer type and all
models, the variance pooling versions outperform their counterparts, with up
to $10.4\%$ relative increase (Deep Sets for COADREAD). We conjecture that the
smaller performance gain in DeepGraphConv comes from “feature smoothing” in
GCNs, which reduces the heterogeneity in patch features. In the Supplemental
information, we present an ablation study with the Cox loss and show that
similar improvements are observed with the addition of variance pooling
module.

The fact that explicitly incorporating intratumoral heterogeneity gives the
greatest benefit in COADREAD and UCEC is interesting, since subsets of these
tumor types have very high mutation rates due to microsatellite instability, a
molecular change with impacts on survival [Levine, 2013, Guinney et al.,
2015]. This might lead to the development of multiple subclones of tumor cells
that have different morphologies, as the effects of those mutations propagate
from DNA to protein.

##  4 Interpretability and biological insights

This section presents two interpretability approaches to probe what biological
signals the trained models learn. We focus on the Attention MeanPool and
VarPool architectures, though similar approaches can be used for other
architectures.

Interpretability via attention scores For each patient, we visualize the
patches with the largest attention scores, i.e., the $a$ in (1) and (3), as
depicted in Figure 2. Note that this interpretability approach is not possible
for the Deep Sets framework, since the attention weights are set to be
uniform.

Interpretability via variance projection scores For variance pooling, we can
examine signals captured by each variance pool projection. Each projection can
be viewed as sorting the patches of a WSI, the biological signal of which we
aim to uncover. We project the patches onto a given projection vector and
visualize the patches along the spectrum. A similar interpretability approach
was taken in Carmichael et al. [2021]. Following (3), for a projection vector
$v$ and $j$th instance, we define the signed attention weighted squared
residuals (SAsqR) as

$\operatorname{}{SAsqR}_{j}=\text{sign}(r_{j})a_{j}r_{j}^{2},\text{ where
}r_{j}:=h_{j}^{T}v-\sum_{\ell=1}^{n}a_{\ell}h_{\ell}^{T}v.$

These values can be used to sort the patches e.g. to identify the patches with
the most extreme negative/positive SAsqR values.

###  4.1 Interpretability visualizations

We illustrate the two interpretability tools discussed in the previous
section. Our analysis outputs a large number of images based on these tools.
The figures below illustrate the trends that were reviewed by two
pathologists.

Figure 2 compares the top attended patches for MeanPool and VarPool in one
COADREAD patient. We observed that MeanPool attends to more compact, dense,
and intact tumor with surrounding stroma, while VarPool attends to fragmented
tumor, often with surrounding extracellular mucin and debris.

(a) Top MeanPool patches.

(b) Top VarPool patches.

Figure 2:  The highest attended patches for the attention MeanPool and VarPool
for a single COADREAD patient. Here, MeanPool typically attends to solid tumor
while VarPool attends to fragmented, more poorly-differentiated tumor. The
scale bar represents $100\,\mu m$.

Figure 3:  Patches along one variance projection for two COADREAD patients.
The patches are ordered from the most negative (left) SAsqR scores to most
positive (right). The leftmost patches contain dense inflammatory infiltrate
comprised mostly of lymphocytes, while the rightmost patches show denser tumor
regions.

Figure 4:  Patches along one variance projection for two UCEC patients. This
projection shows two patterns that are illustrated by these patients. For the
first patient, the leftmost patches display necrosis and the right most
patches display tumor. For the second patient, the leftmost patches display
predominantly muscle and connective tissue and the right most patches display
tumor.

Figure 3 shows patches along a single variance projection for two COADREAD
patients. Figure 4 is similar for UCEC. We have sorted the patches by their
SQsqR values and shown several patches for 11 quantiles: 0, 10, …, 100. We
might expect the rightmost and leftmost patches to represent visual contrast.
We observe dense tumor region on one end of the spectrum and other micro-
environment features, such as lymphocytes, muscle, or necrosis on the other
end.

##  5 Conclusion

We introduce a novel variance pooling module that can be incorporated into
existing multiple instance learning (MIL) frameworks, to encode heterogeneity
information through the second order statistics of low rank projections. We
show that variance pooling can improve survival prediction across several
cancer types. We also provide interpretability tools based on attention and
projection scores to understand the heterogeneous biological signals captured
by our framework. We hope this study will lead to further investigations of
intratumoral heterogeneity – and other tumor microenvironment features – for
computational pathology.

## Acknowledgement

We thank Katherine Hoadley for helpful suggestions. This work was supported in
part by internal funds from BWH Pathology, NIGMS R35GM138216 (F.M.), BWH
President’s Fund, MGH Pathology, BWH Precision Medicine Program, Google Cloud
Research Grant, Nvidia GPU Grant Program and funding from the Fredrick
National Lab. R.C. was additionally supported by the NSF graduate research
fellowship. T.Y.C. was additionally funded by the NIH National Cancer
Institute (NCI) Ruth L. Kirschstein National Service Award, T32CA251062. The
content is solely the responsibility of the authors and does not reflect the
official views of the NIH, NIGMS, NCI, or NSF.

## References

  * Andor et al. [2016] Noemi Andor, Trevor A Graham, Marnix Jansen, Li C Xia, C Athena Aktipis, Claudia Petritsch, Hanlee P Ji, and Carlo C Maley.  Pan-cancer analysis of the extent and consequences of intratumor heterogeneity.  _Nature medicine_ , 22(1):105–113, 2016.
  * Campanella et al. [2019] Gabriele Campanella, Matthew G Hanna, Luke Geneslaw, Allen Miraflor, Vitor Werneck Krauss Silva, Klaus J Busam, Edi Brogi, Victor E Reuter, David S Klimstra, and Thomas J Fuchs.  Clinical-grade computational pathology using weakly supervised deep learning on whole slide images.  _Nature Medicine_ , 2019.
  * Carmichael et al. [2021] Iain Carmichael, Benjamin C Calhoun, Katherine A Hoadley, Melissa A Troester, Joseph Geradts, Heather D Couture, Linnea Olsson, Charles M Perou, Marc Niethammer, Jan Hannig, et al.  Joint and individual analysis of breast cancer histologic images and genomic covariates.  _The Annals of Applied Statistics_ , 15(4):1697–1722, 2021.
  * Chen et al. [2021] Richard J Chen, Ming Y Lu, Muhammad Shaban, Chengkuan Chen, Tiffany Y Chen, Drew FK Williamson, and Faisal Mahmood.  Whole slide images are 2d point clouds: Context-aware survival prediction using patch-based graph convolutional networks.  In _MICCAI_ , pages 339–349. Springer, 2021.
  * Coudray et al. [2018] Nicolas Coudray, Paolo Santiago Ocampo, Theodore Sakellaropoulos, Navneet Narula, Matija Snuderl, David Fenyö, Andre L Moreira, Narges Razavian, and Aristotelis Tsirigos.  Classification and mutation prediction from non–small cell lung cancer histopathology images using deep learning.  _Nature medicine_ , 24(10):1559–1567, 2018.
  * Courtiol et al. [2019] Pierre Courtiol, Charles Maussion, Matahi Moarii, Elodie Pronier, Samuel Pilcer, Meriem Sefta, Pierre Manceron, Sylvain Toldo, Mikhail Zaslavskiy, Nolwenn Le Stang, et al.  Deep learning-based classification of mesothelioma improves prediction of patient outcome.  _Nature medicine_ , 25(10):1519–1525, 2019.
  * Couture et al. [2018] Heather D Couture, Lindsay A Williams, Joseph Geradts, Sarah J Nyante, Ebonee N Butler, JS Marron, Charles M Perou, Melissa A Troester, and Marc Niethammer.  Image analysis with deep learning to predict breast cancer grade, er status, histologic subtype, and intrinsic subtype.  _NPJ breast cancer_ , 4(1):1–8, 2018.
  * Faust et al. [2020] Kevin Faust, Adil Roohi, Alberto J Leon, Emeline Leroux, Anglin Dent, Andrew J Evans, Trevor J Pugh, Sangeetha N Kalimuthu, Ugljesa Djuric, and Phedias Diamandis.  Unsupervised resolution of histomorphologic heterogeneity in renal cell carcinoma using a brain tumor–educated neural network.  _JCO Clinical Cancer Informatics_ , 4:811–821, 2020.
  * Guinney et al. [2015] Justin Guinney, Rodrigo Dienstmann, Xin Wang, Aurélien De Reynies, Andreas Schlicker, Charlotte Soneson, Laetitia Marisa, Paul Roepman, Gift Nyamundanda, Paolo Angelino, et al.  The consensus molecular subtypes of colorectal cancer.  _Nature medicine_ , 21(11):1350–1356, 2015.
  * Gurcan et al. [2009] Metin N Gurcan, Laura E Boucheron, Ali Can, Anant Madabhushi, Nasir M Rajpoot, and Bulent Yener.  Histopathological image analysis: A review.  _IEEE reviews in biomedical engineering_ , 2:147–171, 2009\.
  * Ilse et al. [2018] Maximilian Ilse, Jakub Tomczak, and Max Welling.  Attention-based deep multiple instance learning.  In _Proceedings of the 35th International Conference on Machine Learning_. PMLR, 2018.
  * Levine [2013] Douglas A Levine.  Integrated genomic characterization of endometrial carcinoma.  _Nature_ , 497(7447):67–73, 2013.
  * Liu et al. [2018] Jianfang Liu, Tara Lichtenberg, Katherine A Hoadley, Laila M Poisson, Alexander J Lazar, Andrew D Cherniack, Albert J Kovatich, Christopher C Benz, Douglas A Levine, Adrian V Lee, et al.  An integrated tcga pan-cancer clinical data resource to drive high-quality survival outcome analytics.  _Cell_ , 173(2):400–416, 2018.
  * Lu et al. [2018] Cheng Lu, David Romo-Bucheli, Xiangxue Wang, Andrew Janowczyk, Shridar Ganesan, Hannah Gilmore, David Rimm, and Anant Madabhushi.  Nuclear shape and orientation features from h&e images predict survival in early-stage estrogen receptor-positive breast cancers.  _Laboratory investigation_ , 98(11):1438–1448, 2018.
  * Lu et al. [2021] Ming Y Lu, Drew FK Williamson, Tiffany Y Chen, Richard J Chen, Matteo Barbieri, and Faisal Mahmood.  Data-efficient and weakly supervised computational pathology on whole-slide images.  _Nature biomedical engineering_ , 5(6):555–570, 2021.
  * Luck et al. [2018] Margaux Luck, Tristan Sylvain, Joseph Paul Cohen, Heloise Cardinal, Andrea Lodi, and Yoshua Bengio.  Learning to rank for censored survival data.  _arXiv preprint arXiv:1806.01984_ , 2018.
  * Marusyk et al. [2020] Andriy Marusyk, Michalina Janiszewska, and Kornelia Polyak.  Intratumor heterogeneity: the rosetta stone of therapy resistance.  _Cancer cell_ , 37(4):471–484, 2020.
  * Mobadersany et al. [2018] Pooya Mobadersany, Safoora Yousefi, Mohamed Amgad, David A Gutman, Jill S Barnholtz-Sloan, José E Velázquez Vega, Daniel J Brat, and Lee AD Cooper.  Predicting cancer outcomes from histology and genomics using convolutional networks.  _Proceedings of the National Academy of Sciences_ , 115(13):E2970–E2979, 2018.
  * Raynaud et al. [2018] Franck Raynaud, Marco Mina, Daniele Tavernari, and Giovanni Ciriello.  Pan-cancer inference of intra-tumor heterogeneity reveals associations with different forms of genomic instability.  _PLoS genetics_ , 14(9):e1007669, 2018.
  * Schirris et al. [2022] Yoni Schirris, Efstratios Gavves, Iris Nederlof, Hugo Mark Horlings, and Jonas Teuwen.  Deepsmile: Contrastive self-supervised pre-training benefits msi and hrd classification directly from h&e whole-slide images in colorectal and breast cancer.  _Medical Image Analysis_ , 79:102464, 2022.
  * Steck et al. [2007] Harald Steck, Balaji Krishnapuram, Cary Dehing-Oberije, Philippe Lambin, and Vikas C Raykar.  On ranking in survival analysis: Bounds on the concordance index.  _Advances in neural information processing systems_ , 20, 2007.
  * Wulczyn et al. [2021] Ellery Wulczyn, David F Steiner, Melissa Moran, Markus Plass, Robert Reihs, Fraser Tan, Isabelle Flament-Auvigne, Trissia Brown, Peter Regitnig, Po-Hsuan Cameron Chen, et al.  Interpretable survival prediction for colorectal cancer using deep learning.  _NPJ digital medicine_ , 4(1):1–13, 2021.
  * Zadeh and Schmid [2020] Shekoufeh Gorgi Zadeh and Matthias Schmid.  Bias in cross-entropy-based training of deep survival networks.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 43(9):3126–3137, 2020.
  * Zaheer et al. [2017] Manzil Zaheer, Satwik Kottur, Siamak Ravanbakhsh, Barnabas Poczos, Russ R Salakhutdinov, and Alexander J Smola.  Deep sets.  In _Advances in Neural Information Processing Systems_ , volume 30. Curran Associates, Inc., 2017.
  * Zhao et al. [2020] Yu Zhao, Fan Yang, Yuqi Fang, Hailing Liu, Niyun Zhou, Jun Zhang, Jiarui Sun, Sen Yang, Bjoern Menze, Xinjuan Fan, et al.  Predicting lymph node metastasis using histopathological images based on multiple instance learning with deep graph convolution.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 4837–4846, 2020.
  * Zhu et al. [2017] Xinliang Zhu, Jiawen Yao, Feiyun Zhu, and Junzhou Huang.  Wsisa: Making survival prediction from whole slide histopathological images.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , pages 7234–7242, 2017.

## Appendix

(a) Top MeanPool patches.

(b) Top VarPool patches.

Figure 5:  The highest attended patches for the attention MeanPool and VarPool
for a single GBMLGG patient. MeanPool attends to patches with high tumor
cellularity, while VarPool attends to patches with lower cellularity.

Figure 6:  Patches along one variance projection for two GBMLGG patients. As in the main body, the patches are ordered from the most negative (left) SAsqR scores to most positive (right). The top patient is the patient predicted to have a high risk score, while the bottom patient is predicted to have a low risk score. This projection sorts the patches from lower cellularity (right) to higher cellularity (left) – recall we focus on the right/left extremes.  Table 2: Loss function ablation study. Test set c-Index performance (mean $\pm$ std) on the 100 point scale for cox loss, where the training details were the same as the rank loss experiment. This shows that the variance pooling mechanism is loss-agnostic, equally showing improvement for both rank and cox loss. Frameworks |  BRCA |  GBMLGG |  COADREAD
---|---|---|---
Deep Sets |  $60.34\pm 4.45$ |  $72.76\pm 1.13$ |  $54.01\pm 2.6$
Deep Sets w/ VarPool |  $\mathbf{60.96\pm 5.41}$ |  $\mathbf{73.88\pm 1.26}$ |  $\mathbf{57.95\pm 3.01}$
Attn MeanPool |  $59.89\pm 4.51$ |  $\mathbf{73.25\pm 1.16}$ |  $54.71\pm 2.55$
Attn MeanPool w/ VarPool |  $\mathbf{61.1\pm 4.48}$ |  $73.21\pm 1.18$ |  $\mathbf{58.49\pm 3.12}$
DeepGraphConv |  $57.62\pm 4.95$ |  $\mathbf{71.26\pm 2.3}$ |  $56.03\pm 2.59$
DeepGraphConv w/ VarPool |  $\mathbf{59.78\pm 5.2}$ |  $71.15\pm 2.34$ |  $\mathbf{56.23\pm 3.3}$
|  BLCA |  UCEC |  Overall
Deep Sets |  $56.82\pm 4.31$ |  $52.05\pm 4.35$ |  $59.20$
Deep Sets w/ VarPool |  $\mathbf{57.59\pm 3.84}$ |  $\mathbf{55.14\pm 3.43}$ |  $\mathbf{61.10}$
Attn MeanPool |  $56.95\pm 4.26$ |  $51.16\pm 2.85$ |  $59.19$
Attn MeanPool w/ VarPool |  $\mathbf{57.47\pm 3.43}$ |  $\mathbf{53.53\pm 3.49}$ |  $\mathbf{60.76}$
DeepGraphConv |  $52.35\pm 7.1$ |  $51.15\pm 1.89$ |  $57.68$
DeepGraphConv w/ VarPool |  $\mathbf{52.52\pm 5.96}$ |  $\mathbf{53.2\pm 2.33}$ |  $\mathbf{58.58}$
Table 3: The summary statistics for number of instances per WSI slide in the TCGA dataset. For training, we restricted the number of instances per slide to the $75\%$ quantile for each cancer type. To fit a batch in one tensor, zero-padding and manually zeroing the attention was required for slides with fewer instances. For slides with more instances, we performed random subsampling. These measures were necessary to perform training with minibatches of slides (instead of batch size of 1 slide, as typically done in the literature).  |  BRCA |  GBMLGG |  COADREAD |  BLCA |  UCEC
---|---|---|---|---|---
mean |  $10,826.28$ |  $7,448.28$ |  $8.575.54$ |  $16,599.11$ |  $16,003.92$
median |  $10,388$ |  $4,661$ |  $5,535$ |  $16,742$ |  $16,856$
$75\%$ quantile |  $14,871$ |  $11,803$ |  $14,076$ |  $20,766$ |  $21,624$
min |  $260$ |  $23$ |  $174$ |  $491$ |  $458$
max |  $67,268$ |  $36,996$ |  $36,495$ |  $39,793$ |  $45,482$
std |  $6,999.84$ |  $6,847.36$ |  $6,799.17$ |  $7.118.88$ |  $8,219.33$


# Accelerating variational quantum Monte Carlo using the variational quantum
eigensolver

Ashley Montanaro  Phasecraft Ltd.  Stasja Stanisic  Phasecraft Ltd.

###### Abstract

Variational Monte Carlo (VMC) methods are used to sample classically from
distributions corresponding to quantum states which have an efficient
classical description. VMC methods are based on performing a number of steps
of a Markov chain starting with samples from a simple initial distribution.
Here we propose replacing this initial distribution with samples produced
using a quantum computer, for example using the variational quantum
eigensolver (VQE). We show that, based on the use of initial distributions
generated by numerical simulations and by experiments on quantum hardware,
convergence to the target distribution can be accelerated compared with
classical samples; the energy can be reduced compared with the energy of the
state produced by VQE; and VQE states produced by small quantum computers can
be used to accelerate large instances of VMC. Quantum-enhanced VMC makes
minimal requirements of the quantum computer and offers the prospect of
accelerating classical methods using noisy samples from near-term quantum
computers which are not yet able to accurately represent ground states of
complex quantum systems.

The variational quantum eigensolver (VQE) is a quantum algorithm which
attempts to produce the ground state of a quantum system by minimising the
energy of states of a particular form (“variational ansatz”) peruzzo14 . The
hope is that by optimising over a wider family of states than are accessible
to classical methods, one can produce more accurate approximations to ground
states. VQE may be more suitable for near-term quantum computers than some
other quantum algorithms, because it uses low-depth quantum circuits which are
optimised to make best use of the hardware. However, current quantum hardware
is as yet unable to run VQE on instances beyond the capacity of classical
exact diagonalisation.

In the period before VQE or other algorithms are able to produce ground states
accurately directly, it is of interest to develop ways in which VQE can
improve – and be improved by – existing classical methods. Here we consider
the case of variational Monte Carlo (VMC) methods sorella05 ; yokoyama87 ;
giamarchi91 ; carleo17 . These methods allow the calculation of properties of
quantum states that are too large to represent directly on a classical
computer, but which nevertheless have an efficient classical description, in
the sense that any given amplitude of the state can be computed efficiently
(up to an overall normalisation constant).

The intent behind VMC is to approximate the true ground state
$|\psi_{G}\rangle$ of a quantum system with such a classically efficient state
$|\widetilde{\psi_{G}}\rangle$, and then to compute properties of
$|\widetilde{\psi_{G}}\rangle$ in the hope that these give insight into
properties of $|\psi_{G}\rangle$. Although, from empirical and complexity-
theoretic arguments, it is not expected that all quantum ground states have an
efficient description of this form, in some cases VMC methods have
nevertheless been very effective yokoyama87 ; giamarchi91 ; carleo17 .

Here we propose an approach which we call quantum-enhanced variational Monte
Carlo (QEVMC), which is based on using samples produced by VQE – or some other
quantum algorithm – as input to VMC. We show using numerical simulations and
experimental data from a quantum computer that this hybrid approach can be
advantageous, in three respects:

  1. 1.

Convergence to the target distribution can be accelerated using VQE samples,
compared with a standard classically samplable initial distribution;

  2. 2.

VMC can reduce the energy compared with the energy of the original VQE state;

  3. 3.

VQE samples generated on a small quantum computer can be used within VMC to
accelerate convergence for systems beyond the quantum computer’s capacity.

We validate these three points for two well-studied proof-of-principle
systems: the Fermi-Hubbard model using the Gutzwiller wavefunction
gutzwiller63 , and the Ising model with transverse field using the neural
network quantum state ansatz carleo17 .

The first point above addresses a practical issue faced by VMC: convergence to
the target distribution can be slow choo20 ; zhao21 , which is a particular
challenge given the inherently sequential nature of the MCMC process. One
reason for slow convergence is that the initial distribution in VMC may be far
from the target distribution. VQE is expected to be able to produce initial
distributions which are closer to the target distribution, by better
representing the desired quantum correlations.

Regarding the second point above, note that in a situation where VMC improves
the VQE energy, we could instead have used a classical approach throughout to
achieve this improved energy, without the need for the quantum computer. Thus,
improving the VQE energy using VMC as a goal in itself seems to be primarily
of theoretical, rather than practical, interest. However, connecting to the
first point, using the quantum computer could speed up convergence to this
target energy, and hence make a previously infeasible VMC procedure feasible.

Method |  Type of Hamiltonian |  Needs only offline access |  Type of access needed |  Converges to true ground state
---|---|---|---|---
Quantum-enhanced MCMC layden22  |  Classical |  ✗ |  Sampling |  ✓
QC-QMC huggins22 ; wan23  |  Quantum |  ✓ |  Estimating overlaps |  ✓
CQA/QAEE xu22  |  Quantum |  ✗ |  Estimating overlaps |  ✓
QC-FCIQMC zhang22  |  Quantum |  ✗ |  Estimating elements of a Hamiltonian |  ✓
VQNHE zhang22a  |  Quantum |  ✗ |  Sampling (rel. error) |  ✗
NEM bennewitz22  |  Quantum |  ✓ |  Sampling |  ✗
Data-enhanced VMC czischek22 ; moss23  |  Experimental |  ✓ |  Sampling |  ✗
Quantum-enhanced VMC |  Quantum |  ✓ |  Sampling |  ✗
Table 1: QEVMC compared with related work for producing low-energy states of
quantum/classical systems using a hybrid quantum/classical approach. Methods
which converge to the true ground state in general require exponential time to
achieve this. Using QC-QMC without interactive access to the quantum computer
initially required exponential processing time huggins22 ; this was improved
to polynomial time in wan23 .

A crucial advantage of QEVMC is that it can be implemented using a very simple
protocol. A set of VQE samples, once generated, are stored and used to
generate a distribution which can be sampled from at the start of the VMC
procedure; there is no need to interact with the quantum computer after this
point. In addition, we expect QEVMC to be resilient to errors in the quantum
state. For example, in the simple case where the final state produced by VQE
is a mixed state $\rho_{V}=(1-\epsilon)\psi_{V}+\epsilon\mu$, where
$|\psi_{V}\rangle$ is the intended final state and $\mu$ is the maximally
mixed state, the use of VMC on $\rho_{V}$ is equivalent to using
$|\psi_{V}\rangle$ as input, except for an $\epsilon$ fraction of the samples,
where the uniform distribution is used. The final output of QEVMC should
therefore gracefully interpolate between the output generated with a perfect
input state, and the output generated with the classical uniform distribution
as input.

It is even possible that in some cases, a small amount of noise could be
beneficial to QEVMC, as it could reduce the distance between the initial
distribution generated by VQE and the target distribution.

###  .1 Related work

Quite a number of recent works have explored ways in which near-term quantum
algorithms could accelerate classical approaches for finding low-energy states
of quantum (or classical) systems based on Monte Carlo methods. For a very
recent survey, see mazzola23 . We summarise these and how they compare with
QEVMC in Table 1.

A standard approach for sampling from a desired probability distribution is
the Markov chain Monte Carlo method known as the Metropolis-Hastings
algorithm. This algorithm produces samples from a desired probability
distribution by repeatedly sampling from a proposal distribution that depends
on the current state, and accepting or rejecting the sample based on a certain
rule. Layden et al. layden22  propose replacing these samples with quantum
samples, where $P(x\to y)=|\langle y|U|x\rangle|^{2}$ for some symmetric
unitary $U$. Based on numerical and experimental results, and theoretical
intuition, they argue that this should lead to an algorithm which mixes more
rapidly to the target distribution than classical proposals. Layden et al.
apply their algorithm to producing samples from the Boltzmann distribution of
classical Ising models. Attractive aspects of their approach are that it
requires relatively simple quantum circuits $U$ and should be quite robust to
errors experienced by the quantum computer. However, their method requires
many rounds of communication with the quantum computer, as each new step of
the algorithm requires $U$ to be applied to a different computational basis
state. By contrast, the method proposed here relies only on the use of one set
of samples from the quantum computer, which can even be produced “offline” in
advance.

A central concept in quantum Monte Carlo methods in general is the trial state
$|\psi_{T}\rangle$, which is intended to approximate the ground state
$|\psi_{G}\rangle$ at the end of the algorithm. In VMC, the parameters of
$|\psi_{T}\rangle$ are chosen using an optimisation algorithm, whereas in
projector quantum Monte Carlo, the trial state is intended to approximate the
state corresponding to time-evolving the initial state for various imaginary
times $\tau$.

Huggins et al. huggins22  proposed to use a quantum computer to generate trial
states in projector QMC from an ansatz not accessible classically. Using these
trial states within the QMC framework requires approximately computing
amplitudes using a quantum computer, to small relative error. Achieving a
sufficiently small error may require a large, and even exponential, number of
measurements mazzola22 , although this question is delicate and may depend on
which variant of QMC is used lee22response . The approach of huggins22
initially required either the use of a quantum computer iteratively at each
stage of the algorithm, or exponentially costly classical processing, but the
classical processing cost with offline access to the quantum computer was
improved to polynomial by Wan et al. wan23 . Xu and Li xu22  presented a
related approach, also based on quantum trial states, with an explicit
application to accelerating VMC. Their approach also requires approximate
quantum computation of amplitudes, and needs a new quantum circuit to be
executed for each trial state.

Y. Zhang et al. zhang22  developed an alternative approach to accelerating the
variant of projector quantum Monte Carlo known as FCIQMC. A key issue
afflicting this method is the presence of the sign problem (“non-
stoquasticity” of the Hamiltonian; a stoquastic Hamiltonian is one where all
off-diagonal entries are real and non-positive). Y. Zhang et al. suggest the
use of a quantum circuit $U$ produced using VQE to approximately diagonalise
the Hamiltonian, such that the sign problem is mitigated. Their method
requires the use of $U$ at each step of the QMC algorithm to estimate entries
of the transformed Hamiltonian.

S.-X. Zhang et al. zhang22a  proposed the use of VQE combined with a neural
network, where the neural network is used to effectively implement a
nonunitary operation acting on the state produced via VQE, and expectations of
desired operators are estimated via a similar approach to that used in VMC.
Their method jointly optimises over the VQE parameters and the neural network
parameters, involves modifying the quantum circuit to be executed compared
with the original VQE circuit, and seems to require estimating probabilities
up to small relative error, in order to accurately compute a ratio of terms.

Bennewitz et al. bennewitz22  developed the concept of “neural error
mitigation” (NEM), whereby a neural quantum state (NQS) ansatz is used to
improve the energy produced by a quantum circuit via VQE. Similarly to
quantum-enhanced VMC, their approach requires only samples from measurements
on the quantum state that can be stored offline. Their method is based on
training the weights in the NQS such that the corresponding distribution
matches the distribution produced via VQE, as far as possible, as an initial
step before continuing the VMC optimisation process. This is unlike the
approach proposed here, where the distribution produced by VQE is used
directly as input to an MCMC algorithm. Indeed, in some cases, NQS sampling
can be performed directly, without the need for MCMC.

A similar approach of pre-training NQS weights was used in the context of data
from a Rydberg atom experiment (as opposed to data generated by a quantum
computer running VQE), by Czischek et al. czischek22  and Moss et al. moss23 .
These works show that training NQS weights based on experimental data can
reduce convergence time in VMC optimisation (as opposed to reducing MCMC
convergence time, which is our primary focus here).

Finally, Yang, Lu and Li yang21  have proposed a hybrid quantum Monte Carlo
algorithm for simulating time-dynamics with lower quantum circuit complexity
than standard quantum algorithmic techniques.

##  I Variational Monte Carlo

We begin by reviewing the VMC approach. Imagine that we would like to find the
ground state of a Hamiltonian $H$ on $n$ qubits. The VMC framework assumes
efficient access to the following:

  1. 1.

The ability to sample from an initial distribution $\mathcal{D}_{0}$ on
$\\{0,1\\}^{n}$;

  2. 2.

A mixer (Markov chain on $\\{0,1\\}^{n}$);

  3. 3.

The ability to compute amplitudes $\langle x|\psi_{\bm{\theta}}\rangle$,
perhaps up to an overall normalising constant. Here
$|\psi_{\bm{\theta}}\rangle$ is one of a family of states parametrised by a
vector $\bm{\theta}$.

We can use these tools to sample from the distribution
$\mathcal{D}_{\bm{\theta}}$, where
$\Pr_{\mathcal{D}_{\bm{\theta}}}[x]=|\langle
x|\psi_{\bm{\theta}}\rangle|^{2}=:p_{\bm{\theta}}(x)$, via the Metropolis-
Hastings algorithm. This algorithm proceeds as follows:

  1. 1.

Sample $x\sim\mathcal{D}_{0}$.

  2. 2.

Pick a new configuration $y\in\\{0,1\\}^{n}$ from a _proposal distribution_
$P$.

  3. 3.

Accept the new configuration with probability

$\min\left\\{1,\frac{p_{\bm{\theta}}(y)P(y\to x)}{p_{\bm{\theta}}(x)P(x\to
y)}\right\\}.$

  4. 4.

Go to step 2.

One can prove that (under weak conditions on $P$) eventually this algorithm
produces samples from the distribution $p_{\bm{\theta}}$ up to arbitrary
precision. This procedure has two downsides:

  * •

Samples following the first sample are correlated, with samples taken at close
times experiencing higher levels of correlation. To avoid this issue, one can
discard samples in between the first sample and subsequent ones retained.

  * •

It can take many steps of the random walk to converge to the desired
distribution $\mathcal{D}_{\bm{\theta}}$. This is particularly likely to be an
issue if $\mathcal{D}_{\bm{\theta}}$ is far away from $\mathcal{D}_{0}$.

Note that all we needed to execute the algorithm is the ability to produce
samples from $\mathcal{D}_{0}$, and evaluate probabilities
$p_{\bm{\theta}}(x)$ up to an overall scaling constant. In particular, we did
not need to be able to sample from $\mathcal{D}_{\bm{\theta}}$ (which is what
we were aiming to do), or to evaluate the normalising constant
$\sum_{x}p_{\bm{\theta}}(x)$.

Once we have samples from $\mathcal{D}_{\bm{\theta}}$, we can compute
observables $\langle\psi_{\bm{\theta}}|M|\psi_{\bm{\theta}}\rangle$, using

$\displaystyle\langle\psi_{\bm{\theta}}|M|\psi_{\bm{\theta}}\rangle$
$\displaystyle=\sum_{x\in\\{0,1\\}^{n}}\langle\psi_{\bm{\theta}}|M|x\rangle\langle
x|\psi_{\bm{\theta}}\rangle$
$\displaystyle=\sum_{x\in\\{0,1\\}^{n}}\frac{\langle\psi_{\bm{\theta}}|M|x\rangle}{\langle\psi_{\bm{\theta}}|x\rangle}|\langle
x|\psi_{\bm{\theta}}\rangle|^{2}.$

The quantity
$\frac{\langle\psi_{\bm{\theta}}|M|x\rangle}{\langle\psi_{\bm{\theta}}|x\rangle}$,
called the local energy in the case where $M=H$, can be computed efficiently
for many operators $M$ (e.g. if $M$ is a sum of a small number of Pauli
matrices). Thus, averaging the local energy over samples from
$\mathcal{D}_{\bm{\theta}}$ gives an approximation of
$\langle\psi_{\bm{\theta}}|M|\psi_{\bm{\theta}}\rangle$. Given that this
expression is a ratio of two terms involving $|\psi_{\bm{\theta}}\rangle$,
$|\psi_{\bm{\theta}}\rangle$ does not need to be normalised.

The overall VMC framework then optimises over variational parameters
$\bm{\theta}$ to find a state $|\psi_{\bm{\theta}}\rangle$ whose energy with
respect to $H$ is minimised.

##  II Quantum-enhanced variational Monte Carlo

In VMC, the initial distribution $\mathcal{D}_{0}$ is taken to be a “simple”
distribution which is related to the desired ground state. For example, a
uniform distribution (for a spin system), or a Slater determinant (for a
fermionic system). Here we propose an approach which we call quantum-enhanced
variational Monte Carlo (QEVMC), where these initial samples are replaced with
samples from a quantum state prepared using VQE. The rest of the VMC algorithm
remains unchanged. In particular, these initial samples can be produced by an
experiment using a quantum computer, and then stored offline until they are
used in VMC.

The intuition for why this could be advantageous is as follows. The quantum
state generated using VQE should have a higher fidelity with the true ground
state of $H$ than the initial state used in VMC. Therefore, it should take
fewer steps to reach the ground state distribution (if one could compute its
amplitudes) than starting with the VMC initial state. In VMC, we do not aim to
reach the true ground state, but a (variational) approximation. But if this
approximation is close to the true ground state, achieving high fidelity with
the ground state implies high fidelity with the approximation.

One can make this argument more rigorous as follows. Define the
$\chi^{2}$-divergence between distributions $\nu$ and $\pi$ as

$\chi^{2}(\nu,\pi)=\sum_{x}\pi(x)\left(\frac{\nu(x)}{\pi(x)}-1\right)^{2}=\sum_{x}\frac{\nu(x)^{2}}{\pi(x)}-1.$

Let $C$ be an ergodic Markov chain with stationary distribution $\pi$. Then,
given a sequence of distributions $\nu_{0},\dots,\nu_{n}$, where $\nu_{i+1}$
is produced by applying one step of $C$ to $\nu_{i}$, we have fill91

$4\|\nu_{n}-\pi\|_{1}^{2}\leq\lambda_{2}(M)^{n}\chi^{2}(\nu_{0},\pi),$

where $\lambda_{2}(C)<1$ denotes the second largest eigenvalue of $C$. This
implies that, in order to achieve $\|\nu_{n}-\pi\|_{1}\leq\epsilon$, it is
sufficient to take

$n=O\left(\frac{\log(\chi^{2}(\nu_{0},\pi)/\epsilon)}{\log(1/\lambda_{2}(C))}\right)=O\left(\tau\log(\chi^{2}(\nu_{0},\pi)/\epsilon)\right),$

where $\tau=1/(1-\lambda_{2}(C))$ is the relaxation time of $M$. This is not
straightforward to compute as it depends on the probability of acceptance in
the Metropolis-Hastings algorithm.

The $\chi^{2}$-divergence can be related to the fidelity between the
corresponding quantum states $|\nu\rangle$ and $|\pi\rangle$ obtained by
taking the square roots of the probabilities via Jensen’s inequality, which
shows that

$\frac{1}{\sqrt{\chi^{2}(\nu,\pi)+1}}\leq\sum_{x}\nu(x)\sqrt{\frac{|\pi(x)|}{|\nu(x)|}}=\langle\nu|\pi\rangle$

and hence $\chi^{2}(\nu,\pi)\geq 1/\langle\nu|\pi\rangle^{2}-1$. This is a
lower bound (as opposed to an upper bound), and only holds if $|\nu\rangle$
and $|\pi\rangle$ have the same phase for every amplitude, but it is still
suggestive: in order to obtain a small $\chi^{2}$-divergence, one needs to
have a high fidelity. On the other hand, if the fidelity between the starting
state and the desired final distribution is small (for example, exponentially
small) we will obtain a high $\chi^{2}$-divergence and expect a long mixing
time.

One situation where an initial distribution $\nu_{0}^{VQE}$ produced by VQE
could outperform an initial distribution $\nu_{0}^{C}$ produced classically is
where $\chi^{2}(\nu_{0}^{VQE},\pi)=2^{\alpha n}$, and
$\chi^{2}(\nu_{0}^{C},\pi)=2^{\beta n}$, for some $0<\alpha<\beta$. Then we
would expect a speedup by a factor of $\sim\beta/\alpha$. A more significant
speedup would be obtained if $\chi^{2}(\nu_{0}^{VQE},\pi)=O(1)$ –
corresponding to the VQE state having $\Omega(1)$ fidelity with the VMC state
– where the expected speedup would be a factor of $\Theta(n)$. These
theoretical bounds only provide intuition; the speedups achieved in practice
could be better (or worse).

##  III Results

We apply QEVMC to two well-studied models in condensed matter physics: the
Fermi-Hubbard model, using the Gutzwiller wave function ansatz yokoyama87 ,
and the Ising model with transverse field, using neural network quantum states
carleo17 .

###  III.1 Fermi-Hubbard model

We aim to find the ground state of instances of the Fermi-Hubbard model,

$H=-\sum_{\langle
i,j\rangle,\sigma}a^{\dagger}_{i\sigma}a_{j\sigma}+a_{j\sigma}^{\dagger}a_{i\sigma}+U\sum_{k}n_{k\uparrow}n_{k\downarrow},$

where $\sigma\in\\{\uparrow,\downarrow\\}$ and the sum is over sites $(i,j)$
in a rectangular lattice. This model has been extensively used as a benchmark
for both classical and quantum variational methods yokoyama87 ; Cade2020 ;
Stanisic2021 .

A variational ansatz which is frequently used in VMC is the Gutzwiller wave
function gutzwiller63 ,

$|\psi_{G}\rangle\propto
e^{-c\sum_{k}n_{k\uparrow}n_{k\downarrow}}|\psi_{SD}\rangle,$

where $|\psi_{SD}\rangle$ is a Slater determinant corresponding to the ground
state of the hopping part of $H$, and $c>0$. The intention behind this ansatz
is to penalise basis states with high double occupancy, which will receive a
large energy penalty from the onsite part of $H$. Amplitudes of
$|\psi_{G}\rangle$ can be evaluated efficiently classically, and to produce
this wave function in the VMC framework, it is natural to use as initial state
the Slater determinant $|\psi_{SD}\rangle$, whose corresponding distribution
can be sampled from directly. As a mixer, one can use the random walk which
moves a randomly chosen electron from one site to a neighbouring site
yokoyama87 . Computationally, this corresponds to updating the bitstring by
flipping a neighbouring pair of bits from 01 to 10 (or vice versa) within one
spin sector.

Here we initially applied QEVMC to the Gutzwiller wavefunction for system
sizes that are small enough to be simulated exactly classically, and which are
also accessible to quantum computing experiments. This enabled us to carry out
a full analysis and comparison of the performance of the algorithm.

We applied QEVMC to experimental data that had already been generated in
experiments applying VQE to the Fermi-Hubbard model on quantum hardware
Stanisic2021 . These experiments optimised over the family of quantum circuits
known as the Hamiltonian variational ansatz Wecker2015  to find good circuits
for preparing approximate ground states of the Fermi-Hubbard model. Although
in some cases the fidelity of the state produced with the true ground state
was relatively low, qualitative physical features of the ground state were
nevertheless reproduced Stanisic2021 . In this work we focus on three examples
from Stanisic2021 : the $1\times 4$ Fermi-Hubbard model with two variational
layers; the $1\times 8$ Fermi-Hubbard model with one variational layer; and
the $2\times 4$ Fermi-Hubbard model with one variational layer. These are all
instances where more variational layers would be required to achieve a high-
fidelity approximation to the ground state. In all cases we take $U=4$ and
work at half-filling.

First, we looked at the accuracy of the algorithm as it proceeds, for a
particular choice of Gutzwiller parameter $c$ which approximately minimises
the energy ($c\approx 0.421$ for $1\times 4$, $c\approx 0.431$ for $1\times
8$, $c\approx 0.453$ for $2\times 4$). We measured the accuracy via the total
variation distance from the desired probability distribution, and the error in
the energy.

The initial distributions that we compared were:

  * •

Slater determinant ($U=0$ ground state), the standard starting distribution
for this method;

  * •

Theoretical VQE distributions obtained via classical simulation;

  * •

Experimentally obtained distributions, taking advantage of error mitigation
via postselection on occupation number, and averaging over time-reversal and
reflection symmetries of the Fermi-Hubbard model Stanisic2021 .

The results are shown in Figures 1, 2 and 3. One can see that in all cases,
starting with the VQE distribution does lead to convergence to the desired
distribution. Also, in the cases of the $1\times 4$ experimental results, and
all of the simulated results, the distributions obtained by starting with the
VQE distribution have a lower total variation distance throughout the
algorithm. However, the $1\times 8$ experimental results do not achieve a
lower total variational distance throughout the whole algorithm. In the
$1\times 8$ and $2\times 4$ cases, the final energies obtained by VMC are
lower than the energies achieved by the VQE distribution alone (which are
$-3.37$ and $-3.71$, respectively, for the experimental distributions
Stanisic2021 ). Note that the local energy reported at the start of the VMC
algorithm is not necessarily an accurate reflection of the energy achieved by
the VQE distribution.

We additionally measured the speedup that would be obtained via VQE, by
measuring the number of steps of the Metropolis algorithm that are required to
achieve a given total variation distance level, when starting with a Slater
determinant, vs the VQE initial distribution. The results are shown in terms
of the ratio between the number of steps taken by QEVMC, vs. the number of
steps taken when starting with a Slater determinant. Here we found that the
simulated VQE distributions achieve a significant speedup (e.g. up to
$14\times$ for $1\times 8$). The $1\times 4$ experimental distribution
achieves a speedup up to $6\times$, but there is little speedup obtained in
the $1\times 8$ and $2\times 4$ experiments.

In the above experiments, in order to exactly determine the performance of the
various VMC starting distributions, we directly implemented the corresponding
stochastic matrix multiplication, which is inherently inefficient for large
problem sizes. We additionally implemented the true VMC method as described in
Section I. In addition to providing a check on the performance of the
algorithm in practice, this enabled us to use the VQE distribution within VMC
for problem sizes beyond the capacity of classical exact diagonalisation.

To demonstrate this, we considered $1\times L$ Fermi-Hubbard instances for
$L\in\\{16,24,48\\}$, corresponding to problems on up to 96 qubits. We use the
VQE distribution for $1\times 8$ Fermi-Hubbard to generate initial samples
from this distribution by concatenating $L/8$ independent samples, and aim to
sample from a distribution corresponding to a state with the same Gutzwiller
parameter as in $1\times 8$ ($c\approx 0.431$). The results are shown in
Figure 4. One can see that even for significantly larger problem sizes than
are accessible directly on the quantum computer, use of the VQE initial
distribution enables convergence to the target energy within a reasonable
number of steps, and in some cases apparently faster than using a Slater
determinant. Intriguingly, in the case of the $1\times 48$ system, use of the
experimental VQE distribution seems to give the fastest convergence. However,
given that the total variation distance to the desired distribution is not
accessible for large system sizes, it is unclear whether this effect would
persist across a wide range of other observables.

###  III.2 Ising model with transverse field

Here we aim to find the ground state of the 1D Ising model with transverse
field (TFI),

$H=-J\sum_{i=1}^{N-1}Z_{i}Z_{i+1}-h\sum_{j=1}^{N}X_{j},$

where $N$ is the number of lattice sites, $J=1$, and we assume open boundary
conditions (OBC) unless specified otherwise.

Figure 1: QEVMC applied to a $1\times 4$ Fermi-Hubbard instance at half-
filling. Speedup factor is the ratio between the chain length required to
achieve a given total variation distance from the target distribution when
starting with a Slater determinant, vs. starting with a VQE distribution. VQE
distributions are taken from Stanisic2021 .

Figure 2: QEVMC applied to a $1\times 8$ Fermi-Hubbard instance at half-
filling. VQE distribution taken from Stanisic2021 .

Figure 3: QEVMC applied to a $2\times 4$ Fermi-Hubbard instance at half-
filling. VQE distribution taken from Stanisic2021 .

(a) $1\times 16$

(b) $1\times 24$

(c) $1\times 48$

Figure 4: Using small VQE distributions as starting distributions for larger
VMC runs. Solid lines represent the mean energy obtained over 10,000 runs,
shaded regions are standard error of the mean. Dashed horizontal lines are
estimates of the target energy obtained from a VMC run with 1M steps starting
from a Slater determinant.

A successful ansatz for this model is the neural network quantum state (NQS)
ansatz proposed by Carleo and Troyer carleo17  and given as

$|\psi(\mathcal{W})\rangle\propto\sum_{S}\Psi(S;\mathcal{W})|S\rangle$

where $|S\rangle=|s_{1}\rangle\dots|s_{N}\rangle$ is a computational basis
state, $S$ is the bitstring which denotes that basis state, and the
coefficients are defined as

$\Psi(S;\mathcal{W})=\sum_{\\{h_{i}\\}}\mathrm{exp}\left[\sum_{j=1}^{N}a_{j}\sigma_{j}^{z}+\sum_{i=1}^{M}h_{i}(b_{i}+\sum_{j=1}^{N}W_{ij}\sigma_{j}^{z})\right],$

where $\sigma_{j}^{z}=\langle S|Z_{j}|S\rangle$, and $h_{i}=\\{-1,1\\}$
corresponds to $M$ hidden variables. The set of weights
$\mathcal{W}=\\{a_{i},b_{j},W_{ij}\\}$ are the parameters we optimize to find
the closest NQS to the ground state of some given Hamiltonian. As in the
original paper carleo17 , we are adding only a single hidden layer to the NQS,
therefore we use the traced out expression giving the simplified equation for
any bitstring $S$,

$\Psi(\mathcal{S};\mathcal{W})=\mathrm{exp}\left[\sum_{j=1}^{N}a_{j}\sigma_{j}^{z}\right]\times\prod_{i=1}^{M}F_{i}(S),$

where $F_{i}(S)=2\cosh\left[b_{i}+\sum_{j=1}^{N}W_{ij}\sigma_{j}^{z}\right]$.
While generally the weights for the NQS ansatz are complex-valued, motivated
by the solutions for the TFI model on 40 and 80 qubits in carleo17  which are
real-valued, we take them to be real-valued, halving the parameter search
space on this problem. We examine three regimes. Taking $J=1$, we have
$h\in\\{0.5,1.0,2.0\\}$, with $J/h=1$ being the critical phase-transition
point, and we use the number of nodes in the hidden layer to be $M=\alpha N$
with $\alpha\in\\{1,2,4\\}$. The data presented in the paper here is for $h=1$
unless otherwise specified, however the results for $h=0.5$ and $h=2$ were not
qualitatively different.

To find the NQS that best describes the ground state of TFI model, the
algorithm requires us to first run an optimization over the NQS weights
($\mathcal{W}$), minimizing the expected energy of the NQS with respect to the
TFI Hamiltonian. The standard algorithm used for this optimization is the
stochastic reconfiguration (SR) algorithm Sorella1998 ; Sorella2007 . The
initial NQS for this optimization is chosen by selecting random weights. Then
MCMC is used to acquire the necessary sample distribution and calculate the
expected energy values of the NQS at each iteration of the optimization
algorithm. The initial distribution is taken to be a uniform distribution, and
single spin flip is used as a mixer. This corresponds to updating the
bitstring by flipping a single spin chosen at random from $0$ to $1$ (or vice
versa). Each bitstring is updated independently (this corresponds to running
the MC algorithm many times from different initial bitstrings randomly chosen
from the initial distribution to acquire the bitstring from the final
distribution). Given this final distribution of samples at each iteration of
SR, to calculate the expected energy we calculate the local energy of the
samples. This is given as

$E_{\mathrm{loc}}(S)=\frac{\langle
S|H|\Psi(\mathcal{W})\rangle}{{\Psi(S;\mathcal{W})}},$

where $\Psi(\mathcal{W})$ is the full wavefunction NQS describes, and $S$ is a
given bitstring sample from the final distribution. We then take the average
of these local energy estimates for the final energy estimate.

We investigate here a potential speed-up using an already known VQE
distribution as the initial distribution for the Monte Carlo procedure. The
impact of this could be two-fold. For an already known NQS which describes the
ground state of a Hamiltonian, we could speed up the mixing time from the
initial distribution to the NQS distribution. If the NQS is not known, the
overall number of iterations the SR procedure needs to take to a given error
could offer a speed-up, therefore needing a smaller number of samples overall.

To generate the VQE distribution we use the Hamiltonian Variational ansatz
using an all plus state as the initial state. To generate the Hamiltonian
Variational ansatz, we split up the TFI Hamiltonian into even $ZZ$ terms, odd
$ZZ$ terms, and the the single qubit transverse field terms. In a single layer
of the ansatz we first apply all the even term gates parameterized by a single
parameter, then odd terms parameterized by another parameter, and finally all
the field terms which have a third parameter. Therefore each layer of VQE uses
three parameters, and we find optimal VQE states using the BFGS algorithm for
a set number of layers. We examine solutions up to four layers as we would
like the states produced by VQE not to approximate the ground state too well.

We consider a few different system sizes, similarly to the Fermi-Hubbard
model. We examine various aspects of the speed-up for a $16$ spin system via
exact calculations and we then examine a $24$ spin system with the true VMC
method. Finally, we also look at $L\in\\{40,80\\}$ spin system using the NQS
that was found in the original work carleo17  as approximate ground states of
the TFI model to examine potential speed-up by concatenating $L/20$
independent samples from the VQE distribution for a $20$-spin system.

To examine potential speed-up the VQE distribution gives to the MCMC for a
known NQS, we first find good NQS states using the stochastic reconfiguration
algorithm with a fixed number of iterations. We use 200 iterations for
$h\in\\{0.5,1.0,2.0\\}$. At each iteration, to estimate the energy of the NQS
at the current iteration, we use $5000$ samples from the initial distribution
(here, uniform) and take the last state of the Markov chain of length
$15$/$20$ for system size $16$/$24$. We estimate the final energy by using the
average of local energy estimate from the last five iteration, getting the
relative energy errors wrt to the ground state energy to be less than
$0.066/0.0042/7.4\times 10^{-7}$ for $h\in\\{0.5,1.0,2.0\\}$ demonstrating
good agreement with the ground state energy.

As we now have a good representation of the ground state as an NQS, we can
inspect the difference in mixing time between VQE and the uniform
distribution, using a similar approach to the small Fermi-Hubbard instances
considered above. For the $16$ spin system, we can calculate the exact total
variation distance (TVD) of the proposed distribution to the desired NQS
distribution as a function of the number of steps of the Markov chain. We can
also calculate the length of the Markov chain needed to achieve a specific
relative energy error (wrt the NQS energy) as well as to achieve a specific
TVD (to the NQS distribution). We then compare performance of these different
distributions, by comparing the length of MC chain needed to achieve a
specific TVD taking the uniform distribution as the baseline. We can see some
of the results in Figure 5.

Similarly, for the $24$ spin system, we no longer easily have access to exact
values, but we can look at energy estimates reached at different points in the
Markov chain, and compare the errors of the distributions. For some of the
results, refer to Figure 6 (for this specific result we run the SR algorithm
slightly longer, for $300$ iterations and get a relative energy error with
ground state energy to be 0.0017). We can see that in the case of $16$ spin
system with $h=1$, there is a speed-up of up to $20$ times when using a four
layer VQE state compared to the uniform distribution. For the $24$ spin system
with $h=1$, we see a speed-up of up-to $46$ times when using a one layer VQE
state compared to the uniform distribution.

Figure 5: QEVMC applied to a $16$ spin TFI instance ($h=1$) using an already
trained NQS ($\alpha=1$). The values presented here are calculated using exact
distributions. Initial distributions used marked in different colours,
“Uniform” stands for uniform distribution, “VQE $i$” for VQE distribution
obtained from the Hamiltonian Variational Ansatz with $i$ layers and obtaining
samples from the final state in the computational basis. The second plot has
vertical lines showing the relative energy error of the VQE state, calculated
exactly (color coded as on the legend). We can notice that these energy errors
are higher than the final errors after running VMC using the VQE
distributions.

Figure 6: QEVMC applied to a $24$ spin TFI instance ($h=1$) using an already
trained NQS ($\alpha=1$). Energy estimates were generated using $5000$ samples
for each datapoint. Energy error given is relative error to the estimated NQS
energy. To calculate the NQS energy estimate, a chain of length $500$ was
acquired with $5000$ samples and using the uniform distribution as the
starting distribution; then the mean value of the last $10$ points of the
chain is taken. On the first plot we can see the NQS energy estimate as black
horizontal line, we can also see the exact VQE energy for the four layer state
which is much higher than the energy estimated using VMC.

Next, we can also examine the case of $40$\- and $80$-spin systems in a
similar manner to the $24$ spin system, using the NQS results from carleo17 .
Here we use samples generated by $20$-spin systems with open boundary
conditions (OBC), while the NQS from carleo17  is for systems with periodic
boundary conditions (PBC). While one could consider running VQE for $20$-spin
systems with PBC, we hypothesize that OBC are more likely to give good results
due to how the two distributions are stitched together (specifically, there is
no need for correlation between the first and the 20th spin). Here, the
initial samples were drawn independently from the $20$-spin VQE distribution
and combined to give an initial sample for the larger system. It is an open
question whether alternative ways to piece different distributions together
given the physics of the model could offer better speedups. We look at the
behaviour of the relative energy error over a chain of length $200$ with
$5000$ samples taken for the energy estimate. Results can be seen in Figure 7.
We can see from the Figure that even though the VQE distribution was acquired
on a smaller system, it allows for speedups on larger systems, and the speedup
scales up for a larger system size ($80$ vs $40$ spins).

Figure 7:  QEVMC applied to a $40$ and $80$ spin TFI instance ($h=1$) using an
already trained NQS ($\alpha=4$). Energy error is relative error to estimated
ground state energy using the data from the original paper carleo17 . First
figure shows the relative error as the length of MCMC increase in the $80$
spin instance, second figure shows the needed chain length to achieve a
certain relative error; showing that the VQE as initial distribution allows
certain error rates to be reached quicker; the last two figures show the
speedup factor in terms of the needed length to a certain error when compared
to the uniform distribution as initial distribution in the $40$ and $80$ spin
case respectively. We can notice that the speedup scales up with the increase
in the size of the system.

Finally, we will want to see if it is possible to speed-up running of the SR
algorithm by using VQE distribution instead of uniform distribution as the
starting distribution for the MCMC at each iteration of the SR algorithm. We
will make this comparison by comparing the number of SR iterations we need to
achieve the same relative error using different distributions. We first have
to make sure that the length of the MCMC used before when running the SR
algorithm is indeed suitable for these various distributions. We visually
inspect the length of chain needed for a randomly chosen NQS target
distribution to be approximated given as initial distribution various VQE
generated distributions compared to the uniform distribution. We find that
after $20$ to $30$ steps, there is a qualitative agreement between the
distributions for the settings inspected here ($24$ qubits and
$h\in\\{0.5,1.0,2.0\\}$; in the case of $16$ qubits, the number of steps is
even lower). So fixing as before the MCMC length to $20$ at each iteration of
the SR algorithm (taking $5000$ samples), we find that the VQE distribution as
the initial distribution to mix towards the NQS distribution gives improvement
over uniform distribution in the case of $24$ qubits at various $h$, see
Figure 8, with more VQE layers allowing for higher speedup.

Figure 8: QEVMC applied to a $24$ spin TFI instance during stochastic
reconfiguration algorithm for $h=1$. Colour coding is blue for uniform
distribution, orange/green/red/purple for distribution acquired using 1/2/3/4
layers of VQE. Middle image shows SR iterations vs energy error. Right picture
shows the speed-up factor of the VQE distributions in comparison to uniform.
Energy error is relative error to exact ground state energy.

###  III.3 Origin of the QEVMC speedup

It is a natural question whether one can pin down the origin of the speedup
achieved by QEVMC. Given that the mixer is unchanged compared with the
standard Metropolis-Hastings algorithm, there are two plausible possibilities
for how this speedup is achieved:

  1. 1.

The initial distribution produced by VQE is closer to the target distribution
(so fewer steps are required to reach it);

  2. 2.

The acceptance rate of MCMC steps is higher when starting with the VQE
distribution (so more of the steps make progress towards the target
distribution).

In Figure 9 we present evidence that the second option is not the case, by
showing that the average acceptance rate is lower for QEVMC for both the
Fermi-Hubbard model and TFIM, albeit not substantially. This suggests that the
first option is likely to be the correct explanation for the acceleration
achieved by QEVMC.

(a) $1\times 8$ Fermi-Hubbard

(b) $1\times 16$ TFIM

Figure 9: Acceptance probabilities for QEVMC compared with starting with
classical distributions. Labels are the same as in Figures 1 and 5.

##  IV Discussion

We have shown that QEVMC can outperform standard classical VMC in terms of
convergence speed to the target distribution. Whether this advantage
translates into a genuine speed improvement in practice depends on a number of
factors, including the clock speed of the quantum computer and the behaviour
of the MCMC algorithm. For example, the level of correlation considered to be
acceptable in samples generated following the first sample from the target
distribution affects the gap that must be left between subsequent samples that
are kept, and the level of speedup available to the quantum algorithm. Chains
which require a long “burn in” time to sample from the target distribution are
the most promising candidates to achieve a significant quantum speedup.

An interesting follow-up question in the setting of our NQS results could be
regarding the performance of an adaptive algorithm where the length of the
chain and the number of samples changes as we get closer to the ground state.
As a marker of this proximity to the ground state, one could use the variance
of the local energy estimate. Since the VQE distribution gives an increasing
boost in performance as the NQS becomes closer to the ground state
description, in such an adaptive algorithm, one would expect the VQE
distribution to give further speedups over the uniform distribution.

Several works have proposed the use of experimental quantum states or those
generated via VQE to train weights in a NQS ansatz bennewitz22 ; czischek22 ;
moss23 , and it has been shown that this procedure can reduce the cost of VMC
optimisation czischek22 ; moss23 . It would be interesting to compare the
performance improvements that can be obtained via this approach with the MCMC
convergence speedups achieved here.

An alternative approach to the use of QEVMC is to generate the VMC target
distribution directly on a quantum computer. For example, algorithms are known
for constructing the Gutzwiller wave function murta21 ; seki21  whose quantum
circuit depth is not substantially greater than that of preparing the Slater
determinant $|\psi_{SD}\rangle$. However, these algorithms are based on the
use of postselection, so may require many runs of the corresponding quantum
circuit to prepare the desired state. This state could then be used as initial
distribution for VMC, in a similar way to QEVMC. This could be advantageous in
the realistic scenario that noise in the quantum circuit leads to errors in
the state prepared. Then VMC would enable the generation of samples
arbitrarily close to the target distribution.

The simple structure of the QEVMC algorithm opens up a use case for near-term
quantum computers, where samples from a classically intractable distribution
are produced and stored for later use as input to a classical algorithm. We
expect that other applications of this concept may be found.

### Acknowledgements

We would like to thank M. Schuyler Moss for pointing out references
bennewitz22 ; czischek22 ; moss23 , Jan Lukas Bosse and Lana Mineh for their
work on the VQE implementation used in this project, and the rest of the
Phasecraft team for helpful comments on this work. This project has received
funding from the European Research Council (ERC) under the European Union’s
Horizon 2020 research and innovation programme (grant agreement No. 817581).

## References

  * (1) A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik, and J. L. O’Brien.  A variational eigenvalue solver on a photonic quantum processor.  Nature Communications, 5(1):4213, 2014.
  * (2) S. Sorella.  Wave function optimization in the variational monte carlo method.  Phys. Rev. B, 71:241103, Jun 2005.
  * (3) H. Yokoyama and H. Shiba.  Variational monte-carlo studies of hubbard model. i.  Journal of the Physical Society of Japan, 56(4):1490–1506, 1987\.
  * (4) T. Giamarchi and C. Lhuillier.  Phase diagrams of the two-dimensional hubbard and t-j models by a variational monte carlo method.  Phys. Rev. B, 43:12943–12951, Jun 1991.
  * (5) G. Carleo and M. Troyer.  Solving the quantum many-body problem with artificial neural networks.  Science, 355(6325):602–606, 2017.
  * (6) M. C. Gutzwiller.  Effect of correlation on the ferromagnetism of transition metals.  Phys. Rev. Lett., 10:159–162, Mar 1963.
  * (7) K. Choo, A. Mezzacapo, and G. Carleo.  Fermionic neural-network states for ab-initio electronic structure.  Nature Communications, 11(1):2368, 2020.
  * (8) T. Zhao, S. De, B. Chen, J. Stokes, and S. Veerapaneni.  Overcoming barriers to scalability in variational quantum monte carlo.  In Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis, SC ’21, New York, NY, USA, 2021. Association for Computing Machinery.
  * (9) D. Layden, G. Mazzola, R. V. Mishmash, M. Motta, P. Wocjan, J.-S. Kim, and S. Sheldon.  Quantum-enhanced markov chain monte carlo, 2022.
  * (10) W. J. Huggins, B. A. O’Gorman, N. C. Rubin, D. R. Reichman, R. Babbush, and J. Lee.  Unbiasing fermionic quantum monte carlo with a quantum computer.  Nature, 603(7901):416–420, 2022.
  * (11) K. Wan, W. J. Huggins, J. Lee, and R. Babbush.  Matchgate shadows for fermionic quantum simulation, 2023.
  * (12) X. Xu and Y. Li.  Quantum-assisted monte carlo algorithms for fermions, 2022.
  * (13) Y. Zhang, Y. Huang, J. Sun, D. Lv, and X. Yuan.  Quantum computing quantum monte carlo, 2022.
  * (14) S.-X. Zhang, Z.-Q. Wan, C.-K. Lee, C.-Y. Hsieh, S. Zhang, and H. Yao.  Variational quantum-neural hybrid eigensolver.  Phys. Rev. Lett., 128:120502, Mar 2022.
  * (15) E. R. Bennewitz, F. Hopfmueller, B. Kulchytskyy, J. Carrasquilla, and P. Ronagh.  Neural error mitigation of near-term quantum simulations.  Nature Machine Intelligence, 4(7):618–624, 2022.
  * (16) S. Czischek, M. S. Moss, M. Radzihovsky, E. Merali, and R. G. Melko.  Data-enhanced variational monte carlo simulations for rydberg atom arrays.  Phys. Rev. B, 105:205108, May 2022.
  * (17) M. S. Moss, S. Ebadi, T. T. Wang, G. Semeghini, A. Bohrdt, M. D. Lukin, and R. G. Melko.  Enhancing variational monte carlo using a programmable quantum simulator, 2023.
  * (18) G. Mazzola.  Quantum computing for chemistry and physics applications from a monte carlo perspective, 2023.
  * (19) G. Mazzola and G. Carleo.  Exponential challenges in unbiasing quantum monte carlo algorithms with quantum computers, 2022.
  * (20) J. Lee, D. R. Reichman, R. Babbush, N. C. Rubin, F. D. Malone, B. O’Gorman, and W. J. Huggins.  Response to "exponential challenges in unbiasing quantum monte carlo algorithms with quantum computers", 2022.
  * (21) Y. Yang, B.-N. Lu, and Y. Li.  Accelerated quantum monte carlo with mitigated error on noisy quantum computer.  PRX Quantum, 2:040361, Dec 2021.
  * (22) J. A. Fill.  Eigenvalue bounds on convergence to stationarity for nonreversible Markov chains, with an application to the exclusion process.  The Annals of Applied Probability, 1(1):62–87, 1991.
  * (23) C. Cade, L. Mineh, A. Montanaro, and S. Stanisic.  Strategies for solving the Fermi-Hubbard model on near-term quantum computers.  Physical Review B, 102(23):235122, 2020.
  * (24) S. Stanisic, J. L. Bosse, F. M. Gambetta, R. A. Santos, W. Mruczkiewicz, T. E. O’Brien, E. Ostby, and A. Montanaro.  Observing ground-state properties of the Fermi-Hubbard model using a scalable algorithm on a quantum computer, 2021.
  * (25) D. Wecker, M. B. Hastings, and M. Troyer.  Progress towards practical quantum variational algorithms.  Phys. Rev. A, 92(4):042303, 2015.
  * (26) S. Sorella.  Green function monte carlo with stochastic reconfiguration.  Physical Review Letters, 80(20):4558–4561, 1998.
  * (27) S. Sorella, M. Casula, and D. Rocca.  Weak binding between two aromatic rings: Feeling the van der Waals attraction by quantum Monte Carlo methods.  The Journal of Chemical Physics, 127(1):014105, 07 2007.
  * (28) B. Murta and J. Fernández-Rossier.  Gutzwiller wave function on a digital quantum computer.  Phys. Rev. B, 103:L241113, Jun 2021.
  * (29) K. Seki, Y. Otsuka, and S. Yunoki.  Gutzwiller wave function on a quantum computer using a discrete hubbard-stratonovich transformation.  Phys. Rev. B, 105:155119, Apr 2022.