Himanshu<EMAIL_ADDRESS>1 Soumitri
<EMAIL_ADDRESS>2 Carnegie Mellon University UNC Chapel
Hill AL4FGSBIR

# Active Learning for Fine-Grained Sketch-Based Image Retrieval

###### Abstract

The ability to retrieve a photo by mere free-hand sketching highlights the
immense potential of Fine-grained sketch-based image retrieval (FG-SBIR).
However, its rapid practical adoption, as well as scalability, is limited by
the expense of acquiring faithful sketches for easily available photo
counterparts. A solution to this problem is Active Learning, which could
minimise the need for labeled sketches while maximising performance. Despite
extensive studies in the field, there exists no work that utilises it for
reducing sketching effort in FG-SBIR tasks. To this end, we propose a novel
active learning sampling technique that drastically minimises the need for
drawing photo sketches. Our proposed approach tackles the trade-off between
uncertainty and diversity by utilising the relationship between the existing
photo-sketch pair to a photo that does not have its sketch and augmenting this
relation with its intermediate representations. Since our approach relies only
on the underlying data distribution, it is agnostic of the modelling approach
and hence is applicable to other cross-modal instance-level retrieval tasks as
well. With experimentation over two publicly available fine-grained SBIR
datasets ChairV2 and ShoeV2, we validate our approach and reveal its
superiority over adapted baselines.

$\ast$$\ast$footnotetext: Both authors contributed equally to the paper.

##  1 Introduction

The success of computer vision applications can be largely attributed to deep
learning architectures [Huang et al.(2017)Huang, Liu, Van Der Maaten, and
Weinberger, Szegedy et al.(2015)Szegedy, Liu, Jia, Sermanet, Reed, Anguelov,
Erhan, Vanhoucke, and Rabinovich], which, in turn, have yielded favourable
results due to their access to large-scale labelled databases [Deng et
al.(2009)Deng, Dong, Socher, Li, Li, and Fei-Fei, Lin et al.(2014)Lin, Maire,
Belongie, Hays, Perona, Ramanan, Dollár, and Zitnick] for training. Being in
the age of Big Data, enormous volumes of data is easily available; however,
proper annotation of the same is a painstakingly cumbersome as well as an
expensive process, often requiring specialized qualifications if the task at
hand demands for domain expertise, such as handling medical images [Wu et
al.(2021)Wu, Shen, Huang, Wang, and Yao]. To alleviate this bottleneck,
researchers have proposed various annotation-efficient methods [Beluch et
al.(2018)Beluch, Genewein, Nürnberger, and Köhler, Li et al.(2019)Li, Wu, Liu,
Zhang, and Guan, Morgado et al.(2021)Morgado, Vasconcelos, and Misra, Zhang et
al.(2020)Zhang, Li, Yang, Wang, Zha, and Huang] to standard computer vision
tasks like classification and segmentation. A commonly used technique is
active learning [Sener and Savarese(2018), Sinha et al.(2019)Sinha, Ebrahimi,
and Darrell, Budd et al.(2021)Budd, Robinson, and Kainz], which seeks to find
the most “useful" unlabelled data samples to be annotated for learning, so as
to reduce annotation cost as well as increase overall generalisability on the
supervised learning task to be performed. Apart from conventional visual
tasks, active learning has been applied to other domains such as video
captioning [Chan et al.(2020)Chan, Vijayanarasimhan, Ross, and Canny], hand
pose estimation [Caramalau et al.(2021)Caramalau, Bhattarai, and Kim] and
single-image super-resolution [Wang et al.(2016)Wang, Gao, Zhang, and Li].

Figure 1: Intuition behind our proposed AL framework. The violation index
quantifies the “disturbance” introduced into the learned embedding space by
the incoming photo sample, due to having a greater similarity with a
previously paired sketch present in the latent space. More details are
provided in section 4.

In this paper, we embrace a paradigm shift to tackle the aforementioned
challenges in a domain which is fundamentally very much different from
traditional vision tasks – fine-grained sketch-based image retrieval (FG-SBIR)
[Bhunia et al.(2020)Bhunia, Yang, Hospedales, Xiang, and Song, Bhunia et
al.(2021)Bhunia, Chowdhury, Sain, Yang, Xiang, and Song, Sain et
al.(2021)Sain, Bhunia, Yang, Xiang, and Song], a relatively newer direction of
research from traditional category-level SBIR [Wang et al.(2015)Wang, Kang,
and Li, Collomosse et al.(2017)Collomosse, Bui, Wilber, Fang, and Jin]. As the
name suggests, FG-SBIR aims at exploiting the finer sketch representations for
cross-modal instance-level retrieval, achieved by learning an embedding space
where sketch-photo pairs lie close to each other. The most common approach in
several works [Yu et al.(2016)Yu, Liu, Song, Xiang, Hospedales, and Loy, Pang
et al.(2019)Pang, Li, Yang, Zhang, Hospedales, Xiang, and Song, Bhunia et
al.(2020)Bhunia, Yang, Hospedales, Xiang, and Song] has been to train a
supervised triplet loss-based model [Hoffer and Ailon(2015)] that learns
feature similarities between an image and its corresponding sketch and hence
requires a large number of sketch-photo pairs. However, drawing full sketches
is both time-consuming and difficult, since it requires artistic expertise and
amateurish sketching can only lead to learning degradation. Thus, there is a
need to develop a robust, annotation-efficient pipeline for FG-SBIR. Very few
recent works have been proposed treading on this motivation, such as a
generalisable zero-shot [Pang et al.(2019)Pang, Li, Yang, Zhang, Hospedales,
Xiang, and Song] and semi-supervised learning [Bhunia et al.(2021)Bhunia,
Chowdhury, Sain, Yang, Xiang, and Song] FG-SBIR framework. While the latter
involves training of two networks which makes it computationally expensive,
the performance of [Pang et al.(2019)Pang, Li, Yang, Zhang, Hospedales, Xiang,
and Song] is far from fully-supervised alternatives.

Developing an Active Learning pipeline for FG-SBIR imposes some unique
challenges. Considering a traditional FG-SBIR model which lacks a probability
distribution of the samples in its output, there is an absence of a direct way
to measure the uncertainty of such a model. Hence, it becomes difficult to
select a photo for labelling without having an estimate of its uncertainty
from the model. Moreover, off-the-shelf active learning methods [Sinha et
al.(2019)Sinha, Ebrahimi, and Darrell, Kirsch et al.(2019)Kirsch, Van
Amersfoort, and Gal] that were primarily proposed for classification tasks are
not suitable for FG-SBIR, since for classification the learning mechanism
draws firm discriminatory boundaries among samples, whereas in a cross-modal
instance-level retrieval setup the objective is to draw softer decision
boundaries, as the samples belong to the same category and only differ in
minute fine-grained details. Additionally, the FG-SBIR model holds the photo
and sketch embeddings in a joint space and thus, selection based on only photo
or sketch might yield previously unseen results, compared to selection from a
single modality. Hence, developing a sampling technique for the FG-SBIR
requires the handling of both modalities - photo and its sketch.

The effectiveness of an active learning pipeline depends highly upon the
technique of selecting samples from the unlabelled pool. As a result, a good
technique for sampling photos would lead to the maximum increase in the
model’s performance. To this end, we propose a novel sampling strategy for
active learning that utilises the evolving relations between photos and their
sketches to approximate the influence of a new photo from the unlabelled pool,
on the model’s existing knowledge. Our sampling technique is informed by the
embedding space learnt using the labelled pool of photo-sketch pairs and two
major aspects of an unlabelled photo – its predicted representation and its
approximate potential of influence in the existing embedding space. While the
former is the basis of our sampling technique, the latter helps in tackling
the classic trade-off between uncertainty and diversity-based sampling
techniques. Specifically, to model an unlabelled photo’s potential of
influencing the existing embedding space, we formulate a quantity, violation
index. The violation index of a photo acts as a proxy for measuring the amount
of confusion it could create in the existing sketch-photo pair embedding space
(refer to Figure 1 for an intuitive understanding). Further, for diversity
sampling, we choose $k$-MEANS++ due to its ability to converge faster, as
highlighted in [Ash et al.(2019)Ash, Zhang, Krishnamurthy, Langford, and
Agarwal]. The overall advantage is two-fold: since we adopt a relative method
of approximating the influence using the task network itself, it does not need
an auxiliary learner or source of knowledge. Moreover, our technique relies on
the underlying distribution of the data itself and hence is agnostic of the
model used. The number of samples to be chosen depends on the permissible
budget of annotation; the chosen samples are then queried for their sketches
to be drawn, which are then paired and put into the training examples.

To sum up, the primary contributions of the presented study are as follows:
(1) For the first time, we propose an active learning pipeline for annotation-
efficient fine-grained SBIR; (2) To this end, we formulate a novel sampling
strategy that incorporates uncertainty as well as diversity to quantify the
“usefulness" of an unlabelled sample for querying its label (here, sketch);
(3) With suitable experimentation and comparison with adopted baselines on two
publicly available FG-SBIR datasets, as well as conducting ablations on the
proposed framework, we demonstrate the usefulness of our approach.

##  2 Related Works

Fine-grained SBIR: Although SBIR was originally proposed and studied as a
category-level retrieval task [Wang et al.(2015)Wang, Kang, and Li, Bui et
al.(2018)Bui, Ribeiro, Ponti, and Collomosse, Sain et al.(2022)Sain, Bhunia,
Potlapalli, Chowdhury, Xiang, and Song], recently there has been significant
interest among researchers towards exploiting the _fine-grained_ information
that sketches provide [Bhunia et al.(2020)Bhunia, Yang, Hospedales, Xiang, and
Song, Song et al.(2017)Song, Yu, Song, Xiang, and Hospedales, Pang et
al.(2019)Pang, Li, Yang, Zhang, Hospedales, Xiang, and Song, Bhunia et
al.(2021)Bhunia, Chowdhury, Sain, Yang, Xiang, and Song, Sain et
al.(2023b)Sain, Bhunia, Koley, Chowdhury, Chattopadhyay, Xiang, and Song] for
enhanced cross-domain matching, something other query mediums (e.g. text) fail
to do. The first deep learning-based approach by Yu et al. [Yu et al.(2016)Yu,
Liu, Song, Xiang, Hospedales, and Loy] was further enhanced using cross-domain
generative-discriminative learning [Pang et al.(2017)Pang, Song, Xiang, and
Hospedales] and attention mechanisms [Song et al.(2017)Song, Yu, Song, Xiang,
and Hospedales]. More recent studies include zero-shot-like cross-category FG-
SBIR [Pang et al.(2019)Pang, Li, Yang, Zhang, Hospedales, Xiang, and Song,
Sain et al.(2023a)Sain, Bhunia, Chowdhury, Koley, Xiang, and Song]; cross-
modal co-attention-based hierarchical model [Sain et al.(2020)Sain, Bhunia,
Yang, Xiang, and Song]; on-the-fly SBIR setup for early retrieval [Bhunia et
al.(2020)Bhunia, Yang, Hospedales, Xiang, and Song]; style-agnostic meta-
learning setup [Sain et al.(2021)Sain, Bhunia, Yang, Xiang, and Song] and a
semi-supervised framework [Bhunia et al.(2021)Bhunia, Chowdhury, Sain, Yang,
Xiang, and Song] to tackle data scarcity in FG-SBIR [Sain et al.(2023b)Sain,
Bhunia, Koley, Chowdhury, Chattopadhyay, Xiang, and Song]. However, none of
these works seeks to address the need for a smart data labelling system for
cross-modal instance-level retrieval problems like FG-SBIR so that we can
achieve optimal performance using the minimum possible labelling budget.

Active Learning: Active learning (AL) [Ren et al.(2020)Ren, Xiao, Chang,
Huang, Li, Chen, and Wang, Budd et al.(2021)Budd, Robinson, and Kainz] has
been extensively studied for over two decades, the primary goal being to
develop an effective strategy to reduce annotation effort by “actively"
selecting representative samples and improve learning. Apart from conventional
image classification [Beluch et al.(2018)Beluch, Genewein, Nürnberger, and
Köhler, Sinha et al.(2019)Sinha, Ebrahimi, and Darrell], AL has been widely
used for various computer vision tasks such as medical imaging [Budd et
al.(2021)Budd, Robinson, and Kainz], image and video segmentation [Siddiqui et
al.(2020)Siddiqui, Valentin, and Nießner, Xie et al.(2020)Xie, Feng, Chen,
Sun, Ma, and Song, Belharbi et al.(2021)Belharbi, Ben Ayed, McCaffrey, and
Granger], among others. Broadly speaking, AL approaches may be categorized as:
(1) Uncertainty-based methods [Wang and Shang(2014), Beluch et
al.(2018)Beluch, Genewein, Nürnberger, and Köhler, Yoo and Kweon(2019), Choi
et al.(2021)Choi, Yi, Kim, Choo, Kim, Chang, Gwon, and Chang], which aim to
construct a so-called acquisition function that quantifies the uncertainty of
the model on unlabelled data points, based on which the points are sampled and
queried for annotation; and (2) Representation-based methods [Sener and
Savarese(2018), Sinha et al.(2019)Sinha, Ebrahimi, and Darrell, Ash et
al.(2019)Ash, Zhang, Krishnamurthy, Langford, and Agarwal, Kim et
al.(2021)Kim, Park, Kim, and Chun], which seek to learn a common embedding
space for the labelled and unlabelled data items so as to sample the
unlabelled data points that capture the most diverse regions of the embedding
and thereby better represent the overall data distribution. Existing AL
methods mostly deal with classification and thus cannot be directly adopted
for FG-SBIR where paired photo and sketch need to be aligned in the vicinity
in the joint embedding space. This cross-modal instance-wise matching brings
an exclusively different set of challenges for employing AL in FG-SBIR. We
intend to address and tackle these through this work, which, to the best of
our knowledge, is the first to introduce AL to a cross-modal instance-level
retrieval problem.

##  3 Background and Problem Formulation

Baseline FG-SBIR: Instead of complicated pre-training [Pang et al.(2020)Pang,
Yang, Hospedales, Xiang, and Song] or joint-training [Bhunia et
al.(2021)Bhunia, Chowdhury, Sain, Yang, Xiang, and Song], we use a three-
branch state-of-the-art siamese network [Pang et al.(2019)Pang, Li, Yang,
Zhang, Hospedales, Xiang, and Song] as our baseline retrieval model, which is
considered to be a strong baseline to date [Pang et al.(2019)Pang, Li, Yang,
Zhang, Hospedales, Xiang, and Song, Song et al.(2017)Song, Yu, Song, Xiang,
and Hospedales]. Each branch starts from ImageNet pre-trained VGG-16 [Karen
Simonyan(2015)], sharing equal weights. Given an input image
$I\in\mathbb{R}^{H\times W\times 3}$, we extract the convolutional feature-map
$\mathcal{F}({I})$, which upon global average pooling followed by $l_{2}$
normalisation generates a $d$ dimensional feature embedding. This model has
been trained with an anchor sketch (a), a positive (p) photo, and a negative
(n) photo triplets $\\{a,p,n\\}$ using triplet-loss. Triplet-loss aims at
increasing the distance between anchor sketch and negative photo
$\delta^{-}={||\mathcal{F}(a)-\mathcal{F}(n)||}_{2}$, while simultaneously
decreasing the same between anchor sketch and positive photo
$\delta^{+}={||\mathcal{F}(a)-\mathcal{F}(p)||}_{2}$. Therefore, the triplet-
loss $\mathcal{L}$ with the margin hyperparameter $\mu>0$ can be written as:

$\mathcal{L}=max\\{0,\delta^{+}-\delta^{-}+\mu\\}$ (1)

During inference, given a gallery of $M$ photos $\\{P_{i}\\}_{i=1}^{M}$, we
can compute a list of $d$ dimensional vectors as
$G=\\{\mathcal{F}(P_{i})\\}_{i=1}^{M}$. Now, given a query sketch $S$, and
pair-wise distance metric, we obtain a top-q retrieved list as
$Ret_{q}(F(S),G)$. If the paired (ground-truth) photo appears in the top-q
list, we consider accuracy to be true for that sketch sample.

Active Learning for FG-SBIR: Following the principle of active learning shown
in Equation 2, we aim to minimise the number of rounds $\mathcal{R}$ so that
the fewest possible sketches are needed to be acquired. At the end of each of
round $r$, the performance measures $\mathcal{P}$ are recorded to quantify and
understand the effectiveness of our sampling technique $\mathcal{X}$.

$\min_{\mathcal{R}}\min_{\mathcal{L}}{[\mathcal{X}(\mathcal{{L}}|\mathcal{P}_{0}\subset\cdots\mathcal{P}_{k}\subset\mathcal{P}_{U})]}_{r=1}^{\mathcal{R}}$
(2)

##  4 Proposed Method

Figure 2: Overall workflow of the proposed active learning framework for fine-
grained SBIR. The method starts with a subset of sketch-photo pairs for
learning the FG-SBIR model, following which it is used to compute violation
index of unlabeled photos with respect to the embedding space. Our sampling
technique selects a suitable query set of photos which are passed to the
Oracle (the ground truth sketch provider for queried photos) to obtain their
sketch counterparts. These pairs are subsequently added to the labeled subset.

Overview: We aim to design an active learning framework specially suited for
fine-grained SBIR so as to make it label-efficient. To keep our framework
fairly simple, we have used the baseline triplet network (described in Section
3) as the retrieval model. During the training process, we leverage active
learning to select the most informative examples for querying its labels, the
details of which are given in the subsequent sections. The AL mechanism is
employed in training cycles; at the end of each cycle the sampling technique
selects photos from the unlabelled data and queries for their paired sketches,
which are provided by the Oracle (ground truth annotator) and added to the
labelled training pool.

###  4.1 Introducing Violation Index

In this section, we introduce the core idea of our contribution – the
violation index. Our approach is inspired by the learning process of cross-
modal retrieval frameworks. When an FG-SBIR model learns, it tries to put an
image and its sketch nearby in the embedding space, while pushing away other
images. Formally, consider a model $\mathcal{M}$, the embedding of an image
obtained from final layer of the model as ${E_{I}}$, and that of the
corresponding sketch as ${E_{S}}$. Also, let ${E_{I^{\prime}}}$ be the
embedding of an image other than ${E_{I}}$, $N_{L}$ be the total size of the
labelled dataset $\mathcal{D}_{L}$ and $N_{U}$ be the total size of the
unlabeled dataset $\mathcal{D}_{U}$. To learn the similarities, following
Section 3 let us consider a triplet loss function defined as follows:

${\mathcal{L}}=\operatorname{max}\left({\|E_{I}-E_{S}\|}^{2}-{\|E_{I}-E_{I^{\prime}}\|}^{2}+\mu,0\right)$
(3)

The objective of the learning process is to minimise $\mathcal{L}$ over all
images in the dataset (Eq 4). Hence, an ideal embedding space would comprise
images lying closest to its actual sketch. However, a non-ideal embedding
space is more realistic and introduces some new challenges.

${minimize}\sum_{i=1}^{{}N_{L}}{\mathcal{L}}\left(E_{I}^{(i)},E_{S}^{(i)},E_{I^{\prime}}^{(i)}\right)$
(4)

Now, when a new image $I^{\prime}$ from the unlabelled set $\mathcal{D}_{U}$
is introduced during training, its corresponding sketch might not lie closest
to it in the embedding space, a condition defined in Equation 5. Hence, in the
model’s view, the image now seems closer to one of the existing sketches
rather than its own. This phenomenon leads to perturbation in the existing
embedding space. We quantify the degree of disturbance an image from the
unlabelled pool of images produces in the existing embedding space and call it
its _violation index_.

${\|E_{I^{\prime}}-E_{S}\|}^{2}\leq{\|E_{I}-E_{S}\|}^{2}\ $ (5)

Formally, we define the violation index ($VI$) of the photo embedding
$E_{I^{\prime}}$ as:

$Violation\,Index(E_{I^{\prime}})=\frac{1}{N_{L}}\times\sum_{i=1}^{{}N_{L}}\frac{\|E_{I}^{(i)}-E_{S}^{(i)}\|}{\|E_{I^{\prime}}-E_{S}^{(i)}\|}$
(6)

Thus, the violation index is an improvement over a simpler distance-based
sampling technique since it accounts for the fact that the inherent
imperfections in sketches cause the cross-modal embedding space to be
sensitive to new image-sketch pairs. Moreover, since the metric is built on
relative distances instead, computing the relative similarities between a new
image and existing pairs help surface novel samples that still do not violate
existing sketches or images.

###  4.2 VI-based Active Sampling

Given the violation index of an image shows its average relative distance from
existing image-sketch pairs and hints towards how many of them it violates in
the learned embedding space. Intuitively, images with a low violation index
would be relatively unseen in the training data, whereas ones with a high
violation index might closely resemble one or more images from the training
set.

Images with a low violation index or min violating samples are relatively
unseen in the training data, which means that they may be novel or unique
compared to other images. These images may contain features or characteristics
that are not commonly found in the training set. On the other hand, images
with a high violation index or max violating samples are more likely to
closely resemble one or more images from the training set. This suggests that
these images may contain familiar features that are present in the training
data.

When considering the selection strategy in Active Learning, a naive solution
could be to select the samples with minimum violation indices,i.e., the
sampling technique $X$ gives us the set
$\\{x\in{D_{U}}:\mid{\emph{VISet}_{D_{U}}}\cap(-\infty,\emph{VI}(E_{x})\mid<K\\}$
where $\emph{VISet}_{D_{U}}$ is the set of violation indices for the images in
the unlabelled set and $\emph{VI}(E_{x})$ is the violation index for image $x$
(Eq 6). Although with this approach the perturbation in the existing embedding
space is minimized, it loses out on the opportunity to reduce uncertainty by
learning through closely resembling images. As such, there exists a tradeoff
between reducing existing uncertainty and learning novel instances. Hence, a
better selection strategy would be an ensemble of minimum and maximum
violating image samples, i.e, a better $X$ gives us the set

$\displaystyle\\{x\in{D_{U}}:\mid{\emph{VISet}_{D_{U}}}\cap(-\infty,\emph{VI}(E_{x}))\mid<p\\}\,\,\,\,\cup$
$\displaystyle\\{x\in{D_{U}}:\mid{\emph{VISet}_{D_{U}}}\cap(\emph{VI}(E_{x}),\infty)\mid>(K-p)\\}$

where $p\in I$, is a hyper-parameter such that $0\leq p\leq N_{U}$.

Since the violation index does not inherently capture the diversity of images,
a further enhancement of our selection strategy includes diversity sampling to
maximize novelty in the selected subset. To do this, we adopt a $kmeans++$
based clustering of the image embeddings followed by violation index-based
selection inside each cluster. We select $kmeans++$ due to its ability to
converge faster [Ash et al.(2019)Ash, Zhang, Krishnamurthy, Langford, and
Agarwal]. The overall workflow has been depicted in Figure 2.

##  5 Experiments and Results

Datasets: We use QMUL-Shoe-V2 [Pang et al.(2019)Pang, Li, Yang, Zhang,
Hospedales, Xiang, and Song, Riaz Muhammad et al.(2018)Riaz Muhammad, Yang,
Song, Xiang, and Hospedales] and QMUL-Chair-V2 [Song et al.(2018)Song, Pang,
Song, Xiang, and Hospedales] datasets that have been specifically designed for
FG-SBIR. QMUL-Shoe-V2 contains a total of 6,730 sketches and 2,000 photos, of
which we use 6,051 and 1,800 respectively for training and the rest for
testing. For QMUL-Chair-V2, we split it as 1,275/725 sketches and 300/100
photos for training/testing respectively. For each photo, one of its possible
sketches is selected randomly, and is considered as its label. Initially, we
consider 300 photo-sketch pairs as the labelled set and the rest as unlabelled
(i.e. absence of their corresponding sketches).

Implementation: We implemented our framework in PyTorch [Paszke et
al.(2017)Paszke, Gross, Chintala, Chanan, Yang, DeVito, Lin, Desmaison,
Antiga, and Lerer] accelerated by an 11 GB Nvidia RTX 2080-Ti GPU. ImageNet
[Russakovsky et al.(2015)Russakovsky, Deng, Su, Krause, Satheesh, Ma, Huang,
Karpathy, Khosla, Bernstein, et al.] pre-trained VGG-16 [Karen Simonyan(2015)]
network (embedding dimension $D=256$) is used as the backbone network for both
sketch and photo branches. In all experiments, we use Adam optimizer [Kingma
and Ba(2014)] with learning rate of $1e-4$, batch size $16$ and train the base
model with a triplet objective. In the active learning setup, we conduct 5
cycles of complete training, where at the end of each round, we employ our
sampling technique to add $K$ samples to the labelled pool after the provision
of its actual sketch.

In each active learning round, we obtain the embeddings for the photos and
sketches from the labelled set and the photos from the unlabelled set.
Following this, our sampling technique utilises these embeddings to select $K$
photos from the unlabelled set, whose corresponding sketches from the dataset
are then used to label them. Once labelled, these photos and their sketches
are added to the labelled set.

Evaluation Metrics: Following the standard FG-SBIR setting [Bhunia et
al.(2020)Bhunia, Yang, Hospedales, Xiang, and Song, Bhunia et al.(2021)Bhunia,
Chowdhury, Sain, Yang, Xiang, and Song], we quantify the performance of the
sampling technique using acc.@$q$ metric ($q=1,10$), i.e. the percentage of
sketches having true-match photos appearing in the top-$q$ list, after each AL
round.

###  5.1 Baselines

To the best of our knowledge, there has been no previous work on active
learning for fine-grained sketch-based image retrievals. Thus, we compare ours
with a few SoTA active learning techniques adapted suitably for fine-grained
SBIR. We choose three widely used baseline sampling techniques, namely: random
sampling, kmeans sampling, and coreset sampling. As evident from our choice of
baseline methods, we do not consider uncertainty-based sampling techniques due
to the limitations discussed in Section 3.3. Random sampling is a widely used
baseline for evaluating active sampling techniques. Here, we randomly select
photos for labelling based on our-predefined budget. In K-means sampling
[Zhdanov(2019)], we cluster photos into as many clusters as the labeling
budget and select the photos closest to the cluster centroid. In case of a
tie, we randomly select any one of the closest photos. Finally, Coresets
sampling [Sener and Savarese(2018)] is another diversity sampling technique in
which we utilise "coresets" to find the most representative photos for
sampling. To find the coresets, we use a farthest-first approach to select
maximally distant photos in the embedding space.

###  5.2 Performance Analysis

$\displaystyle{20}$$\displaystyle{25}$$\displaystyle{30}$$\displaystyle{35}$$\displaystyle{40}$$\displaystyle{45}$$\displaystyle{50}$%
of Labelled
Data$\displaystyle{20.0}$$\displaystyle{22.5}$$\displaystyle{25.0}$$\displaystyle{27.5}$$\displaystyle{30.0}$Top-1
Accuracy (%)QMUL-Shoe-V2
Dataset$\displaystyle{20}$$\displaystyle{25}$$\displaystyle{30}$$\displaystyle{35}$$\displaystyle{40}$$\displaystyle{45}$$\displaystyle{50}$%
of Labelled Data$\displaystyle{30}$$\displaystyle{40}$$\displaystyle{50}$QMUL-
Chair-V2 DatasetSampling TechniqueRandomK-meansViolation Index (ours)Core-
setViolation Index (ours)Core-setK-meansRandom

Figure 3: Comparing VI-based active sampling (ours) with SoTA AL baselines for
fine-grained SBIR. All results are reported as mean of 5 runs and an error
band of $\sigma=1.0$.

To compare our approach with baselines, we report the mean and one standard
deviation of acc.@$1$ values across 5 different runs. We also report the
initial accuracy of the model. Figure 3 compares our violation index-based
approach and baseline techniques. Before any active learning, the acc.@$1$
after training it on 8% labelled data were 13% (49.8%) and 11.7% (42.8%) on
the Shoe and Chair datasets respectively.

Our method outperforms the baselines on both QMUL-Shoe-V2 and QMUL-Chair-V2
datasets and obtains consistently higher mean accuracies. We consider the mean
accuracy differences between our approach and the baselines to further
substantiate the results shown in Figure 3. On the QMUL-Chair-V2 dataset, we
achieve a mean gain of 2.6% acc.@$1$ compared to the 3 baseline techniques. On
the QMUL-Shoe-V2 dataset, we observe a mean increase of 1.1% acc.@$1$.

Our proposed approach achieves comparable performance by utilizing only 40-50%
of the dataset. Compared to the state-of-the-art acc.@$1$ of 36.47% on 100% of
the QMUL-Shoe-V2 dataset [Pang et al.(2019)Pang, Li, Yang, Zhang, Hospedales,
Xiang, and Song], we obtain a mean acc.@$1$ of 29.8% by only utilizing 40% of
the dataset. This is achieved when we use $\alpha$=0 as our hyperparameter. We
also see consistent results on the QMUL-Chair-V2 dataset (using $\alpha$=0.7),
where we obtain 49.6% acc.@$1$ by only using 50% of the dataset.

###  5.3 Ablation Study

To better understand the effects of violation index-based sampling, we perform
ablations on violation index-based selection and the choice of hyper-
parameter.

Significance of violation index and diversity-based sampling: For this, we
consider selecting only the unlabelled images that have the smallest and
largest violation indices in active learning rounds performed on the QMUL-
Shoe-V2 dataset.

As per Figure 4, an interesting observation is made: In early active learning
cycles, selecting minimum VI yields higher acc.@$1$ compared to maximum VI. As
the training dataset grows, this relationship is inverted, with higher
acc.@$1$ achieved through Maximum VI-based selection. Empirical results
support our hypothesis on violation indices, indicating that with less
training data, increased retrieval accuracy primarily stems from selecting
diverse instances. Conversely, with larger training data, reducing model
uncertainty on similar instances contributes more to acc.@$1$. Thus, selecting
maximum VI samples yields optimal results.

Figure 4: Ablations on diversity-based sampling and violation index. We
compare only VI-based sampling, diversity-based sampling, and a combination of
both. Figure 5: Ablations on value of hyper-parameter $\alpha$ on the ShoeV2
and ChairV2 datasets.

Another significant aspect of our experiment involves the need for diversity-
based sampling, not just violation index-based sampling. As shown in Figure 4,
violation index-based selection within diverse clusters obtained from kmeans++
consistently outperforms vanilla violation-based approach and kmeans++. Thus,
the violation index captures semantic relations between labeled and unlabeled
datasets but does not explicitly utilize relations between unlabeled images
for diverse selection.

Sensitivity to hyper-parameter:  We analyze the sensitivity of our method to
the hyperparameter $\alpha$, the results shown in Figure 5. We vary the value
of $\alpha$ from 0 to 1 with a gap of 0.1 and report the mean acc.@$1$ on both
datasets. On the QMUL-Shoe-V2 dataset, we observe a steady increase in
acc.@$1$ as we increase $\alpha$. On this dataset, $\alpha$=0 produces the
best results. On the QMUL-Chair-V2 dataset, we observe a steady increase in
initial and final acc.@$1$ as we increase $\alpha$. As per Figure 6 and on
this dataset, $\alpha$=0.7 produces the best results. On both datasets, there
is a sudden drop in performance as we increase $\alpha$ from 0.0 to 0.1 and
0.2. We also present the same results compared with the baseline method in
Figure 6. We believe that studying this interesting behavior is an open
research avenue.

Figure 6: Ablations on value of hyper-parameter $\alpha$ and its comparison
with baselines on the QMUL-Shoe-V2 and QMUL-Chair-V2 datasets.

##  6 Conclusion

We have proposed an active learning framework to tackle the annotation
bottleneck in fine-grained SBIR. To this end, we brought forth a quantifiable
metric, violation index, that measures the latent space displacements due to
addition of a new instance. With suitable experiments and ablations, we have
shown the robustness of our model compared to classification-specific AL
works, especially in the low-data regimes. Our model is modality agnostic and
thus can be leveraged for any cross-modal retrieval task. In future, we plan
to extend our studies by investigating its applicability to such other tasks.

## References

  * [Ash et al.(2019)Ash, Zhang, Krishnamurthy, Langford, and Agarwal] Jordan T Ash, Chicheng Zhang, Akshay Krishnamurthy, John Langford, and Alekh Agarwal.  Deep batch active learning by diverse, uncertain gradient lower bounds.  In _ICLR_ , 2019.
  * [Belharbi et al.(2021)Belharbi, Ben Ayed, McCaffrey, and Granger] Soufiane Belharbi, Ismail Ben Ayed, Luke McCaffrey, and Eric Granger.  Deep active learning for joint classification & segmentation with weak annotator.  In _WACV_ , 2021.
  * [Beluch et al.(2018)Beluch, Genewein, Nürnberger, and Köhler] William H Beluch, Tim Genewein, Andreas Nürnberger, and Jan M Köhler.  The power of ensembles for active learning in image classification.  In _CVPR_ , 2018.
  * [Bhunia et al.(2020)Bhunia, Yang, Hospedales, Xiang, and Song] Ayan Kumar Bhunia, Yongxin Yang, Timothy M Hospedales, Tao Xiang, and Yi-Zhe Song.  Sketch less for more: On-the-fly fine-grained sketch-based image retrieval.  In _CVPR_ , 2020.
  * [Bhunia et al.(2021)Bhunia, Chowdhury, Sain, Yang, Xiang, and Song] Ayan Kumar Bhunia, Pinaki Nath Chowdhury, Aneeshan Sain, Yongxin Yang, Tao Xiang, and Yi-Zhe Song.  More photos are all you need: Semi-supervised learning for fine-grained sketch based image retrieval.  In _CVPR_ , 2021.
  * [Budd et al.(2021)Budd, Robinson, and Kainz] Samuel Budd, Emma C Robinson, and Bernhard Kainz.  A survey on active learning and human-in-the-loop deep learning for medical image analysis.  _Medical Image Analysis_ , 2021.
  * [Bui et al.(2018)Bui, Ribeiro, Ponti, and Collomosse] Tu Bui, Leonardo Ribeiro, Moacir Ponti, and John Collomosse.  Deep manifold alignment for mid-grain sketch based image retrieval.  In _ACCV_ , 2018.
  * [Caramalau et al.(2021)Caramalau, Bhattarai, and Kim] Razvan Caramalau, Binod Bhattarai, and Tae-Kyun Kim.  Active learning for bayesian 3d hand pose estimation.  In _WACV_ , 2021.
  * [Chan et al.(2020)Chan, Vijayanarasimhan, Ross, and Canny] David M Chan, Sudheendra Vijayanarasimhan, David A Ross, and John F Canny.  Active learning for video description with cluster-regularized ensemble ranking.  In _ACCV_ , 2020.
  * [Choi et al.(2021)Choi, Yi, Kim, Choo, Kim, Chang, Gwon, and Chang] Jongwon Choi, Kwang Moo Yi, Jihoon Kim, Jinho Choo, Byoungjip Kim, Jinyeop Chang, Youngjune Gwon, and Hyung Jin Chang.  VaB-AL: Incorporating class imbalance and difficulty with variational bayes for active learning.  In _CVPR_ , 2021.
  * [Collomosse et al.(2017)Collomosse, Bui, Wilber, Fang, and Jin] John Collomosse, Tu Bui, Michael J Wilber, Chen Fang, and Hailin Jin.  Sketching with style: Visual search with sketches and aesthetic context.  In _ICCV_ , 2017.
  * [Deng et al.(2009)Deng, Dong, Socher, Li, Li, and Fei-Fei] Jia Deng, Wei Dong, Richard Socher, Li-Jia Li, Kai Li, and Li Fei-Fei.  Imagenet: A large-scale hierarchical image database.  In _CVPR_ , 2009.
  * [Hoffer and Ailon(2015)] Elad Hoffer and Nir Ailon.  Deep metric learning using triplet network.  In _SIMBAD_ , 2015.
  * [Huang et al.(2017)Huang, Liu, Van Der Maaten, and Weinberger] Gao Huang, Zhuang Liu, Laurens Van Der Maaten, and Kilian Q Weinberger.  Densely connected convolutional networks.  In _CVPR_ , 2017.
  * [Karen Simonyan(2015)] Andrew Zisserman Karen Simonyan.  Very deep convolutional networks for large-scale image recognition.  In _ICLR_ , 2015.
  * [Kim et al.(2021)Kim, Park, Kim, and Chun] Kwanyoung Kim, Dongwon Park, Kwang In Kim, and Se Young Chun.  Task-aware variational adversarial active learning.  In _CVPR_ , 2021.
  * [Kingma and Ba(2014)] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.
  * [Kirsch et al.(2019)Kirsch, Van Amersfoort, and Gal] Andreas Kirsch, Joost Van Amersfoort, and Yarin Gal.  Batchbald: Efficient and diverse batch acquisition for deep bayesian active learning.  In _NeurIPS_ , 2019.
  * [Li et al.(2019)Li, Wu, Liu, Zhang, and Guan] Qimai Li, Xiao-Ming Wu, Han Liu, Xiaotong Zhang, and Zhichao Guan.  Label efficient semi-supervised learning via graph filtering.  In _CVPR_ , 2019.
  * [Lin et al.(2014)Lin, Maire, Belongie, Hays, Perona, Ramanan, Dollár, and Zitnick] Tsung-Yi Lin, Michael Maire, Serge Belongie, James Hays, Pietro Perona, Deva Ramanan, Piotr Dollár, and C Lawrence Zitnick.  Microsoft coco: Common objects in context.  In _ECCV_ , 2014.
  * [Morgado et al.(2021)Morgado, Vasconcelos, and Misra] Pedro Morgado, Nuno Vasconcelos, and Ishan Misra.  Audio-visual instance discrimination with cross-modal agreement.  In _CVPR_ , 2021.
  * [Pang et al.(2017)Pang, Song, Xiang, and Hospedales] Kaiyue Pang, Yi-Zhe Song, Tony Xiang, and Timothy M Hospedales.  Cross-domain generative learning for fine-grained sketch-based image retrieval.  In _BMVC_ , 2017.
  * [Pang et al.(2019)Pang, Li, Yang, Zhang, Hospedales, Xiang, and Song] Kaiyue Pang, Ke Li, Yongxin Yang, Honggang Zhang, Timothy M Hospedales, Tao Xiang, and Yi-Zhe Song.  Generalising fine-grained sketch-based image retrieval.  In _CVPR_ , 2019.
  * [Pang et al.(2020)Pang, Yang, Hospedales, Xiang, and Song] Kaiyue Pang, Yongxin Yang, Timothy M Hospedales, Tao Xiang, and Yi-Zhe Song.  Solving mixed-modal jigsaw puzzle for fine-grained sketch-based image retrieval.  In _CVPR_ , 2020.
  * [Paszke et al.(2017)Paszke, Gross, Chintala, Chanan, Yang, DeVito, Lin, Desmaison, Antiga, and Lerer] Adam Paszke, Sam Gross, Soumith Chintala, Gregory Chanan, Edward Yang, Zachary DeVito, Zeming Lin, Alban Desmaison, Luca Antiga, and Adam Lerer.  Automatic differentiation in pytorch.  2017\.
  * [Ren et al.(2020)Ren, Xiao, Chang, Huang, Li, Chen, and Wang] Pengzhen Ren, Yun Xiao, Xiaojun Chang, Po-Yao Huang, Zhihui Li, Xiaojiang Chen, and Xin Wang.  A survey of deep active learning.  _arXiv preprint arXiv:2009.00236_ , 2020.
  * [Riaz Muhammad et al.(2018)Riaz Muhammad, Yang, Song, Xiang, and Hospedales] Umar Riaz Muhammad, Yongxin Yang, Yi-Zhe Song, Tao Xiang, and Timothy M Hospedales.  Learning deep sketch abstraction.  In _CVPR_ , 2018.
  * [Russakovsky et al.(2015)Russakovsky, Deng, Su, Krause, Satheesh, Ma, Huang, Karpathy, Khosla, Bernstein, et al.] Olga Russakovsky, Jia Deng, Hao Su, Jonathan Krause, Sanjeev Satheesh, Sean Ma, Zhiheng Huang, Andrej Karpathy, Aditya Khosla, Michael Bernstein, et al.  Imagenet large scale visual recognition challenge.  _International journal of computer vision_ , 115:211–252, 2015.
  * [Sain et al.(2020)Sain, Bhunia, Yang, Xiang, and Song] Aneeshan Sain, Ayan Kumar Bhunia, Yongxin Yang, Tao Xiang, and Yi-Zhe Song.  Cross-modal hierarchical modelling for fine-grained sketch based image retrieval.  In _BMVC_ , 2020.
  * [Sain et al.(2021)Sain, Bhunia, Yang, Xiang, and Song] Aneeshan Sain, Ayan Kumar Bhunia, Yongxin Yang, Tao Xiang, and Yi-Zhe Song.  Stylemeup: Towards style-agnostic sketch-based image retrieval.  In _CVPR_ , 2021.
  * [Sain et al.(2022)Sain, Bhunia, Potlapalli, Chowdhury, Xiang, and Song] Aneeshan Sain, Ayan Kumar Bhunia, Vaishnav Potlapalli, Pinaki Nath Chowdhury, Tao Xiang, and Yi-Zhe Song.  Sketch3t: Test-time training for zero-shot sbir.  In _CVPR_ , 2022.
  * [Sain et al.(2023a)Sain, Bhunia, Chowdhury, Koley, Xiang, and Song] Aneeshan Sain, Ayan Kumar Bhunia, Pinaki Nath Chowdhury, Subhadeep Koley, Tao Xiang, and Yi-Zhe Song.  Clip for all things zero-shot sketch-based image retrieval, fine-grained or not.  In _CVPR_ , 2023a.
  * [Sain et al.(2023b)Sain, Bhunia, Koley, Chowdhury, Chattopadhyay, Xiang, and Song] Aneeshan Sain, Ayan Kumar Bhunia, Subhadeep Koley, Pinaki Nath Chowdhury, Soumitri Chattopadhyay, Tao Xiang, and Yi-Zhe Song.  Exploiting unlabelled photos for stronger fine-grained sbir.  In _CVPR_ , 2023b.
  * [Sener and Savarese(2018)] Ozan Sener and Silvio Savarese.  Active learning for convolutional neural networks: A core-set approach.  In _ICLR_ , 2018.
  * [Siddiqui et al.(2020)Siddiqui, Valentin, and Nießner] Yawar Siddiqui, Julien Valentin, and Matthias Nießner.  Viewal: Active learning with viewpoint entropy for semantic segmentation.  In _CVPR_ , 2020.
  * [Sinha et al.(2019)Sinha, Ebrahimi, and Darrell] Samarth Sinha, Sayna Ebrahimi, and Trevor Darrell.  Variational adversarial active learning.  In _ICCV_ , 2019.
  * [Song et al.(2017)Song, Yu, Song, Xiang, and Hospedales] Jifei Song, Qian Yu, Yi-Zhe Song, Tao Xiang, and Timothy M Hospedales.  Deep spatial-semantic attention for fine-grained sketch-based image retrieval.  In _ICCV_ , 2017.
  * [Song et al.(2018)Song, Pang, Song, Xiang, and Hospedales] Jifei Song, Kaiyue Pang, Yi-Zhe Song, Tao Xiang, and Timothy M Hospedales.  Learning to sketch with shortcut cycle consistency.  In _CVPR_ , 2018.
  * [Szegedy et al.(2015)Szegedy, Liu, Jia, Sermanet, Reed, Anguelov, Erhan, Vanhoucke, and Rabinovich] Christian Szegedy, Wei Liu, Yangqing Jia, Pierre Sermanet, Scott Reed, Dragomir Anguelov, Dumitru Erhan, Vincent Vanhoucke, and Andrew Rabinovich.  Going deeper with convolutions.  In _CVPR_ , 2015.
  * [Wang and Shang(2014)] Dan Wang and Yi Shang.  A new active labeling method for deep learning.  In _IJCNN_ , 2014.
  * [Wang et al.(2015)Wang, Kang, and Li] Fang Wang, Le Kang, and Yi Li.  Sketch-based 3d shape retrieval using convolutional neural networks.  In _CVPR_ , 2015.
  * [Wang et al.(2016)Wang, Gao, Zhang, and Li] Haijun Wang, Xinbo Gao, Kaibing Zhang, and Jie Li.  Single-image super-resolution using active-sampling gaussian process regression.  _IEEE TIP_ , 2016.
  * [Wu et al.(2021)Wu, Shen, Huang, Wang, and Yao] Zihan Wu, Rongbo Shen, Junzhou Huang, Liansheng Wang, and Jianhua Yao.  Strongly supervised mitosis detection in breast histopathology images using weak labels.  In _ISBI_ , 2021.
  * [Xie et al.(2020)Xie, Feng, Chen, Sun, Ma, and Song] Shuai Xie, Zunlei Feng, Ying Chen, Songtao Sun, Chao Ma, and Mingli Song.  Deal: Difficulty-aware active learning for semantic segmentation.  In _ACCV_ , 2020.
  * [Yoo and Kweon(2019)] Donggeun Yoo and In So Kweon.  Learning loss for active learning.  In _CVPR_ , 2019.
  * [Yu et al.(2016)Yu, Liu, Song, Xiang, Hospedales, and Loy] Qian Yu, Feng Liu, Yi-Zhe Song, Tao Xiang, Timothy M Hospedales, and Chen-Change Loy.  Sketch me that shoe.  In _CVPR_ , 2016.
  * [Zhang et al.(2020)Zhang, Li, Yang, Wang, Zha, and Huang] Beichen Zhang, Liang Li, Shijie Yang, Shuhui Wang, Zheng-Jun Zha, and Qingming Huang.  State-relabeling adversarial active learning.  In _CVPR_ , 2020.
  * [Zhdanov(2019)] Fedor Zhdanov.  Diverse mini-batch active learning.  _arXiv preprint arXiv:1901.05954_ , 2019.


# Efficient algorithm for generating Pauli coordinates for an arbitrary linear
operator

Daniel Gunlycke<EMAIL_ADDRESS>U.S. Naval Research Laboratory,
Washington, DC 20375, United States  Mark C. Palenik  U.S. Naval Research
Laboratory, Washington, DC 20375, United States  Alex R. Emmert  SEAP Student,
U.S. Naval Research Laboratory, Washington, DC 20375, United States  Sean A.
Fischer  U.S. Naval Research Laboratory, Washington, DC 20375, United States

###### Abstract

Several linear algebra routines for quantum computing use a basis of tensor
products of identity and Pauli operators to describe linear operators, and
obtaining the coordinates for any given linear operator from its matrix
representation requires a basis transformation, which for an
$\mathrm{N}\times\mathrm{N}$ matrix generally involves
$\mathcal{O}(\mathrm{N}^{4})$ arithmetic operations. Herein, we present an
efficient algorithm that for our particular basis transformation only involves
$\mathcal{O}(\mathrm{N}^{2}\log_{2}\mathrm{N})$ operations. Because this
algorithm requires fewer than $\mathcal{O}(\mathrm{N}^{3})$ operations, for
large $\mathrm{N}$, it could be used as a preprocessing step for quantum
computing algorithms for certain applications. As a demonstration, we apply
our algorithm to a Hamiltonian describing a system of relativistic interacting
spin-zero bosons and calculate the ground-state energy using the variational
quantum eigensolver algorithm on a quantum computer.

##  I Introduction

It has long been recognized that quantum computing offers inherent advantages
over classical computing Benioff (1980); Feynman (1982), allowing quantum
computers to solve certain mathematical tasks otherwise unfeasible Shor
(1994); Grover (1996); Lloyd (1996). Recent advances in fabricating quantum
computers with quantum registers containing tens of qubits have stimulated
renewed efforts in making practical quantum algorithms to run on these quantum
processors, offering polynomial, if not exponential speedup over corresponding
classical calculations Arute  _et al._ (2019). Such algorithms are now
increasingly being used in various application areas, including quantum
simulation of quantum chemical systems Kandala  _et al._ (2017, 2019); Fischer
and Gunlycke (2019); Rubin  _et al._ (2020).

One significant challenge in the implementation of existing algorithms, and
potentially the development of new ones is the loading of information into the
quantum computer Aaronson (2015). Furthermore, the simulation of Hamiltonian
dynamics Lloyd (1996), the variational quantum eigensolver (VQE) algorithm
Peruzzo  _et al._ (2014), and the HHL algorithm for linear systems of
equations Harrow  _et al._ (2009) all assume that the linear operator of
interest be described in what we herein refer to as the Pauli representation.
This representation is a vector representation in a unique basis composed of
tensor products of identity and Pauli operators. Such a basis is particularly
amenable to quantum computing, whether for the evaluation of expectation
values or the representation of unitary operators used in the quantum logic
gates. While there are efficient methods for certain classes of Hamiltonians
Ortiz  _et al._ (2001); Bravyi and Kitaev (2002), a general transformation
that could be applied to arbitrary operators is desirable to expand the reach
of these algorithms and create the opportunity for new algorithms.

In this article, we present an algorithm that derives the Pauli coordinates in
the Pauli representation of an arbitrary linear operator on a finite-
dimensional vector space from its matrix representation. As a bijective map
from one $\mathrm{N}\times\mathrm{N}$ matrix to another is described by
$\mathrm{N}^{4}$ coefficients, an algorithm performing a general
transformation of linear operators described by $\mathrm{N}\times\mathrm{N}$
matrices would involve at least $\mathcal{O}(\mathrm{N}^{4})$ arithmetic
operations. As most linear algebra tasks, including eigenvalue decomposition,
can be accomplished in $\mathcal{O}(\mathrm{N}^{3})$ operations, having an
$\mathcal{O}(\mathrm{N}^{4})$ preprocessing step for any subsequent quantum
algorithm would render the overall method inefficient. Herein however, we show
that the particular transformation from a matrix representation of a linear
operator to the corresponding Pauli representation can be accomplished using
an algorithm requiring only $\mathcal{O}(\mathrm{N}^{2}\log_{2}\mathrm{N})$
operations, or less. This efficient preprocessing algorithm opens the
possibility for quantum algorithms of a variety of linear algebra tasks to
display an overall quantum advantage.

In the next section, Sec. II, we present an efficient method for generating
the sought Pauli coordinates in the Pauli representation of an arbitrary
linear operator. This method is the basis for the algorithm presented and
analyzed in Sec. III. As an illustration, Sec. IV presents calculations that
produce the ground-state energy for a system of relativistic interacting spin-
zero bosons using the transformation algorithm herein as a preprocessing step.
This section covers the derivation of input matrix elements of the Hamiltonian
on a restricted Fock space, the transformation itself, and subsequent ground-
state energy calculations performed using exact diagonalization and the VQE
algorithm with the quantum portion executed on an IBM Q Simulator and the IBM
Q Santiago quantum computer. Lastly, Sec. V summarizes the main conclusions.

##  II Method

Consider an arbitrary finite linear system and let the state space for this
system be the $\mathrm{n}$-dimensional vector space $\mathcal{V}_{\mathrm{n}}$
over the real or complex field $\mathbb{F}$. The system can then be described
by a linear operator $\hat{A}$ in the set of endomorphisms
$\operatorname{End}(\mathcal{V}_{\mathrm{n}})$. Herein, we define $\hat{A}$ by
its matrix representation

$A=\begin{pmatrix}a_{0,0}&\ldots&a_{0,j}&\ldots&a_{0,n-1}\\\
\vdots&\ddots&&&\vdots\\\ a_{i,0}&&a_{i,j}&&a_{i,n-1}\\\
\vdots&&&\ddots&\vdots\\\
a_{n-1,0}&\ldots&a_{n-1,j}&\ldots&a_{n-1,n-1}\end{pmatrix},$ (1)

and assume that all matrix elements $a_{i,j}\in\mathbb{F}$ are given.

Our objective is to represent the linear operator $\hat{\mathrm{A}}$ on the
$\mathrm{N}$-dimensional Hilbert space $\mathcal{H}^{\mathrm{Q}}$ for a
quantum register comprising $\mathrm{Q}$ identical qubits using the basis
composed of tensor products of identity and Pauli operators, where
$\mathcal{H}^{\mathrm{Q}}$ denotes the $\mathrm{Q}^{\mathrm{th}}$ tensor power
of the two-dimensional Hilbert space $\mathcal{H}$. To keep track of the
qubits in this register, we introduce the index set
$\mathcal{Q}=\\{0,1,...,\mathrm{Q}-1\\}$ and label their spaces
$\mathcal{H}_{q}$, where $\mathcal{H}_{q}=\mathcal{H}$, for all
$q\in\mathcal{Q}$.

Our first task is to inject $\mathcal{V}_{\mathrm{n}}$ into
$\mathcal{V}_{\mathrm{N}}\cong\mathcal{H}^{\mathrm{Q}}$, which requires that
$\mathrm{n}\leq\mathrm{N}$, where $\mathrm{N}=2^{\mathrm{Q}}$. We recommend
letting $\mathrm{Q}=\lceil\log_{2}\mathrm{n}\rceil$, but any larger
$\mathrm{Q}$ will work as well. In any case, we define our injection such that
the elements

$\mathrm{a}_{i,j}=\left\\{\begin{array}[]{ll}a_{i,j},&\mathrm{for}~{}i,j\in\\{0,1,...,\mathrm{n}-1\\},\\\\[5.0pt]
\Delta\delta_{i,j},&\mathrm{for}~{}i,j\in\\{\mathrm{n},\mathrm{n}+1,...,\mathrm{N}-1\\},\\\\[5.0pt]
0,&\mathrm{otherwise},\end{array}\right.$ (2)

form the matrix representation for
$\hat{\mathrm{A}}\in\operatorname{End}(\mathcal{V}_{\mathrm{N}})$, where
$\Delta\in\mathbb{F}$ is a constant and $\delta$ is the Kronecker delta.
Although $\Delta$ can take any value, including zero, depending on how the
final produced representation of $\hat{\mathrm{A}}$ is going to be used, it
can be advantageous to choose $\Delta$ such that its absolute value is large
compared to those of all other eigenvalues of $\hat{\mathrm{A}}$.

Next, we derive an equivalent representation of the linear operator
$\hat{\mathrm{A}}\in\operatorname{End}(\mathcal{H}^{\mathrm{Q}})$ that takes
advantage of the tensor product form of the full transformation monoid

$\operatorname{End}(\mathcal{H}^{\mathrm{Q}})=\operatorname{End}(\mathcal{H}_{1})\otimes\operatorname{End}(\mathcal{H}_{2})\otimes\ldots\otimes\operatorname{End}(\mathcal{H}_{\mathrm{Q}}).$
(3)

Consider the operators $\hat{\tau}_{0}$, $\hat{\tau}_{1}$, $\hat{\tau}_{2}$,
and $\hat{\tau}_{3}$ defined by their respective matrix representations

$\begin{array}[]{ccc}\tau_{0}\\!=\\!\left(\\!\begin{array}[]{cc}1&0\\\
0&0\end{array}\\!\right)\\!,~{}\tau_{1}\\!=\\!\left(\\!\begin{array}[]{cc}0&1\\\
0&0\end{array}\\!\right)\\!,~{}\tau_{2}\\!=\\!\left(\\!\begin{array}[]{cc}0&0\\\
1&0\end{array}\\!\right)\\!,~{}\tau_{3}\\!=\\!\left(\\!\begin{array}[]{cc}0&0\\\
0&1\end{array}\\!\right)\\!.\end{array}$ (4)

It is clear that the set
$\\{\hat{\tau}_{0},\hat{\tau}_{1},\hat{\tau}_{2},\hat{\tau}_{3}\\}$ forms a
basis for $\operatorname{End}(\mathcal{H}_{q})$, for each $q\in\mathcal{Q}$.
Using the operators in
$\\{\hat{\tau}_{0},\hat{\tau}_{1},\hat{\tau}_{2},\hat{\tau}_{3}\\}$, we can
construct tensor products of the form

$\hat{\mathrm{T}}_{r_{\mathcal{Q}}}=\hat{\tau}_{r_{\mathrm{Q}-1}}\otimes\hat{\tau}_{r_{\mathrm{Q}-2}}\otimes\ldots\otimes\hat{\tau}_{r_{0}},$
(5)

where $r_{\mathcal{Q}}=(r_{q})_{q\in\mathcal{Q}}$ are families of elements in
$\mathcal{R}=\\{0,1,2,3\\}$ indexed by $\mathcal{Q}$. As the set
$\\{\hat{\mathrm{T}}_{r_{\mathcal{Q}}}\\}$ of all such tensor products forms a
basis for $\operatorname{End}(\mathcal{H}^{\mathrm{Q}})$, we have the linear
combination

$\hat{\mathrm{A}}=\sum_{r_{\mathcal{Q}}}c^{(0)}_{r_{\mathcal{Q}}}\,\hat{\mathrm{T}}_{r_{\mathcal{Q}}},$
(6)

where the coordinates $c^{(0)}_{r_{\mathcal{Q}}}\in\mathbb{F}$ uniquely
specifies our linear operator $\hat{\mathrm{A}}$. In other words, there is an
isomorphism between the matrix representation
$\mathbb{F}^{\mathrm{N}\times\mathrm{N}}$ and the vector representation
$\mathbb{F}^{4^{\mathrm{Q}}}$ of $\hat{\mathrm{A}}$ in the basis
$\\{\hat{\mathrm{T}}_{r_{\mathcal{Q}}}\\}$. This isomorphism can be described
by a one-to-one correspondence between the set of matrix elements
$\\{\mathrm{a}_{i,j}\\}$ from Eq. (2) and the coordinate set
$\\{\mathrm{c}^{(0)}_{r_{\mathcal{Q}}}\\}$ that copies the elements

$\mathrm{a}_{i,j}\mapsto\mathrm{c}^{(0)}_{r_{\mathcal{Q}}}=\mathrm{a}_{i,j},$
(7)

for all $i,j\in\mathcal{N}=\\{0,1,...,\mathrm{N}\\}$, in accordance with a
separate bijection between the sets $\\{(i,j)\\}$ and $\\{r_{\mathcal{Q}}\\}$.

Before defining this latter bijection, let us first associate the indices $i$
and $j$ with the families $i_{\mathcal{Q}}=(i_{q})_{q\in\mathcal{Q}}$ and
$j_{\mathcal{Q}}=(j_{q})_{q\in\mathcal{Q}}$, respectively, of elements in
$\mathcal{B}=\\{0,1\\}$ indexed by $\mathcal{Q}$. We generate these families
exploiting the isomorphisms $i_{\mathcal{Q}}\cong(i)_{2}$ and
$j_{\mathcal{Q}}\cong(j)_{2}$, where $(i)_{2}$ and $(j)_{2}$ denote the binary
representations of $i$ and $j$, respectively. Formally, we identify the
elements $i_{q}$ and $j_{q}$ with the $q^{\mathrm{th}}$ digits in $(i)_{2}$
and $(j)_{2}$, respectively, for all $q\in\mathcal{Q}$. Lastly, for each
generated pair of families $i_{\mathcal{Q}}$ and $j_{\mathcal{Q}}$, we produce
the family $r_{\mathcal{Q}}$ using the mapping
$\mathcal{B}\times\mathcal{B}\rightarrow\mathcal{R}$ defined by

$(i_{q},j_{q})\mapsto r_{q}=2*i_{q}+j_{q},$ (8)

again for all $q\in\mathcal{Q}$. See Table 1 below. With a well-defined
bijection from $\\{(i,j)\\}$ to $\\{r_{\mathcal{Q}}\\}$, we can now obtain all
the coordinates $\mathrm{c}^{(0)}_{r_{\mathcal{Q}}}$ in Eq. (6) using Eq. (7).

Equipped with a vector representation of the linear operator
$\hat{\mathrm{A}}\in\operatorname{End}(\mathcal{H}^{\mathrm{Q}})$, we next
formulate a basis transformation of $\operatorname{End}(\mathcal{H}_{q})$ to
the basis composed of operators in
$\\{\hat{\mathrm{I}},\hat{\sigma}_{x},\hat{\sigma}_{y},\hat{\sigma}_{z}\\}$,
where $\hat{\mathrm{I}}$ is the two-dimensional identity operator and
$\hat{\sigma}_{x}$, $\hat{\sigma}_{y}$, and $\hat{\sigma}_{z}$ are the three
Pauli operators. However, rather than using the basis
$\\{\hat{\mathrm{I}},\hat{\sigma}_{x},\hat{\sigma}_{y},\hat{\sigma}_{z}\\}$
for $\mathcal{H}$ directly in the construction of tensor products for
$\mathcal{H}^{\mathrm{Q}}$, let us first introduce the closely related basis
$\\{\hat{\sigma}_{0},\hat{\sigma}_{1},\hat{\sigma}_{2},\hat{\sigma}_{3}\\}$
defined by the operators $\hat{\sigma}_{0}=\hat{\mathrm{I}}$,
$\hat{\sigma}_{1}=\hat{\sigma}_{x}$, $\hat{\sigma}_{2}=-i\hat{\sigma}_{y}$,
and $\hat{\sigma}_{3}=\hat{\sigma}_{z}$, where $i$ is the imaginary unit. This
basis has the advantage that all nonzero elements of the matrix
representations of its basis operators,

$\begin{array}[]{ccc}\sigma_{0}\\!=\\!\left(\\!\begin{array}[]{cc}1&0\\\
0&1\end{array}\\!\right)\\!,\,\sigma_{1}\\!=\\!\left(\\!\begin{array}[]{cc}0&1\\\
1&0\end{array}\\!\right)\\!,\,\sigma_{2}\\!=\\!\left(\\!\begin{array}[]{cc}0&\\!-1\\\
1&0\end{array}\\!\\!\right)\\!,\,\sigma_{3}\\!=\\!\left(\\!\begin{array}[]{cc}1&0\\\
0&\\!-1\end{array}\\!\\!\right)\\!,\end{array}$ (9)

respectively, are contained in $\mathcal{S}=\\{\pm 1\\}$, which in some
instances might result in a faster implementation of the algorithm below.

The tensor products of the operators in
$\\{\hat{\mathrm{I}},\hat{\sigma}_{x},\hat{\sigma}_{y},\hat{\sigma}_{z}\\}$
can then be expressed on the form

$\hat{\mathrm{S}}_{r_{\mathcal{Q}}}=\Theta_{r_{\mathcal{Q}}}\,\hat{\sigma}_{r_{\mathrm{Q}-1}}\otimes\hat{\sigma}_{r_{\mathrm{Q}-2}}\otimes\ldots\otimes\hat{\sigma}_{r_{0}},$
(10)

where the phase factor

$\Theta_{r_{\mathcal{Q}}}=i\,^{\sum_{q}\\!\delta_{r_{q},2}}$ (11)

ensures that we form the desired basis
$\\{\hat{\mathrm{S}}_{r_{\mathcal{Q}}}\\}$ for
$\operatorname{End}(\mathcal{H}^{\mathrm{Q}})$. In this basis, the operator
$\hat{\mathrm{A}}\in\operatorname{End}(\mathcal{H}^{\mathrm{Q}})$ is given by
the linear combination

$\hat{\mathrm{A}}=\sum_{r_{\mathcal{Q}}}\mathrm{c}_{r_{\mathcal{Q}}}\,\hat{\mathrm{S}}_{r_{\mathcal{Q}}},$
(12)

where $\mathrm{c}_{r_{\mathcal{Q}}}\in\mathbb{F}$ are the $4^{\mathrm{Q}}$
Pauli coordinates we seek.

The basis transformation from $\\{\hat{\mathrm{T}}_{r_{\mathcal{Q}}}\\}$ to
$\\{\hat{\mathrm{S}}_{r_{\mathcal{Q}}}\\}$ is given by

$\hat{\mathrm{S}}_{r_{\mathcal{Q}}}=\Theta_{r_{\mathcal{Q}}}\sum_{r_{\mathcal{Q}}^{\prime}}\mathrm{M}_{r_{\mathcal{Q}},r_{\mathcal{Q}}^{\prime}}\,\hat{\mathrm{T}}_{r_{\mathcal{Q}}^{\prime}},$
(13)

for all families $r_{\mathcal{Q}}\in\mathcal{R}^{\mathrm{Q}}$, where
$\mathrm{M}_{r_{\mathcal{Q}},r_{\mathcal{Q}}^{\prime}}$ are
$4^{\mathrm{Q}}\times 4^{\mathrm{Q}}=\mathrm{N}^{4}$ coefficients describing
the basis transformation operator
$\widehat{\mathrm{M}}\in\operatorname{End}(\operatorname{End}(\mathcal{H}^{\mathrm{Q}}))$.
The transformation operator is of the form

$\widehat{\mathrm{M}}=\bigotimes_{q}\widehat{\mathrm{m}}_{q},$ (14)

and it follows from the definitions in Eqs. (4) and (9) that the matrix
representations of the operators $\widehat{\mathrm{m}}_{q}$ are identical for
all $q\in\mathcal{Q}$ with $\mathrm{m}_{q}=\mathrm{m}$, where

$\mathrm{m}=\begin{pmatrix}1&0&0&1\\\ 0&1&1&0\\\ 0&1&-1&0\\\
1&0&0&-1\end{pmatrix}.$ (15)

This transformation matrix is involutory up to a structure constant,
satisfying $\mathrm{m}^{2}=2\mathrm{I}$, where $\mathrm{I}$ is a unit matrix.
Thus, the corresponding inverse transformation matrix is given by
$\mathrm{m}^{-1}=\mathrm{m}/2$.

Putting everything together, the sought coordinates can now be obtained from
the bijective mapping
$\mathbb{F}^{4^{\mathrm{Q}}}\rightarrow\mathbb{F}^{4^{\mathrm{Q}}}$ defined by
$\mathrm{c}^{(0)}_{r_{\mathcal{Q}}}\mapsto\mathrm{c}_{r_{\mathcal{Q}}}$, for
all $r_{\mathcal{Q}}\in\mathcal{R}^{\mathrm{Q}}$, where we have

$\mathrm{c}_{r_{\mathcal{Q}}}=2^{-\mathrm{Q}}\,\Theta_{r_{\mathcal{Q}}}^{-1}\sum_{r_{\mathcal{Q}}^{\prime}}\prod_{q\in\mathcal{Q}}\mathrm{m}_{r_{q},r_{q}^{\prime}}\,\mathrm{c}^{(0)}_{r_{\mathcal{Q}}^{\prime}},$
(16)

from Eqs. (5–6, 10, 12–14) and the relation $\mathrm{m}^{-1}=\mathrm{m}/2$.
Because the basis transformation in Eq. (13) is separable, as Eq. (14) shows,
we can perform the transformation in $\mathrm{Q}$ independent steps.
Specifically, we have

$\mathrm{c}_{r_{\mathcal{Q}}}=2^{-\mathrm{Q}}\,\Theta_{r_{\mathcal{Q}}}^{-1}\,\mathrm{c}^{(\mathrm{Q})}_{r_{\mathcal{Q}}},$
(17)

where $\mathrm{c}^{(\mathrm{Q})}_{r_{\mathcal{Q}}}$ is obtained after
iterating

$\mathrm{c}^{(q+1)}_{r_{\mathcal{Q}}}=\prod_{q^{\prime}\neq
q}\delta_{r_{q^{\prime}},r_{q^{\prime}}^{\prime}}\sum_{r_{q}^{\prime}}\mathrm{m}_{r_{q},r_{q}^{\prime}}\,\mathrm{c}^{(q)}_{r_{\mathcal{Q}}^{\prime}},$
(18)

for all $r_{\mathcal{Q}}\in\mathcal{R}^{\mathrm{Q}}$, over all
$q\in\mathcal{Q}$. The presence of the delta functions is a manifestation of
the separability of the transformation. Note that it follows from Eq. (15)
that for our particular basis transformation, there are exactly two nonzero
terms in the sum in Eq. (18). Thus, we can express this equation as the linear
combination

$\mathrm{c}^{(q+1)}_{r_{\mathcal{Q}}}=\theta_{r_{q}}\,\mathrm{c}^{(q)}_{r_{\mathcal{Q}}}+\theta_{r_{q}^{*}}\,\mathrm{c}^{(q)}_{r_{\mathcal{Q}}^{*}},$
(19)

for all $q\in\mathcal{Q}$, where
$r_{\mathcal{Q}}^{*}\in\mathcal{R}^{\mathrm{Q}}$ is defined by

$r_{q^{\prime}}^{*}=\left\\{\begin{array}[]{ll}3-r_{q^{\prime}},&\mathrm{for}~{}q^{\prime}=q,\\\\[5.0pt]
r_{q^{\prime}},&\mathrm{for}~{}q^{\prime}\neq q,\end{array}\right.$ (20)

for all $q^{\prime}\in\mathcal{Q}$, and where the two coefficients
$\theta_{r_{q}},\theta_{r_{q}^{*}}\in\mathcal{S}$ follows from Eq. (15) and
are given in Table 1. With these definitions, we can now obtain all the
$4^{\mathrm{Q}}$ Pauli coordinates $\mathrm{c}_{r_{\mathcal{Q}}}$ using Eqs.
(11, 17, 19).

Table 1: Isomorphism between the $q^{\mathrm{th}}$ elements in the families $i_{\mathcal{Q}}$, $j_{\mathcal{Q}}$, and $r_{\mathcal{Q}}$, describing the positions in the matrix and vector representations of the linear operator $\hat{\mathrm{A}}\in\operatorname{End}(\mathcal{H}^{\mathrm{Q}})$. The coefficients $\theta_{r_{q}}$ and $\theta_{r_{q}^{*}}$ describe the coordinate relationship in Eq. (19). $\quad(i_{q},j_{q})\quad$ | $\quad r_{q}\quad$ | $\quad\theta_{r_{q}}\quad$ | $\quad(i_{q}^{*},j_{q}^{*})\quad$ | $\quad r_{q}^{*}\quad$ | $\quad\theta_{r_{q}^{*}}\quad$
---|---|---|---|---|---
$(0,0)$ | $0$ | $+1$ | $(1,1)$ | $3$ | $+1$
$(0,1)$ | $1$ | $+1$ | $(1,0)$ | $2$ | $+1$
$(1,0)$ | $2$ | $-1$ | $(0,1)$ | $1$ | $+1$
$(1,1)$ | $3$ | $-1$ | $(0,0)$ | $0$ | $+1$

##  III Algorithm

###  III.1 General matrix representations

To implement the method described above, we use the following algorithm:

  1. 1.

Identify the $\mathrm{n}\times\mathrm{n}$ matrix representing the linear
operator $\hat{A}\in\operatorname{End}(\mathcal{V}_{\mathrm{n}})$ over the
real or complex field $\mathbb{F}$.

  2. 2.

Using Eq. (2), define the matrix elements $\mathrm{a}_{i,j}$ describing
$\hat{\mathrm{A}}\in\operatorname{End}(\mathcal{V}_{\mathrm{N}})$, indexed by
$\mathcal{N}\times\mathcal{N}$, by the matrix elements of $\hat{A}$ and copies
of the free parameter $\Delta\in\mathbb{F}$.

  3. 3.

Associate each matrix index pair $(i,j)$ with a vector index $r_{\mathcal{Q}}$
using Eq. (8), for all $q\in\mathcal{Q}$, where $i_{q}$ and $j_{q}$ are the
$q^{\mathrm{th}}$ digits in binary representations of $i$ and $j$,
respectively.

  4. 4.

Using Eq. (7), map the elements $\mathrm{a}_{i,j}$ to the initial coordinates
$\mathrm{c}^{(0)}_{r_{\mathcal{Q}}}$ describing
$\hat{\mathrm{A}}\in\operatorname{End}(\mathcal{H}^{\mathrm{Q}})$ in the basis
$\\{\hat{\mathrm{T}}_{r_{\mathcal{Q}}}\\}$.

  5. 5.

Iterate over $q\in\mathcal{Q}$, and for each such iteration calculate a new
coordinate set $\\{\mathrm{c}^{(q+1)}_{r_{\mathcal{Q}}}\\}$ using Eq. (19),
where $r_{\mathcal{Q}}^{*}$ and the signs
$\theta_{r_{q}},\theta_{r_{q}^{*}}\in\mathcal{S}$ are defined by Eq. (20) and
Table 1, respectively.

  6. 6.

Generate the final coordinates $\mathrm{c}_{r_{\mathcal{Q}}}$ describing
$\hat{\mathrm{A}}\in\operatorname{End}(\mathcal{H}^{\mathrm{Q}})$ in the
desired basis $\\{\hat{\mathrm{S}}_{r_{\mathcal{Q}}}\\}$ using Eq. (17) with
$\Theta_{r_{\mathcal{Q}}}$ given by Eq. (11).

As the transformations in step 5 requires $\mathrm{Q}$, or equivalently
$\log_{2}\mathrm{N}$ iterations of Eq. (19), each involving the calculation of
$\mathrm{N}^{2}$ coordinates, the total number of coordinates that need to be
calculated is

$\mathrm{L}_{\mathrm{dense}}=\mathrm{N}^{2}\log_{2}\mathrm{N}.$ (21)

From this result, we conclude that the number of needed arithmetic operations
in a fast implementation scales as
$\mathcal{O}(\mathrm{N}^{2}\log_{2}\mathrm{N})$.

The algorithm herein is thus remarkably fast, considering that the
$\mathrm{N}^{2}\times\mathrm{N}^{2}$ matrix $\mathrm{M}$ describing the basis
transformation operator $\widehat{\mathrm{M}}$ in Eq. (14) has
$\mathrm{N}^{4}$ elements and ordinarily would require at least
$\mathcal{O}(\mathrm{N}^{4})$ operations.

Figure 1: Coordinate couplings for $\mathrm{Q}=3$. In each iteration
$q\in\mathcal{Q}$, each of the $\mathrm{N}^{2}=64$ old coordinates
$c^{(q)}_{r_{\mathcal{Q}}}$ with index $r_{\mathcal{Q}}$ (depicted by black
dots) contributes to two new coordinates
$\mathrm{c}^{(q+1)}_{r_{\mathcal{Q}}}$ and
$\mathrm{c}^{(q+1)}_{r_{\mathcal{Q}}^{*}}$, where $r_{\mathcal{Q}}^{*}$ is
given by Eq. (20). The matrix representation of the full transformation
operator therefore contains $\mathrm{N}^{3}$ nonzero elements and cannot be
calculated explicitly in our $\mathcal{O}(\mathrm{N}^{2}\log_{2}\mathrm{N})$
algorithm.

It is also worth considering what would happen if we were to evaluate Eq. (16)
directly without taking advantage of the separability of the transformation.
As each coordinate in Eq. (19) is a linear combination of two previously
obtained coordinates, one finds that after $\mathrm{Q}$ iterations, each final
coordinate $\mathrm{c}_{r_{\mathcal{Q}}}$ is a linear combination of
$2^{\mathrm{Q}}$ distinct initial coordinates
$\mathrm{c}^{(0)}_{r_{\mathcal{Q}}}$. Conversely, because of the
transformation is involutory, each initial coordinate contributes to exactly
$2^{\mathrm{Q}}$ final coordinates. See Fig. 1. This inverse relationship is
useful to track coordinates when the linear operator $\hat{A}$ is represented
by a sparse matrix.

In the general case, there are exactly $\mathrm{N}$ nonzero terms in the sum
in Eq. (16) and a direct evaluation of the $\mathrm{N}^{2}$ coordinates using
Eq. (16) would require at least $\mathcal{O}(\mathrm{N}^{3})$ operations.
Remarkably, we could not even compute explicitly, let alone store the
$\mathrm{N}^{3}$ nonzero matrix elements describing the transformation
$\widehat{\mathrm{M}}$, while retaining the
$\mathcal{O}(\mathrm{N}^{2}\log_{2}\mathrm{N})$ scaling of our algorithm.

Figure 2: Index mappings and couplings in different base representations.
Owing to the qubit-wise mapping in Eq. (8), the binary representations
$(i)_{2}$ and $(j)_{2}$ of the matrix indices $i\cong(i)_{10}$ and
$j\cong(j)_{10}$, respectively, are interlaced in $(r)_{2}$ representing the
coordinate index $r_{\mathcal{Q}}$. Moreover, the coordinate with index
$r_{\mathcal{Q}}$ couples to itself and the coordinate with index
$r_{\mathcal{Q}}^{*}$ given by Eq. (20). The mapping $r_{q}\mapsto
r_{q}^{*}=3-r_{q}$ can be executed by flipping the two bits in the
$q^{\mathrm{th}}$ crumb (gray) of the binary representation of
$r_{\mathcal{Q}}$ to generate $r_{\mathcal{Q}}^{*}$.

Owing to the digit-wise mapping in Eq. (8), each coordinate index element
$r_{q}\in\mathcal{R}$ can be stored in a unit composed of exactly two
bits—also known as a crumb. Collectively, the most significant bits in each
crumb form a binary representation of the row index $i$, and the least
significant bits form a binary representation of the column index $j$. Thus,
we obtain $r_{\mathcal{Q}}$ for each matrix element index pair $(i,j)$ by
“interlacing” $i$ and $j$. For instance, the element at row six, column three
with $i\cong(6)_{10}\cong(110)_{2}$ and $j\cong(3)_{10}\cong(011)_{2}$
immediately maps to the coordinate at
$r_{\mathcal{Q}}\cong(101101)_{2}\cong(231)_{4}$, as shown in the left-hand
side of Fig. 2.

Other operations can also easily be implemented in a bit representation.
First, we generate the coordinate index $r_{\mathcal{Q}}^{*}$ by copying
$r_{\mathcal{Q}}$ and flipping the two bits in the $q^{\mathrm{th}}$ crumb.
Second, we set the sign $\theta_{r_{q}}$ to $+1$ ($-1$), when the
$(2q)^{\mathrm{th}}$ bit is $0$ ($1$). We always have $\theta_{r_{q}^{*}}=+1$.
Lastly, we calculate the exponent in $\Theta_{r_{\mathcal{Q}}}$ by counting
the number of crumbs that equal $(10)_{2}$.

For parallel implementations of the algorithm, it is useful to note that one
can take advantage of the pairing of all the coordinates in
$\\{\mathrm{c}^{(q)}_{r_{\mathcal{Q}}}\\}$. If all
$\mathrm{c}^{(q)}_{r_{\mathcal{Q}}}$ and
$\mathrm{c}^{(q)}_{r_{\mathcal{Q}}^{*}}$ are kept together in the distribution
of the coordinate array, one can update each coordinate pair using

$\begin{pmatrix}\mathrm{c}^{(q+1)}_{r_{\mathcal{Q}}}\\\
\mathrm{c}^{(q+1)}_{r_{\mathcal{Q}}^{*}}\end{pmatrix}=\begin{pmatrix}(-1)^{i_{q}}&1\\\
1&(-1)^{i_{q}^{*}}\end{pmatrix}\begin{pmatrix}\mathrm{c}^{(q)}_{r_{\mathcal{Q}}}\\\
\mathrm{c}^{(q)}_{r_{\mathcal{Q}}^{*}}\end{pmatrix}.$ (22)

This allows the coordinate data to be stored in a single $\mathrm{N}^{2}$
array and be locally updated. Efficient data storage is important, as we have
found that our $\mathrm{N}^{2}\log_{2}\mathrm{N}$ algorithm is fast enough in
our parallel implementation that data storage, which scales as
$\mathrm{N}^{2}$, is the real limiting factor.

###  III.2 Sparse matrix representations

The algorithm for sparse matrix representations of $\hat{A}$ is the same as
that for the general case above, except that we track and only operate on
nonzero coordinates in steps 4–6.

Let the number of initial nonzero coordinates be $l$. In the sparse limit,
nonzero coordinates do not couple to each other within our transformation. The
number of nonzero coordinates then doubles with each iteration, so that the
number of coordinates after iteration $q$ is $l_{q}=2^{q+1}l$. When the number
of initial nonzero coordinates exceeds $\mathrm{N}$, at most $\mathrm{Q}^{*}$
iterations can be performed before $2^{q+1}l$ exceeds the maximum
$\mathrm{N}^{2}$ coordinates. Assuming that the coordinate doubling continues
up to this point, we have
$\mathrm{Q}^{*}=\min\\{\lfloor-\log_{2}\ell\rfloor,\mathrm{Q}\\}$, where
$\ell=l/\mathrm{N}^{2}$ is the initial coordinate density. The produced
coordinate density during iteration $q$ can then be expressed as

$\ell_{q}=\left\\{\begin{array}[]{ll}2^{q+1}\ell,&\mathrm{when}~{}q\leq\mathrm{Q}^{*}-1,\\\\[5.0pt]
1,&\mathrm{when}~{}q>\mathrm{Q}^{*}-1,\end{array}\right.$ (23)

for all $q\in\mathcal{Q}$. This represents the worst case scenario, as in
reality, some nonzero coordinates could couple to each other, which would lead
to a lower density, in particular when $q\lesssim\mathrm{Q}^{*}$.

To determine how the number of operations scales at large $\mathrm{N}$, we
estimate the total number of nonzero coordinates $\mathrm{L}$ that needs to be
calculated. In the worst case scenario mentioned above, the total number of
nonzero coordinates produced in all $\mathrm{Q}$ iterations is

$\mathrm{L}=\mathrm{N}^{2}\sum_{q\in\mathcal{Q}}\ell_{q},$ (24)

where the coordinate density $\ell_{q}$ is given in Eq. (23).

Figure 3: Scaling of the number of nonzero coordinates in our algorithm. The
diamond data points, representing the actual number of calculated coordinates
for random initial coordinate sets, match exactly the values from the
corresponding curves obtained from Eqs. (25) and (21) for $l=1$ and
$l=\mathrm{N}^{2}$, respectively. In the case $l=\mathrm{N}$, the numerical
data points are all below the corresponding curve obtained from either Eqs.
(25) or (26), owing to nonzero coordinate couplings. The dashed lines
represent linear ($\mathrm{k}=1$), quadratic ($\mathrm{k}=2$), and cubic
($\mathrm{k}=3$) scaling for $\mathrm{N}$ large. Even in the worst case
scenario, our algorithm scaling is well below the classical cubic scaling for
conventional linear algebra routines, as well as any implementation using the
full transformation matrix ($\mathrm{k}=4$).

The expression for the total number of coordinates takes two different forms
depending on whether $\mathrm{Q}^{*}$ is equal to or less than $\mathrm{Q}$.
The critical point is when the initial number of nonzero coordinates is
$l=\mathrm{N}$—or equivalently $\ell=1/\mathrm{N}$. For sparse matrix
representations with $l\leq\mathrm{N}$, we have $\mathrm{Q}^{*}=\mathrm{Q}$,
for which the upper bound for the total number of nonzero coordinates that we
need to calculate is

$\mathrm{L}_{\mathrm{sparse}}=2(\mathrm{N}-1)l.$ (25)

This expression shows that in the very sparse region, where $l$ is a fixed
number independent of $\mathrm{N}$, any fast implementation of our algorithm
scales linearly with $\mathrm{N}$. Evidence of this scaling can also be seen
in the blue curve in Fig. 3, representing $l=1$. Specifically, we note that
this curve for large $\mathrm{Q}$ approaches $\mathrm{k}=1$, which corresponds
to $\mathcal{O}(\mathrm{N})$ scaling.

In the case $l$ is a fraction of $\mathrm{N}$, the scaling becomes quadratic.
This scaling is also immediately evident from Eq. (25). This case is
illustrated by the green curve in Fig. 3, representing $l=\mathrm{N}$, which
as expected approach $\mathrm{k}=2$ asymptotically, for large $\mathrm{Q}$,
which corresponds to $\mathcal{O}(\mathrm{N}^{2})$ scaling.

When the initial number of nonzero coordinates $l\geq\mathrm{N}$, the number
of coordinates becomes saturated. The sum in Eq. (24) then splits into two
types of terms. Using $2^{\mathrm{Q}^{*}}\ell\approx 1$, we arrive at the
approximate expression

$\mathrm{L}_{\mathrm{intermediate}}=\mathrm{N}^{2}\left[\log_{2}\mathrm{N}\ell+2(1-\ell)\right].$
(26)

For $\ell=1/\mathrm{N}$, the first term is zero and the right-hand side
coincides with that in Eq. (25) with $l=\mathrm{N}$. Thus, at this critical
point, we have
$\mathrm{L}_{\mathrm{intermediate}}=\mathrm{L}_{\mathrm{sparse}}=2\mathrm{N}(\mathrm{N}-1)$,
which as already mentioned scales quadratically.

Increasing the initial nonzero coordinate density $\ell$ eventually lands us
at the dense limit $\ell=1$. At this point, we have
$\mathrm{L}_{\mathrm{intermediate}}=\mathrm{L}_{\mathrm{dense}}=\mathrm{N}^{2}\log_{2}\mathrm{N}$
as expected. The upper bound for the number of nonzero coordinates that we
need to calculate is shown by the red curve representing $l=\mathrm{N}^{2}$ in
Fig. 3. Because the difference between $\mathcal{O}(\mathrm{N}^{2})$ and
$\mathcal{O}(\mathrm{N}^{2}\log_{2}\mathrm{N})$ is relatively small when
$\mathrm{N}$ is not too large, one needs to make a judgement whether the
overhead of tracking nonzero coordinates is worthwhile on a case by case
basis. In any case, while not quite quadratic, the central point of this work
is that this curve is well below the well known $\mathrm{k}=3$ scaling for
conventional linear algebra routines.

To see how our total-number-of-nonzero-coordinate expressions above hold up in
practice, we have performed coordinate tracking within our algorithm and
applied our implementation to initial coordinate sets with $l$ nonzeros at
random indices $r_{\mathcal{Q}}$. The results are shown as diamonds in Fig. 3
for $l=1,\mathrm{N},\mathrm{N}^{2}$. In the first case ($l=1$) when there is
only a single nonzero element, no nonzero coordinate couplings are possible,
and thus the calculated $\mathrm{L}$ is a lower limit that matches exactly
that given by Eq. (25) for $l=1$. In the third case ($l=\mathrm{N}^{2}$), the
initial coordinate set is dense with $\ell=1$, which means that all
coordinates are generally nonzero. This is the upper limit, which also does
not permit variability. Consequently, the calculated $\mathrm{L}$ is exactly
that of Eq. (21). In between these two limits, we observe variability in the
calculated $\mathrm{L}$ caused by nonzero coordinate couplings. See how the
data points do not exactly match the curve for $l=\mathrm{N}$. Also note that
all data points are below the curve, which is expected as the curve as
mentioned above, is at this critical point not only given by the approximate
Eq. (26) but also the upper bound in Eq. (25) for $l=\mathrm{N}$.

This observation confirms that a finite number of coordinate couplings only
reduce the amount of needed calculations and never increase it. The most
extreme case of this effect occurs, for instance, when the initial coordinate
set represents a diagonal operator. In this case, $\ell_{q}=\mathrm{N}^{-1}$,
for all $q\in\mathcal{Q}$, resulting in
$\mathrm{L}=\mathrm{N}\log_{2}\mathrm{N}$ instead of
$\mathrm{L}=2\mathrm{N}(\mathrm{N}-1)$.

##  IV Simulating relativistic interacting spin-zero bosons

The Jordan-Wigner transform Jordan and Wigner (1928) allows for the efficient
mapping of a Hamiltonian of a fermionic system to a Hamiltonian of a quantum
register by representing the fermionic creation and annihilation operators as
tensor products of identity and Pauli operators Ortiz  _et al._ (2001).
However, no corresponding transformation exists for the bosonic operators that
satisfy the bosonic commutation relations. As a result, an alternative method
is needed to transform bosonic systems.

To illustrate how the method herein could be used, we solve below the ground-
state energy of an interacting system of relativistic spin-zero bosons. As is
customary in quantum field theory, we describe this system by a Hamiltonian

$\hat{H}=\hat{H}_{0}+\hat{V}$ (27)

that separates the terms $\hat{H}_{0}$ describing the corresponding free non-
interacting system and the terms $\hat{V}$ describing the boson interactions.
Free relativistic spin-zero bosons are described by the Klein–Gordon
Hamiltonian

$\hat{H}_{0}=\frac{1}{2}\int\Big{\\{}\pi^{2}(\vec{r})+\big{[}\hbar
c\nabla\phi(\vec{r})\big{]}^{2}+\big{[}mc^{2}\phi(\vec{r})\big{]}^{2}\Big{\\}}\,\mathrm{d}^{3}r,$
(28)

where $\hbar$, $c$, and $m$ are the reduced Planck constant, the speed of
light, and the boson mass, respectively, and $\phi(\vec{r})$ and
$\pi(\vec{r})$ are quantum fields satisfying the commutation relations
$[\phi(\vec{r}),\pi(\vec{r}^{\prime})]=i\delta(\vec{r}-\vec{r}^{\prime})$ and
$[\phi(\vec{r}),\phi(\vec{r}^{\prime})]=[\pi(\vec{r}),\pi(\vec{r}^{\prime})]=0$.
The interacting terms are herein modeled by the quartic interaction

$\hat{V}=\frac{\lambda}{4!}\big{(}\hbar
c\big{)}^{3}\int\phi^{4}(\vec{r})\,\mathrm{d}^{3}r,$ (29)

where $\lambda$ is a dimensionless coupling constant.

We assume that the bosons are contained in a three-dimensional box with
Cartesian dimensions $L_{x}\\!\times\\!L_{y}\\!\times\\!L_{z}$ centered at the
Cartesian coordinate $(L_{x},L_{y},L_{z})/2$. Moreover, we assume that this
box has $M=M_{x}\\!\times\\!M_{y}\\!\times\\!M_{z}$ modes described by the
wave vectors $\vec{k}_{\vec{\mu}}=(k_{\mu_{x}},k_{\mu_{y}},k_{\mu_{z}})$,
where $\vec{\mu}=(\mu_{x},\mu_{y},\mu_{z})$ and
$k_{\mu_{\alpha}}L_{\alpha}=\mu_{\alpha}\pi$, for all integers
$\mu_{\alpha}\in\mathcal{M}_{\alpha}=\\{1,2,...,M_{\alpha}\\}$, for all
Cartesian components $\alpha\in\\{x,y,z\\}$.

With these modes, we can define the fields

$\displaystyle\phi(\vec{r})$
$\displaystyle=\sum_{\vec{\mu}}\frac{1}{\sqrt{2\hbar\omega_{\vec{\mu}}}}\left(\hat{a}_{\vec{\mu}}+\hat{a}_{\vec{\mu}}^{\dagger}\right)\prod_{\alpha}\sqrt{\frac{2}{L_{\alpha}}}\sin
k_{\mu_{\alpha}}\\!r_{\alpha},$ (30) $\displaystyle\pi(\vec{r})$
$\displaystyle=\sum_{\vec{\mu}}(-i)\sqrt{\frac{\hbar\omega_{\vec{\mu}}}{2}}\left(\hat{a}_{\vec{\mu}}-\hat{a}_{\vec{\mu}}^{\dagger}\right)\prod_{\alpha}\sqrt{\frac{2}{L_{\alpha}}}\sin
k_{\mu_{\alpha}}\\!r_{\alpha},$ (31)

where $\vec{r}=(r_{x},r_{y},r_{z})$ with $r_{\alpha}=\alpha$, for all
components $\alpha$, and $\hat{a}_{\vec{\mu}}$ and
$\hat{a}_{\vec{\mu}}^{\dagger}$ are annihilation and creation operators for
bosons in mode $\vec{\mu}$ with energy

$\hbar\omega_{\vec{\mu}}=\sqrt{\big{|}\hbar
c\vec{k}_{\vec{\mu}}\big{|}^{2}+\big{(}mc^{2}\big{)}^{2}}.$ (32)

Using the bosonic commutation relations
$[\hat{a}_{\vec{\mu}},\hat{a}_{\vec{\nu}}^{\dagger}]=\delta_{\vec{\mu},\vec{\nu}}$
and
$[\hat{a}_{\vec{\mu}},\hat{a}_{\vec{\nu}}]=[\hat{a}_{\vec{\mu}}^{\dagger},\hat{a}_{\vec{\nu}}^{\dagger}]=0$,
for all
$\vec{\mu},\vec{\nu}\in\mathcal{M}_{x}\times\,\mathcal{M}_{y}\times\,\mathcal{M}_{z}$,
we find as expected that this choice of fields diagonalizes the Hamiltonian of
the free non-interacting system, which becomes

$\hat{H}_{0}=\sum_{\vec{\mu}}\hbar\omega_{\vec{\mu}}\,\hat{a}_{\vec{\mu}}^{\dagger}\hat{a}_{\vec{\mu}},$
(33)

after the the zero-point energy has been dropped. The energy of the
$n_{\vec{\mu}}$ particles in mode $\vec{\mu}$ of the non-interacting system is
then $E_{\vec{\mu}}=\hbar\omega_{\vec{\mu}}n_{\vec{\mu}}$. As neither the mode
frequencies $\hbar\omega_{\vec{\mu}}$ nor the particle numbers $n_{\vec{\mu}}$
can be negative, the vacuum state $|0\rangle$, for which $n_{\vec{\mu}}=0$,
for all $\mu$, is the ground state of the non-interacting system with energy
$E_{0}=0$.

In this representation, the Hamiltonian describing the boson interactions
becomes

$\displaystyle\hat{V}=\sum_{\vec{\mu}\vec{\nu}\vec{\xi}\vec{o}}V_{\vec{\mu}\vec{\nu}\vec{\xi}\vec{o}}\Big{(}$
$\displaystyle
3\delta_{\vec{\mu}\vec{\nu}}\delta_{\vec{\xi}\vec{o}}+6\delta_{\vec{\mu}\vec{\nu}}\hat{a}_{\vec{\xi}}\hat{a}_{\vec{o}}+12\delta_{\vec{\mu}\vec{\nu}}\hat{a}^{\dagger}_{\vec{\xi}}\hat{a}_{\vec{o}}$
$\displaystyle\\!+6\delta_{\vec{\mu}\vec{\nu}}\hat{a}^{\dagger}_{\vec{\xi}}\hat{a}^{\dagger}_{\vec{o}}+\hat{a}_{\vec{\mu}}\hat{a}_{\vec{\nu}}\hat{a}_{\vec{\xi}}\hat{a}_{\vec{o}}+4\hat{a}^{\dagger}_{\vec{\mu}}\hat{a}_{\vec{\nu}}\hat{a}_{\vec{\xi}}\hat{a}_{\vec{o}}$
$\displaystyle\\!+6\hat{a}^{\dagger}_{\vec{\mu}}\hat{a}^{\dagger}_{\vec{\nu}}\hat{a}_{\vec{\xi}}\hat{a}_{\vec{o}}+4\hat{a}^{\dagger}_{\vec{\mu}}\hat{a}^{\dagger}_{\vec{\nu}}\hat{a}^{\dagger}_{\vec{\xi}}\hat{a}_{\vec{o}}+\hat{a}^{\dagger}_{\vec{\mu}}\hat{a}^{\dagger}_{\vec{\nu}}\hat{a}^{\dagger}_{\vec{\xi}}\hat{a}^{\dagger}_{\vec{o}}\Big{)},$
(34)

with the coefficients

$\displaystyle V_{\vec{\mu}\vec{\nu}\vec{\xi}\vec{o}}$
$\displaystyle=\frac{\lambda}{4!}\frac{1}{\Omega}\bigg{(}\frac{\hbar
c}{4}\bigg{)}^{3}\frac{1}{\big{(}2\hbar\big{)}^{2}\\!\\!\sqrt{\vphantom{\hbar}\omega_{\vec{\mu}}\,\omega_{\vec{\nu}}\,\omega_{\vec{\xi}}\,\omega_{\vec{o}}}}$
$\displaystyle\times$
$\displaystyle\prod_{\alpha}\\!\sum_{s_{\mu}s_{\nu}s_{\xi}s_{o}}\\!\\!\\!(-1)^{s_{\mu}+s_{\nu}+s_{\xi}+s_{o}}\delta_{s_{\mu}\mu_{\alpha}+s_{\nu}\nu_{\alpha}+s_{\xi}\xi_{\alpha}+s_{o}o_{\alpha},0},$
(35)

where $\Omega=L_{x}L_{y}L_{z}$ is the volume of the box and the sum is over
all signs $s_{\mu},s_{\nu},s_{\xi},s_{o}\in\mathcal{S}$. Note that the
Kronecker delta provides wave vector conservation for each Cartesian component
$\alpha$.

Because any number of bosons can occupy a single mode, we must limit the
number of particles, as well as the number of modes in the box. Denote the
maximum number of particles $N$. We can then form a basis from all possible
Fock states

$|n_{1},n_{2},...,n_{M}\rangle=\prod_{\vec{\mu}}\frac{1}{\sqrt{n_{\vec{\mu}}\,!}}\big{(}a^{\dagger}_{\vec{\mu}}\big{)}^{n_{\vec{\mu}}}|0\rangle,$
(36)

where the total number of particles $\sum_{\vec{\mu}}n_{\vec{\mu}}\leq N$. The
number of such Fock states, and hence the dimension of our Hilbert space is

$\mathrm{n}=\frac{(M+N)!}{M!\,N!}.$ (37)

Letting the arbitrary linear operator $\hat{A}$ be the Hamiltonian $\hat{H}$,
the matrix elements describing our operator are given by

$a_{i,j}=\langle
n_{1}^{(i)},n_{2}^{(i)},...,n_{M}^{(i)}|\hat{H}|n_{1}^{(j)},n_{2}^{(j)},...,n_{M}^{(j)}\rangle,$
(38)

where $i,j\in\\{0,1,...,\mathrm{n}-1\\}$ label the $\mathrm{n}$ Fock states.
Using the algorithm presented herein, we then define the elements
$\mathrm{a}_{i,j}$ for the matrix representation of the Hamiltonian
$\hat{\mathrm{H}}$ and calculate the coefficients
$\mathrm{c}_{r_{\mathcal{Q}}}$, which we subsequently use as input in the VQE
algorithm. This algorithm requires that the operator is of the form of Eq.
(12).

As a demonstration of this approach, we performed a set of transformations and
VQE calculations for a relativistic system containing four massless
interacting bosons in a box of dimension $2\times 1\times 1$ in units of $L$.
We restricted the number of modes to two; specifically we let $M_{x}=2$ and
$M_{y}=M_{z}=1$. As $N=4$ and $M=2$, the dimension given by Eq. (37) of the
Hilbert space $\mathcal{V}_{\mathrm{n}}$ for our chosen system is
$\mathrm{n}=15$. The Hamiltonian $\hat{H}$ operating on
$\mathcal{V}_{\mathrm{n}}$ has $57$ nonzero matrix elements, which we obtained
from Eq. (38) for various coupling constants $\lambda$. Next, we injected
$\mathcal{V}_{\mathrm{n}}$ into the $16$-dimensional Hilbert space for a
quantum register $\mathcal{H}^{\mathrm{Q}}$ with $\mathrm{Q}=4$ qubits. For
this injection, we chose $\Delta=100$ arb. units, which is in principle large
enough to prevent the fictitious state introduced by the injection from: (1)
becoming the ground state associated with the Hamiltonian $\hat{\mathrm{H}}$
on $\mathcal{H}^{\mathrm{Q}}$, and (2) appreciably affecting the physical
states of our system. In reality, choosing a large $\Delta$ has the drawback
that the noise present in quantum computing calculations, which introduces
weights into fictitious states, results to an artificial increase in the
calculated ground-state energies.

We applied the above transformation to produce the Pauli coordinates
describing $\hat{H}$ for various values of $\lambda$, which we subsequently
used as input to the VQE algorithm. We then executed this algorithm using the
Aer, Aqua, and Ignis application programming interfaces included in IBM Qiskit
Aleksandrowicz  _et al._ (2019). We chose a variational form of the quantum
algorithm with a quantum circuit consisting of RY single-qubit gates and CNOT
entangling gates repeated with the depth of four Fischer and Gunlycke (2019).
We optimized the RY rotations using the StatevectorSimulator and the Powell
optimizer in Qiskit with the convergence criteria that the relative error in
the solution is less than $10^{-4}$. Lastly, we calculated the final ground-
state energies independently on the simulator and the IBM Q Santiago quantum
computer. For each circuit, we averaged the output over 8,192 identical runs
to reduce statistical noise.

Figure 4: Scaled ground-state energy $EL$ in natural units of four massless
relativistic bosons in a two-mode box of dimensions $2\times 1\times 1$ in
units of length $L$ as a function of the dimensionless interaction coupling
constant $\lambda$. The energy has been obtained from exact diagonalization of
the Hamiltonian, as well as from the VQE algorithm executed on an IBM Q
Simulator and the IBM Q Santiago quantum computer. The dotted lines show the
approximate scaled energy $E_{\mathrm{approx}}L=E_{0}L+\beta\lambda L$ for
weak interactions with the slope $\beta\approx 5.374\times 10^{-3}$ obtained
from perturbation theory.

Figure 4 shows the final ground-state energy for different interaction
strengths calculated independently on the simulator and the IBM Q Santiago
quantum computer. Also shown is the ground-state energy obtained from exact
diagonalization of $\hat{\mathrm{H}}$. Though we observe appreciable absolute
errors in the quantum computing calculations, the overall shape of the curve
is similar to that of the exact diagonalization curve, suggesting that the
quantum computing calculations at least produce a reasonable approximation of
the ground state.

Let us now explore the ground-state energy for weak interactions. The ground-
state energy $E$ can then be approximated by the first-order perturbation
energy $E^{(1)}=\langle 0|\hat{V}|0\rangle$, which we can express as

$E^{(1)}=\frac{3\lambda}{4!}\frac{(\hbar
c)^{3}}{4\Omega}\sum_{\vec{\mu}\vec{\xi}}\frac{1}{\hbar\omega_{\vec{\mu}}\,\hbar\omega_{\vec{\xi}}}\prod_{\alpha}\bigg{(}1+\frac{\delta_{\mu_{\alpha},\xi_{\alpha}}}{2}\bigg{)}.$
(39)

Note that this expression is independent of $N$, suggesting that the
interacting ground-state energy depends only weakly on the maximum allowed
number of particles in the box. Irrespective of $N$, there are exactly two
modes $(1,1,1)$ and $(2,1,1)$ in our considered system, which we label $1$ and
$2$, respectively, for short. The ground-state energy for our two-mode
interacting system is then approximately

$E^{(1)}=\frac{3\lambda}{4!}\frac{(\hbar
c)^{3}}{64L}\bigg{[}\frac{27}{(\hbar\omega_{1}L)^{2}}+\frac{27}{(\hbar\omega_{2}L)^{2}}+\frac{16}{(\hbar\omega_{1}L)(\hbar\omega_{2}L)}\bigg{]},$
(40)

where the scaled mode energies from Eq. (32) are $\hbar\omega_{1}L=3\pi\hbar
c/2\approx 4.712$ and $\hbar\omega_{2}L=\sqrt{3}\pi\hbar c\approx 5.441$, in
natural units. The approximate scaled ground-state energy for weak
interactions can thus be expressed $E^{(1)}L=\beta\lambda$, where the slope
$\beta\approx 5.374\times 10^{-3}$. As shown in Fig. 4, $\beta\lambda$
provides an excellent approximation to the exact ground-state energy for weak
interactions. The corresponding estimated slope from our first two data points
calculated on IBM Q Santiago is $\beta_{\mathrm{qc}}\approx 8.192\times
10^{-3}$. Although by no means a perfect estimate, the comparison shows that
the quantum computer we used already has sufficient accuracy to provide
ballpark prediction.

##  V Conclusions

One of the challenges facing quantum computing applications is the
conformation of classical input data to the input requirements of quantum
algorithms such as the VQE and the HHL algorithms. As the preprocessing
algorithm present herein requires no more than
$\mathcal{O}(\mathrm{N}^{2}\log_{2}\mathrm{N})$ arithmetic operations, it
could in conjunction with a quantum or hybrid quantum/classical algorithm
offer an overall speedup over purely classical linear algebra algorithms
requiring $\mathcal{O}(\mathrm{N}^{3})$ arithmetic operations. Further speedup
is still possible for specific linear operators describing systems exhibiting
some form of symmetry. However, the real limitation we ran into for large
$\mathrm{N}$ was not compute time, but rather available memory. The number of
matrix elements or coordinates needed to describe a general linear operator is
$\mathrm{N}^{2}$. We therefore expect that in the foreseeable future, quantum
computing applications will continue to be most powerful for problems
requiring large calculations with limited amount of input data. Even so, our
hope is that the preprocessing algorithm herein will nevertheless make a
broader class of problems accessible to quantum computing.

###### Acknowledgements.

This work has been supported by the Office of Naval Research (ONR) through the
U.S. Naval Research Laboratory (NRL) and by NRL through a Naval Innovative
Science and Engineering (NISE) program. We acknowledge quantum computing
resources from IBM through a collaboration with the Air Force Research
Laboratory (AFRL).

## References

  * Benioff (1980) P. Benioff, J. Stat. Phys. 22, 563 (1980).
  * Feynman (1982) R. P. Feynman, Int. J. Theor. Phys. 21, 467 (1982).
  * Shor (1994) P. W. Shor, in  _Proceedings of the 35th Annual Symposium on Foundations of Computer Science_ (IEEE Computer Society Press, Los Alamitos, CA, 1994) pp. 124–134.
  * Grover (1996) L. K. Grover, in  _Proceedings of the $28^{\mathrm{th}}$ annual ACM symposium on Theory of Computing_ (ACM Publications, New York, NY, 1996) pp. 212–219.
  * Lloyd (1996) S. Lloyd, Science 273, 1073 (1996).
  * Arute  _et al._ (2019) F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, R. Biswas, S. Boixo, F. G. S. L. Brandao, D. A. Buell, B. Burkett, Y. Chen, Z. Chen, B. Chiaro, R. Collins, W. Courtney, A. Dunsworth, E. Farhi, B. Foxen, A. Fowler, C. Gidney, M. Giustina, R. Graff, K. Guerin, S. Habegger, M. P. Harrigan, M. J. Hartmann, A. Ho, M. Hoffmann, T. Huang, T. S. Humble, S. V. Isakov, E. Jeffrey, Z. Jiang, D. Kafri, K. Kechedzhi, J. Kelly, P. V. Klimov, S. Knysh, A. Korotkov, F. Kostritsa, D. Landhuis, M. Lindmark, E. Lucero, D. Lyakh, S. Mandrà, J. R. McClean, M. McEwen, A. Megrant, X. Mi, K. Michielsen, M. Mohseni, J. Mutus, O. Naaman, M. Neeley, C. Neill, M. Y. Niu, E. Ostby, A. Petukhov, J. C. Platt, C. Quintana, E. G. Rieffel, P. Roushan, N. C. Rubin, D. Sank, K. J. Satzinger, V. Smelyanskiy, K. J. Sung, M. D. Trevithick, A. Vainsencher, B. Villalonga, T. White, Z. J. Yao, P. Yeh, A. Zalcman, H. Neven, and J. M. Martinis, Nature 574, 505 (2019).
  * Kandala  _et al._ (2017) A. Kandala, A. Mezzacapo, K. Temme, M. Takita, M. Brink, J. M. Chow, and J. M. Gambetta, Nature 549, 242 (2017).
  * Kandala  _et al._ (2019) A. Kandala, K. Temme, A. D. Córcoles, A. Mezzacapo, J. M. Chow, and J. M. Gambetta, Nature 567, 491 (2019).
  * Fischer and Gunlycke (2019) S. A. Fischer and D. Gunlycke, “Symmetry configuration mapping for representing quantum systems on quantum computers,” (2019), arXiv:1907.01493 [quant-ph] .
  * Rubin  _et al._ (2020) N. C. Rubin, F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, S. Boixo, M. Broughton, B. B. Buckley, D. A. Buell, B. Burkett, N. Bushnell, Y. Chen, Z. Chen, B. Chiaro, R. Collins, W. Courtney, S. Demura, A. Dunsworth, E. Farhi, A. Fowler, B. Foxen, C. Gidney, M. Giustina, R. Graff, S. Habegger, M. P. Harrigan, A. Ho, S. Hong, T. Huang, W. J. Huggins, L. Ioffe, S. V. Isakov, E. Jeffrey, Z. Jiang, C. Jones, D. Kafri, K. Kechedzhi, J. Kelly, S. Kim, P. V. Klimov, A. Korotkov, F. Kostritsa, D. Landhuis, P. Laptev, M. Lindmark, E. Lucero, O. Martin, J. M. Martinis, J. R. McClean, M. McEwen, A. Megrant, X. Mi, M. Mohseni, W. Mruczkiewicz, J. Mutus, O. Naaman, M. Neeley, C. Neill, H. Neven, M. Y. Niu, T. E. O’Brien, E. Ostby, A. Petukhov, H. Putterman, C. Quintana, P. Roushan, D. Sank, K. J. Satzinger, V. Smelyanskiy, D. Strain, K. J. Sung, M. Szalay, T. Y. Takeshita, A. Vainsencher, T. White, N. Wiebe, Z. J. Yao, P. Yeh, and A. Zalcman, Science 369, 1084 (2020).
  * Aaronson (2015) S. Aaronson, Nat. Phys. 11, 291 (2015).
  * Peruzzo  _et al._ (2014) A. Peruzzo, J. R. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik, and J. L. O’Brien, Nat. Commun. 5, 4213 (2014).
  * Harrow  _et al._ (2009) A. W. Harrow, A. Hassidim, and S. Lloyd, Phys. Rev. Lett. 103, 150502 (2009).
  * Ortiz  _et al._ (2001) G. Ortiz, J. E. Gubernatis, E. Knill, and R. Laflamme, Phys. Rev. A 64, 022319 (2001).
  * Bravyi and Kitaev (2002) S. Bravyi and A. Kitaev, Ann. Phys. 298, 210 (2002).
  * Jordan and Wigner (1928) P. Jordan and E. Wigner, Z. Phys. 47, 631 (1928).
  * Aleksandrowicz  _et al._ (2019) G. Aleksandrowicz, T. Alexander, P. Barkoutsos, L. Bello, Y. Ben-Haim, D. Bucher, F. J. Cabrera-Hernádez, J. Carballo-Franquis, A. Chen, C.-F. Chen, J. M. Chow, A. D. Córcoles-Gonzales, A. J. Cross, A. Cross, J. Cruz-Benito, C. Culver, S. D. L. P. González, E. D. L. Torre, D. Ding, E. Dumitrescu, I. Duran, P. Eendebak, M. Everitt, I. F. Sertage, A. Frisch, A. Fuhrer, J. Gambetta, B. G. Gago, J. Gomez-Mosquera, D. Greenberg, I. Hamamura, V. Havlicek, J. Hellmers, Ł. Herok, H. Horii, S. Hu, T. Imamichi, T. Itoko, A. Javadi-Abhari, N. Kanazawa, A. Karazeev, K. Krsulich, P. Liu, Y. Luh, Y. Maeng, M. Marques, F. J. Martín-Fernández, D. T. McClure, D. McKay, S. Meesala, A. Mezzacapo, N. Moll, D. M. Rodríguez, G. Nannicini, P. Nation, P. Ollitrault, L. J. O’Riordan, H. Paik, J. Pérez, A. Phan, M. Pistoia, V. Prutyanov, M. Reuter, J. Rice, A. R. Davila, R. H. P. Rudy, M. Ryu, N. Sathaye, C. Schnabel, E. Schoute, K. Setia, Y. Shi, A. Silva, Y. Siraichi, S. Sivarajah, J. A. Smolin, M. Soeken, H. Takahashi, I. Tavernelli, C. Taylor, P. Taylour, K. Trabing, M. Treinish, W. Turner, D. Vogt-Lee, C. Vuillot, J. A. Wildstrom, J. Wilson, E. Winston, C. Wood, S. Wood, S. Wörner, I. Y. Akhalwaya, and C. Zoufal, “Qiskit: An open-source framework for quantum computing,”  (2019).